Skip to content

dev/emc: Release FV3 Solver Updates (cherry-pick) #266

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
laurenchilutti merged 2 commits into from
Apr 28, 2023
Merged

dev/emc: Release FV3 Solver Updates (cherry-pick) #266

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
10c6c2f
Release FV3 solver updates (#261)
laurenchilutti Apr 20, 2023
5ba14f2
adding fast_tau_w_sec to openmp statement within nh_utils.F90
laurenchilutti Apr 25, 2023
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
6 changes: 4 additions & 2 deletions model/dyn_core.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -632,7 +632,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
#endif
ptop, phis, omga, ptc, &
q_con, delpc, gz, pkc, ws3, flagstruct%p_fac, &
flagstruct%a_imp, flagstruct%scale_z )
flagstruct%a_imp, flagstruct%scale_z, pfull, &
flagstruct%fast_tau_w_sec, flagstruct%rf_cutoff )
call timing_off('Riem_Solver')

if (gridstruct%nested) then
Expand Down Expand Up @@ -1050,7 +1051,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
ptop, zs, q_con, w, delz, pt, delp, zh, &
pe, pkc, pk3, pk, peln, ws, &
flagstruct%scale_z, flagstruct%p_fac, flagstruct%a_imp, &
flagstruct%use_logp, remap_step, beta<-0.1)
flagstruct%use_logp, remap_step, beta<-0.1, &
flagstruct%fast_tau_w_sec)
call timing_off('Riem_Solver')

call timing_on('COMM_TOTAL')
Expand Down
3 changes: 3 additions & 0 deletions model/fv_arrays.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -698,6 +698,9 @@ module fv_arrays_mod
!< considered; and for non-hydrostatic models values of 10 or less should be
!< considered, with smaller values for higher-resolution.
real :: rf_cutoff = 30.E2 !< Pressure below which no Rayleigh damping is applied if tau > 0.
real :: fast_tau_w_sec = 0.0 !< Time scale (seconds) for Rayleigh damping applied to vertical velocity only.
!< Values of 0.2 are very effective at eliminating spurious vertical motion in
!< the stratosphere. Default is 0.0, which disables this.
logical :: filter_phys = .false.
logical :: dwind_2d = .false. !< Whether to use a simpler & faster algorithm for interpolating
!< the A-grid (cell-centered) wind tendencies computed from the physics
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 @@ -331,6 +331,7 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split,
real , pointer :: consv_te
real , pointer :: tau
real , pointer :: tau_w
real , pointer :: fast_tau_w_sec
real , pointer :: rf_cutoff
logical , pointer :: filter_phys
logical , pointer :: dwind_2d
Expand Down Expand Up @@ -887,6 +888,7 @@ subroutine set_namelist_pointers(Atm)
consv_te => Atm%flagstruct%consv_te
tau => Atm%flagstruct%tau
tau_w => Atm%flagstruct%tau_w
fast_tau_w_sec => Atm%flagstruct%fast_tau_w_sec
rf_cutoff => Atm%flagstruct%rf_cutoff
filter_phys => Atm%flagstruct%filter_phys
dwind_2d => Atm%flagstruct%dwind_2d
Expand Down Expand Up @@ -1055,7 +1057,7 @@ subroutine read_namelist_fv_core_nml(Atm)
dry_mass, grid_type, do_Held_Suarez, do_reed_physics, reed_cond_only, &
consv_te, fill, filter_phys, fill_dp, fill_wz, fill_gfs, consv_am, RF_fast, &
range_warn, dwind_2d, inline_q, z_tracer, reproduce_sum, adiabatic, do_vort_damp, no_dycore, &
tau, tau_w, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, &
tau, tau_w, fast_tau_w_sec, tau_h2o, rf_cutoff, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, &
na_init, nudge_dz, hybrid_z, Make_NH, n_zs_filter, nord_zs_filter, full_zs_filter, reset_eta, &
pnats, dnats, dnrts, a2b_ord, remap_t, p_ref, d2_bg_k1, d2_bg_k2, &
c2l_ord, dx_const, dy_const, umax, deglat, &
Expand Down
12 changes: 6 additions & 6 deletions model/nh_core.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, &
ptop, zs, q_con, w, delz, pt, &
delp, zh, pe, ppe, pk3, pk, peln, &
ws, scale_m, p_fac, a_imp, &
use_logp, last_call, fp_out)
use_logp, last_call, fp_out, fast_tau_w_sec)
!--------------------------------------------
! !OUTPUT PARAMETERS
! Ouput: gz: grav*height at edges
Expand All @@ -82,7 +82,7 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, &
integer, intent(in):: ms, is, ie, js, je, km, ng
integer, intent(in):: isd, ied, jsd, jed
real, intent(in):: dt !< the BIG horizontal Lagrangian time step
real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m
real, intent(in):: akap, cp, ptop, p_fac, a_imp, scale_m, fast_tau_w_sec
real, intent(in):: zs(isd:ied,jsd:jed)
logical, intent(in):: last_call, use_logp, fp_out
real, intent(in):: ws(is:ie,js:je)
Expand Down Expand Up @@ -116,10 +116,10 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, &
!$OMP parallel do default(none) shared(is,ie,js,je,km,delp,ptop,peln1,pk3,ptk,akap,rgrav,zh,pt, &
!$OMP w,a_imp,dt,gama,ws,p_fac,scale_m,ms,delz,last_call, &
#ifdef MULTI_GASES
!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad ) &
!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,kapad,fast_tau_w_sec ) &
!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2,kapad2)
#else
!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con ) &
!$OMP peln,pk,fp_out,ppe,use_logp,zs,pe,cappa,q_con,fast_tau_w_sec ) &
!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2)
#endif
do 2000 j=js, je
Expand Down Expand Up @@ -206,15 +206,15 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, &
kapad2, &
#endif
pe2, dm, &
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac)
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, fast_tau_w_sec)
else
call SIM_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
#ifdef MULTI_GASES
kapad2, &
#endif
pe2, dm, &
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), &
a_imp, p_fac, scale_m)
a_imp, p_fac, scale_m, fast_tau_w_sec)
endif


Expand Down
82 changes: 62 additions & 20 deletions model/nh_utils.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,11 @@ module nh_utils_mod
#else
use constants_mod, only: rdgas, cp_air, grav
#endif
use constants_mod, only: pi_8
use tp_core_mod, only: fv_tp_2d
use sw_core_mod, only: fill_4corners, del6_vt_flux
use fv_arrays_mod, only: fv_grid_bounds_type, fv_grid_type,fv_nest_BC_type_3d
use mpp_mod, only: mpp_pe
#ifdef MULTI_GASES
use multi_gases_mod, only: vicpqd, vicvqd
#endif
Expand All @@ -71,6 +73,10 @@ module nh_utils_mod

real, parameter:: r3 = 1./3.

real, allocatable :: rff(:)
logical :: RFw_initialized = .false.
integer :: k_rf = 0

CONTAINS

subroutine update_dz_c(is, ie, js, je, km, ng, dt, dp0, zs, area, ut, vt, gz, ws, &
Expand Down Expand Up @@ -346,11 +352,12 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
kapad, &
#endif
ptop, hs, w3, pt, q_con, &
delp, gz, pef, ws, p_fac, a_imp, scale_m)
delp, gz, pef, ws, p_fac, a_imp, scale_m, &
pfull, fast_tau_w_sec, rf_cutoff)

integer, intent(in):: is, ie, js, je, ng, km
integer, intent(in):: ms
real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m
real, intent(in):: dt, akap, cp, ptop, p_fac, a_imp, scale_m, fast_tau_w_sec, rf_cutoff
real, intent(in):: ws(is-ng:ie+ng,js-ng:je+ng)
real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: pt, delp
real, intent(in), dimension(is-ng:,js-ng:,1:):: q_con, cappa
Expand All @@ -359,6 +366,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
#endif
real, intent(in):: hs(is-ng:ie+ng,js-ng:je+ng)
real, intent(in), dimension(is-ng:ie+ng,js-ng:je+ng,km):: w3
real, intent(in) :: pfull(km)
! OUTPUT PARAMETERS
real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: gz
real, intent( out), dimension(is-ng:ie+ng,js-ng:je+ng,km+1):: pef
Expand All @@ -369,6 +377,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
real, dimension(is-1:ie+1,km ):: kapad2
#endif
real gama, rgrav
real(kind=8) :: rff_temp
integer i, j, k
integer is1, ie1

Expand All @@ -378,12 +387,26 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
is1 = is - 1
ie1 = ie + 1

!Set up rayleigh damping
if (fast_tau_w_sec > 1.e-5 .and. .not. RFw_initialized) then
allocate(rff(km))
RFw_initialized = .true.
do k=1,km
if (pfull(k) > rf_cutoff) exit
k_rf = k
rff_temp = real(dt/fast_tau_w_sec,kind=8) &
* sin(0.5d0*pi_8*log(real(rf_cutoff/pfull(k),kind=8))/log(real(rf_cutoff/ptop, kind=8)))**2
rff(k) = 1.0d0 / ( 1.0d0+rff_temp )
enddo
endif


!$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, &
#ifdef MULTI_GASES
!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad) &
!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,kapad,fast_tau_w_sec) &
!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg, kapad2)
#else
!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa) &
!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,fast_tau_w_sec) &
!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg)
#endif
do 2000 j=js-1, je+1
Expand Down Expand Up @@ -455,7 +478,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
kapad2, &
#endif
pe2, &
dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac)
dm, pm2, pem, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), p_fac, fast_tau_w_sec)
endif

do k=2,km+1
Expand Down Expand Up @@ -622,15 +645,15 @@ subroutine Riem_Solver3test(ms, dt, is, ie, js, je, km, ng, &
kapad2, &
#endif
pe2, dm, &
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac)
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), p_fac, -1.)
else
call SIM_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, &
#ifdef MULTI_GASES
kapad2, &
#endif
pe2, dm, &
pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), &
a_imp, p_fac, scale_m)
a_imp, p_fac, scale_m, -1.)
endif

do k=1, km
Expand Down Expand Up @@ -1385,9 +1408,9 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
kapad2, &
#endif
pe, dm2, &
pm2, pem, w2, dz2, pt2, ws, p_fac)
pm2, pem, w2, dz2, pt2, ws, p_fac, fast_tau_w_sec)
integer, intent(in):: is, ie, km
real, intent(in):: dt, rgas, gama, kappa, p_fac
real, intent(in):: dt, rgas, gama, kappa, p_fac, fast_tau_w_sec
real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
real, intent(in ):: ws(is:ie)
real, intent(in ), dimension(is:ie,km+1):: pem
Expand Down Expand Up @@ -1464,14 +1487,14 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
do k=2, km
do i=is, ie
#ifdef MOIST_CAPPA
aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
aa(i,k) = t1g*0.5*(gm2(i,k-1)+gm2(i,k))/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k))
#else
#ifdef MULTI_GASES
gamax = 1./(1.-kapad2(i,k))
t1gx = gamax * 2.*dt*dt
aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
aa(i,k) = t1gx/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k))
#else
aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k)+pp(i,k))
aa(i,k) = t1g/(dz2(i,k-1)+dz2(i,k)) * (pem(i,k))
#endif
#endif
enddo
Expand All @@ -1489,14 +1512,14 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
enddo
do i=is, ie
#ifdef MOIST_CAPPA
p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
p1(i) = t1g*gm2(i,km)/dz2(i,km)*(pem(i,km+1))
#else
#ifdef MULTI_GASES
gamax = 1./(1.-kapad2(i,km))
t1gx = gamax * 2.*dt*dt
p1(i) = t1gx/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
p1(i) = t1gx/dz2(i,km)*(pem(i,km+1))
#else
p1(i) = t1g/dz2(i,km)*(pem(i,km+1)+pp(i,km+1))
p1(i) = t1g/dz2(i,km)*(pem(i,km+1))
#endif
#endif
gam(i,km) = aa(i,km) / bet(i)
Expand All @@ -1509,6 +1532,16 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
enddo
enddo

!!! Try Rayleigh damping of w
if (fast_tau_w_sec > 1.e-5) then
!currently not damping to heat
do k=1,k_rf
do i=is,ie
w2(i,k) = w2(i,k)*rff(k)
enddo
enddo
endif

do i=is, ie
pe(i,1) = 0.
enddo
Expand Down Expand Up @@ -1549,16 +1582,16 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
enddo
enddo

end subroutine SIM1_solver
end subroutine SIM1_solver

subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
#ifdef MULTI_GASES
kapad2, &
#endif
pe2, dm2, &
pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m, fast_tau_w_sec)
integer, intent(in):: is, ie, km
real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m
real, intent(in):: dt, rgas, gama, kappa, p_fac, alpha, scale_m, fast_tau_w_sec
real, intent(in), dimension(is:ie,km):: dm2, pt2, pm2, gm2, cp2
real, intent(in ):: ws(is:ie)
real, intent(in ), dimension(is:ie,km+1):: pem
Expand Down Expand Up @@ -1637,8 +1670,7 @@ subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &

do k=1, km+1
do i=is, ie
! pe2 is Full p
pe2(i,k) = pem(i,k) + pp(i,k)
pe2(i,k) = pem(i,k)
enddo
enddo

Expand Down Expand Up @@ -1697,6 +1729,16 @@ subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, &
enddo
enddo

!!! Try Rayleigh damping of w
if (fast_tau_w_sec > 1.e-5) then
!currently not damping to heat
do k=1,k_rf
do i=is,ie
w2(i,k) = w2(i,k)*rff(k)
enddo
enddo
endif

do i=is, ie
pe2(i,1) = 0.
enddo
Expand Down