Merge pull request NOAA-GFDL#67 from GEOS-ESM/feature/sdrabenh/hwt03_…

…changes_into_v11.0.0

Merge hwt_v03 changes into v11.0.0-rc1

model/dyn_core.F90

@@ -557,7 +557,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
  enddo
  endif
  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, &
+ call update_dz_c(is, ie, js, je, npz, ng, dt2, flagstruct%dz_min, 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)
  call timing_off('UPDATE_DZ_C')
@@ -695,16 +695,18 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
  if ( flagstruct%do_vort_damp ) then
 ! damping on delp and vorticity:
  nord_v(k)=0; 
+ damp_vt(k) = 0.5*d2_divg
  endif
  d_con_k = 0.
- elseif ( (k<=MAX(2,flagstruct%n_sponge-1)).and. flagstruct%d2_bg_k2>0.01 ) then
+ elseif ( k<=MAX(2,flagstruct%n_sponge-1) .and. flagstruct%d2_bg_k2>0.01 ) then
  nord_k=0; d2_divg = max(flagstruct%d2_bg, flagstruct%d2_bg_k2)
  nord_w=0; damp_w = d2_divg
  if ( flagstruct%do_vort_damp ) then
  nord_v(k)=0; 
+ damp_vt(k) = 0.5*d2_divg
  endif
  d_con_k = 0.
- elseif ( (k<=MAX(3,flagstruct%n_sponge))  .and. flagstruct%d2_bg_k2>0.05 ) then
+ elseif ( k<=MAX(3,flagstruct%n_sponge) .and. flagstruct%d2_bg_k2>0.05 ) then
  nord_k=0; d2_divg = max(flagstruct%d2_bg, 0.2*flagstruct%d2_bg_k2)
  nord_w=0; damp_w = d2_divg
  d_con_k = 0.
@@ -843,7 +845,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
 #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, delz, ws, rdt, gridstruct, bd, flagstruct%lim_fac)
+ gridstruct%rarea, dp_ref, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, flagstruct%dz_min, gridstruct, bd, flagstruct%lim_fac)
  call timing_off('UPDATE_DZ')
 
  if (idiag%id_ws>0 .and. last_step) then
@@ -1165,12 +1167,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_split, zvir, cp, akap,
  do k=1,n_con
  delt = abs(bdt*flagstruct%delt_max)
 ! Sponge layers:
- if ( flagstruct%n_sponge == 0) then 
- if ( k == 1 ) delt = 0.1*delt
- if ( k == 2 ) delt = 0.5*delt
- else
- delt = delt*MIN(1.0,FLOAT(k)/FLOAT(flagstruct%n_sponge))
- endif
+ if ( k == 1 ) delt = 0.1*delt
+ if ( k == 2 ) delt = 0.5*delt
  do j=js,je
  do i=is,ie
 #ifdef MOIST_CAPPA
@@ -1801,7 +1799,7 @@ subroutine one_grad_p(u, v, pk, gz, divg2, delp, dt, ng, gridstruct, bd, npx, np
  !$OMP parallel do default(none) shared(is,ie,js,je,wk2,divg2)
  do j=js,je+1
  do i=is,ie
- wk2(i,j) = divg2(i,j)-divg2(i+1,j)
+ wk2(i,j) = (divg2(i,j)-divg2(i+1,j))
  enddo
  enddo
 

model/fv_arrays.F90

@@ -614,6 +614,8 @@ module fv_arrays_mod
  !< and 1. The default value is 0.75. Only used if 
  !< 'hydrostatic' = .false.
 
+ real :: dz_min = 2.0 !< Controls minimum thickness in NH solver
+
  integer :: n_split = 0 !< The number of small dynamics (acoustic) time steps between
  !< vertical remapping. 0 by default, in which case the model 
  !< produces a good first guess by examining the resolution, 

model/fv_control.F90

@@ -222,6 +222,7 @@ module fv_control_mod
  logical , pointer :: reset_eta 
  real , pointer :: p_fac
  real , pointer :: a_imp
+ real , pointer :: dz_min
  integer , pointer :: n_split 
  ! Default 
  integer , pointer :: m_split 
@@ -656,7 +657,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split)
 
  namelist /fv_grid_nml/ grid_name, grid_file
  namelist /fv_core_nml/npx, npy, ntiles, npz, npz_rst, layout, io_layout, ncnst, nwat, &
- use_logp, p_fac, a_imp, k_split, n_split, m_split, q_split, print_freq, write_3d_diags, do_schmidt, &
+ use_logp, p_fac, a_imp, dz_min, k_split, n_split, m_split, q_split, print_freq, write_3d_diags, do_schmidt, &
  hord_mt, hord_vt, hord_tm, hord_dp, hord_tr, shift_fac, stretch_fac, target_lat, target_lon, &
  kord_mt, kord_wz, kord_tm, kord_tr, fv_debug, fv_land, nudge, do_sat_adj, do_f3d, &
  external_ic, read_increment, ncep_ic, nggps_ic, ecmwf_ic, use_new_ncep, use_ncep_phy, fv_diag_ic, &
@@ -916,6 +917,7 @@ subroutine run_setup(Atm, dt_atmos, grids_on_this_pe, p_split)
  if(is_master()) then
  write(*,*) 'Off center implicit scheme param=', a_imp
  write(*,*) ' p_fac=', p_fac
+ write(*,*) ' dz_min=', dz_min
  endif
  endif
 
@@ -1220,6 +1222,7 @@ subroutine setup_pointers(Atm)
  reset_eta => Atm%flagstruct%reset_eta
  p_fac => Atm%flagstruct%p_fac
  a_imp => Atm%flagstruct%a_imp
+ dz_min => Atm%flagstruct%dz_min
  n_split => Atm%flagstruct%n_split
  m_split => Atm%flagstruct%m_split
  k_split => Atm%flagstruct%k_split

model/fv_dynamics.F90

@@ -922,7 +922,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
  call range_check('VA_dyn', ua, is, ie, js, je, ng, npz, gridstruct%agrid, &
  -280., 280., bad_range)
  call range_check('TA_dyn', pt, is, ie, js, je, ng, npz, gridstruct%agrid, &
- 130., 335., bad_range)
+ 150., 333., bad_range)
  if ( .not. hydrostatic ) then
  call range_check('W_dyn', w, is, ie, js, je, ng, npz, gridstruct%agrid, &
  -100., 100., bad_range)

model/nh_utils.F90

@@ -63,18 +63,17 @@ module nh_utils_mod
  public sim3p0_solver, rim_2d
  public Riem_Solver_c
 
- real, parameter:: dz_min = 2.
  real, parameter:: r3 = 1./3.
 
 CONTAINS 
 
- subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, &
+ subroutine update_dz_c(is, ie, js, je, km, ng, dt, dz_min, dp0, zs, area, ut, vt, gz, ws, &
  npx, npy, sw_corner, se_corner, ne_corner, nw_corner, bd, grid_type)
 ! !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
@@ -193,7 +192,7 @@ 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,dz_min,km)
  do j=js1, je1
  do i=is1, ie1
  ws(i,j) = ( zs(i,j) - gz(i,j,km+1) ) * rdt
@@ -209,12 +208,12 @@ 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, delz, ws, rdt, gridstruct, bd, lim_fac)
+ dp0, zs, zh, crx, cry, xfx, yfx, delz, ws, rdt, dz_min, gridstruct, bd, lim_fac)
 
  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)
@@ -308,7 +307,7 @@ 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 i=is,ie
  ws(i,j) = ( zs(i,j) - zh(i,j,km+1) ) * rdt