Skip to content
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

Jump to bottom

Support for NSSL 2-Moment Microphysics #33

Open
wants to merge 23 commits into
base: dev/emc
Choose a base branch
from
Open

Support for NSSL 2-Moment Microphysics #33

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
14b111c
Update dev/emc from gsd/develop 2020/06/30 (#30)
climbfuji Jul 2, 2020
0c5684d
Fix for multi_gases to 32 bit compiling (#19)
XiaqiongZhou-NOAA Jul 7, 2020
00a3798
Merge branch 'dev/emc' of https://github.com/NOAA-EMC/GFDL_atmos_cube…
tsupinie Jul 20, 2020
55559f6
Add NSSL 2-Moment Microphysics
tsupinie May 22, 2020
14a08c4
Add files to .gitignore
tsupinie Jul 8, 2020
cdb17af
Butterfly effect (from @DusanJovic-NOAA) (#21)
climbfuji Jun 5, 2020
388c65a
Merge GFDL new dynamic core to EMC fork/branch (#15)
XiaqiongZhou-NOAA Jun 9, 2020
beab6a8
Fixes for GNU compilation issues (#22)
climbfuji Jun 22, 2020
f4de207
model/fv_regional_bc.F90: bugfix, use correct MPI variable type in ex…
climbfuji Jun 24, 2020
e012ce1
model/fv_regional_bc.F90: allocate bufr1 to bufr4 to required size fo…
climbfuji Jun 25, 2020
3110292
model/fv_regional_bc.F90: bugfix, use correct message size
climbfuji Jun 25, 2020
531a0d6
model/fv_regional_bc.F90: add comments and a sanity check for buffer …
climbfuji Jun 25, 2020
cdb30ea
Update dev/emc from gsd/develop 2020/06/30 (#30)
climbfuji Jul 2, 2020
546f8df
Fix for multi_gases to 32 bit compiling (#19)
XiaqiongZhou-NOAA Jul 7, 2020
ed6814d
Add a missing hailwat variable
tsupinie Jul 14, 2020
7c0e7b5
Fix some dumb typos in fv_nesting.F90
tsupinie Jul 20, 2020
2838249
Merge branch 'dev/emc' into nssl-microphysics
tsupinie Jul 20, 2020
0f6710a
Clean up merge issues
tsupinie Sep 22, 2020
7fac711
Merge branch 'dev/emc' of https://github.com/NOAA-EMC/GFDL_atmos_cube…
tsupinie Sep 22, 2020
e57c0cd
Merge branch 'dev/emc' into nssl-microphysics
tsupinie Sep 22, 2020
685a21b
Fix some memory errors in NSSL microphysics
tsupinie Oct 16, 2020
cc34b28
Merge branch 'dev/emc' of https://github.com/NOAA-EMC/GFDL_atmos_cube…
tsupinie Feb 19, 2021
907571a
Merge branch 'dev/emc' into nssl-microphysics
tsupinie Feb 19, 2021
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
7 changes: 7 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
*.mod
*.i90
*.o

depend
driver/fvGFS/*.inc
libfv3core.a
2 changes: 1 addition & 1 deletion model/fv_arrays.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ module fv_arrays_mod

logical :: initialized = .false.
real sphum, liq_wat, ice_wat ! GFDL physics
real rainwat, snowwat, graupel
real rainwat, snowwat, graupel, hailwat

real :: efx(max_step), efx_sum, efx_nest(max_step), efx_sum_nest, mtq(max_step), mtq_sum
integer :: steps
Expand Down
64 changes: 57 additions & 7 deletions model/fv_dynamics.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,10 @@ module fv_dynamics_mod
! <td>neg_adj2</td>
! </tr>
! <tr>
! <td>fv_sg_mod</td>
! <td>neg_adj4</td>
! </tr>
! <tr>
! <td>fv_timing_mod</td>
! <td>timing_on, timing_off</td>
! </tr>
Expand All @@ -116,6 +120,10 @@ module fv_dynamics_mod
! <td>neg_adj3</td>
! </tr>
! <tr>
! <td>fv_sg_mod</td>
! <td>neg_adj4</td>
! </tr>
! <tr>
! <td>tracer_manager_mod</td>
! <td>get_tracer_index</td>
! </tr>
Expand All @@ -137,7 +145,7 @@ module fv_dynamics_mod
use mpp_mod, only: mpp_pe
use field_manager_mod, only: MODEL_ATMOS
use tracer_manager_mod, only: get_tracer_index
use fv_sg_mod, only: neg_adj3, neg_adj2
use fv_sg_mod, only: neg_adj3, neg_adj2, neg_adj4
use fv_nesting_mod, only: setup_nested_grid_BCs
use fv_regional_mod, only: regional_boundary_update, set_regional_BCs
use fv_regional_mod, only: dump_field, H_STAGGER, U_STAGGER, V_STAGGER
Expand Down Expand Up @@ -270,7 +278,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
integer:: kord_tracer(ncnst)
integer :: i,j,k, n, iq, n_map, nq, nr, nwat, k_split
integer :: sphum, liq_wat = -999, ice_wat = -999 ! GFDL physics
integer :: rainwat = -999, snowwat = -999, graupel = -999, cld_amt = -999
integer :: rainwat = -999, snowwat = -999, graupel = -999, hailwat = -999, cld_amt = -999
integer :: theta_d = -999
logical used, do_omega
#ifdef MULTI_GASES
Expand Down Expand Up @@ -382,6 +390,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat')
snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat')
graupel = get_tracer_index (MODEL_ATMOS, 'graupel')
hailwat = get_tracer_index (MODEL_ATMOS, 'hailwat')
cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt')
endif

Expand All @@ -407,12 +416,12 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
#else
!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,nwat,q,q_con,sphum,liq_wat, &
#endif
!$OMP rainwat,ice_wat,snowwat,graupel) private(cvm,i,j,k)
!$OMP rainwat,ice_wat,snowwat,graupel,hailwat) private(cvm,i,j,k)
do k=1,npz
do j=js,je
#ifdef USE_COND
call moist_cp(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, q_con(is:ie,j,k), cvm)
ice_wat, snowwat, graupel, hailwat, q, q_con(is:ie,j,k), cvm)
#endif
do i=is,ie
dp1(i,j,k) = zvir*q(i,j,k,sphum)
Expand All @@ -425,7 +434,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
#else
!$OMP parallel do default(none) shared(is,ie,js,je,isd,ied,jsd,jed,npz,dp1,zvir,q,q_con,sphum,liq_wat, &
#endif
!$OMP rainwat,ice_wat,snowwat,graupel,pkz,flagstruct, &
!$OMP rainwat,ice_wat,snowwat,graupel,hailwat,pkz,flagstruct, &
#ifdef MULTI_GASES
!$OMP kapad, &
#endif
Expand All @@ -440,7 +449,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
do j=js,je
#ifdef MOIST_CAPPA
call moist_cv(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, q_con(is:ie,j,k), cvm)
ice_wat, snowwat, graupel, hailwat, q, q_con(is:ie,j,k), cvm)
#endif
do i=is,ie
#ifdef MULTI_GASES
Expand Down Expand Up @@ -505,7 +514,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
gridstruct%rsin2, gridstruct%cosa_s, &
zvir, cp_air, rdgas, hlv, te_2d, ua, va, teq, &
flagstruct%moist_phys, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, hydrostatic, idiag%id_te)
ice_wat, snowwat, graupel, hailwat, hydrostatic, idiag%id_te)
if( idiag%id_te>0 ) then
used = send_data(idiag%id_te, teq, fv_time)
! te_den=1.E-9*g_sum(teq, is, ie, js, je, ng, area, 0)/(grav*4.*pi*radius**2)
Expand Down Expand Up @@ -689,6 +698,9 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,snowwat), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)
if ( graupel > 0 ) &
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,graupel), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)
if ( hailwat > 0 ) &
call fill2D(is, ie, js, je, ng, npz, q(isd,jsd,1,hailwat), delp, gridstruct%area, domain, gridstruct%bounded_domain, npx, npy)

call timing_off('Fill2D')
endif
#endif
Expand Down Expand Up @@ -799,6 +811,44 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
used = send_data(idiag%id_mdt, dtdt_m, fv_time)
endif

if( nwat==7 ) then
if (cld_amt > 0) then
call neg_adj4(is, ie, js, je, ng, npz, &
flagstruct%hydrostatic, &
peln, delz, &
pt, delp, q(isd,jsd,1,sphum), &
q(isd,jsd,1,liq_wat), &
q(isd,jsd,1,rainwat), &
q(isd,jsd,1,ice_wat), &
q(isd,jsd,1,snowwat), &
q(isd,jsd,1,graupel), &
q(isd,jsd,1,hailwat), &
q(isd,jsd,1,cld_amt), flagstruct%check_negative)
else
call neg_adj4(is, ie, js, je, ng, npz, &
flagstruct%hydrostatic, &
peln, delz, &
pt, delp, q(isd,jsd,1,sphum), &
q(isd,jsd,1,liq_wat), &
q(isd,jsd,1,rainwat), &
q(isd,jsd,1,ice_wat), &
q(isd,jsd,1,snowwat), &
q(isd,jsd,1,graupel), &
q(isd,jsd,1,hailwat), check_negative=flagstruct%check_negative)
endif
if ( flagstruct%fv_debug ) then
call prt_mxm('T_dyn_a3', pt, is, ie, js, je, ng, npz, 1., gridstruct%area_64, domain)
call prt_mxm('SPHUM_dyn', q(isd,jsd,1,sphum ), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('liq_wat_dyn', q(isd,jsd,1,liq_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('rainwat_dyn', q(isd,jsd,1,rainwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
call prt_mxm('hailwat_dyn', q(isd,jsd,1,hailwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
endif
endif


if( nwat == 6 ) then
if (cld_amt > 0) then
call neg_adj3(is, ie, js, je, ng, npz, &
Expand Down
67 changes: 51 additions & 16 deletions model/fv_mapz.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
real rcp, rg, rrg, bkh, dtmp, k1k
integer:: i,j,k
integer:: kdelz
integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, ccn_cm3, iq, n, kp, k_next
integer:: nt, liq_wat, ice_wat, rainwat, snowwat, cld_amt, graupel, hailwat, ccn_cm3, iq, n, kp, k_next
integer :: ierr

ccpp_associate: associate( fast_mp_consv => CCPP_interstitial%fast_mp_consv, &
Expand All @@ -241,6 +241,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
rainwat = get_tracer_index (MODEL_ATMOS, 'rainwat')
snowwat = get_tracer_index (MODEL_ATMOS, 'snowwat')
graupel = get_tracer_index (MODEL_ATMOS, 'graupel')
hailwat = get_tracer_index (MODEL_ATMOS, 'hailwat')
cld_amt = get_tracer_index (MODEL_ATMOS, 'cld_amt')
ccn_cm3 = get_tracer_index (MODEL_ATMOS, 'ccn_cm3')

Expand All @@ -250,7 +251,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &

!$OMP parallel do default(none) shared(is,ie,js,je,km,pe,ptop,kord_tm,hydrostatic, &
!$OMP pt,pk,rg,peln,q,nwat,liq_wat,rainwat,ice_wat,snowwat, &
!$OMP graupel,q_con,sphum,cappa,r_vir,rcp,k1k,delp, &
!$OMP graupel,hailwat,q_con,sphum,cappa,r_vir,rcp,k1k,delp, &
!$OMP delz,akap,pkz,te,u,v,ps, gridstruct, last_step, &
!$OMP ak,bk,nq,isd,ied,jsd,jed,kord_tr,fill, adiabatic, &
#ifdef MULTI_GASES
Expand Down Expand Up @@ -293,7 +294,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
do k=1,km
#ifdef MOIST_CAPPA
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, gz, cvm)
ice_wat, snowwat, graupel, hailwat, q, gz, cvm)
do i=is,ie
q_con(i,j,k) = gz(i)
#ifdef MULTI_GASES
Expand Down Expand Up @@ -490,7 +491,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
do k=1,km
#ifdef MOIST_CAPPA
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, gz, cvm)
ice_wat, snowwat, graupel, hailwat, q, gz, cvm)
do i=is,ie
q_con(i,j,k) = gz(i)
#ifdef MULTI_GASES
Expand Down Expand Up @@ -612,7 +613,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
!$OMP parallel default(none) shared(is,ie,js,je,km,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, &
!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP graupel,hailwat,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
!$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, &
Expand All @@ -627,7 +628,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
!$OMP parallel default(none) shared(is,ie,js,je,km,kmp,ptop,u,v,pe,ua,va,isd,ied,jsd,jed,kord_mt, &
!$OMP te_2d,te,delp,hydrostatic,hs,rg,pt,peln, adiabatic, &
!$OMP cp,delz,nwat,rainwat,liq_wat,ice_wat,snowwat, &
!$OMP graupel,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP graupel,hailwat,q_con,r_vir,sphum,w,pk,pkz,last_step,consv, &
!$OMP do_adiabatic_init,zsum1,zsum0,te0_2d,domain, &
!$OMP ng,gridstruct,E_Flux,pdt,dtmp,reproduce_sum,q, &
!$OMP mdt,cld_amt,cappa,dtdt,out_dt,rrg,akap,do_sat_adj, &
Expand Down Expand Up @@ -689,7 +690,7 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
do k=1,km
#ifdef MOIST_CAPPA
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, gz, cvm)
ice_wat, snowwat, graupel, hailwat, q, gz, cvm)
do i=is,ie
! KE using 3D winds:
q_con(i,j,k) = gz(i)
Expand Down Expand Up @@ -816,11 +817,20 @@ subroutine Lagrangian_to_Eulerian(last_step, consv, ps, pe, delp, pkz, pk, &
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/ virq(q(i,j,k,1:num_gas))
#else
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i)))
#endif
enddo
elseif ( nwat==7 ) then
do i=is,ie
gz(i) = q(i,j,k,liq_wat)+q(i,j,k,rainwat)+q(i,j,k,ice_wat)+q(i,j,k,snowwat)+q(i,j,k,graupel)+q(i,j,k,hailwat)
#ifdef MULTI_GASES
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/ virq(q(i,j,k,1:num_gas))
#else
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k))/((1.+r_vir*q(i,j,k,sphum))*(1.-gz(i)))
#endif
enddo
else
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, gz, cvm)
ice_wat, snowwat, graupel, hailwat, q, gz, cvm)
do i=is,ie
#ifdef MULTI_GASES
pt(i,j,k) = (pt(i,j,k)+dtmp*pkz(i,j,k)) / virq(q(i,j,k,1:num_gas))
Expand Down Expand Up @@ -867,13 +877,13 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
u, v, w, delz, pt, delp, q, qc, pe, peln, hs, &
rsin2_l, cosa_s_l, &
r_vir, cp, rg, hlv, te_2d, ua, va, teq, &
moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hydrostatic, id_te)
moist_phys, nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, hydrostatic, id_te)
!------------------------------------------------------
! Compute vertically integrated total energy per column
!------------------------------------------------------
! !INPUT PARAMETERS:
integer, intent(in):: km, is, ie, js, je, isd, ied, jsd, jed, id_te
integer, intent(in):: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, nwat
integer, intent(in):: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, nwat
real, intent(inout), dimension(isd:ied,jsd:jed,km):: ua, va
real, intent(in), dimension(isd:ied,jsd:jed,km):: pt, delp
real, intent(in), dimension(isd:ied,jsd:jed,km,*):: q
Expand Down Expand Up @@ -908,7 +918,7 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
#ifdef MULTI_GASES
!$OMP num_gas, &
#endif
!$OMP q,nwat,liq_wat,rainwat,ice_wat,snowwat,graupel,sphum) &
!$OMP q,nwat,liq_wat,rainwat,ice_wat,snowwat,graupel,hailwat,sphum) &
!$OMP private(phiz, tv, cvm, qd)
do j=js,je

Expand Down Expand Up @@ -954,7 +964,7 @@ subroutine compute_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
do k=1,km
#ifdef MOIST_CAPPA
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, qd, cvm)
ice_wat, snowwat, graupel, hailwat, q, qd, cvm)
#endif
do i=is,ie
#ifdef MOIST_CAPPA
Expand Down Expand Up @@ -3325,9 +3335,9 @@ end subroutine mappm
!! including the heating capacity of water vapor and condensates.
!>@details See \cite emanuel1994atmospheric for information on variable heat capacities.
subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, qd, cvm, t1)
ice_wat, snowwat, graupel, hailwat, q, qd, cvm, t1)
integer, intent(in):: is, ie, isd,ied, jsd,jed, km, nwat, j, k
integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat
#ifdef MULTI_GASES
real, intent(in), dimension(isd:ied,jsd:jed,km,num_gas):: q
#else
Expand Down Expand Up @@ -3416,6 +3426,18 @@ subroutine moist_cv(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai
cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
#else
cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
#endif
enddo
case(7)
do i=is,ie
qv(i) = q(i,j,k,sphum)
ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat) + q(i,j,k,graupel) + q(i,j,k,hailwat)
qd(i) = ql(i) + qs(i)
#ifdef MULTI_GASES
cvm(i) = (1.-(qv(i)+qd(i)))*cv_air*vicvqd(q(i,j,k,1:num_gas)) + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
#else
cvm(i) = (1.-(qv(i)+qd(i)))*cv_air + qv(i)*cv_vap + ql(i)*c_liq + qs(i)*c_ice
#endif
enddo
case default
Expand All @@ -3435,10 +3457,10 @@ end subroutine moist_cv
!>@brief The subroutine 'moist_cp' computes the FV3-consistent moist heat capacity under constant pressure,
!! including the heating capacity of water vapor and condensates.
subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
ice_wat, snowwat, graupel, q, qd, cpm, t1)
ice_wat, snowwat, graupel, hailwat, q, qd, cpm, t1)

integer, intent(in):: is, ie, isd,ied, jsd,jed, km, nwat, j, k
integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
integer, intent(in):: sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat
real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q
real, intent(out), dimension(is:ie):: cpm, qd
real, intent(in), optional:: t1(is:ie)
Expand Down Expand Up @@ -3529,6 +3551,19 @@ subroutine moist_cp(is,ie, isd,ied, jsd,jed, km, j, k, nwat, sphum, liq_wat, rai
cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
#endif
enddo
case(7)
do i=is,ie
qv(i) = q(i,j,k,sphum)
ql(i) = q(i,j,k,liq_wat) + q(i,j,k,rainwat)
qs(i) = q(i,j,k,ice_wat) + q(i,j,k,snowwat) + q(i,j,k,graupel) + q(i,j,k,hailwat)
qd(i) = ql(i) + qs(i)
#ifdef MULTI_GASES
cpm(i) = (1.-(qv(i)+qd(i)))*cp_air*vicpqd(q(i,j,k,:)) + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
#else
cpm(i) = (1.-(qv(i)+qd(i)))*cp_air + qv(i)*cp_vapor + ql(i)*c_liq + qs(i)*c_ice
#endif
enddo

case default
call mpp_error (NOTE, 'fv_mapz::moist_cp - using default cp_air')
do i=is,ie
Expand Down
Loading