From 43f7ed39fcd302e8404a152011f7a02c5d76ddc9 Mon Sep 17 00:00:00 2001 From: Ted Mansell <37668594+MicroTed@users.noreply.github.com> Date: Fri, 4 Mar 2022 10:37:02 -0600 Subject: [PATCH] Support for cloud microphysics hail species (#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 <52886575+LarissaReames-NOAA@users.noreply.github.com> Co-authored-by: LarissaReames-NOAA --- model/fv_arrays.F90 | 2 +- model/fv_dynamics.F90 | 66 +++- model/fv_mapz.F90 | 68 +++- model/fv_nesting.F90 | 37 ++- model/fv_regional_bc.F90 | 13 +- model/fv_sg.F90 | 446 ++++++++++++++++++++++++++- model/fv_update_phys.F90 | 10 +- tools/coarse_grained_diagnostics.F90 | 27 +- tools/external_ic.F90 | 42 ++- tools/fv_diagnostics.F90 | 48 ++- 10 files changed, 676 insertions(+), 83 deletions(-) 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) )