diff --git a/cmake/Functions/MPAS_Functions.cmake b/cmake/Functions/MPAS_Functions.cmake index 15d9f63fc1..6ba923273e 100644 --- a/cmake/Functions/MPAS_Functions.cmake +++ b/cmake/Functions/MPAS_Functions.cmake @@ -92,6 +92,7 @@ function(mpas_fortran_target target) list(APPEND MPAS_FORTRAN_TARGET_COMPILE_DEFINITIONS _MPI=1) # Enable OpenMP support if(MPAS_OPENMP) + list(APPEND MPAS_FORTRAN_TARGET_COMPILE_DEFINITIONS MPAS_OPENMP=1) target_link_libraries(${target} PUBLIC OpenMP::OpenMP_Fortran) endif() diff --git a/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F b/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F index d21061b0fb..4a1f2283f3 100644 --- a/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F +++ b/src/core_atmosphere/diagnostics/mpas_pv_diagnostics.F @@ -1029,7 +1029,7 @@ subroutine calc_gradxu_cell(gradxu, addEarthVort, & du_dz = calc_vertDeriv_one(val0, valm, z0-zm) !v val0 = dotProduct(velCell0, xyzLocal(:,2),3) - valm = dotProduct(velCellp, xyzLocal(:,2),3) + valm = dotProduct(velCellm, xyzLocal(:,2),3) dv_dz = calc_vertDeriv_one(val0, valm, z0-zm) else !u diff --git a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F index 227fbde862..4d3bd972d9 100644 --- a/src/core_atmosphere/dynamics/mpas_atm_time_integration.F +++ b/src/core_atmosphere/dynamics/mpas_atm_time_integration.F @@ -16,6 +16,10 @@ module atm_time_integration +#ifdef MPAS_OPENMP + use omp_lib +#endif + use mpas_pool_routines use mpas_kind_types use mpas_constants @@ -3135,9 +3139,18 @@ subroutine atm_rk_integration_setup( state, diag, nVertLevels, num_scalars, & call mpas_pool_get_array(state, 'scalars', scalars_2, 2) - !$acc kernels - theta_m_2(:,cellEnd+1) = 0.0_RKIND - !$acc end kernels +#ifdef MPAS_OPENMP + ! Assign 0 to cellEnd+1 element only on the last thread, because on + ! all other threads cellEnd+1 is the same memory location + ! as cellStart on the following threads + if (omp_get_thread_num() == omp_get_num_threads() - 1) then +#endif + !$acc kernels + theta_m_2(:,cellEnd+1) = 0.0_RKIND + !$acc end kernels +#ifdef MPAS_OPENMP + end if +#endif !$acc parallel default(present) !$acc loop gang worker @@ -3422,6 +3435,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, rcv = rgas/(cp-rgas) c2 = cp*rcv +!$omp single !$acc parallel default(present) !$acc loop gang worker ! MGD bad to have all threads setting this variable? @@ -3429,6 +3443,7 @@ subroutine atm_compute_vert_imp_coefs_work(nCells, moist_start, moist_end, dts, cofrz(k) = rdzw(k) end do !$acc end parallel +!$omp end single !$acc parallel default(present) !$acc loop gang worker private(b_tri,c_tri) @@ -4384,12 +4399,14 @@ subroutine atm_recover_large_step_variables_work(nCells, nEdges, nCellsSolve, nE ! Avoid FP errors caused by a potential division by zero below by ! initializing the "garbage cell" of rho_zz to a non-zero value +!$omp single !$acc parallel default(present) !$acc loop gang vector do k=1,nVertLevels rho_zz(k,nCells+1) = 1.0 end do !$acc end parallel +!$omp end single ! compute new density everywhere so we can compute u from ru. ! we will also need it to compute theta_m below @@ -5294,10 +5311,12 @@ subroutine atm_advance_scalars_mono_work(field_name, block, state, nCells, nEdge end do #ifndef MPAS_OPENACC +!$omp single do k=1,nVertLevels scalar_old(k,nCells+1) = 0.0_RKIND scalar_new(k,nCells+1) = 0.0_RKIND end do +!$omp end single #endif !$acc end parallel @@ -7724,11 +7743,20 @@ subroutine atm_rk_dynamics_substep_finish( state, diag, nVertLevels, dynamics_su call mpas_pool_get_array(state, 'rho_zz', rho_zz_1, 1) call mpas_pool_get_array(state, 'rho_zz', rho_zz_2, 2) - ! Interim fix for the atm_compute_dyn_tend_work subroutine accessing uninitialized values - ! in garbage cells of theta_m - !$acc kernels - theta_m_1(:,cellEnd+1) = 0.0_RKIND - !$acc end kernels +#ifdef MPAS_OPENMP + ! Assign 0 to cellEnd+1 element only on the last thread, because on + ! all other threads cellEnd+1 is the same memory location + ! as cellStart on the following threads + if (omp_get_thread_num() == omp_get_num_threads() - 1) then +#endif + ! Interim fix for the atm_compute_dyn_tend_work subroutine accessing uninitialized values + ! in garbage cells of theta_m + !$acc kernels + theta_m_1(:,cellEnd+1) = 0.0_RKIND + !$acc end kernels +#ifdef MPAS_OPENMP + end if +#endif inv_dynamics_split = 1.0_RKIND / real(dynamics_split) diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F index 0116dcf561..d26ad7c714 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F @@ -357,9 +357,11 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) !In Registry.xml, dzs is a function of nCells. In the Noah lsm scheme, dzs is independent !of cell locations: +!$omp single do n = 1,num_soils dzs_p(n) = dzs(n,its) enddo +!$omp end single do j = jts,jte do n = 1,num_soils diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F index dfa327fea8..9091105a87 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_lw.F @@ -469,9 +469,11 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi if(config_o3climatology) then !ozone mixing ratio: +!$omp single do k = 1, num_oznLevels pin_p(k) = pin(k) enddo +!$omp end single do j = jts,jte do k = 1, num_oznLevels do i = its,ite @@ -491,7 +493,7 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi p2d(i,k) = pres_hyd_p(i,k,j) / 100._RKIND enddo enddo - call vinterp_ozn(1,ncols,ncols,nlevs,p2d,pin_p,num_oznlevels,o3clim_p(1,1,j),o32d) + call vinterp_ozn(1,ncols,ncols,nlevs,p2d,pin_p,num_oznlevels,o3clim_p(its:ite,1:num_oznLevels,j),o32d) do i = its,ite do k = kts,kte o3vmr(k,i) = o32d(i,k) @@ -501,9 +503,11 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi if(allocated(p2d)) deallocate(p2d) if(allocated(o32d)) deallocate(o32d) else +!$omp single do k = 1, num_oznLevels pin_p(k) = 0.0_RKIND enddo +!$omp end single do j = jts,jte do k = 1, num_oznLevels do i = its,ite @@ -571,9 +575,11 @@ subroutine radiation_lw_from_MPAS(xtime_s,configs,mesh,state,time_lev,diag_physi endif !ozone mixing ratio: +!$omp single do k = 1, num_oznlevels pin_p(k) = pin(k) enddo +!$omp end single do n = 1, num_months do j = jts,jte do k = 1, num_oznlevels diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F index 3a14275826..0498d15651 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_radiation_sw.F @@ -598,9 +598,11 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i !ozone volume mixing ratio: if(config_o3climatology) then +!$omp single do k = 1, num_oznLevels pin_p(k) = pin(k) enddo +!$omp end single do j = jts,jte do k = 1, num_oznLevels do i = its,ite @@ -609,9 +611,11 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i enddo enddo else +!$omp single do k = 1, num_oznLevels pin_p(k) = 0.0_RKIND enddo +!$omp end single do j = jts,jte do k = 1, num_oznLevels do i = its,ite @@ -672,9 +676,11 @@ subroutine radiation_sw_from_MPAS(configs,mesh,state,time_lev,diag_physics,atm_i enddo endif !ozone mixing ratio: +!$omp single do k = 1, num_oznlevels pin_p(k) = pin(k) enddo +!$omp end single do n = 1, num_months do j = jts,jte do k = 1, num_oznlevels diff --git a/src/core_atmosphere/physics/mpas_atmphys_sfc_diagnostics.F b/src/core_atmosphere/physics/mpas_atmphys_sfc_diagnostics.F index 6e3a84c394..0360157390 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_sfc_diagnostics.F +++ b/src/core_atmosphere/physics/mpas_atmphys_sfc_diagnostics.F @@ -41,7 +41,6 @@ subroutine atmphys_sfc_diagnostics(configs,mesh,diag,diag_physics,sfc_input,outp !local variables and pointers: character(len=StrKIND),pointer:: lsm_scheme - integer,pointer:: nCellsSolve integer:: i real(kind=RKIND),dimension(:),pointer:: psfc @@ -58,8 +57,6 @@ subroutine atmphys_sfc_diagnostics(configs,mesh,diag,diag_physics,sfc_input,outp call mpas_pool_get_config(configs,'config_lsm_scheme',lsm_scheme) - call mpas_pool_get_dimension(mesh,'nCellsSolve',nCellsSolve) - call mpas_pool_get_array(diag,'surface_pressure',psfc) call mpas_pool_get_array(diag_physics,'chs2',chs2) @@ -77,7 +74,7 @@ subroutine atmphys_sfc_diagnostics(configs,mesh,diag,diag_physics,sfc_input,outp sf_select: select case(trim(lsm_scheme)) case("sf_noah") - do i = 1,nCellsSolve + do i = its,ite rho = psfc(i)/(R_d*tsk(i)) if(cqs2(i) .lt. 1.e-5) then q2(i) = qsfc(i) @@ -95,7 +92,7 @@ subroutine atmphys_sfc_diagnostics(configs,mesh,diag,diag_physics,sfc_input,outp case("sf_noahmp") call mpas_pool_get_array(output_noahmp,'q2mxy',q2mxy) call mpas_pool_get_array(output_noahmp,'t2mxy',t2mxy) - do i = 1,nCellsSolve + do i = its,ite rho = psfc(i)/(R_d*tsk(i)) if((xland(i)-1.5 .gt. 0._RKIND) .or. (xland(i)-1.5.le.0._RKIND .and. xice(i).ge.xice_threshold)) then if(cqs2(i) .lt. 1.e-5) then diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index 134009f537..caee0479ef 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -382,6 +382,7 @@ module mpas_atmphys_vars !do not run urban physics, idiff is set to zero. integer:: ysu_pblmix +!$omp threadprivate(ysu_pblmix) integer,dimension(:,:),allocatable:: & kpbl_p !index of PBL top [-] @@ -670,6 +671,8 @@ module mpas_atmphys_vars integer,parameter:: taer_asy_opt = 3!input option for nwfa, nifa aerosol asymmetry factor. integer:: aer_opt !=[0,3] : 0 for no aerosols, 3 for "water-" and "ice-friendly" aerosols. +!$omp threadprivate(aer_opt) + integer,dimension(:,:),allocatable:: & taer_type_p !=[1,2,3]: 1 for rural, 2 is urban and 3 is maritime in WRF. In MPAS, !aer_type is initialized as a function of landmask (=1 over land; =2 over diff --git a/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F b/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F index 8e24340501..2e569ea114 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F +++ b/src/core_atmosphere/physics/physics_wrf/module_mp_thompson.F @@ -110,6 +110,7 @@ MODULE module_mp_thompson ! REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6 REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6 REAL, PRIVATE:: Nt_c + !$omp threadprivate(Nt_c) !..Declaration of constants for assumed CCN/IN aerosols when none in !.. the input data. Look inside the init routine for modifications @@ -127,6 +128,7 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: mu_g = 0.0 REAL, PARAMETER, PRIVATE:: mu_i = 0.0 REAL, PRIVATE:: mu_c + !$omp threadprivate(mu_c) !..Sum of two gamma distrib for snow (Field et al. 2005). !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3) diff --git a/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F b/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F index ed41320ac7..950746df1d 100644 --- a/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F +++ b/src/core_atmosphere/physics/physics_wrf/module_sf_mynn.F @@ -226,6 +226,8 @@ subroutine sfclay_mynn( & real(kind=RKIND),dimension(its:ite):: u10_hv,v10_hv,th2_hv,t2_hv,q2_hv,wstar_hv,qstar_hv real(kind=RKIND),dimension(its:ite):: cd_hv,cda_hv,ck_hv,cka_hv,ustm_hv + real(kind=RKIND),dimension(its:ite,1:kte):: dz8w_hv,u3d_hv,v3d_hv,p3d_hv,t3d_hv,rho3d_hv,qv3d_hv,qc3d_hv,pattern_spp_pbl_hv + !----------------------------------------------------------------------------------------------------------------- f_spp = .false. @@ -244,9 +246,60 @@ subroutine sfclay_mynn( & qstar_hv(i) = qstar(i,j) enddo - call sf_mynn_pre_run(its,ite,kte,itimestep,dz8w,u3d,v3d,p3d,t3d,rho3d,qv3d,qc3d,f_spp, & - pattern_spp_pbl,ust_hv,mol_hv,qsfc_hv,qstar_hv,dz8w1d,u1d,v1d,p1d,t1d,rho1d, & - qv1d,qc1d,rstoch1d,dz2w1d,u1d2,v1d2,errmsg,errflg) +#if 1 + dz8w_hv = dz8w(its:ite,1:kte,j) + u3d_hv = u3d(its:ite,1:kte,j) + v3d_hv = v3d(its:ite,1:kte,j) + p3d_hv = p3d(its:ite,1:kte,j) + t3d_hv = t3d(its:ite,1:kte,j) + rho3d_hv = rho3d(its:ite,1:kte,j) + qv3d_hv = qv3d(its:ite,1:kte,j) + qc3d_hv = qc3d(its:ite,1:kte,j) + if(f_spp) then + pattern_spp_pbl_hv = pattern_spp_pbl(its:ite,1:kte,j) + end if + + call sf_mynn_pre_run(its,ite,kte,itimestep, & + dz8w_hv,u3d_hv,v3d_hv,p3d_hv,t3d_hv,rho3d_hv,qv3d_hv,qc3d_hv, & + f_spp,pattern_spp_pbl_hv, & + ust_hv,mol_hv,qsfc_hv,qstar_hv, & + dz8w1d,u1d,v1d,p1d,t1d,rho1d,qv1d,qc1d,rstoch1d, & + dz2w1d,u1d2,v1d2,errmsg,errflg) +#else + do i = its,ite + dz8w1d(i) = dz8w(i,1,j) + u1d(i) = u3d(i,1,j) + v1d(i) = v3d(i,1,j) + qv1d(i) = qv3d(i,1,j) + qc1d(i) = qc3d(i,1,j) + p1d(i) = p3d(i,1,j) + t1d(i) = t3d(i,1,j) + rho1d(i) = rho3d(i,1,j) + !--- 2nd model level winds - for diags with high-resolution grids: + dz2w1d(i) = dz8w(i,2,j) + u1d2(i) = u3d(i,2,j) + v1d2(i) = v3d(i,2,j) + enddo + + if(f_spp) then + do i = its,ite + rstoch1d(i) = pattern_spp_pbl(i,1,j) + enddo + else + do i = its,ite + rstoch1d(i)=0._RKIND + enddo + endif + + if(itimestep == 1) then + do i = its,ite + ust_hv(i) = max(0.04*sqrt(u1d(i)*u1d(i) + v1d(i)*v1d(i)),0.001) + mol_hv(i) = 0._RKIND + qsfc_hv(i) = qv1d(i)/(1.+qv1d(i)) + qstar_hv(i) = 0._RKIND + enddo + endif +#endif !input arguments: do i = its,ite