zm_diff.patch

diff --git a/components/cam/bld/build-namelist b/components/cam/bld/build-namelist
index fb28d0411e..c6f2e47187 100755
--- a/components/cam/bld/build-namelist
+++ b/components/cam/bld/build-namelist
@@ -3366,6 +3366,7 @@ add_default($nl,'fix_g1_err_ndrop','val'=>'.false.');
 
 # Deep convection scheme
 add_default($nl, 'deep_scheme');
+add_default($nl, 'deep_num_steps');
 
 # Radiation scheme
 add_default($nl, 'radiation_scheme');
diff --git a/components/cam/bld/namelist_files/namelist_defaults_cam.xml b/components/cam/bld/namelist_files/namelist_defaults_cam.xml
index d2282e9ff3..8611b71f79 100644
--- a/components/cam/bld/namelist_files/namelist_defaults_cam.xml
+++ b/components/cam/bld/namelist_files/namelist_defaults_cam.xml
@@ -1016,6 +1016,8 @@
 <deep_scheme clubb_do_deep="1" >CLUBB_SGS</deep_scheme>
 <deep_scheme unicon="1"        >UNICON   </deep_scheme>
 
+<deep_num_steps>1</deep_num_steps>
+
 <!-- PBL scheme -->
 <eddy_scheme pbl="uw"             >diag_TKE </eddy_scheme>
 <eddy_scheme pbl="hb"             >HB       </eddy_scheme>
diff --git a/components/cam/bld/namelist_files/namelist_definition.xml b/components/cam/bld/namelist_files/namelist_definition.xml
index 12029f66be..d0adb19a9d 100644
--- a/components/cam/bld/namelist_files/namelist_definition.xml
+++ b/components/cam/bld/namelist_files/namelist_definition.xml
@@ -2076,6 +2076,14 @@ This default logical is set in cospsimulator_intr.F90
 Default: FALSE
 </entry>
 
+<!-- Deep convection -->
+
+<entry id="deep_num_steps" type="integer" category="conv"
+       group="phys_ctl_nl" valid_values="" >
+Number of deep convection substeps.
+Default: 1
+</entry>
+
 <!-- Cloud Macro/Micro-physics -->
 
 <entry id="cld_macmic_num_steps" type="integer" category="conv"
diff --git a/components/cam/src/physics/cam/phys_control.F90 b/components/cam/src/physics/cam/phys_control.F90
index bb61a6d152..ef24bc161c 100644
--- a/components/cam/src/physics/cam/phys_control.F90
+++ b/components/cam/src/physics/cam/phys_control.F90
@@ -112,6 +112,8 @@ logical, public, protected :: print_fixer_message  = .false.     ! switch on err
 integer, public, protected :: ieflx_opt = 0
 logical, public, protected :: l_ieflx_fix = .false.
 
+! Deep convection substeps
+integer           :: deep_num_steps = 1
 ! Macro/micro-physics co-substeps
 integer           :: cld_macmic_num_steps = 1
 
@@ -184,7 +186,7 @@ subroutine phys_ctl_readnl(nlfile)
       use_qqflx_fixer, & 
       print_fixer_message, & 
       use_hetfrz_classnuc, use_gw_oro, use_gw_front, use_gw_convect, &
-      cld_macmic_num_steps, micro_do_icesupersat, &
+      deep_num_steps, cld_macmic_num_steps, micro_do_icesupersat, &
       fix_g1_err_ndrop, ssalt_tuning, resus_fix, convproc_do_aer, &
       convproc_do_gas, convproc_method_activate, liqcf_fix, regen_fix, demott_ice_nuc, pergro_mods, pergro_test_active, &
       mam_amicphys_optaa, n_so4_monolayers_pcage,micro_mg_accre_enhan_fac, &
@@ -270,6 +272,7 @@ subroutine phys_ctl_readnl(nlfile)
    call mpibcast(l_st_mac,                        1 , mpilog,  0, mpicom)
    call mpibcast(l_st_mic,                        1 , mpilog,  0, mpicom)
    call mpibcast(l_rad,                           1 , mpilog,  0, mpicom)
+   call mpibcast(deep_num_steps,                  1 , mpiint,  0, mpicom)
    call mpibcast(cld_macmic_num_steps,            1 , mpiint,  0, mpicom)
    call mpibcast(prc_coef1,                       1 , mpir8,   0, mpicom)
    call mpibcast(prc_exp,                         1 , mpir8,   0, mpicom)
@@ -338,6 +341,14 @@ subroutine phys_ctl_readnl(nlfile)
       endif
    endif
 
+   ! Deep convection substepping support.
+   if (deep_num_steps > 1) then
+      if (deep_scheme /= "ZM") then
+         call endrun ("Setting deep_num_steps > 1 is only supported with the &
+              &ZM deep convection scheme.")
+      end if
+   end if
+
    ! Macro/micro co-substepping support.
    if (cld_macmic_num_steps > 1) then
       if (microp_scheme /= "MG" .or. (macrop_scheme /= "park" .and. macrop_scheme /= "CLUBB_SGS")) then
@@ -415,7 +426,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
                         use_mass_borrower_out, & 
                         l_ieflx_fix_out, & 
                         use_qqflx_fixer_out, & 
-                        print_fixer_message_out, & 
+                        print_fixer_message_out, deep_num_steps_out, &
                         cld_macmic_num_steps_out, micro_do_icesupersat_out, &
                         fix_g1_err_ndrop_out, ssalt_tuning_out,resus_fix_out,convproc_do_aer_out,  &
                         convproc_do_gas_out, convproc_method_activate_out, mam_amicphys_optaa_out, n_so4_monolayers_pcage_out, &
@@ -492,6 +503,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
    logical,           intent(out), optional :: l_rad_out
    logical,           intent(out), optional :: mg_prc_coeff_fix_out
    logical,           intent(out), optional :: rrtmg_temp_fix_out
+   integer,           intent(out), optional :: deep_num_steps_out
    integer,           intent(out), optional :: cld_macmic_num_steps_out
    real(r8),          intent(out), optional :: prc_coef1_out
    real(r8),          intent(out), optional :: prc_exp_out
@@ -555,6 +567,7 @@ subroutine phys_getopts(deep_scheme_out, shallow_scheme_out, eddy_scheme_out, &
    if ( present(l_st_mac_out            ) ) l_st_mac_out          = l_st_mac
    if ( present(l_st_mic_out            ) ) l_st_mic_out          = l_st_mic
    if ( present(l_rad_out               ) ) l_rad_out             = l_rad
+   if ( present(deep_num_steps_out      ) ) deep_num_steps_out    = deep_num_steps
    if ( present(cld_macmic_num_steps_out) ) cld_macmic_num_steps_out = cld_macmic_num_steps
    if ( present(prc_coef1_out           ) ) prc_coef1_out            = prc_coef1
    if ( present(prc_exp_out             ) ) prc_exp_out              = prc_exp
diff --git a/components/cam/src/physics/cam/physpkg.F90 b/components/cam/src/physics/cam/physpkg.F90
index 0e2baa9157..a686915216 100644
--- a/components/cam/src/physics/cam/physpkg.F90
+++ b/components/cam/src/physics/cam/physpkg.F90
@@ -66,6 +66,9 @@ module physpkg
   integer ::  prec_sh_idx        = 0
   integer ::  snow_sh_idx        = 0
   integer ::  rice2_idx          = 0
+  integer ::  rprddp_idx         = 0
+  integer ::  nevapr_dpcu_idx    = 0
+  integer ::  ttend_dp_idx       = 0
   integer :: species_class(pcnst)  = -1 !BSINGH: Moved from modal_aero_data.F90 as it is being used in second call to zm deep convection scheme (convect_deep_tend_2)
 
   save
@@ -83,6 +86,7 @@ module physpkg
   character(len=16) :: shallow_scheme
   character(len=16) :: macrop_scheme
   character(len=16) :: microp_scheme 
+  integer           :: deep_num_steps          ! Number of deep convection substeps
   integer           :: cld_macmic_num_steps    ! Number of macro/micro substeps
   logical           :: do_clubb_sgs
   logical           :: use_subcol_microp   ! if true, use subcolumns in microphysics
@@ -164,6 +168,7 @@ subroutine phys_register
     call phys_getopts(shallow_scheme_out       = shallow_scheme, &
                       macrop_scheme_out        = macrop_scheme,   &
                       microp_scheme_out        = microp_scheme,   &
+                      deep_num_steps_out       = deep_num_steps, &
                       cld_macmic_num_steps_out = cld_macmic_num_steps, &
                       do_clubb_sgs_out         = do_clubb_sgs,     &
                       do_aerocom_ind3_out      = do_aerocom_ind3,  &
@@ -880,6 +885,10 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out )
     prec_sh_idx  = pbuf_get_index('PREC_SH')
     snow_sh_idx  = pbuf_get_index('SNOW_SH')
 
+    rprddp_idx  = pbuf_get_index('RPRDDP')
+    nevapr_dpcu_idx  = pbuf_get_index('NEVAPR_DPCU')
+    ttend_dp_idx  = pbuf_get_index('TTEND_DP')
+
     if (shallow_scheme .eq. 'UNICON') then
         rice2_idx    = pbuf_get_index('rice2')
     endif
@@ -1851,6 +1860,8 @@ subroutine tphysbc (ztodt,               &
     use subcol,          only: subcol_gen, subcol_ptend_avg
     use subcol_utils,    only: subcol_ptend_copy, is_subcol_on
     use phys_control,    only: use_qqflx_fixer, use_mass_borrower
+    use zm_average_gathered, only: zm_gath_avg, zm_gath_avg_init, &
+                                   zm_gath_avg_accum, zm_gath_avg_output
 
     implicit none
 
@@ -1906,6 +1917,9 @@ subroutine tphysbc (ztodt,               &
 
     integer  i,k,m,ihist                       ! Longitude, level, constituent indices
     integer :: ixcldice, ixcldliq              ! constituent indices for cloud liquid and ice water.
+    ! for ZM substepping
+    integer :: deep_it                         ! iteration variables
+    real(r8) :: deep_ztodt                     ! modified timestep
     ! for macro/micro co-substepping
     integer :: macmic_it                       ! iteration variables
     real(r8) :: cld_macmic_ztodt               ! modified timestep
@@ -1939,6 +1953,11 @@ subroutine tphysbc (ztodt,               &
 
     real(r8),pointer :: rice2(:)                  ! reserved ice from UNICON [m/s]
 
+    ! Variables from deep convection used by wet deposition or gravity waves
+    real(r8),pointer :: rprddp(:,:)
+    real(r8),pointer :: nevapr_dpcu(:,:)
+    real(r8),pointer :: ttend_dp(:,:)
+
     ! carma precipitation variables
     real(r8) :: prec_sed_carma(pcols)          ! total precip from cloud sedimentation (CARMA)
     real(r8) :: snow_sed_carma(pcols)          ! snow from cloud ice sedimentation (CARMA)
@@ -1961,6 +1980,21 @@ subroutine tphysbc (ztodt,               &
     real(r8) :: prec_sed_macmic(pcols)
     real(r8) :: snow_sed_macmic(pcols)
 
+    real(r8) :: cmfmc_accum(pcols, pverp)
+    real(r8) :: cmfcme_accum(pcols, pver)
+    real(r8) :: dlf_accum(pcols, pver)
+    real(r8) :: pflx_accum(pcols, pverp)
+    real(r8) :: zdu_accum(pcols, pver)
+
+    real(r8) :: prec_dp_accum(pcols)
+    real(r8) :: snow_dp_accum(pcols)
+    real(r8) :: rprddp_accum(pcols,pver)
+    real(r8) :: nevapr_dpcu_accum(pcols,pver)
+    real(r8) :: ttend_dp_accum(pcols,pver)
+    real(r8) :: rliq_accum(pcols)
+
+    type(zm_gath_avg) :: zga
+
     ! energy checking variables
     real(r8) :: zero(pcols)                    ! array of zeros
     real(r8) :: zero_sc(pcols*psubcols)        ! array of zeros
@@ -2279,21 +2313,6 @@ end if
     ! Moist convection
     !===================================================
     call t_startf('moist_convection')
-    !
-    ! Since the PBL doesn't pass constituent perturbations, they
-    ! are zeroed here for input to the moist convection routine
-    !
-    call t_startf ('convect_deep_tend')
-    call convect_deep_tend(  &
-         cmfmc,      cmfcme,             &
-         dlf,        pflx,    zdu,       &
-         rliq,    &
-         ztodt,   &
-         state,   ptend, cam_in%landfrac, pbuf, mu, eu, du, md, ed, dp,   &
-         dsubcld, jt, maxg, ideep, lengath) 
-    call t_stopf('convect_deep_tend')
-
-    call physics_update(state, ptend, ztodt, tend)
 
     call pbuf_get_field(pbuf, prec_dp_idx, prec_dp )
     call pbuf_get_field(pbuf, snow_dp_idx, snow_dp )
@@ -2307,13 +2326,88 @@ end if
     call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw )
 
     if (use_subcol_microp) then
-      call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol)
-      call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol)
+       call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol)
+       call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol)
     end if
 
-    ! Check energy integrals, including "reserved liquid"
-    flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol)
-    call check_energy_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero)
+    call pbuf_get_field(pbuf, rprddp_idx, rprddp)
+    call pbuf_get_field(pbuf, nevapr_dpcu_idx, nevapr_dpcu)
+    call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)
+
+    cmfmc_accum = 0._r8
+    cmfcme_accum = 0._r8
+    dlf_accum = 0._r8
+    pflx_accum = 0._r8
+    zdu_accum = 0._r8
+
+    prec_dp_accum = 0._r8
+    snow_dp_accum = 0._r8
+    rprddp_accum = 0._r8
+    nevapr_dpcu_accum = 0._r8
+    ttend_dp_accum = 0._r8
+    rliq_accum = 0._r8
+
+    deep_ztodt = ztodt / deep_num_steps
+    call zm_gath_avg_init(state%pdel, zga)
+
+    do deep_it = 1, deep_num_steps
+       !
+       ! Since the PBL doesn't pass constituent perturbations, they
+       ! are zeroed here for input to the moist convection routine
+       !
+       call t_startf ('convect_deep_tend')
+       call convect_deep_tend(  &
+            cmfmc,      cmfcme,             &
+            dlf,        pflx,    zdu,       &
+            rliq,    &
+            deep_ztodt,   &
+            state,   ptend, cam_in%landfrac, pbuf, mu, eu, du, md, ed, dp,   &
+            dsubcld, jt, maxg, ideep, lengath) 
+       call t_stopf('convect_deep_tend')
+
+       call physics_ptend_scale(ptend, 1._r8/deep_num_steps, ncol)
+       call physics_update(state, ptend, ztodt, tend)
+
+       ! Check energy integrals, including "reserved liquid"
+       flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol)
+       call check_energy_chng(state, tend, "convect_deep", nstep, ztodt, zero, &
+            flx_cnd/deep_num_steps, snow_dp/deep_num_steps, zero)
+
+       call zm_gath_avg_accum( &
+            lengath, ideep, mu, eu, du, &
+            md, ed, dp, dsubcld, jt, &
+            maxg, zga)
+
+       cmfmc_accum(:ncol,:) = cmfmc_accum(:ncol,:) + cmfmc(:ncol,:)
+       cmfcme_accum(:ncol,:) = cmfcme_accum(:ncol,:) + cmfcme(:ncol,:)
+       dlf_accum(:ncol,:) = dlf_accum(:ncol,:) + dlf(:ncol,:)
+       pflx_accum(:ncol,:) = pflx_accum(:ncol,:) + pflx(:ncol,:)
+       zdu_accum(:ncol,:) = zdu_accum(:ncol,:) + zdu(:ncol,:)
+
+       prec_dp_accum(:ncol) = prec_dp_accum(:ncol) + prec_dp(:ncol)
+       snow_dp_accum(:ncol) = snow_dp_accum(:ncol) + snow_dp(:ncol)
+       rprddp_accum(:ncol,:) = rprddp_accum(:ncol,:) + rprddp(:ncol,:)
+       nevapr_dpcu_accum(:ncol,:) = nevapr_dpcu_accum(:ncol,:) + nevapr_dpcu(:ncol,:)
+       ttend_dp_accum(:ncol,:) = ttend_dp_accum(:ncol,:) + ttend_dp(:ncol,:)
+       rliq_accum(:ncol) = rliq_accum(:ncol) + rliq(:ncol)
+    end do
+    call zm_gath_avg_output( &
+         zga, lengath, ideep, mu, eu, &
+         du, md, ed, dp, dsubcld, &
+         jt, maxg)
+
+    cmfmc(:ncol,:) = cmfmc_accum(:ncol,:) / deep_num_steps
+    cmfcme(:ncol,:) = cmfcme_accum(:ncol,:) / deep_num_steps
+    dlf(:ncol,:) = dlf_accum(:ncol,:) / deep_num_steps
+    pflx(:ncol,:) = pflx_accum(:ncol,:) / deep_num_steps
+    zdu(:ncol,:) = zdu_accum(:ncol,:) / deep_num_steps
+
+    prec_dp(:ncol) = prec_dp_accum(:ncol) / deep_num_steps
+    snow_dp(:ncol) = snow_dp_accum(:ncol) / deep_num_steps
+    rprddp(:ncol,:) = rprddp_accum(:ncol,:) / deep_num_steps
+    nevapr_dpcu(:ncol,:) = nevapr_dpcu_accum(:ncol,:) / deep_num_steps
+    ttend_dp(:ncol,:) = ttend_dp_accum(:ncol,:) / deep_num_steps
+    rliq(:ncol) = rliq_accum(:ncol) / deep_num_steps
 
     !
     ! Call Hack (1994) convection scheme to deal with shallow/mid-level convection
diff --git a/components/cam/src/physics/cam/zm_average_gathered.F90 b/components/cam/src/physics/cam/zm_average_gathered.F90
new file mode 100644
index 0000000000..b9a7427cbc
--- /dev/null
+++ b/components/cam/src/physics/cam/zm_average_gathered.F90
@@ -0,0 +1,162 @@
+module zm_average_gathered
+
+  use shr_kind_mod, only: r8 => shr_kind_r8
+#ifndef TEST_MODE
+  use ppgrid, only: pcols, pver, pverp
+  use shr_assert_mod, only: shr_assert, shr_assert_in_domain
+#endif
+
+  implicit none
+  private
+
+  public zm_gath_avg
+  public zm_gath_avg_init
+  public zm_gath_avg_accum
+  public zm_gath_avg_output
+
+#ifdef TEST_MODE
+  ! For testing
+  integer, parameter, public :: pcols = 8
+  integer, parameter, public :: pver = 8
+  integer, parameter, public :: pverp = pver + 1
+#endif
+
+  type zm_gath_avg
+     ! How often has each column had convection?
+     integer :: con_count(pcols)
+     ! How many substeps have occurred?
+     integer :: step_count
+     ! Un-gathered fields to average and gather
+     real(r8) :: mu(pcols, pver)
+     real(r8) :: eu(pcols, pver)
+     real(r8) :: du(pcols, pver)
+     real(r8) :: md(pcols, pver)
+     real(r8) :: ed(pcols, pver)
+     ! layer thickness
+     real(r8) :: dp(pcols, pver)
+     ! top and bottom of convective region
+     integer :: jt(pcols)
+     integer :: maxg(pcols)
+  end type zm_gath_avg
+
+contains
+
+  subroutine zm_gath_avg_init(pdel, zga)
+    real(r8), intent(in) :: pdel(pcols, pver)
+    type(zm_gath_avg), intent(out) :: zga
+
+    zga%con_count = 0
+    zga%step_count = 0
+
+    zga%mu = 0.
+    zga%eu = 0.
+    zga%du = 0.
+    zga%md = 0.
+    zga%ed = 0.
+
+    zga%dp = 0.01 * pdel
+
+    zga%jt = pver
+    zga%maxg = 1
+  end subroutine zm_gath_avg_init
+
+  subroutine zm_gath_avg_accum( &
+       lengath, ideep, mu, eu, du, &
+       md, ed, dp, dsubcld, jt, &
+       maxg, zga)
+    integer, intent(in) :: lengath
+    integer, intent(in) :: ideep(pcols)
+    real(r8), intent(in) :: mu(pcols, pver)
+    real(r8), intent(in) :: eu(pcols, pver)
+    real(r8), intent(in) :: du(pcols, pver)
+    real(r8), intent(in) :: md(pcols, pver)
+    real(r8), intent(in) :: ed(pcols, pver)
+    real(r8), intent(in) :: dp(pcols, pver)
+    real(r8), intent(in) :: dsubcld(pcols)
+    integer, intent(in) :: jt(pcols)
+    integer, intent(in) :: maxg(pcols)
+
+    type(zm_gath_avg), intent(inout) :: zga
+
+    integer :: i, k
+
+    zga%step_count = zga%step_count + 1
+
+    do i = 1, lengath
+       zga%con_count(ideep(i)) = zga%con_count(ideep(i)) + 1
+
+       zga%mu(ideep(i),:) = zga%mu(ideep(i),:) + mu(i,:)
+       zga%eu(ideep(i),:) = zga%eu(ideep(i),:) + eu(i,:)
+       zga%du(ideep(i),:) = zga%du(ideep(i),:) + du(i,:)
+       zga%md(ideep(i),:) = zga%md(ideep(i),:) + md(i,:)
+       zga%ed(ideep(i),:) = zga%ed(ideep(i),:) + ed(i,:)
+
+#ifdef TEST_MODE
+       if (any(zga%dp(ideep(i),:) /= dp(i,:))) then
+          stop 1
+       end if
+       if (dsubcld(i) /= sum(zga%dp(ideep(i),maxg(i):pver))) then
+          stop 2
+       end if
+#else
+       ! Just to test that I have correctly traced this during development.
+       ! Do not put this in production.
+       do k = 1, pver
+          call shr_assert_in_domain(abs(zga%dp(ideep(i), k) - dp(i,k)), &
+               lt=1.e-2_r8, varname="dp error")
+       end do
+       call shr_assert_in_domain(abs(dsubcld(i) - sum(zga%dp(ideep(i),maxg(i):pver))), &
+            lt=1.e-2_r8, varname="dsubcld error")
+#endif
+
+       zga%jt(ideep(i)) = min(zga%jt(ideep(i)), jt(i))
+       zga%maxg(ideep(i)) = max(zga%maxg(ideep(i)), maxg(i))
+    end do
+
+  end subroutine zm_gath_avg_accum
+
+  subroutine zm_gath_avg_output( &
+       zga, lengath, ideep, mu, eu, &
+       du, md, ed, dp, dsubcld, &
+       jt, maxg)
+    type(zm_gath_avg), intent(in) :: zga
+
+    integer, intent(out) :: lengath
+    integer, intent(out) :: ideep(pcols)
+    real(r8), intent(out) :: mu(pcols, pver)
+    real(r8), intent(out) :: eu(pcols, pver)
+    real(r8), intent(out) :: du(pcols, pver)
+    real(r8), intent(out) :: md(pcols, pver)
+    real(r8), intent(out) :: ed(pcols, pver)
+    real(r8), intent(out) :: dp(pcols, pver)
+    real(r8), intent(out) :: dsubcld(pcols)
+    integer, intent(out) :: jt(pcols)
+    integer, intent(out) :: maxg(pcols)
+
+    integer :: i
+
+    lengath = 0
+    do i = 1, pcols
+       if (zga%con_count(i) > 0) then
+          lengath = lengath + 1
+          ideep(lengath) = i
+       end if
+    end do
+
+    do i = 1, lengath
+       mu(i,:) = zga%mu(ideep(i),:) / zga%step_count
+       eu(i,:) = zga%eu(ideep(i),:) / zga%step_count
+       du(i,:) = zga%du(ideep(i),:) / zga%step_count
+       md(i,:) = zga%md(ideep(i),:) / zga%step_count
+       ed(i,:) = zga%ed(ideep(i),:) / zga%step_count
+
+       jt(i) = zga%jt(ideep(i))
+       maxg(i) = zga%maxg(ideep(i))
+
+       dp(i,:) = zga%dp(ideep(i),:)
+       dsubcld(i) = sum(dp(i,maxg(i):pver))
+    end do
+
+  end subroutine zm_gath_avg_output
+
+end module zm_average_gathered
diff --git a/components/cam/test/zm_average_gathered_test.F90 b/components/cam/test/zm_average_gathered_test.F90
new file mode 100644
index 0000000000..651b0ca492
--- /dev/null
+++ b/components/cam/test/zm_average_gathered_test.F90
@@ -0,0 +1,241 @@
+program zm_average_gathered_test
+
+  use shr_kind_mod, only: r8 => shr_kind_r8
+
+  use zm_average_gathered, only: &
+       zm_gath_avg, &
+       zm_gath_avg_init, &
+       zm_gath_avg_accum, &
+       zm_gath_avg_output, &
+       pcols, &
+       pver, &
+       pverp
+
+  implicit none
+
+  real(r8) :: pdel(pcols, pver)
+
+  real(r8) :: mu(pcols, pver)
+  real(r8) :: eu(pcols, pver)
+  real(r8) :: du(pcols, pver)
+  real(r8) :: md(pcols, pver)
+  real(r8) :: ed(pcols, pver)
+  real(r8) :: dp(pcols, pver)
+  real(r8) :: dsubcld(pcols)
+  integer :: jt(pcols)
+  integer :: maxg(pcols)
+  integer :: ideep(pcols)
+  integer :: lengath
+
+  type(zm_gath_avg) :: zga
+  
+  real(r8) :: con_count(pcols, pver)
+
+  integer :: i, k
+
+  do k = 1, pver
+     do i = 1, pcols
+        pdel(i,k) = k + i*0.5
+     end do
+  end do
+
+  call zm_gath_avg_init(pdel, zga)
+
+  lengath = 3
+  ideep(1) = 2
+  ideep(2) = 4
+  ideep(3) = 5
+
+  jt(1) = 3
+  maxg(1) = 6
+
+  mu(1, :) = 2.
+  eu(1, :) = 2. - 1.
+  du(1, :) = 2. - 2.
+  md(1, :) = 2. - 3.
+  ed(1, :) = 2. - 4.
+
+  dp(1, :) = 0.01 * pdel(2, :)
+
+  jt(2) = 3
+  maxg(2) = 6
+
+  mu(2, :) = 1.
+  eu(2, :) = 1. - 1.
+  du(2, :) = 1. - 2.
+  md(2, :) = 1. - 3.
+  ed(2, :) = 1. - 4.
+
+  dp(2, :) = 0.01 * pdel(4, :)
+
+  jt(3) = 3
+  maxg(3) = 6
+
+  mu(3, :) = 3.
+  eu(3, :) = 3. - 1.
+  du(3, :) = 3. - 2.
+  md(3, :) = 3. - 3.
+  ed(3, :) = 3. - 4.
+
+  dp(3, :) = 0.01 * pdel(5, :)
+
+  dsubcld = 0.
+  do k = 1, pver
+     do i = 1, lengath
+        if (k >= maxg(i)) then
+           dsubcld(i) = dsubcld(i) + dp(i,k)
+        end if
+     end do
+  end do
+
+  call zm_gath_avg_accum( &
+       lengath, ideep, mu, eu, du, &
+       md, ed, dp, dsubcld, jt, &
+       maxg, zga)
+
+  lengath = 3
+  ideep(1) = 4
+  ideep(2) = 5
+  ideep(3) = 7
+
+  jt(1) = 4
+  maxg(1) = 5
+
+  mu(1, :) = 2.
+  eu(1, :) = 2. - 1.
+  du(1, :) = 2. - 2.
+  md(1, :) = 2. - 3.
+  ed(1, :) = 2. - 4.
+
+  dp(1, :) = 0.01 * pdel(4, :)
+
+  jt(2) = 2
+  maxg(2) = 7
+
+  mu(2, :) = 1.
+  eu(2, :) = 1. - 1.
+  du(2, :) = 1. - 2.
+  md(2, :) = 1. - 3.
+  ed(2, :) = 1. - 4.
+
+  dp(2, :) = 0.01 * pdel(5, :)
+
+  jt(3) = 3
+  maxg(3) = 6
+
+  mu(3, :) = 4.
+  eu(3, :) = 4. - 1.
+  du(3, :) = 4. - 2.
+  md(3, :) = 4. - 3.
+  ed(3, :) = 4. - 4.
+
+  dp(3, :) = 0.01 * pdel(7, :)
+
+  dsubcld = 0.
+  do k = 1, pver
+     do i = 1, lengath
+        if (k >= maxg(i)) then
+           dsubcld(i) = dsubcld(i) + dp(i,k)
+        end if
+     end do
+  end do
+
+  call zm_gath_avg_accum( &
+       lengath, ideep, mu, eu, du, &
+       md, ed, dp, dsubcld, jt, &
+       maxg, zga)
+
+  dp = 0.
+
+  call zm_gath_avg_output( &
+       zga, lengath, ideep, mu, eu, &
+       du, md, ed, dp, dsubcld, &
+       jt, maxg)
+
+  if (lengath /= 4) then
+     print *, "Wrong lengath: ", lengath
+  end if
+
+  if (any(ideep(1:4) /= [2, 4, 5, 7])) then
+     print *, "Wrong ideep: ", ideep
+  end if
+
+  if (any(jt(1:4) /= [3, 3, 2, 3])) then
+     print *, "Wrong jt: ", jt
+  end if
+
+  if (any(maxg(1:4) /= [6, 6, 7, 6])) then
+     print *, "Wrong maxg: ", maxg
+  end if
+
+  if (any(mu(1, :) /= 1.)) then
+     print *, "Wrong mu1: ", mu(1, :)
+  end if
+
+  if (any(mu(2, :) /= 1.5)) then
+     print *, "Wrong mu2: ", mu(2, :)
+  end if
+
+  if (any(mu(3, :) /= 2.)) then
+     print *, "Wrong mu3: ", mu(3, :)
+  end if
+
+  if (any(mu(4, :) /= 2.)) then
+     print *, "Wrong mu4: ", mu(4, :)
+  end if
+
+  con_count = 0.
+  con_count(1,:) = 0.5
+  con_count(2,:) = 1
+  con_count(3,:) = 1.
+  con_count(4,:) = 0.5
+
+  if (any(mu(1:4,:) /= eu(1:4,:) + 1.*con_count(1:4,:))) then
+     print *, "Wrong eu: ", mu(1:4,:), eu(1:4,:) + 1.*con_count(1:4,:)
+  end if
+
+  if (any(mu(1:4,:) /= du(1:4,:) + 2.*con_count(1:4,:))) then
+     print *, "Wrong du: ", mu(1:4,:), du(1:4,:)
+  end if
+
+  if (any(mu(1:4,:) /= md(1:4,:) + 3.*con_count(1:4,:))) then
+     print *, "Wrong md: ", mu(1:4,:), md(1:4,:)
+  end if
+
+  if (any(mu(1:4,:) /= ed(1:4,:) + 4.*con_count(1:4,:))) then
+     print *, "Wrong ed: ", mu(1:4,:), ed(1:4,:)
+  end if
+
+  if (any(dp(1,:) /= 0.01 * pdel(2,:))) then
+     print *, "Wrong dp1: ", dp(1,:)
+  end if
+
+  if (any(dp(2,:) /= 0.01 * pdel(4,:))) then
+     print *, "Wrong dp2: ", dp(2,:)
+  end if
+
+  if (any(dp(3,:) /= 0.01 * pdel(5,:))) then
+     print *, "Wrong dp3: ", dp(3,:)
+  end if
+
+  if (any(dp(4,:) /= 0.01 * pdel(7,:))) then
+     print *, "Wrong dp4: ", dp(4,:)
+  end if
+
+  if (dsubcld(1) /= 0.01 * pdel(2, 6) + 0.01 * pdel(2, 7) + 0.01 * pdel(2, 8) ) then
+     print *, "Wrong dsubcld1: ", dsubcld(1)
+  end if
+
+  if (dsubcld(2) /= 0.01 * pdel(4, 6) + 0.01 * pdel(4, 7) + 0.01 * pdel(4, 8) ) then
+     print *, "Wrong dsubcld2: ", dsubcld(2)
+  end if
+
+  if (dsubcld(3) /= 0.01 * pdel(5, 7) + 0.01 * pdel(5, 8) ) then
+     print *, "Wrong dsubcld3: ", dsubcld(3)
+  end if
+
+  if (dsubcld(4) /= 0.01 * pdel(7, 6) + 0.01 * pdel(7, 7) + 0.01 * pdel(7, 8) ) then
+     print *, "Wrong dsubcld4: ", dsubcld(4)
+  end if
+
+end program zm_average_gathered_test