Skip to content

Release FV3 solver updates #261

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 5 commits into from
Apr 20, 2023
Merged

Release FV3 solver updates #261

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0ac3a21
Bringing over internal updates to nh_core and nh_utils which includes…
laurenchilutti Apr 19, 2023
f3c456f
resolve inconsistent use of single and double precision in rayleigh d…
laurenchilutti Apr 19, 2023
ec8bc8c
rename tau_w
laurenchilutti Apr 19, 2023
c6d02f6
Modify the rff calculation to be fully double precision. Fixed a mis…
laurenchilutti Apr 19, 2023
6b4b89a
fully implement the name change of tau_w to fast_tau_w_sec
laurenchilutti Apr 20, 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 @@ -525,7 +525,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
call Riem_Solver_C( ms, dt2, is, ie, js, je, npz, ng, &
akap, cappa, cp, 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 @@ -922,7 +923,8 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp,
akap, cappa, 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 @@ -692,6 +692,9 @@ module fv_arrays_mod
real :: rf_cutoff = 30.E2 !< Pressure below which no Rayleigh damping is applied if tau > 0.
real :: te_err = 1.e-5 !< 64bit: 1.e-14, 32bit: 1.e-7; turn off to save computer time
real :: tw_err = 1.e-8 !< 64bit: 1.e-14, 32bit: 1.e-7; turn off to save computer time
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
6 changes: 5 additions & 1 deletion model/fv_control.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -232,6 +232,8 @@ subroutine fv_control_init(Atm, dt_atmos, this_grid, grids_on_this_pe, p_split)
real , pointer :: ke_bg
real , pointer :: consv_te
real , pointer :: tau
real , pointer :: tau_w
real , pointer :: fast_tau_w_sec
real , pointer :: rf_cutoff
real , pointer :: te_err
real , pointer :: tw_err
Expand Down Expand Up @@ -763,6 +765,8 @@ subroutine set_namelist_pointers(Atm)
ke_bg => Atm%flagstruct%ke_bg
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
te_err => Atm%flagstruct%te_err
tw_err => Atm%flagstruct%tw_err
Expand Down Expand Up @@ -930,7 +934,7 @@ subroutine read_namelist_fv_core_nml(Atm)
dry_mass, grid_type, do_Held_Suarez, &
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_h2o, rf_cutoff, te_err, tw_err, nf_omega, hydrostatic, fv_sg_adj, sg_cutoff, breed_vortex_inline, &
tau, fast_tau_w_sec, tau_h2o, rf_cutoff, te_err, tw_err, 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
10 changes: 5 additions & 5 deletions model/nh_core.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,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, tau_w)
!--------------------------------------------
! !OUTPUT PARAMETERS
! Ouput: gz: grav*height at edges
Expand All @@ -54,7 +54,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, tau_w
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 All @@ -81,7 +81,7 @@ 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, &
!$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,tau_w ) &
!$OMP private(cp2, gm2, dm, dz2, pm2, pem, peg, pelng, pe2, peln2, w2)
do 2000 j=js, je

Expand Down Expand Up @@ -146,11 +146,11 @@ subroutine Riem_Solver3(ms, dt, is, ie, js, je, km, ng, &
dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.)
elseif ( a_imp > 0.999 ) then
call SIM1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, 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, tau_w)
else
call SIM_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, 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, tau_w)
endif

do k=1, km
Expand Down
75 changes: 57 additions & 18 deletions model/nh_utils.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,11 @@ module nh_utils_mod
! To do list:
! include moisture effect in pt
!------------------------------
use constants_mod, only: rdgas, cp_air, grav
use constants_mod, only: rdgas, cp_air, grav, pi
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

implicit none
private
Expand All @@ -44,6 +45,10 @@ module nh_utils_mod
#endif
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 @@ -312,16 +317,18 @@ end subroutine update_dz_d

subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
akap, cappa, cp, 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, tau_w, 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, tau_w, 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
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 @@ -338,8 +345,21 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
is1 = is - 1
ie1 = ie + 1

!Set up rayleigh damping
if (tau_w > 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(k) = dt/tau_w * sin(0.5d0*pi*log(rf_cutoff/pfull(k))/log(rf_cutoff/ptop))**2d0
rff(k) = 1.0d0 / ( 1.0d0+rff(k) )
enddo
endif


Comment on lines +348 to +362

@bensonr bensonr Apr 19, 2023

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In some places here we are forcing constants to double precision, others are single, and one has no designation. I think we need to ensure consistency.

@laurenchilutti laurenchilutti Apr 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have a recommendation of either all double precision or all single?

!$OMP parallel do default(none) shared(js,je,is1,ie1,km,delp,pef,ptop,gz,rgrav,w3,pt, &
!$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,tau_w) &
!$OMP private(cp2,gm2, dm, dz2, w2, pm2, pe2, pem, peg)
do 2000 j=js-1, je+1

Expand Down Expand Up @@ -395,7 +415,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, &
dm, pm2, w2, dz2, pt(is1:ie1,j,1:km), ws(is1,j), .true.)
else
call SIM1_solver(dt, is1, ie1, km, rdgas, gama, gm2, cp2, akap, 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, tau_w)
endif

do k=2,km+1
Expand Down Expand Up @@ -529,11 +549,11 @@ subroutine Riem_Solver3test(ms, dt, is, ie, js, je, km, ng, &
dm, pm2, w2, dz2, pt(is:ie,j,1:km), ws(is,j), .false.)
elseif ( a_imp > 0.999 ) then
call SIM1_solver(dt, is, ie, km, rdgas, gama, gm2, cp2, akap, 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, 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 @@ -1191,9 +1211,9 @@ end subroutine SIM3p0_solver


subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2, &
pm2, pem, w2, dz2, pt2, ws, p_fac)
pm2, pem, w2, dz2, pt2, ws, p_fac, tau_w)
integer, intent(in):: is, ie, km
real, intent(in):: dt, rgas, gama, kappa, p_fac
real, intent(in):: dt, rgas, gama, kappa, p_fac, tau_w
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 @@ -1259,9 +1279,9 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2,
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
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
enddo
enddo
Expand All @@ -1278,9 +1298,9 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2,
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
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
gam(i,km) = aa(i,km) / bet(i)
bet(i) = dm2(i,km) - (aa(i,km)+p1(i) + aa(i,km)*gam(i,km))
Expand All @@ -1292,6 +1312,16 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2,
enddo
enddo

!!! Try Rayleigh damping of w
if (tau_w > 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 @@ -1321,12 +1351,12 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2,
enddo
enddo

end subroutine SIM1_solver
end subroutine SIM1_solver

subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe2, dm2, &
pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m)
pm2, pem, w2, dz2, pt2, ws, alpha, p_fac, scale_m, tau_w)
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, tau_w
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 @@ -1394,8 +1424,7 @@ subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe2, dm2,

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 @@ -1442,6 +1471,16 @@ subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe2, dm2,
enddo
enddo

!!! Try Rayleigh damping of w
if (tau_w > 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