Skip to content

new option psm_bc for edge velocity profiles for height advection #77

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bensonr merged 1 commit into from
Mar 18, 2021
Merged

new option psm_bc for edge velocity profiles for height advection #77

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion model/dyn_core.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions model/fv_arrays.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion model/fv_control.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down
142 changes: 133 additions & 9 deletions model/nh_utils.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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
Expand Down