Merge pull request NOAA-GFDL#62 from GEOS-ESM/geos/develop

GitFlow: Merge geos/develop into geos/main
climbfuji · Mar 23, 2023 · a83b965 · a83b965
2 parents ae54f82 + a3440a1
commit a83b965
Show file tree

Hide file tree

Showing 3 changed files with 71 additions and 12 deletions.
diff --git a/model/dyn_core.F90 b/model/dyn_core.F90
@@ -297,16 +297,16 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
     iep1 = ie + 1
     jep1 = je + 1
 
+!$OMP parallel do default(none) shared(npz,dp_ref,ak,bk)
+       do k=1,npz
+          dp_ref(k) = (ak(k+1)-ak(k)) + (bk(k+1)-bk(k))*1.E5
+       enddo
+
     if ( .not.hydrostatic ) then
 
          rgrav = 1.0/grav
            k1k =  akap / (1.-akap)    ! rg/Cv=0.4
 
-!$OMP parallel do default(none) shared(npz,dp_ref,ak,bk)
-       do k=1,npz
-          dp_ref(k) = (ak(k+1)-ak(k)) + (bk(k+1)-bk(k))*1.E5  
-       enddo
-
 !$OMP parallel do default(none) shared(isd,ied,jsd,jed,zs,phis,rgrav)
        do j=jsd,jed
           do i=isd,ied

diff --git a/model/fv_control.F90 b/model/fv_control.F90
@@ -767,13 +767,13 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split)
       ! Define n_split if not in namelist
       if (ntiles==6) then
 #ifdef MAPL_MODE
-         dimx = ceiling(stretch_fac)*4.0*(npx-1)
-         if (.not. hydrostatic ) then
+         dimx = stretch_fac*4.0*(npx-1)
+         if (.not. hydrostatic) then
                                 ns0 = 6
             if ( dimx >=  360 ) ns0 = 7
             if ( dimx >= 1440 ) ns0 = 8
-            if ( stretch_fac > 1.0 ) ns0=ns0+1
-        endif
+         endif
+         if (stretch_fac > 1.0) ns0 = 7
 #else
          dimx = 4.0*(npx-1)
          if ( hydrostatic ) then

diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90
@@ -108,6 +108,10 @@ module fv_mapz_mod
   real, parameter:: cp_vap = cp_vapor        !< 1846.
   real, parameter:: tice = 273.16
 
+  logical, parameter :: w_limiter = .true.
+  real, parameter :: w_max = 90.
+  real, parameter :: w_min = -60.
+
   real(kind=8) :: E_Flux = 0.
   private
 
@@ -207,7 +211,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 !-----------------------------------------------------------------------
   real(kind=8), dimension(is:ie,js:je):: te_2d, zsum0, zsum1
   real, dimension(is:ie,js:je):: dpln
-  real, dimension(is:ie,km)  :: q2, dp2
+  real, dimension(is:ie,km)  :: q2, dp2, w2
   real, dimension(is:ie,km+1):: pe1, pe2, pk1, pk2, pn2, phis
   real, dimension(is:ie+1,km+1):: pe0, pe3
   real, dimension(is:ie):: gz, cvm
@@ -332,7 +336,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
 !$OMP                                  ak,bk,nq,isd,ied,jsd,jed,kord_tr,fill, adiabatic, &
 !$OMP                                  hs,w,ws,kord_wz,rrg,kord_mt,consv,remap_option,gmao_remap)    &
 !$OMP                          private(gz,cvm,bkh,dp2,   &
-!$OMP                                  pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,dpln,dlnp)
+!$OMP                                  pe0,pe1,pe2,pe3,pk1,pk2,pn2,phis,q2,w2,dpln,dlnp)
   do 1000 j=js,je+1
 
      do k=1,km+1
@@ -572,7 +576,62 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk,   &
               delz(i,j,k) = -delz(i,j,k)*dp2(i,k)
            enddo
         enddo
-   endif
+
+         !Fix excessive w - momentum conserving --- sjl
+         ! gz(:) used here as a temporary array
+         if ( w_limiter ) then
+            do k=1,km
+               do i=is,ie
+                  w2(i,k) = w(i,j,k)
+               enddo
+            enddo
+            do k=1, km-1
+               do i=is,ie
+                  if ( w2(i,k) > w_max ) then
+                     gz(i) = (w2(i,k)-w_max) * dp2(i,k)
+                     w2(i,k  ) = w_max
+                     w2(i,k+1) = w2(i,k+1) + gz(i)/dp2(i,k+1)
+                     !print*, ' W_LIMITER down: ', i,j,k, w2(i,k:k+1), w(i,j,k:k+1)
+                  elseif ( w2(i,k) < w_min ) then
+                     gz(i) = (w2(i,k)-w_min) * dp2(i,k)
+                     w2(i,k  ) = w_min
+                     w2(i,k+1) = w2(i,k+1) + gz(i)/dp2(i,k+1)
+                     !print*, ' W_LIMITER down: ', i,j,k, w2(i,k:k+1), w(i,j,k:k+1)
+                  endif
+               enddo
+            enddo
+            do k=km, 2, -1
+               do i=is,ie
+                  if ( w2(i,k) > w_max ) then
+                     gz(i) = (w2(i,k)-w_max) * dp2(i,k)
+                     w2(i,k  ) = w_max
+                     w2(i,k-1) = w2(i,k-1) + gz(i)/dp2(i,k-1)
+                     !print*, ' W_LIMITER up: ', i,j,k, w2(i,k-1:k), w(i,j,k-1:k)
+                  elseif ( w2(i,k) < w_min ) then
+                     gz(i) = (w2(i,k)-w_min) * dp2(i,k)
+                     w2(i,k  ) = w_min
+                     w2(i,k-1) = w2(i,k-1) + gz(i)/dp2(i,k-1)
+                     !print*, ' W_LIMITER up: ', i,j,k, w2(i,k-1:k), w(i,j,k-1:k)
+                  endif
+               enddo
+            enddo
+            do i=is,ie
+               if (w2(i,1) > w_max*2. ) then
+                  w2(i,1) = w_max*2 ! sink out of the top of the domain
+                  !print*, ' W_LIMITER top limited: ', i,j,1, w2(i,1), w(i,j,1)
+               elseif (w2(i,1) < w_min*2. ) then
+                  w2(i,1) = w_min*2.
+                  !print*, ' W_LIMITER top limited: ', i,j,1, w2(i,1), w(i,j,1)
+               endif
+            enddo
+            do k=1,km
+               do i=is,ie
+                  w(i,j,k) = w2(i,k)
+               enddo
+            enddo
+         endif
+
+    endif
 
   endif !(j < je+1)