diff --git a/model/fv_arrays.F90 b/model/fv_arrays.F90
index a353e0a65..27f44c04b 100644
--- a/model/fv_arrays.F90
+++ b/model/fv_arrays.F90
@@ -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
diff --git a/model/fv_dynamics.F90 b/model/fv_dynamics.F90
index 0523ccbe5..1143324b7 100644
--- a/model/fv_dynamics.F90
+++ b/model/fv_dynamics.F90
@@ -95,6 +95,10 @@ module fv_dynamics_mod
! | neg_adj2 |
!
!
+! | fv_sg_mod |
+! neg_adj4 |
+!
+!
! | fv_timing_mod |
! timing_on, timing_off |
!
@@ -115,6 +119,10 @@ module fv_dynamics_mod
! neg_adj3 |
!
!
+! | fv_sg_mod |
+! neg_adj4 |
+!
+!
! | tracer_manager_mod |
! get_tracer_index |
!
@@ -136,7 +144,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
@@ -279,7 +287,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
@@ -395,6 +403,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
@@ -420,12 +429,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)
@@ -438,7 +447,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
@@ -453,7 +462,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
@@ -518,7 +527,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)
@@ -724,6 +733,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
@@ -782,6 +794,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)
@@ -849,7 +862,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, &
diff --git a/model/fv_mapz.F90 b/model/fv_mapz.F90
index d3313a444..add25e28b 100644
--- a/model/fv_mapz.F90
+++ b/model/fv_mapz.F90
@@ -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, &
@@ -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')
@@ -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
@@ -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
@@ -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
@@ -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, &
@@ -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, &
@@ -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)
@@ -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))
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
diff --git a/model/fv_nesting.F90 b/model/fv_nesting.F90
index 92124963f..c41dd3e9c 100644
--- a/model/fv_nesting.F90
+++ b/model/fv_nesting.F90
@@ -1416,15 +1416,15 @@ subroutine setup_pt_NH_BC(pt_BC, delp_BC, delz_BC, sphum_BC, q_BC, nq, &
real, parameter:: c_ice = 1972. !< heat capacity of ice at 0C: c=c_ice+7.3*(T-Tice)
real, parameter:: cv_vap = cp_vapor - rvgas !< 1384.5
- real, dimension(:,:,:), pointer :: liq_watBC_west, ice_watBC_west, rainwatBC_west, snowwatBC_west, graupelBC_west
- real, dimension(:,:,:), pointer :: liq_watBC_east, ice_watBC_east, rainwatBC_east, snowwatBC_east, graupelBC_east
- real, dimension(:,:,:), pointer :: liq_watBC_north, ice_watBC_north, rainwatBC_north, snowwatBC_north, graupelBC_north
- real, dimension(:,:,:), pointer :: liq_watBC_south, ice_watBC_south, rainwatBC_south, snowwatBC_south, graupelBC_south
+ real, dimension(:,:,:), pointer :: liq_watBC_west, ice_watBC_west, rainwatBC_west, snowwatBC_west, graupelBC_west, hailwatBC_west
+ real, dimension(:,:,:), pointer :: liq_watBC_east, ice_watBC_east, rainwatBC_east, snowwatBC_east, graupelBC_east, hailwatBC_east
+ real, dimension(:,:,:), pointer :: liq_watBC_north, ice_watBC_north, rainwatBC_north, snowwatBC_north, graupelBC_north, hailwatBC_north
+ real, dimension(:,:,:), pointer :: liq_watBC_south, ice_watBC_south, rainwatBC_south, snowwatBC_south, graupelBC_south, hailwatBC_south
real :: dp1, q_liq, q_sol, q_con = 0., cvm, pkz, rdg, cv_air
integer :: i,j,k, istart, iend
- integer :: liq_wat, ice_wat, rainwat, snowwat, graupel
+ integer :: liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat
real, parameter:: tice = 273.16 !< For GFS Partitioning
real, parameter:: t_i0 = 15.
@@ -1557,11 +1557,22 @@ subroutine setup_pt_NH_BC(pt_BC, delp_BC, delz_BC, sphum_BC, q_BC, nq, &
graupelBC_north => dum_north
graupelBC_south => dum_south
endif
+ if (hailwat > 0) then
+ hailwatBC_west => q_BC(hailwat)%west_t1
+ hailwatBC_east => q_BC(hailwat)%east_t1
+ hailwatBC_north => q_BC(hailwat)%north_t1
+ hailwatBC_south => q_BC(hailwat)%south_t1
+ else
+ hailwatBC_west => dum_west
+ hailwatBC_east => dum_east
+ hailwatBC_north => dum_north
+ hailwatBC_south => dum_south
+ endif
if (is == 1) then
call setup_pt_NH_BC_k(pt_BC%west_t1, sphum_BC%west_t1, delp_BC%west_t1, delz_BC%west_t1, &
- liq_watBC_west, rainwatBC_west, ice_watBC_west, snowwatBC_west, graupelBC_west, &
+ liq_watBC_west, rainwatBC_west, ice_watBC_west, snowwatBC_west, graupelBC_west, hailwatBC_west, &
#ifdef USE_COND
q_con_BC%west_t1, &
#ifdef MOIST_CAPPA
@@ -1585,7 +1596,7 @@ subroutine setup_pt_NH_BC(pt_BC, delp_BC, delz_BC, sphum_BC, q_BC, nq, &
end if
call setup_pt_NH_BC_k(pt_BC%south_t1, sphum_BC%south_t1, delp_BC%south_t1, delz_BC%south_t1, &
- liq_watBC_south, rainwatBC_south, ice_watBC_south, snowwatBC_south, graupelBC_south, &
+ liq_watBC_south, rainwatBC_south, ice_watBC_south, snowwatBC_south, graupelBC_south, hailwatBC_south, &
#ifdef USE_COND
q_con_BC%south_t1, &
#ifdef MOIST_CAPPA
@@ -1599,7 +1610,7 @@ subroutine setup_pt_NH_BC(pt_BC, delp_BC, delz_BC, sphum_BC, q_BC, nq, &
if (ie == npx-1) then
call setup_pt_NH_BC_k(pt_BC%east_t1, sphum_BC%east_t1, delp_BC%east_t1, delz_BC%east_t1, &
- liq_watBC_east, rainwatBC_east, ice_watBC_east, snowwatBC_east, graupelBC_east, &
+ liq_watBC_east, rainwatBC_east, ice_watBC_east, snowwatBC_east, graupelBC_east, hailwatBC_east, &
#ifdef USE_COND
q_con_BC%east_t1, &
#ifdef MOIST_CAPPA
@@ -1622,7 +1633,7 @@ subroutine setup_pt_NH_BC(pt_BC, delp_BC, delz_BC, sphum_BC, q_BC, nq, &
end if
call setup_pt_NH_BC_k(pt_BC%north_t1, sphum_BC%north_t1, delp_BC%north_t1, delz_BC%north_t1, &
- liq_watBC_north, rainwatBC_north, ice_watBC_north, snowwatBC_north, graupelBC_north, &
+ liq_watBC_north, rainwatBC_north, ice_watBC_north, snowwatBC_north, graupelBC_north, hailwatBC_north, &
#ifdef USE_COND
q_con_BC%north_t1, &
#ifdef MOIST_CAPPA
@@ -1636,7 +1647,7 @@ end subroutine setup_pt_NH_BC
subroutine setup_pt_NH_BC_k(ptBC,sphumBC,delpBC,delzBC, &
- liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC, &
+ liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC,hailwatBC, &
#ifdef USE_COND
q_conBC, &
#ifdef MOIST_CAPPA
@@ -1648,7 +1659,7 @@ subroutine setup_pt_NH_BC_k(ptBC,sphumBC,delpBC,delzBC, &
integer, intent(IN) :: isd_BC, ied_BC, istart, iend, jstart, jend, npz
real, intent(OUT), dimension(isd_BC:ied_BC,jstart:jend,npz) :: ptBC
real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: sphumBC, delpBC, delzBC
- real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC
+ real, intent(IN), dimension(isd_BC:ied_BC,jstart:jend,npz) :: liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC,hailwatBC
#ifdef USE_COND
real, intent(OUT), dimension(isd_BC:ied_BC,jstart:jend,npz) :: q_conBC
#ifdef MOIST_CAPPA
@@ -1669,7 +1680,7 @@ subroutine setup_pt_NH_BC_k(ptBC,sphumBC,delpBC,delzBC, &
rdg = -rdgas / grav
cv_air = cp_air - rdgas
-!$OMP parallel do default(none) shared(istart,iend,jstart,jend,npz,zvir,ptBC,sphumBC,delpBC,delzBC,liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC, &
+!$OMP parallel do default(none) shared(istart,iend,jstart,jend,npz,zvir,ptBC,sphumBC,delpBC,delzBC,liq_watBC,rainwatBC,ice_watBC,snowwatBC,graupelBC,hailwatBC, &
#ifdef USE_COND
!$OMP q_conBC, &
#ifdef MOIST_CAPPA
@@ -1684,7 +1695,7 @@ subroutine setup_pt_NH_BC_k(ptBC,sphumBC,delpBC,delzBC, &
dp1 = zvir*sphumBC(i,j,k)
#ifdef USE_COND
q_liq = liq_watBC(i,j,k) + rainwatBC(i,j,k)
- q_sol = ice_watBC(i,j,k) + snowwatBC(i,j,k) + graupelBC(i,j,k)
+ q_sol = ice_watBC(i,j,k) + snowwatBC(i,j,k) + graupelBC(i,j,k) + hailwatBC(i,j,k)
q_con = q_liq + q_sol
q_conBC(i,j,k) = q_con
#ifdef MOIST_CAPPA
diff --git a/model/fv_regional_bc.F90 b/model/fv_regional_bc.F90
index 1d83443b8..06a5158b4 100644
--- a/model/fv_regional_bc.F90
+++ b/model/fv_regional_bc.F90
@@ -137,6 +137,7 @@ module fv_regional_mod
!
integer,save :: cld_amt_index & !<--
,graupel_index & !
+ ,hailwat_index & !
,ice_water_index & ! Locations of
,liq_water_index & ! tracer vbls
,o3mr_index & ! in the tracers
@@ -735,6 +736,7 @@ subroutine setup_regional_BC(Atm &
rain_water_index = get_tracer_index(MODEL_ATMOS, 'rainwat')
snow_water_index = get_tracer_index(MODEL_ATMOS, 'snowwat')
graupel_index = get_tracer_index(MODEL_ATMOS, 'graupel')
+ hailwat_index = get_tracer_index(MODEL_ATMOS, 'hailwat')
cld_amt_index = get_tracer_index(MODEL_ATMOS, 'cld_amt')
o3mr_index = get_tracer_index(MODEL_ATMOS, 'o3mr')
! write(0,*)' setup_regional_bc'
@@ -3479,7 +3481,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
real(kind=R_GRID):: pst
!!! High-precision
integer i,ie,is,j,je,js,k,l,m, k2,iq
- integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt
+ integer sphum, o3mr, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, cld_amt
!
!---------------------------------------------------------------------------------
!
@@ -3489,6 +3491,7 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
rainwat = rain_water_index
snowwat = snow_water_index
graupel = graupel_index
+ hailwat = hailwat_index
cld_amt = cld_amt_index
o3mr = o3mr_index
@@ -3752,7 +3755,10 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
! and may not provide a very good result
!
if ( .not. data_source_fv3gfs ) then
- if ( Atm%flagstruct%nwat .eq. 6 ) then
+ if ( Atm%flagstruct%nwat .eq. 6 .or. Atm%flagstruct%nwat .eq. 7 ) then
+ if ( hailwat > 0 ) then
+ BC_side%q_BC(is:ie,j,1:npz,hailwat) = 0.
+ endif
do k=1,npz
do i=is,ie
qn1(i,k) = BC_side%q_BC(i,j,k,liq_wat)
@@ -3795,7 +3801,8 @@ subroutine remap_scalar_nggps_regional_bc(Atm &
enddo
enddo
endif
- endif ! data source /= FV3GFS GAUSSIAN NEMSIO/NETCDF and GRIB2 FILE
+
+ endif ! data source /= FV3GFS GAUSSIAN NEMSIO FILE
!
! For GFS spectral input, omega in pa/sec is stored as w in the input data so actual w(m/s) is calculated
! For GFS nemsio input, omega is 0, so best not to use for input since boundary data will not exist for w
diff --git a/model/fv_sg.F90 b/model/fv_sg.F90
index 06b401dbf..2446b1144 100644
--- a/model/fv_sg.F90
+++ b/model/fv_sg.F90
@@ -68,7 +68,7 @@ module fv_sg_mod
implicit none
private
-public fv_subgrid_z, qsmith, neg_adj3, neg_adj2
+public fv_subgrid_z, qsmith, neg_adj3, neg_adj2, neg_adj4
real, parameter:: esl = 0.621971831
real, parameter:: tice = 273.16
@@ -147,7 +147,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
integer i, j, k, kk, n, m, iq, km1, im, kbot, l
real, parameter:: ustar2 = 1.E-4
real:: cv_air, xvir
- integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
+ integer :: sphum, liq_wat, rainwat, snowwat, graupel, hailwat, ice_wat, cld_amt
cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
rk = cp_air/rdgas + 1.
@@ -188,6 +188,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
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
@@ -200,7 +201,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
!$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
!$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
-!$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra, t_max, t_min, &
+!$OMP snowwat,cv_air,m,graupel,hailwat,pkz,rk,rz,fra, t_max, t_min, &
#ifdef MULTI_GASES
!$OMP u_dt,rdt,v_dt,xvir,nwat,km) &
#else
@@ -308,7 +309,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
- elseif ( nwat==5 ) then
+ elseif ( nwat == 5 ) then
do i=is,ie
q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat)
@@ -320,7 +321,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
- else
+ elseif ( nwat == 6 ) then
do i=is,ie
q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
@@ -330,6 +331,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
#else
cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#endif
+ enddo
+ elseif ( nwat == 7 ) then
+ do i=is,ie
+ q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
+ q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel) + q0(i,k,hailwat)
+#ifdef MULTI_GASES
+ cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#else
+ cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
endif
@@ -456,9 +469,12 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
elseif ( nwat==5 ) then ! K_warm_rain scheme with fake ice
qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
q0(i,km1,snowwat) + q0(i,km1,rainwat)
- else
+ elseif ( nwat==6 ) then
qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
+ elseif ( nwat==7 ) then
+ qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
+ q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel) + q0(i,km1,hailwat)
endif
! u:
h0 = mc*(u0(i,k)-u0(i,k-1))
@@ -595,7 +611,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
- else
+ elseif ( nwat == 6 ) THEN
do i=is,ie
q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
@@ -605,6 +621,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
#else
cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#endif
+ enddo
+ elseif ( nwat == 7 ) THEN
+ do i=is,ie
+ q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
+ q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel) + q0(i,kk,hailwat)
+#ifdef MULTI_GASES
+ cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#else
+ cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
endif
@@ -715,7 +743,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
integer i, j, k, kk, n, m, iq, km1, im, kbot
real, parameter:: ustar2 = 1.E-4
real:: cv_air, xvir
- integer :: sphum, liq_wat, rainwat, snowwat, graupel, ice_wat, cld_amt
+ integer :: sphum, liq_wat, rainwat, snowwat, graupel, hailwat, ice_wat, cld_amt
cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
rk = cp_air/rdgas + 1.
@@ -745,6 +773,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
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
@@ -757,7 +786,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
!$OMP parallel do default(none) shared(im,is,ie,js,je,nq,kbot,qa,ta,sphum,ua,va,delp,peln, &
!$OMP hydrostatic,pe,delz,g2,w,liq_wat,rainwat,ice_wat, &
-!$OMP snowwat,cv_air,m,graupel,pkz,rk,rz,fra,cld_amt, &
+!$OMP snowwat,cv_air,m,graupel,hailwat,pkz,rk,rz,fra,cld_amt, &
!$OMP u_dt,rdt,v_dt,xvir,nwat) &
!$OMP private(kk,lcp2,icp2,tcp3,dh,dq,den,qs,qsw,dqsdt,qcon,q0, &
!$OMP t0,u0,v0,w0,h0,pm,gzh,tvm,tmp,cpm,cvm, q_liq,q_sol,&
@@ -874,7 +903,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
- else
+ elseif ( nwat == 6 ) THEN
do i=is,ie
q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
@@ -884,6 +913,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
#else
cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#endif
+ enddo
+ elseif ( nwat == 7 ) THEN
+ do i=is,ie
+ q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
+ q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel) + q0(i,k,hailwat)
+#ifdef MULTI_GASES
+ cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#else
+ cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
endif
@@ -1005,9 +1046,12 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
elseif ( nwat==5 ) then
qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
q0(i,km1,snowwat) + q0(i,km1,rainwat)
- else
+ elseif ( nwat==6 ) then
qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel)
+ elseif ( nwat==7 ) then
+ qcon(i,km1) = q0(i,km1,liq_wat) + q0(i,km1,ice_wat) + &
+ q0(i,km1,snowwat) + q0(i,km1,rainwat) + q0(i,km1,graupel) + q0(i,km1,hailwat)
endif
! u:
h0 = mc*(u0(i,k)-u0(i,k-1))
@@ -1143,7 +1187,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
- else
+ elseif ( nwat == 6 ) THEN
do i=is,ie
q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel)
@@ -1153,6 +1197,18 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
#else
cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#endif
+ enddo
+ elseif ( nwat == 7 ) THEN
+ do i=is,ie
+ q_liq = q0(i,kk,liq_wat) + q0(i,kk,rainwat)
+ q_sol = q0(i,kk,ice_wat) + q0(i,kk,snowwat) + q0(i,kk,graupel) + q0(i,kk,hailwat)
+#ifdef MULTI_GASES
+ cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,kk,:)) + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,kk,:)) + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
+#else
+ cpm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cp_air + q0(i,kk,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ cvm(i) = (1.-(q0(i,kk,sphum)+q_liq+q_sol))*cv_air + q0(i,kk,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#endif
enddo
endif
@@ -1199,7 +1255,7 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
! Saturation adjustment
!----------------------
#ifndef GFS_PHYS
- if ( nwat > 5 ) then
+ if ( nwat == 6 ) then
do k=1, kbot
if ( hydrostatic ) then
do i=is, ie
@@ -1211,6 +1267,9 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
#endif
q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
+! IF ( nwat == 7 ) THEN
+! q_sol = q_sol + q0(i,k,hailwat)
+! ENDIF
#ifdef MULTI_GASES
cpm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cp_air*vicpqd(q0(i,k,:)) + q0(i,k,sphum)*cp_vapor + q_liq*c_liq + q_sol*c_ice
#else
@@ -1224,6 +1283,9 @@ subroutine fv_subgrid_z( isd, ied, jsd, jed, is, ie, js, je, km, nq, dt, &
den(i,k) = -delp(i,j,k)/(grav*delz(i,j,k))
q_liq = q0(i,k,liq_wat) + q0(i,k,rainwat)
q_sol = q0(i,k,ice_wat) + q0(i,k,snowwat) + q0(i,k,graupel)
+! IF ( nwat == 7 ) THEN
+! q_sol = q_sol + q0(i,k,hailwat)
+! ENDIF
#ifdef MULTI_GASES
cvm(i) = (1.-(q0(i,k,sphum)+q_liq+q_sol))*cv_air*vicvqd(q0(i,k,:)) + q0(i,k,sphum)*cv_vap + q_liq*c_liq + q_sol*c_ice
#else
@@ -1878,6 +1940,364 @@ subroutine neg_adj3(is, ie, js, je, ng, kbot, hydrostatic, peln, delz, pt, dp,
end subroutine neg_adj3
+! TAS: XXX check to make sure this doesn't need any modifications vs neg_adj2 or neg_adj3
+ subroutine neg_adj4(is, ie, js, je, ng, kbot, hydrostatic, &
+ peln, delz, pt, dp, qv, ql, qr, qi, qs, qg, qh, qa, check_negative)
+
+! This is designed for 7-class micro-physics schemes
+ integer, intent(in):: is, ie, js, je, ng, kbot
+ logical, intent(in):: hydrostatic
+ real, intent(in):: dp(is-ng:ie+ng,js-ng:je+ng,kbot) ! total delp-p
+ real, intent(in):: delz(is:,js:,1:)
+ real, intent(in):: peln(is:ie,kbot+1,js:je) ! ln(pe)
+ logical, intent(in), OPTIONAL :: check_negative
+ real, intent(inout), dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: &
+ pt, qv, ql, qr, qi, qs, qg, qh
+ real, intent(inout), OPTIONAL, dimension(is-ng:ie+ng,js-ng:je+ng,kbot):: qa
+! Local:
+ logical:: sat_adj = .false.
+ real, parameter :: t48 = tice - 48.
+ real, dimension(is:ie,kbot):: dpk, q2
+real, dimension(is:ie,js:je):: pt2, qv2, ql2, qi2, qs2, qr2, qg2, qh2, dp2, p2, icpk, lcpk
+ real:: cv_air
+ real:: dq, qsum, dq1, q_liq, q_sol, cpm, sink, qsw, dwsdt
+ integer i, j, k
+
+ cv_air = cp_air - rdgas ! = rdgas * (7/2-1) = 2.5*rdgas=717.68
+
+ if ( present(check_negative) ) then
+ if ( check_negative ) then
+ call prt_negative('Temperature', pt, is, ie, js, je, ng, kbot, 165.)
+ call prt_negative('sphum', qv, is, ie, js, je, ng, kbot, -1.e-8)
+ call prt_negative('liq_wat', ql, is, ie, js, je, ng, kbot, -1.e-7)
+ call prt_negative('rainwat', qr, is, ie, js, je, ng, kbot, -1.e-7)
+ call prt_negative('ice_wat', qi, is, ie, js, je, ng, kbot, -1.e-7)
+ call prt_negative('snowwat', qs, is, ie, js, je, ng, kbot, -1.e-7)
+ call prt_negative('graupel', qg, is, ie, js, je, ng, kbot, -1.e-7)
+ call prt_negative('hailwat', qh, is, ie, js, je, ng, kbot, -1.e-7)
+ endif
+ endif
+
+ if ( hydrostatic ) then
+ d0_vap = cp_vapor - c_liq
+ else
+ d0_vap = cv_vap - c_liq
+ endif
+ lv00 = hlv0 - d0_vap*t_ice
+
+!$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,ql,qi,qs,qr,qg,qh,dp,pt, &
+!$OMP lv00, d0_vap,hydrostatic,peln,delz,cv_air,sat_adj) &
+!$OMP private(dq,dq1,qsum,dp2,p2,pt2,qv2,ql2,qi2,qs2,qg2,qh2,qr2, &
+!$OMP lcpk,icpk,qsw,dwsdt,sink,q_liq,q_sol,cpm)
+ do k=1, kbot
+ do j=js, je
+ do i=is, ie
+ qv2(i,j) = qv(i,j,k)
+ ql2(i,j) = ql(i,j,k)
+ qi2(i,j) = qi(i,j,k)
+ qs2(i,j) = qs(i,j,k)
+ qr2(i,j) = qr(i,j,k)
+ qg2(i,j) = qg(i,j,k)
+ qh2(i,j) = qh(i,j,k)
+ dp2(i,j) = dp(i,j,k)
+ pt2(i,j) = pt(i,j,k)
+ enddo
+ enddo
+
+ if ( hydrostatic ) then
+ do j=js, je
+ do i=is, ie
+ p2(i,j) = dp2(i,j)/(peln(i,k+1,j)-peln(i,k,j))
+ q_liq = max(0., ql2(i,j) + qr2(i,j))
+ q_sol = max(0., qi2(i,j) + qs2(i,j))
+ cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cp_air + qv2(i,j)*cp_vapor + q_liq*c_liq + q_sol*c_ice
+ lcpk(i,j) = hlv / cpm
+ icpk(i,j) = hlf / cpm
+ enddo
+ enddo
+ else
+ do j=js, je
+ do i=is, ie
+ p2(i,j) = -dp2(i,j)/(grav*delz(i,j,k))*rdgas*pt2(i,j)*(1.+zvir*qv2(i,j))
+ q_liq = max(0., ql2(i,j) + qr2(i,j))
+ q_sol = max(0., qi2(i,j) + qs2(i,j))
+ cpm = (1.-(qv2(i,j)+q_liq+q_sol))*cv_air + qv2(i,j)*cv_vap + q_liq*c_liq + q_sol*c_ice
+ lcpk(i,j) = (lv00+d0_vap*pt2(i,j)) / cpm
+ icpk(i,j) = (Li0+dc_ice*pt2(i,j)) / cpm
+ enddo
+ enddo
+ endif
+
+! Fix the negatives:
+!-----------
+! Ice-phase:
+!-----------
+ do j=js, je
+ do i=is, ie
+ qsum = qi2(i,j) + qs2(i,j)
+ if ( qsum > 0. ) then
+ if ( qi2(i,j) < 0. ) then
+ qi2(i,j) = 0.
+ qs2(i,j) = qsum
+ elseif ( qs2(i,j) < 0. ) then
+ qs2(i,j) = 0.
+ qi2(i,j) = qsum
+ endif
+ else
+! borrow froom graupel
+ qi2(i,j) = 0.
+ qs2(i,j) = 0.
+ qg2(i,j) = qg2(i,j) + qsum
+ endif
+
+! At this stage qi and qs should be positive definite
+! If graupel < 0 then borrow from qs then qi
+ if ( qg2(i,j) < 0. ) then
+ dq = min( qs2(i,j), -qg2(i,j) )
+ qs2(i,j) = qs2(i,j) - dq
+ qg2(i,j) = qg2(i,j) + dq
+ if ( qg2(i,j) < 0. ) then
+! if qg is still negative
+ dq = min( qi2(i,j), -qg2(i,j) )
+ qi2(i,j) = qi2(i,j) - dq
+ qg2(i,j) = qg2(i,j) + dq
+ endif
+ endif
+
+! If qg is still negative then borrow from rain water: phase change
+ if ( qg2(i,j)<0. .and. qr2(i,j)>0. ) then
+ dq = min( qr2(i,j), -qg2(i,j) )
+ qg2(i,j) = qg2(i,j) + dq
+ qr2(i,j) = qr2(i,j) - dq
+ pt2(i,j) = pt2(i,j) + dq*icpk(i,j) ! conserve total energy
+ endif
+! If qg is still negative then borrow from cloud water: phase change
+ if ( qg2(i,j)<0. .and. ql2(i,j)>0. ) then
+ dq = min( ql2(i,j), -qg2(i,j) )
+ qg2(i,j) = qg2(i,j) + dq
+ ql2(i,j) = ql2(i,j) - dq
+ pt2(i,j) = pt2(i,j) + dq*icpk(i,j)
+ endif
+! Last resort; borrow from water vapor
+ if ( qg2(i,j)<0. .and. qv2(i,j)>0. ) then
+ dq = min( 0.999*qv2(i,j), -qg2(i,j) )
+ qg2(i,j) = qg2(i,j) + dq
+ qv2(i,j) = qv2(i,j) - dq
+ pt2(i,j) = pt2(i,j) + dq*(icpk(i,j)+lcpk(i,j))
+ endif
+
+!--------------
+! Liquid phase:
+!--------------
+ qsum = ql2(i,j) + qr2(i,j)
+ if ( qsum > 0. ) then
+ if ( qr2(i,j) < 0. ) then
+ qr2(i,j) = 0.
+ ql2(i,j) = qsum
+ elseif ( ql2(i,j) < 0. ) then
+ ql2(i,j) = 0.
+ qr2(i,j) = qsum
+ endif
+ else
+ ql2(i,j) = 0.
+ qr2(i,j) = qsum ! rain water is still negative
+! fill negative rain with qg first
+ dq = min( max(0.0, qg2(i,j)), -qr2(i,j) )
+ qr2(i,j) = qr2(i,j) + dq
+ qg2(i,j) = qg2(i,j) - dq
+ pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
+ if ( qr(i,j,k) < 0. ) then
+! fill negative rain with available qi & qs (cooling)
+ dq = min( qi2(i,j)+qs2(i,j), -qr2(i,j) )
+ qr2(i,j) = qr2(i,j) + dq
+ dq1 = min( dq, qs2(i,j) )
+ qs2(i,j) = qs2(i,j) - dq1
+ qi2(i,j) = qi2(i,j) + dq1 - dq
+ pt2(i,j) = pt2(i,j) - dq*icpk(i,j)
+ endif
+! fix negative rain water with available vapor
+ if ( qr2(i,j)<0. .and. qv2(i,j)>0. ) then
+ dq = min( 0.999*qv2(i,j), -qr2(i,j) )
+ qv2(i,j) = qv2(i,j) - dq
+ qr2(i,j) = qr2(i,j) + dq
+ pt2(i,j) = pt2(i,j) + dq*lcpk(i,j)
+ endif
+ endif
+ enddo
+ enddo
+
+!******************************************
+! Fast moist physics: Saturation adjustment
+!******************************************
+#ifndef GFS_PHYS
+ if ( sat_adj ) then
+
+ do j=js, je
+ do i=is, ie
+! Melting of cloud ice into cloud water ********
+ if ( qi2(i,j)>1.e-8 .and. pt2(i,j) > tice ) then
+ sink = min( qi2(i,j), (pt2(i,j)-tice)/icpk(i,j) )
+ ql2(i,j) = ql2(i,j) + sink
+ qi2(i,j) = qi2(i,j) - sink
+ pt2(i,j) = pt2(i,j) - sink*icpk(i,j)
+ endif
+
+! vapor <---> liquid water --------------------------------
+ qsw = wqsat2_moist(pt2(i,j), qv2(i,j), p2(i,j), dwsdt)
+ sink = min( ql2(i,j), (qsw-qv2(i,j))/(1.+lcpk(i,j)*dwsdt) )
+ qv2(i,j) = qv2(i,j) + sink
+ ql2(i,j) = ql2(i,j) - sink
+ pt2(i,j) = pt2(i,j) - sink*lcpk(i,j)
+!-----------------------------------------------------------
+
+! freezing of cloud water ********
+ if( ql2(i,j)>1.e-8 .and. pt2(i,j) < t48 ) then
+! Enforce complete freezing below t_00 (-48 C)
+ sink = min( ql2(i,j), (t48-pt2(i,j))/icpk(i,j) )
+ ql2(i,j) = ql2(i,j) - sink
+ qi2(i,j) = qi2(i,j) + sink
+ pt2(i,j) = pt2(i,j) + sink*icpk(i,j)
+ endif ! significant ql existed
+ enddo
+ enddo
+ endif
+#endif
+
+!----------------------------------------------------------------
+! Update fields:
+ do j=js, je
+ do i=is, ie
+ qv(i,j,k) = qv2(i,j)
+ ql(i,j,k) = ql2(i,j)
+ qi(i,j,k) = qi2(i,j)
+ qs(i,j,k) = qs2(i,j)
+ qr(i,j,k) = qr2(i,j)
+ qg(i,j,k) = qg2(i,j)
+ pt(i,j,k) = pt2(i,j)
+ enddo
+ enddo
+
+ enddo
+
+!$OMP parallel do default(none) shared(is,ie,js,je,kbot,dp,qg,qr) &
+!$OMP private(dpk, q2)
+ do j=js, je
+! Graupel:
+ do k=1,kbot
+ do i=is,ie
+ dpk(i,k) = dp(i,j,k)
+ q2(i,k) = qg(i,j,k)
+ enddo
+ enddo
+ call fillq(ie-is+1, kbot, q2, dpk)
+ do k=1,kbot
+ do i=is,ie
+ qg(i,j,k) = q2(i,k)
+ enddo
+ enddo
+! Rain water:
+ do k=1,kbot
+ do i=is,ie
+ q2(i,k) = qr(i,j,k)
+ enddo
+ enddo
+ call fillq(ie-is+1, kbot, q2, dpk)
+ do k=1,kbot
+ do i=is,ie
+ qr(i,j,k) = q2(i,k)
+ enddo
+ enddo
+ enddo
+
+!-----------------------------------
+! Fix water vapor
+!-----------------------------------
+! Top layer: borrow from below
+ k = 1
+!$OMP parallel do default(none) shared(is,ie,js,je,k,qv,dp)
+ do j=js, je
+ do i=is, ie
+ if( qv(i,j,k) < 0. ) then
+ qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
+ qv(i,j,k ) = 0.
+ endif
+ enddo
+ enddo
+
+! this OpenMP do-loop cannot be parallelized with recursion on k/k-1
+!$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) &
+!$OMP private(dq)
+ do j=js, je
+ do k=2,kbot-1
+ do i=is, ie
+ if( qv(i,j,k) < 0. .and. qv(i,j,k-1) > 0. ) then
+ dq = min(-qv(i,j,k)*dp(i,j,k), qv(i,j,k-1)*dp(i,j,k-1))
+ qv(i,j,k-1) = qv(i,j,k-1) - dq/dp(i,j,k-1)
+ qv(i,j,k ) = qv(i,j,k ) + dq/dp(i,j,k )
+ endif
+ if( qv(i,j,k) < 0. ) then
+ qv(i,j,k+1) = qv(i,j,k+1) + qv(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
+ qv(i,j,k ) = 0.
+ endif
+ enddo
+ enddo
+ enddo
+
+! Bottom layer; Borrow from above
+!$OMP parallel do default(none) shared(is,ie,js,je,kbot,qv,dp) private(dq)
+ do j=js, je
+ do i=is, ie
+ if( qv(i,j,kbot) < 0. ) then
+ do k=kbot-1,1,-1
+ if ( qv(i,j,kbot)>=0. ) goto 123
+ if ( qv(i,j,k) > 0. ) then
+ dq = min(-qv(i,j,kbot)*dp(i,j,kbot), qv(i,j,k)*dp(i,j,k))
+ qv(i,j,k ) = qv(i,j,k ) - dq/dp(i,j,k)
+ qv(i,j,kbot) = qv(i,j,kbot) + dq/dp(i,j,kbot)
+ endif
+ enddo ! k-loop
+123 continue
+ endif
+ enddo ! i-loop
+ enddo ! j-loop
+
+
+ if (present(qa)) then
+!-----------------------------------
+! Fix negative cloud fraction
+!-----------------------------------
+! this OpenMP do-loop cannot be parallelized by the recursion on k/k+1
+!$OMP parallel do default(none) shared(is,ie,js,je,kbot,qa,dp)
+ do j=js, je
+ do k=1,kbot-1
+ do i=is, ie
+ if( qa(i,j,k) < 0. ) then
+ qa(i,j,k+1) = qa(i,j,k+1) + qa(i,j,k)*dp(i,j,k)/dp(i,j,k+1)
+ qa(i,j,k ) = 0.
+ endif
+ enddo
+ enddo
+ enddo
+
+! Bottom layer; Borrow from above
+!$OMP parallel do default(none) shared(is,ie,js,je,qa,kbot,dp) &
+!$OMP private(dq)
+ do j=js, je
+ do i=is, ie
+ if( qa(i,j,kbot) < 0. .and. qa(i,j,kbot-1)>0.) then
+ dq = min(-qa(i,j,kbot)*dp(i,j,kbot), qa(i,j,kbot-1)*dp(i,j,kbot-1))
+ qa(i,j,kbot-1) = qa(i,j,kbot-1) - dq/dp(i,j,kbot-1)
+ qa(i,j,kbot ) = qa(i,j,kbot ) + dq/dp(i,j,kbot )
+ endif
+! if qa is still < 0
+ qa(i,j,kbot) = max(0., qa(i,j,kbot))
+ enddo
+ enddo
+
+ endif
+
+ end subroutine neg_adj4
+
subroutine neg_adj2(is, ie, js, je, ng, kbot, hydrostatic, &
peln, delz, pt, dp, qv, ql, qr, qi, qs, qa, check_negative)
diff --git a/model/fv_update_phys.F90 b/model/fv_update_phys.F90
index f9979f560..31d7649bb 100644
--- a/model/fv_update_phys.F90
+++ b/model/fv_update_phys.F90
@@ -223,6 +223,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
integer i, j, k, m, n, nwat
integer sphum, liq_wat, ice_wat, cld_amt ! GFDL AM physics
integer rainwat, snowwat, graupel ! GFDL Cloud Microphysics
+ integer hailwat ! NSSL Cloud Microphysics
integer w_diff ! w-tracer for PBL diffusion
real:: qstar, dbk, rdg, zvir, p_fac, cv_air, gama_dt, tbad
logical :: bad_range
@@ -265,6 +266,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
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')
if ( .not. hydrostatic ) then
@@ -324,7 +326,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
!$OMP parallel do default(none) &
!$OMP shared(is,ie,js,je,npz,flagstruct,pfull,q_dt,sphum,q,qdiag, &
!$OMP nq,w_diff,dt,nwat,liq_wat,rainwat,ice_wat,snowwat, &
-!$OMP graupel,delp,cld_amt,hydrostatic,pt,t_dt,delz,adj_vmr,&
+!$OMP graupel,hailwat,delp,cld_amt,hydrostatic,pt,t_dt,delz,adj_vmr,&
!$OMP gama_dt,cv_air,ua,u_dt,va,v_dt,isd,ied,jsd,jed, &
#ifdef MULTI_GASES
!$OMP nn, nm, &
@@ -435,7 +437,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
if ( hydrostatic ) then
do j=js,je
call moist_cp(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
- ice_wat, snowwat, graupel, q, qc, cvm, pt(is:ie,j,k) )
+ ice_wat, snowwat, graupel, hailwat, q, qc, cvm, pt(is:ie,j,k) )
do i=is,ie
pt(i,j,k) = pt(i,j,k) + t_dt(i,j,k)*dt*con_cp/cvm(i)
enddo
@@ -445,7 +447,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
! Constant pressure
do j=js,je
call moist_cp(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
- ice_wat, snowwat, graupel, q, qc, cvm, pt(is:ie,j,k) )
+ ice_wat, snowwat, graupel, hailwat, q, qc, cvm, pt(is:ie,j,k) )
do i=is,ie
delz(i,j,k) = delz(i,j,k) / pt(i,j,k)
pt(i,j,k) = pt(i,j,k) + t_dt(i,j,k)*dt*con_cp/cvm(i)
@@ -463,7 +465,7 @@ subroutine fv_update_phys ( dt, is, ie, js, je, isd, ied, jsd, jed, ng, nq,
else
do j=js,je
call moist_cv(is,ie,isd,ied,jsd,jed, npz, j, k, nwat, sphum, liq_wat, rainwat, &
- ice_wat, snowwat, graupel, q, qc, cvm, pt(is:ie,j,k))
+ ice_wat, snowwat, graupel, hailwat, q, qc, cvm, pt(is:ie,j,k))
do i=is,ie
pt(i,j,k) = pt(i,j,k) + t_dt(i,j,k)*dt*con_cp/cvm(i)
enddo
diff --git a/tools/coarse_grained_diagnostics.F90 b/tools/coarse_grained_diagnostics.F90
index 8f910abdd..b0f02e734 100644
--- a/tools/coarse_grained_diagnostics.F90
+++ b/tools/coarse_grained_diagnostics.F90
@@ -228,6 +228,14 @@ subroutine populate_coarse_diag_type(Atm, coarse_diagnostics)
coarse_diagnostics(index)%units = 'kg/kg/s'
coarse_diagnostics(index)%reduction_method = MASS_WEIGHTED
+ index = index + 1
+ coarse_diagnostics(index)%axes = 3
+ coarse_diagnostics(index)%module_name = DYNAMICS
+ coarse_diagnostics(index)%name = 'qh_dt_phys_coarse'
+ coarse_diagnostics(index)%description = 'coarse-grained hail tracer tendency from physics'
+ coarse_diagnostics(index)%units = 'kg/kg/s'
+ coarse_diagnostics(index)%reduction_method = MASS_WEIGHTED
+
index = index + 1
coarse_diagnostics(index)%axes = 3
coarse_diagnostics(index)%module_name = DYNAMICS
@@ -366,6 +374,15 @@ subroutine populate_coarse_diag_type(Atm, coarse_diagnostics)
coarse_diagnostics(index)%vertically_integrated = .true.
coarse_diagnostics(index)%reduction_method = AREA_WEIGHTED
+ index = index + 1
+ coarse_diagnostics(index)%axes = 2
+ coarse_diagnostics(index)%module_name = DYNAMICS
+ coarse_diagnostics(index)%name = 'int_qh_dt_phys_coarse'
+ coarse_diagnostics(index)%description = 'coarse-grained vertically integrated hail tracer tendency from physics'
+ coarse_diagnostics(index)%units = 'kg/m**2/s'
+ coarse_diagnostics(index)%vertically_integrated = .true.
+ coarse_diagnostics(index)%reduction_method = AREA_WEIGHTED
+
index = index + 1
coarse_diagnostics(index)%axes = 2
coarse_diagnostics(index)%module_name = DYNAMICS
@@ -1214,17 +1231,18 @@ subroutine compute_cvm(q, pt, isc, iec, jsc, jec, npz, isd, ied, jsd, jed, nwat,
real, dimension(isd:ied,jsd:jed,1:npz), intent(in) :: pt
real, dimension(isc:iec,jsc:jec,1:npz), intent(out) :: cvm
real, dimension(isc:iec) :: qc, cvm_tmp
- integer :: j, k, sphum, liq_wat, ice_wat, rainwat, snowwat, graupel
+ integer :: j, k, sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat
sphum = get_tracer_index (MODEL_ATMOS, 'sphum')
liq_wat = get_tracer_index (MODEL_ATMOS, 'liq_wat')
ice_wat = get_tracer_index (MODEL_ATMOS, 'ice_wat')
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')
do j = jsc, jec
do k = 1, npz
call moist_cv(isc, iec, isd, ied, jsd, jed, npz, j, k, nwat, sphum, &
- liq_wat, rainwat, ice_wat, snowwat, graupel, &
+ liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, &
q, qc, cvm_tmp, pt(isc:iec,j,k))
cvm(isc:iec,j,k) = cvm_tmp
enddo
@@ -1237,17 +1255,18 @@ subroutine compute_cpm(q, pt, isc, iec, jsc, jec, npz, isd, ied, jsd, jed, nwat,
real, dimension(isd:ied,jsd:jed,1:npz), intent(in) :: pt
real, dimension(isc:iec,jsc:jec,1:npz), intent(out) :: cpm
real, dimension(isc:iec) :: qc, cpm_tmp
- integer :: j, k, sphum, liq_wat, ice_wat, rainwat, snowwat, graupel
+ integer :: j, k, sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat
sphum = get_tracer_index (MODEL_ATMOS, 'sphum')
liq_wat = get_tracer_index (MODEL_ATMOS, 'liq_wat')
ice_wat = get_tracer_index (MODEL_ATMOS, 'ice_wat')
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')
do j = jsc, jec
do k = 1, npz
call moist_cp(isc, iec, isd, ied, jsd, jed, npz, j, k, nwat, sphum, &
- liq_wat, rainwat, ice_wat, snowwat, graupel, &
+ liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, &
q, qc, cpm_tmp, pt(isc:iec,j,k))
cpm(isc:iec,j,k) = cpm_tmp
enddo
diff --git a/tools/external_ic.F90 b/tools/external_ic.F90
index 3fe232bd3..56cd554f7 100644
--- a/tools/external_ic.F90
+++ b/tools/external_ic.F90
@@ -220,7 +220,7 @@ subroutine get_external_ic( Atm, fv_domain, cold_start )
integer :: is, ie, js, je
integer :: isd, ied, jsd, jed, ng
- integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, sgs_tke, cld_amt
+ integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, sgs_tke, cld_amt
integer :: liq_aero, ice_aero
#ifdef MULTI_GASES
integer :: spo, spo2, spo3
@@ -318,6 +318,7 @@ subroutine get_external_ic( Atm, fv_domain, cold_start )
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')
#ifdef MULTI_GASES
spo = get_tracer_index(MODEL_ATMOS, 'spo')
spo2 = get_tracer_index(MODEL_ATMOS, 'spo2')
@@ -340,6 +341,8 @@ subroutine get_external_ic( Atm, fv_domain, cold_start )
call prt_maxmin('snowwat', Atm%q(:,:,:,snowwat), is, ie, js, je, ng, Atm%npz, 1.)
if ( graupel > 0 ) &
call prt_maxmin('graupel', Atm%q(:,:,:,graupel), is, ie, js, je, ng, Atm%npz, 1.)
+ if ( hailwat > 0 ) &
+ call prt_maxmin('hailwat', Atm%q(:,:,:,hailwat), is, ie, js, je, ng, Atm%npz, 1.)
#ifdef MULTI_GASES
if ( spo > 0 ) &
call prt_maxmin('SPO', Atm%q(:,:,:,spo), is, ie, js, je, ng, Atm%npz, 1.)
@@ -489,7 +492,7 @@ subroutine get_nggps_ic (Atm, fv_domain)
real(kind=R_GRID), dimension(2):: p1, p2, p3
real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
integer:: i,j,k,nts, ks, naxis_dims
- integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, tke, ntclamt
+ integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, tke, ntclamt
namelist /external_ic_nml/ filtered_terrain, levp, gfs_dwinds, &
checker_tr, nt_checker
@@ -788,6 +791,7 @@ subroutine get_nggps_ic (Atm, fv_domain)
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')
ntclamt = get_tracer_index(MODEL_ATMOS, 'cld_amt')
if (data_source_fv3gfs) then
do k=1,npz
@@ -937,8 +941,11 @@ subroutine read_gfs_ic()
if (data_source_fv3gfs) call register_restart_field(GFS_restart, 't', temp, dim_names_3d3, is_optional=.true.)
! prognostic tracers
+
do nt = 1, ntracers
- q(:,:,:,nt) = -999.99
+
+ q(:,:,:,nt) = -999.99
+
call get_tracer_names(MODEL_ATMOS, nt, tracer_name)
call register_restart_field(GFS_restart, trim(tracer_name), q(:,:,:,nt), dim_names_3d3, is_optional=.true.)
enddo
@@ -1018,7 +1025,7 @@ subroutine get_hrrr_ic (Atm, fv_domain)
real(kind=R_GRID), dimension(2):: p1, p2, p3
real(kind=R_GRID), dimension(3):: e1, e2, ex, ey
integer:: i,j,k,nts, ks
- integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, tke, ntclamt
+ integer:: liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, tke, ntclamt
namelist /external_ic_nml/ filtered_terrain, levp, gfs_dwinds, &
checker_tr, nt_checker
! variables for reading the dimension from the hrrr_ctrl
@@ -1949,7 +1956,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
logical:: found
integer :: is, ie, js, je
integer :: isd, ied, jsd, jed
- integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, sgs_tke, cld_amt
+ integer :: sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, sgs_tke, cld_amt
#ifdef MULTI_GASES
integer :: spo, spo2, spo3
#else
@@ -1999,6 +2006,7 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
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')
#ifdef MULTI_GASES
spo = get_tracer_index(MODEL_ATMOS, 'spo')
spo2 = get_tracer_index(MODEL_ATMOS, 'spo2')
@@ -2012,11 +2020,14 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
if (is_master()) then
print *, 'sphum = ', sphum
print *, 'liq_wat = ', liq_wat
- if ( Atm%flagstruct%nwat .eq. 6 ) then
+ if ( Atm%flagstruct%nwat .ge. 6 ) then
print *, 'rainwat = ', rainwat
print *, 'iec_wat = ', ice_wat
print *, 'snowwat = ', snowwat
print *, 'graupel = ', graupel
+ IF ( Atm%flagstruct%nwat == 7 ) then
+ print *, 'hailwat = ', hailwat
+ ENDIF
endif
#ifdef MULTI_GASES
print *, ' spo3 = ', spo3
@@ -2363,7 +2374,8 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
enddo
enddo
- qc(:,:,:,graupel) = 0. ! note Graupel must be tracer #6
+ qc(:,:,:,graupel) = 0. ! note Graupel must be tracer #6 (hail assumed not present,
+ ! otherwise qc needs have dimension 7)
deallocate ( qec )
if(is_master()) write(*,*) 'done interpolate tracers (qec) into cubic (qc)'
@@ -2562,12 +2574,12 @@ subroutine get_ecmwf_ic( Atm, fv_domain )
wt = Atm%delp(i,j,k)
if ( Atm%flagstruct%nwat .eq. 2 ) then
qt = wt*(1.+Atm%q(i,j,k,liq_wat))
- elseif ( Atm%flagstruct%nwat .eq. 6 ) then
+ elseif ( Atm%flagstruct%nwat .eq. 6 .or. Atm%flagstruct%nwat .eq. 7 ) then
qt = wt*(1. + Atm%q(i,j,k,liq_wat) + &
Atm%q(i,j,k,ice_wat) + &
Atm%q(i,j,k,rainwat) + &
Atm%q(i,j,k,snowwat) + &
- Atm%q(i,j,k,graupel))
+ Atm%q(i,j,k,graupel)) ! assume hailwat is zero if nwat=7
endif
m_fac = wt / qt
do iq=1,ntracers
@@ -2953,7 +2965,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
real(kind=R_GRID):: pst
!!! High-precision
integer i,j,k,l,m, k2,iq
- integer sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, cld_amt, sgs_tke
+ integer sphum, liq_wat, ice_wat, rainwat, snowwat, graupel, hailwat, cld_amt, sgs_tke
#ifdef MULTI_GASES
integer spo, spo2, spo3
#else
@@ -2972,6 +2984,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
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')
#ifdef MULTI_GASES
spo = get_tracer_index(MODEL_ATMOS, 'spo')
@@ -2988,6 +3001,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
print *, 'nwat = ', Atm%flagstruct%nwat
print *, 'sphum = ', sphum
print *, 'clwmr = ', liq_wat
+ print *, 'liq_wat = ', liq_wat
#ifdef MULTI_GASES
print *, 'spo = ', spo
print *, 'spo2 = ', spo2
@@ -2995,11 +3009,12 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
#else
print *, ' o3mr = ', o3mr
#endif
- if ( Atm%flagstruct%nwat .eq. 6 ) then
+ if ( Atm%flagstruct%nwat .ge. 6 ) then
print *, 'rainwat = ', rainwat
print *, 'ice_wat = ', ice_wat
print *, 'snowwat = ', snowwat
print *, 'graupel = ', graupel
+ IF ( Atm%flagstruct%nwat == 7 ) print *, 'hailwat = ', hailwat
endif
print *, 'sgs_tke = ', sgs_tke
print *, 'cld_amt = ', cld_amt
@@ -3016,7 +3031,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
#endif
!$OMP parallel do default(none) &
-!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,data_source_fv3gfs,&
+!$OMP shared(sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,hailwat,data_source_fv3gfs,&
!$OMP cld_amt,ncnst,npz,is,ie,js,je,km,k2,ak0,bk0,psc,zh,omga,qa,Atm,z500,t_in) &
!$OMP private(l,m,pst,pn,gz,pe0,pn0,pe1,pn1,dp2,qp,qn1,gz_fv)
@@ -3234,7 +3249,7 @@ subroutine remap_scalar(Atm, km, npz, ncnst, ak0, bk0, psc, qa, zh, omga, t_in)
endif
#endif
- if (Atm%flagstruct%nwat .eq. 6 ) then
+ if (Atm%flagstruct%nwat .eq. 6) then ! no need to check for nwat=7 (hail) since only nwat=3,6 treated here
Atm%q(i,j,k,rainwat) = 0.
Atm%q(i,j,k,snowwat) = 0.
Atm%q(i,j,k,graupel) = 0.
@@ -4368,3 +4383,4 @@ end subroutine get_staggered_grid
end module external_ic_mod
+
diff --git a/tools/fv_diagnostics.F90 b/tools/fv_diagnostics.F90
index 35d59d4b6..7c6d89420 100644
--- a/tools/fv_diagnostics.F90
+++ b/tools/fv_diagnostics.F90
@@ -172,7 +172,7 @@ module fv_diagnostics_mod
logical :: prt_minmax =.false.
logical :: m_calendar
integer sphum, liq_wat, ice_wat, cld_amt ! GFDL physics
- integer rainwat, snowwat, graupel
+ integer rainwat, snowwat, graupel, hailwat
#ifdef MULTI_GASES
integer spo, spo2, spo3
#else
@@ -284,6 +284,7 @@ subroutine fv_diag_init(Atm, axes, Time, npx, npy, npz, p_ref)
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')
#ifdef MULTI_GASES
spo = get_tracer_index (MODEL_ATMOS, 'spo')
spo2 = get_tracer_index (MODEL_ATMOS, 'spo2')
@@ -1691,7 +1692,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
call nh_total_energy(isc, iec, jsc, jec, isd, ied, jsd, jed, npz, &
Atm(n)%w, Atm(n)%delz, Atm(n)%pt, Atm(n)%delp, &
Atm(n)%q, Atm(n)%phis, Atm(n)%gridstruct%area, Atm(n)%domain, &
- sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, Atm(n)%flagstruct%nwat, &
+ sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat, Atm(n)%flagstruct%nwat, &
Atm(n)%ua, Atm(n)%va, Atm(n)%flagstruct%moist_phys, a2)
#endif
call prt_maxmin('UA_top', Atm(n)%ua(isc:iec,jsc:jec,1), &
@@ -2725,6 +2726,17 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
enddo
enddo
endif
+ if (hailwat > 0) then
+ do k=1,npz
+ do j=jsc,jec
+ do i=isc,iec
+ a2(i,j) = a2(i,j) + Atm(n)%delp(i,j,k) * &
+ Atm(n)%q(i,j,k,hailwat)
+ enddo
+ enddo
+ enddo
+ endif
+
used = send_data(id_iw, a2*ginv, Time)
endif
if ( id_lw>0 ) then
@@ -2833,7 +2845,8 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
endif
! Cloud top temperature & cloud top press:
- if ( (id_ctt>0 .or. id_ctp>0 .or. id_ctz>0).and. Atm(n)%flagstruct%nwat==6) then
+ if ( (id_ctt>0 .or. id_ctp>0 .or. id_ctz>0) &
+ .and. (Atm(n)%flagstruct%nwat==6 .or. Atm(n)%flagstruct%nwat==7)) then
allocate ( var1(isc:iec,jsc:jec) )
allocate ( var2(isc:iec,jsc:jec) )
!$OMP parallel do default(shared) private(tmp)
@@ -2842,6 +2855,9 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
do k=2,npz
tmp = atm(n)%q(i,j,k,liq_wat)+atm(n)%q(i,j,k,rainwat)+atm(n)%q(i,j,k,ice_wat)+ &
atm(n)%q(i,j,k,snowwat)+atm(n)%q(i,j,k,graupel)
+ IF ( Atm(n)%flagstruct%nwat==7) THEN
+ tmp = tmp + atm(n)%q(i,j,k,hailwat)
+ ENDIF
if( tmp>5.e-6 ) then
a2(i,j) = Atm(n)%pt(i,j,k)
var1(i,j) = 0.01*Atm(n)%pe(i,k,j)
@@ -2955,6 +2971,16 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
enddo
enddo
endif
+ if (hailwat > 0) then
+!$OMP parallel do default(shared)
+ do k=1,npz
+ do j=jsc,jec
+ do i=isc,iec
+ wk(i,j,k) = wk(i,j,k) + Atm(n)%q(i,j,k,hailwat)*Atm(n)%delp(i,j,k)
+ enddo
+ enddo
+ enddo
+ endif
used = send_data(id_qp, wk, Time)
endif
@@ -3005,7 +3031,7 @@ subroutine fv_diag(Atm, zvir, Time, print_freq)
do j=jsc,jec
#ifdef USE_COND
call moist_cv(isc,iec,isd,ied,jsd,jed,npz,j,k,Atm(n)%flagstruct%nwat,sphum,liq_wat,rainwat, &
- ice_wat,snowwat,graupel,Atm(n)%q,Atm(n)%q_con(isc:iec,j,k),cvm)
+ ice_wat,snowwat,graupel,hailwat,Atm(n)%q,Atm(n)%q_con(isc:iec,j,k),cvm)
do i=isc,iec
a3(i,j,k) = Atm(n)%pt(i,j,k)*cvm(i)*wk(i,j,k)
enddo
@@ -4375,6 +4401,7 @@ subroutine prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain
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')
if ( nwat==0 ) then
psmo = g_sum(domain, ps(is:ie,js:je), is, ie, js, je, n_g, area, 1)
@@ -4400,6 +4427,8 @@ subroutine prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain
call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,snowwat), psq(is,js,snowwat))
if (graupel > 0) &
call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,graupel), psq(is,js,graupel))
+ if (hailwat > 0) &
+ call z_sum(is, ie, js, je, km, n_g, delp, q(is-n_g,js-n_g,1,hailwat), psq(is,js,hailwat))
! Mean water vapor in the "stratosphere" (75 mb and above):
@@ -4443,6 +4472,9 @@ subroutine prt_mass(km, nq, is, ie, js, je, n_g, nwat, ps, delp, q, area, domain
write(*,*) 'Total snow ', trim(gn), '=', qtot(snowwat)*ginv
if (graupel > 0) &
write(*,*) 'Total graupel ', trim(gn), '=', qtot(graupel)*ginv
+ if (hailwat > 0) &
+ write(*,*) 'Total hailwat ', trim(gn), '=', qtot(hailwat)*ginv
+
write(*,*) '---------------------------------------------'
elseif ( nwat==2 ) then
write(*,*) 'GFS condensate (kg/m^2)', trim(gn), '=', qtot(liq_wat)*ginv
@@ -5774,13 +5806,13 @@ end subroutine eqv_pot
subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
w, delz, pt, delp, q, hs, area, domain, &
sphum, liq_wat, rainwat, ice_wat, &
- snowwat, graupel, nwat, ua, va, moist_phys, te)
+ snowwat, graupel, hailwat, nwat, ua, va, moist_phys, te)
!------------------------------------------------------
! Compute vertically integrated total energy per column
!------------------------------------------------------
! !INPUT PARAMETERS:
integer, intent(in):: km, is, ie, js, je, isd, ied, jsd, jed
- integer, intent(in):: nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel
+ integer, intent(in):: nwat, sphum, liq_wat, rainwat, ice_wat, snowwat, graupel, hailwat
real, intent(in), dimension(isd:ied,jsd:jed,km):: ua, va, pt, delp, w
real, intent(in), dimension(is:ie,js:je,km) :: delz
real, intent(in), dimension(isd:ied,jsd:jed,km,nwat):: q
@@ -5804,7 +5836,7 @@ subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
#ifdef MULTI_GASES
!$OMP num_gas, &
#endif
-!$OMP w,q,pt,delp,delz,hs,cv_air,moist_phys,sphum,liq_wat,rainwat,ice_wat,snowwat,graupel) &
+!$OMP w,q,pt,delp,delz,hs,cv_air,moist_phys,sphum,liq_wat,rainwat,ice_wat,snowwat,graupel,hailwat) &
!$OMP private(phiz,cvm, qc)
do j=js,je
@@ -5822,7 +5854,7 @@ subroutine nh_total_energy(is, ie, js, je, isd, ied, jsd, jed, km, &
if ( moist_phys ) then
do k=1,km
call moist_cv(is,ie,isd,ied,jsd,jed, km, j, k, nwat, sphum, liq_wat, rainwat, &
- ice_wat, snowwat, graupel, q, qc, cvm)
+ ice_wat, snowwat, graupel, hailwat, q, qc, cvm)
do i=is,ie
te(i,j) = te(i,j) + delp(i,j,k)*( cvm(i)*pt(i,j,k) + hlv*q(i,j,k,sphum) + &
0.5*(phiz(i,k)+phiz(i,k+1)+ua(i,j,k)**2+va(i,j,k)**2+w(i,j,k)**2) )