Skip to content

Commit

Permalink
Update to the latest GFDL dev/emc (#70)
Browse files Browse the repository at this point in the history 
* Add get_nth_domain_info subroutine

* Deallocate local arrays in fv_dyn_bundle_setup

* adding back a line that was mistakenly deleted in a previous commit (NOAA-GFDL#166)

* Do not print debug messages to stderr

* fix for 4diau with iau_filter_increments=T (NOAA-GFDL#167)

* fix for 4diau with iau_filter_increments=T

* fix time interval for 3DIAU

* fix typo in comment

* fix bug found in review by @MingjingTong-NOAA

* change tnext to integer variable itnext

Co-authored-by: jswhit2 <[email protected]>

* Use mpp_error instead of write statements in model/fv_regional_bc.F90

* Attempt at integrating fixes on top of dev/emc branch. (NOAA-GFDL#173)

* Support for cloud microphysics hail species (NOAA-GFDL#171)

* Incorporated Tim Supinie's updates to handle a hail category (nwat=7)

* Added print for hail max/min (non-debug section)

* Add hallwat to !OMP shared variables in main loop

Addresses compile failure when OMP is turned on (e.g., on Hera)

* Update to states as used in SFE2021
(Removed SFE code for now)

* Cleanup after rebase

* Redo removal of alt. namelist read

* Removed tools/module_diag_hailcast.F90

* Fixes from Jili Dong

* Add hailwat to tools/external_ic.F90

* Set logic in fv_regional_bc.F90 to match current code; Change logic in external_ic.F90 to check value of nwat

* Recommended changes

* Check actual index (hailwat) instead of assuming a positive value based on 'nwat' value

* Split nwat=7 bits into separate loops (see if it affects nwat=6)

Co-authored-by: Larissa Reames <[email protected]>
Co-authored-by: LarissaReames-NOAA <[email protected]>

Co-authored-by: Dusan Jovic <[email protected]>
Co-authored-by: laurenchilutti <[email protected]>
Co-authored-by: Jeff Whitaker <[email protected]>
Co-authored-by: jswhit2 <[email protected]>
Co-authored-by: MatthewPyle-NOAA <[email protected]>
Co-authored-by: Ted Mansell <[email protected]>
Co-authored-by: Larissa Reames <[email protected]>
Co-authored-by: LarissaReames-NOAA <[email protected]>
  • Loading branch information
@MinsukJi-NOAA @DusanJovic-NOAA
@laurenchilutti @jswhit @jswhit2 @MatthewPyle-NOAA @MicroTed @LarissaReames-NOAA
9 people authored Mar 30, 2022
1 parent c0d57b9 commit 19792ad
Show file tree
Hide file tree
Showing 10 changed files with 676 additions and 83 deletions.
2 changes: 1 addition & 1 deletion model/fv_arrays.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ module fv_arrays_mod
integer :: id_ws, id_te, id_amdt, id_mdt, id_divg, id_aam
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
66 changes: 58 additions & 8 deletions model/fv_dynamics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,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 @@ -115,6 +119,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 @@ -140,7 +148,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 @@ -283,7 +291,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
integer, parameter :: max_packs=13
Expand Down Expand Up @@ -399,6 +407,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 @@ -424,12 +433,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 @@ -442,7 +451,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 @@ -457,7 +466,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 @@ -522,7 +531,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 @@ -728,6 +737,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 @@ -786,6 +798,7 @@ subroutine fv_dynamics(npx, npy, npz, nq_tot, ng, bdt, consv_te, fill,
if (ice_wat > 0) call prt_mxm('ice_wat_dyn', q(isd,jsd,1,ice_wat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if (snowwat > 0) call prt_mxm('snowwat_dyn', q(isd,jsd,1,snowwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if (graupel > 0) call prt_mxm('graupel_dyn', q(isd,jsd,1,graupel), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
if (hailwat > 0) call prt_mxm('hailwat_dyn', q(isd,jsd,1,hailwat), is, ie, js, je, ng, npz, 1.,gridstruct%area_64, domain)
endif
#ifdef AVEC_TIMERS
call avec_timer_stop(6)
Expand Down Expand Up @@ -853,7 +866,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==6 ) then
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)
IF ( hailwat > 0 ) 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, &
flagstruct%hydrostatic, &
Expand Down
68 changes: 52 additions & 16 deletions model/fv_mapz.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -226,7 +226,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, kmp, kp, k_next
integer :: ierr

ccpp_associate: associate( fast_mp_consv => CCPP_interstitial%fast_mp_consv, &
Expand All @@ -242,6 +242,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 @@ -251,7 +252,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 @@ -294,7 +295,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 @@ -545,7 +546,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 @@ -667,7 +668,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 @@ -682,7 +683,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 @@ -744,7 +745,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 @@ -874,11 +875,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 @@ -925,13 +935,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 @@ -966,7 +976,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 @@ -1012,7 +1022,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 @@ -3408,9 +3418,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 @@ -3499,6 +3509,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 @@ -3518,10 +3540,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
#ifdef MULTI_GASES
real, intent(in), dimension(isd:ied,jsd:jed,km,num_gas):: q
#else
Expand Down Expand Up @@ -3617,6 +3639,20 @@ 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

0 comments on commit 19792ad

Please sign in to comment.