diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 72d39c096..95653c83d 100644 --- a/model/dyn_core.F90 +++ b/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 diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 1319a7c1d..af210d9f3 100644 --- a/model/fv_arrays.F90 +++ b/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 diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 35061db41..8f5fe42fe 100644 --- a/model/fv_control.F90 +++ b/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, & diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 543c4c3f5..a200c2932 100644 --- a/model/nh_utils.F90 +++ b/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