Fix OpenMP reproducibility

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()
 

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

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,13 +3435,15 @@ 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?
       do k=1,nVertLevels
          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)
 

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

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

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

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

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

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)

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