update model/nh_utils.F90 with dz_min change (#41)

model/nh_utils.F90

@@ -66,7 +66,7 @@ module nh_utils_mod
    public sim3p0_solver, rim_2d
    public Riem_Solver_c
 
-   real, parameter:: dz_min = 2.
+   real, parameter:: dz_min = 6.
    real, parameter:: r3 = 1./3.
 
 CONTAINS 
@@ -198,14 +198,14 @@ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws
 ! Enforce monotonicity of height to prevent blowup
 !$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km)
   do j=js1, je1
-     do i=is1, ie1
-        ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
-     enddo
-     do k=km, 1, -1
+     do k=2, km+1
         do i=is1, ie1
-           gz(i,j,k) = max( gz(i,j,k), gz(i,j,k+1) + dz_min )
+           gz(i,j,k) = min( gz(i,j,k), gz(i,j,k-1) - dz_min )
         enddo
      enddo
+     do i=is1, ie1
+        ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
+     enddo
   enddo
 
   end subroutine update_dz_c
@@ -314,15 +314,15 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area,
 
 !$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt)
   do j=js, je
-     do i=is,ie
-        ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
-     enddo
-     do k=km, 1, -1
+     do k=2, km+1
         do i=is, ie
 ! Enforce monotonicity of height to prevent blowup
-           zh(i,j,k) = max( zh(i,j,k), zh(i,j,k+1) + dz_min )
+           zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min )
         enddo
      enddo
+     do i=is,ie
+        ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
+     enddo
   enddo
 
   end subroutine update_dz_d