Support for cloud microphysics hail species (NOAA-EMC#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 <larissa.reames@noaa.gov>

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

model/fv_dynamics.F90

@@ -95,6 +95,10 @@ module fv_dynamics_mod
 !     <td>neg_adj2</td>
 !   </tr>
 !   <tr>
+!     <td>fv_sg_mod</td>
+!     <td>neg_adj4</td>
+!   </tr>
+!   <tr>
 !     <td>fv_timing_mod</td>
 !     <td>timing_on, timing_off</td>
 !   </tr>
@@ -115,6 +119,10 @@ module fv_dynamics_mod
 !     <td>neg_adj3</td>
 !   </tr>
 !   <tr>
+!     <td>fv_sg_mod</td>
+!     <td>neg_adj4</td>
+!   </tr>
+!   <tr>
 !     <td>tracer_manager_mod</td>
 !     <td>get_tracer_index</td>
 !   </tr>
@@ -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,         &

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