diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 95653c83d..ee22ef355 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -1024,7 +1024,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, 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, & - flagstruct%dz_min) + flagstruct%dz_min, flagstruct%psm_bc) 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 af210d9f3..82e7023d4 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -866,6 +866,10 @@ module fv_arrays_mod real :: dz_min = 2 !< Minimum thickness depth to to enforce monotonicity of height to prevent blowup. !< 2 by default + integer :: psm_bc = 0 !< Option to use origional BCs (0) or zero-gradient BCs (1) + !< to reconstruct interface u/v with the Parabolic Spline Method + !< for the advection of height. 0 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 8f5fe42fe..4bb7af5a4 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -373,6 +373,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split) real, pointer :: add_noise logical , pointer :: butterfly_effect real, pointer :: dz_min + integer, pointer :: psm_bc integer , pointer :: a2b_ord integer , pointer :: c2l_ord @@ -936,6 +937,7 @@ subroutine set_namelist_pointers(Atm) add_noise => Atm%flagstruct%add_noise butterfly_effect => Atm%flagstruct%butterfly_effect dz_min => Atm%flagstruct%dz_min + psm_bc => Atm%flagstruct%psm_bc a2b_ord => Atm%flagstruct%a2b_ord c2l_ord => Atm%flagstruct%c2l_ord ndims => Atm%flagstruct%ndims @@ -1445,7 +1447,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, & - dz_min, nested, twowaynest, nudge_qv, & + dz_min, psm_bc, 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 a200c2932..e29f3a727 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -210,11 +210,11 @@ 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, dz_min) + dp0, zs, zh, crx, cry, xfx, yfx, ws, rdt, gridstruct, bd, lim_fac, dz_min, psm_bc) type(fv_grid_bounds_type), intent(IN) :: bd integer, intent(in):: is, ie, js, je, ng, km, npx, npy - integer, intent(in):: hord + integer, intent(in):: hord, psm_bc real, intent(in) :: rdt, dz_min real, intent(in) :: dp0(km) real, intent(in) :: area(is-ng:ie+ng,js-ng:je+ng) @@ -251,15 +251,27 @@ subroutine update_dz_d(ndif, damp, hord, is, ie, js, je, km, ng, npx, npy, area, isd = is - ng; ied = ie + ng jsd = js - ng; jed = je + ng + if (psm_bc == 0 ) then !$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, & !$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv) - do j=jsd,jed - call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & - dp0, uniform_grid, 0) - if ( j<=je+1 .and. j>=js ) & - call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & - dp0, uniform_grid, 0) - enddo + do j=jsd,jed + call edge_profile(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & + dp0, uniform_grid, 0) + if ( j<=je+1 .and. j>=js ) & + call edge_profile(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & + dp0, uniform_grid, 0) + enddo + else +!$OMP parallel do default(none) shared(jsd,jed,crx,xfx,crx_adv,xfx_adv,is,ie,isd,ied, & +!$OMP km,dp0,uniform_grid,js,je,cry,yfx,cry_adv,yfx_adv) + do j=jsd,jed + call edge_profile_0grad(crx, xfx, crx_adv, xfx_adv, is, ie+1, jsd, jed, j, km, & + dp0, uniform_grid, 0) + if ( j<=je+1 .and. j>=js ) & + call edge_profile_0grad(cry, yfx, cry_adv, yfx_adv, isd, ied, js, je+1, j, km, & + dp0, uniform_grid, 0) + enddo + endif !$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,km,area,xfx_adv,yfx_adv, & !$OMP damp,zh,crx_adv,cry_adv,npx,npy,hord,gridstruct,bd, & @@ -1889,6 +1901,118 @@ subroutine edge_profile(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_gr end subroutine edge_profile + subroutine edge_profile_0grad(q1, q2, q1e, q2e, i1, i2, j1, j2, j, km, dp0, uniform_grid, limiter) +! Optimized for wind profile reconstruction: +! Added this option by Henry Juang and Xiaqiong Zhou 1/21/2021 + integer, intent(in):: i1, i2, j1, j2 + integer, intent(in):: j, km + integer, intent(in):: limiter + logical, intent(in):: uniform_grid + real, intent(in):: dp0(km) + real, intent(in), dimension(i1:i2,j1:j2,km):: q1, q2 + real, intent(out), dimension(i1:i2,j1:j2,km+1):: q1e, q2e +!----------------------------------------------------------------------- + real, dimension(i1:i2,km+1):: qe1, qe2, gam ! edge values + real gak(km) + real bet, r2o3, r4o3 + real g0, gk, xt1, xt2, a_bot + integer i, k + + if ( uniform_grid ) then +!------------------------------------------------ +! Optimized coding for uniform grid: SJL Apr 2007 +!------------------------------------------------ + r2o3 = 2./3. + r4o3 = 4./3. + do i=i1,i2 + qe1(i,1) = r4o3*q1(i,j,1) + r2o3*q1(i,j,2) + qe2(i,1) = r4o3*q2(i,j,1) + r2o3*q2(i,j,2) + enddo + + gak(1) = 7./3. + do k=2,km + gak(k) = 1. / (4. - gak(k-1)) + do i=i1,i2 + qe1(i,k) = (3.*(q1(i,j,k-1) + q1(i,j,k)) - qe1(i,k-1)) * gak(k) + qe2(i,k) = (3.*(q2(i,j,k-1) + q2(i,j,k)) - qe2(i,k-1)) * gak(k) + enddo + enddo + + bet = 1. / (1.5 - 3.5*gak(km)) + do i=i1,i2 + qe1(i,km+1) = (4.*q1(i,j,km) + q1(i,j,km-1) - 3.5*qe1(i,km)) * bet + qe2(i,km+1) = (4.*q2(i,j,km) + q2(i,j,km-1) - 3.5*qe2(i,km)) * bet + enddo + + do k=km,1,-1 + do i=i1,i2 + qe1(i,k) = qe1(i,k) - gak(k)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gak(k)*qe2(i,k+1) + enddo + enddo + else +! Assuming grid varying in vertical only + g0 = dp0(1) / dp0(2) + bet = 1.5 + 2.*g0 + do i=i1,i2 + qe1(i,2) = 3.*( 0.5*q1(i,j,1) + g0*q1(i,j,2) ) / bet + qe2(i,2) = 3.*( 0.5*q2(i,j,1) + g0*q2(i,j,2) ) / bet + gam(i,1) = g0/bet + enddo + + +! for k=2,km + do k=2,km-1 + gk = dp0(k) / dp0(k+1) + do i=i1,i2 + bet = 2. + 2.*gk - gam(i,k-1) + qe1(i,k+1) = ( 3.*(q1(i,j,k)+gk*q1(i,j,k+1)) - qe1(i,k) ) / bet + qe2(i,k+1) = ( 3.*(q2(i,j,k)+gk*q2(i,j,k+1)) - qe2(i,k) ) / bet + gam(i,k) = gk / bet + enddo + enddo + +!km+1 + do i=i1,i2 + bet = 2.- gam(i,km-1) + qe1(i,km+1) = ( 3.*q1(i,j,km) - qe1(i,km) ) / bet + qe2(i,km+1) = ( 3.*q2(i,j,km) - qe2(i,km) ) / bet + enddo + + do i=i1,i2 + do k=km,2,-1 + qe1(i,k) = qe1(i,k) - gam(i,k-1)*qe1(i,k+1) + qe2(i,k) = qe2(i,k) - gam(i,k-1)*qe2(i,k+1) + enddo + qe1(i,1)=1.5*q1(i,j,1)-0.5*qe1(i,2) + qe2(i,1)=1.5*q2(i,j,1)-0.5*qe2(i,2) + enddo + + endif + +!------------------ +! Apply constraints +!------------------ + if ( limiter/=0 ) then ! limit the top & bottom winds + do i=i1,i2 +! Top + if ( q1(i,j,1)*qe1(i,1) < 0. ) qe1(i,1) = 0. + if ( q2(i,j,1)*qe2(i,1) < 0. ) qe2(i,1) = 0. +! Surface: + if ( q1(i,j,km)*qe1(i,km+1) < 0. ) qe1(i,km+1) = 0. + if ( q2(i,j,km)*qe2(i,km+1) < 0. ) qe2(i,km+1) = 0. + enddo + endif + + do k=1,km+1 + do i=i1,i2 + q1e(i,j,k) = qe1(i,k) + q2e(i,j,k) = qe2(i,k) + enddo + enddo + + end subroutine edge_profile_0grad + !TODO LMH 25may18: do not need delz defined on full compute domain; pass appropriate BCs instead subroutine nh_bc(ptop, grav, kappa, cp, delp, delzBC, pt, phis, & #ifdef MULTI_GASES