From 0ac3a2137bba5786f40e1beb80fc18c380c5be0e Mon Sep 17 00:00:00 2001 From: Lauren Chilutti Date: Wed, 19 Apr 2023 09:09:49 -0400 Subject: [PATCH 1/5] Bringing over internal updates to nh_core and nh_utils which includes a revision to the semi-implicit solver to linearize vertical sound wave propagation about the hydrostatic state. --- model/dyn_core.F90 | 10 +++--- model/fv_arrays.F90 | 3 ++ model/fv_control.F90 | 4 ++- model/nh_core.F90 | 10 +++--- model/nh_utils.F90 | 75 +++++++++++++++++++++++++++++++++----------- 5 files changed, 74 insertions(+), 28 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index 7f7bb4688..a31591f68 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -525,8 +525,9 @@ 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 ) - call timing_off('Riem_Solver') + flagstruct%a_imp, flagstruct%scale_z, pfull, & + flagstruct%tau_w, flagstruct%rf_cutoff ) + call timing_off('Riem_Solver') if (gridstruct%nested) then call nh_bc(ptop, grav, akap, cp, delpc, neststruct%delz_BC, ptc, phis, & @@ -922,8 +923,9 @@ 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) - call timing_off('Riem_Solver') + flagstruct%use_logp, remap_step, beta<-0.1, & + flagstruct%tau_w) + call timing_off('Riem_Solver') call timing_on('COMM_TOTAL') if ( gridstruct%square_domain ) then diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index ed56e0243..868feeebf 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -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 :: tau_w = 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 diff --git a/model/fv_control.F90 b/model/fv_control.F90 index a8842c64e..7d7648cc6 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -232,6 +232,7 @@ 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 :: rf_cutoff real , pointer :: te_err real , pointer :: tw_err @@ -763,6 +764,7 @@ 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 rf_cutoff => Atm%flagstruct%rf_cutoff te_err => Atm%flagstruct%te_err tw_err => Atm%flagstruct%tw_err @@ -930,7 +932,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, tau_w, 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, & diff --git a/model/nh_core.F90 b/model/nh_core.F90 index 00055a928..d1d1771d3 100644 --- a/model/nh_core.F90 +++ b/model/nh_core.F90 @@ -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 @@ -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) @@ -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 @@ -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 diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index adea41188..cfc5aa189 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -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 @@ -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, & @@ -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 @@ -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.5*pi*log(rf_cutoff/pfull(k))/log(rf_cutoff/ptop))**2 + rff(k) = 1.0d0 / ( 1.0d0+rff(k) ) + enddo + endif + + !$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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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 From f3c456fcf95ad42d8e53af66a7da6a2fc3b59570 Mon Sep 17 00:00:00 2001 From: Lauren Chilutti Date: Wed, 19 Apr 2023 14:09:39 -0400 Subject: [PATCH 2/5] resolve inconsistent use of single and double precision in rayleigh damping. Address spacing change. --- model/dyn_core.F90 | 4 ++-- model/nh_utils.F90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index a31591f68..a9c466572 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -527,7 +527,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, q_con, delpc, gz, pkc, ws3, flagstruct%p_fac, & flagstruct%a_imp, flagstruct%scale_z, pfull, & flagstruct%tau_w, flagstruct%rf_cutoff ) - call timing_off('Riem_Solver') + call timing_off('Riem_Solver') if (gridstruct%nested) then call nh_bc(ptop, grav, akap, cp, delpc, neststruct%delz_BC, ptc, phis, & @@ -925,7 +925,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, flagstruct%scale_z, flagstruct%p_fac, flagstruct%a_imp, & flagstruct%use_logp, remap_step, beta<-0.1, & flagstruct%tau_w) - call timing_off('Riem_Solver') + call timing_off('Riem_Solver') call timing_on('COMM_TOTAL') if ( gridstruct%square_domain ) then diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index cfc5aa189..dc8e7beca 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -352,7 +352,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & do k=1,km if (pfull(k) > rf_cutoff) exit k_rf = k - rff(k) = dt/tau_w * sin(0.5*pi*log(rf_cutoff/pfull(k))/log(rf_cutoff/ptop))**2 + 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 From ec8bc8cf2ad2fb02c313141e4a543b3ef3d272fa Mon Sep 17 00:00:00 2001 From: Lauren Chilutti Date: Wed, 19 Apr 2023 12:42:28 -0400 Subject: [PATCH 3/5] rename tau_w --- model/dyn_core.F90 | 4 ++-- model/fv_arrays.F90 | 6 +++--- model/fv_control.F90 | 4 +++- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/model/dyn_core.F90 b/model/dyn_core.F90 index a9c466572..8ae8bbb70 100644 --- a/model/dyn_core.F90 +++ b/model/dyn_core.F90 @@ -526,7 +526,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, akap, cappa, cp, ptop, phis, omga, ptc, & q_con, delpc, gz, pkc, ws3, flagstruct%p_fac, & flagstruct%a_imp, flagstruct%scale_z, pfull, & - flagstruct%tau_w, flagstruct%rf_cutoff ) + flagstruct%fast_tau_w_sec, flagstruct%rf_cutoff ) call timing_off('Riem_Solver') if (gridstruct%nested) then @@ -924,7 +924,7 @@ subroutine dyn_core(npx, npy, npz, ng, sphum, nq, bdt, n_map, n_split, zvir, cp, pe, pkc, pk3, pk, peln, ws, & flagstruct%scale_z, flagstruct%p_fac, flagstruct%a_imp, & flagstruct%use_logp, remap_step, beta<-0.1, & - flagstruct%tau_w) + flagstruct%fast_tau_w_sec) call timing_off('Riem_Solver') call timing_on('COMM_TOTAL') diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90 index 868feeebf..856efa183 100644 --- a/model/fv_arrays.F90 +++ b/model/fv_arrays.F90 @@ -692,9 +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 :: tau_w = 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. + 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 diff --git a/model/fv_control.F90 b/model/fv_control.F90 index 7d7648cc6..cce56fd1e 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -233,6 +233,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 real , pointer :: te_err real , pointer :: tw_err @@ -765,6 +766,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 te_err => Atm%flagstruct%te_err tw_err => Atm%flagstruct%tw_err @@ -932,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_w, 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, & From c6d02f65405de6c327c943a74cbe237e1c568551 Mon Sep 17 00:00:00 2001 From: Lauren Chilutti Date: Wed, 19 Apr 2023 15:08:43 -0400 Subject: [PATCH 4/5] Modify the rff calculation to be fully double precision. Fixed a mistake in the previous cherry-pick commit --- model/fv_control.F90 | 2 -- model/nh_utils.F90 | 7 ++++--- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/model/fv_control.F90 b/model/fv_control.F90 index cce56fd1e..f9d171dd9 100644 --- a/model/fv_control.F90 +++ b/model/fv_control.F90 @@ -232,7 +232,6 @@ 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 @@ -765,7 +764,6 @@ 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 diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index dc8e7beca..9171b951f 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -24,7 +24,7 @@ module nh_utils_mod ! To do list: ! include moisture effect in pt !------------------------------ - use constants_mod, only: rdgas, cp_air, grav, pi + use constants_mod, only: rdgas, cp_air, grav, 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 @@ -336,6 +336,7 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & real, dimension(is-1:ie+1,km ):: dm, dz2, w2, pm2, gm2, cp2 real, dimension(is-1:ie+1,km+1):: pem, pe2, peg real gama, rgrav + real(kind=8) :: rff_temp integer i, j, k integer is1, ie1 @@ -352,8 +353,8 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & 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) ) + rff_temp = real(dt/tau_w,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 From 6b4b89afa59879b3c7b88e3d7cfa6ef47249d122 Mon Sep 17 00:00:00 2001 From: Lauren Chilutti Date: Thu, 20 Apr 2023 08:32:43 -0400 Subject: [PATCH 5/5] fully implement the name change of tau_w to fast_tau_w_sec --- model/nh_core.F90 | 10 +++++----- model/nh_utils.F90 | 25 +++++++++++++------------ 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/model/nh_core.F90 b/model/nh_core.F90 index d1d1771d3..df7802cdc 100644 --- a/model/nh_core.F90 +++ b/model/nh_core.F90 @@ -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, tau_w) + use_logp, last_call, fp_out, fast_tau_w_sec) !-------------------------------------------- ! !OUTPUT PARAMETERS ! Ouput: gz: grav*height at edges @@ -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, tau_w + 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) @@ -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,tau_w ) & +!$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) do 2000 j=js, je @@ -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, tau_w) + 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, pe2, dm, & pm2, pem, w2, dz2, pt(is:ie,j,1:km), ws(is,j), & - a_imp, p_fac, scale_m, tau_w) + a_imp, p_fac, scale_m, fast_tau_w_sec) endif do k=1, km diff --git a/model/nh_utils.F90 b/model/nh_utils.F90 index 9171b951f..6cd65fede 100644 --- a/model/nh_utils.F90 +++ b/model/nh_utils.F90 @@ -318,11 +318,11 @@ 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, & - pfull, tau_w, rf_cutoff) + 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, tau_w, rf_cutoff + 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 @@ -347,20 +347,21 @@ subroutine Riem_Solver_c(ms, dt, is, ie, js, je, km, ng, & ie1 = ie + 1 !Set up rayleigh damping - if (tau_w > 1.e-5 .and. .not. RFw_initialized) then + 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/tau_w,kind=8) * sin(0.5d0*pi_8*log(real(rf_cutoff/pfull(k),kind=8))/log(real(rf_cutoff/ptop, kind=8)))**2 + 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, & -!$OMP a_imp,dt,gama,akap,ws,p_fac,scale_m,ms,hs,q_con,cappa,tau_w) & +!$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) do 2000 j=js-1, je+1 @@ -416,7 +417,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, tau_w) + 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 @@ -1212,9 +1213,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, tau_w) + 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, tau_w + 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 @@ -1314,7 +1315,7 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2, enddo !!! Try Rayleigh damping of w - if (tau_w > 1.e-5) then + if (fast_tau_w_sec > 1.e-5) then !currently not damping to heat do k=1,k_rf do i=is,ie @@ -1355,9 +1356,9 @@ subroutine SIM1_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe, dm2, 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, tau_w) + 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, tau_w + 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 @@ -1473,7 +1474,7 @@ subroutine SIM_solver(dt, is, ie, km, rgas, gama, gm2, cp2, kappa, pe2, dm2, enddo !!! Try Rayleigh damping of w - if (tau_w > 1.e-5) then + if (fast_tau_w_sec > 1.e-5) then !currently not damping to heat do k=1,k_rf do i=is,ie