Merge pull request NOAA-EMC#76 from bensonr/emc_master_sync

changes dz_min from defined constant to a namelist option

model/dyn_core.F90

@@ -619,7 +619,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
                                             call timing_on('UPDATE_DZ_C')
          call update_dz_c(is, ie, js, je, npz, ng, dt2, dp_ref, zs, gridstruct%area, ut, vt, gz, ws3, &
              npx, npy, gridstruct%sw_corner, gridstruct%se_corner, &
-             gridstruct%ne_corner, gridstruct%nw_corner, bd, gridstruct%grid_type)
+             gridstruct%ne_corner, gridstruct%nw_corner, bd, gridstruct%grid_type, flagstruct%dz_min)
                                             call timing_off('UPDATE_DZ_C')
 
                                                call timing_on('Riem_Solver')
@@ -1023,7 +1023,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
 #ifndef SW_DYNAMICS
                                             call timing_on('UPDATE_DZ')
         call update_dz_d(nord_v, damp_vt, flagstruct%hord_tm, is, ie, js, je, npz, ng, npx, npy, gridstruct%area,  &
-                         gridstruct%rarea, dp_ref, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, flagstruct%lim_fac)
+                         gridstruct%rarea, dp_ref, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, flagstruct%lim_fac, &
+                         flagstruct%dz_min)
                                             call timing_off('UPDATE_DZ')
     if ( flagstruct%fv_debug ) then
          if ( .not. flagstruct%hydrostatic )    then

model/fv_arrays.F90

@@ -863,6 +863,9 @@ module fv_arrays_mod
                                            !< at the center of the domain (the center of tile 1), if set to .true.
                                            !< The default value is .false.
 
+   real :: dz_min = 2        !< Minimum thickness depth to  to enforce monotonicity of height to prevent blowup.
+                             !< 2 by default
+
    integer :: a2b_ord = 4   !< Order of interpolation used by the pressure gradient force
                             !< to interpolate cell-centered (A-grid) values to the grid corners.
                             !< The default value is 4 (recommended), which uses fourth-order

model/fv_control.F90

@@ -372,6 +372,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split)
      logical , pointer :: nudge_qv
      real,     pointer :: add_noise
      logical , pointer :: butterfly_effect
+     real,     pointer :: dz_min
 
      integer , pointer :: a2b_ord
      integer , pointer :: c2l_ord
@@ -934,6 +935,7 @@ subroutine set_namelist_pointers(Atm)
        nudge_qv                      => Atm%flagstruct%nudge_qv
        add_noise                     => Atm%flagstruct%add_noise
        butterfly_effect              => Atm%flagstruct%butterfly_effect
+       dz_min                        => Atm%flagstruct%dz_min
        a2b_ord                       => Atm%flagstruct%a2b_ord
        c2l_ord                       => Atm%flagstruct%c2l_ord
        ndims                         => Atm%flagstruct%ndims
@@ -1443,7 +1445,7 @@ subroutine read_namelist_fv_core_nml(Atm)
             c2l_ord, dx_const, dy_const, umax, deglat,      &
             deglon_start, deglon_stop, deglat_start, deglat_stop, &
             phys_hydrostatic, use_hydro_pressure, make_hybrid_z, old_divg_damp, add_noise, butterfly_effect, &
-            nested, twowaynest, nudge_qv, &
+            dz_min, nested, twowaynest, nudge_qv, &
             nestbctype, nestupdate, nsponge, s_weight, &
             check_negative, nudge_ic, halo_update_type, gfs_phil, agrid_vel_rst,     &
             do_uni_zfull, adj_mass_vmr, fac_n_spl, fhouri, update_blend, regional, bc_update_interval,  &

model/nh_utils.F90

@@ -65,18 +65,17 @@ module nh_utils_mod
    public sim3p0_solver, rim_2d
    public Riem_Solver_c
 
-   real, parameter:: dz_min = 6.
    real, parameter:: r3 = 1./3.
 
 CONTAINS
 
   subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, &
-       npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type)
+       npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type, dz_min)
 ! !INPUT PARAMETERS:
   type(fv_grid_bounds_type), intent(IN) :: bd
   integer, intent(in):: is, ie, js, je, ng, km, npx, npy, grid_type
   logical, intent(IN):: sw_corner, se_corner, ne_corner, nw_corner
-  real, intent(in):: dt
+  real, intent(in):: dt, dz_min
   real, intent(in):: dp0(km)
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: ut, vt
   real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng):: area
@@ -195,28 +194,28 @@ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws
 6000 continue
 
 ! Enforce monotonicity of height to prevent blowup
-!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km)
+!$OMP parallel do default(none) shared(is1,ie1,js1,je1,ws,zs,gz,rdt,km,dz_min)
   do j=js1, je1
-     do k=2, km+1
-        do i=is1, ie1
-           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
+     do k=km, 1, -1
+        do i=is1, ie1
+           gz(i,j,k) = max( gz(i,j,k), gz(i,j,k+1) + dz_min )
+        enddo
+     enddo
   enddo
 
   end subroutine update_dz_c
 
 
   subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, rarea,   &
-                         dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac)
+                         dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac, dz_min)
 
   type(fv_grid_bounds_type), intent(IN) :: bd
   integer, intent(in):: is, ie, js, je, ng, km, npx, npy
   integer, intent(in):: hord
-  real, intent(in)   :: rdt
+  real, intent(in)   :: rdt, dz_min
   real, intent(in)   :: dp0(km)
   real, intent(in)   :: area(is-ng:ie+ng,js-ng:je+ng)
   real, intent(in)   :: rarea(is-ng:ie+ng,js-ng:je+ng)
@@ -309,17 +308,17 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area,
 
   enddo
 
-!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt)
+!$OMP parallel do default(none) shared(is,ie,js,je,km,ws,zs,zh,rdt,dz_min)
   do j=js, je
-     do k=2, km+1
+     do i=is,ie
+        ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt
+     enddo
+     do k=km, 1, -1
         do i=is, ie
 ! Enforce monotonicity of height to prevent blowup
-           zh(i,j,k) = min( zh(i,j,k), zh(i,j,k-1) - dz_min )
+           zh(i,j,k) = max( 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