diff --git a/FV3/gfsphysics/GFS_layer/GFS_diagnostics.F90 b/FV3/gfsphysics/GFS_layer/GFS_diagnostics.F90
index 08df54e43..7703e3237 100644
--- a/FV3/gfsphysics/GFS_layer/GFS_diagnostics.F90
+++ b/FV3/gfsphysics/GFS_layer/GFS_diagnostics.F90
@@ -1348,6 +1348,92 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%cnvprcp(:)
     enddo
 
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'shallow_convective_precipitation_rate_average'
+    ExtDiag(idx)%desc = 'averaged surface shallow convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%cnvfac = cn_th
+    ExtDiag(idx)%time_avg = .TRUE.
+    ExtDiag(idx)%time_avg_kind = 'full'
+    ExtDiag(idx)%intpl_method = 'bilinear'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%shallow_cnvprcp(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'shallow_convective_precipitation_rate_bucket_average'
+    ExtDiag(idx)%desc = 'averaged bucket surface shallow convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%cnvfac = cn_th
+    ExtDiag(idx)%time_avg = .TRUE.
+    ExtDiag(idx)%intpl_method = 'bilinear'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%shallow_cnvprcpb(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'shallow_convective_precipitation_rate'
+    ExtDiag(idx)%desc = 'surface shallow convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%shallow_cnvprcp(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'deep_convective_precipitation_rate_average'
+    ExtDiag(idx)%desc = 'averaged surface deep convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%cnvfac = cn_th
+    ExtDiag(idx)%time_avg = .TRUE.
+    ExtDiag(idx)%time_avg_kind = 'full'
+    ExtDiag(idx)%intpl_method = 'bilinear'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%deep_cnvprcp(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'deep_convective_precipitation_rate_bucket_average'
+    ExtDiag(idx)%desc = 'averaged bucket surface deep convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%cnvfac = cn_th
+    ExtDiag(idx)%time_avg = .TRUE.
+    ExtDiag(idx)%intpl_method = 'bilinear'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%deep_cnvprcpb(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'deep_convective_precipitation_rate'
+    ExtDiag(idx)%desc = 'surface deep convective precipitation rate'
+    ExtDiag(idx)%unit = 'kg/m**2/s'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    ExtDiag(idx)%coarse_graining_method = 'area_weighted'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%deep_cnvprcp(:)
+    enddo
+
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'spfhmin2m'
diff --git a/FV3/gfsphysics/GFS_layer/GFS_physics_driver.F90 b/FV3/gfsphysics/GFS_layer/GFS_physics_driver.F90
index b5700764c..f16370055 100644
--- a/FV3/gfsphysics/GFS_layer/GFS_physics_driver.F90
+++ b/FV3/gfsphysics/GFS_layer/GFS_physics_driver.F90
@@ -3829,6 +3829,8 @@ subroutine GFS_physics_driver                                 &
 !
       do i=1,im
         Diag%rainc(i) = frain * rain1(i)
+        Diag%deep_cnvprcp(i) = Diag%deep_cnvprcp(i) + Diag%rainc(i)
+        Diag%deep_cnvprcpb(i) = Diag%deep_cnvprcpb(i) + Diag%rainc(i)
       enddo
 !
       if (Model%lssav) then
@@ -4102,6 +4104,8 @@ subroutine GFS_physics_driver                                 &
 
             do i=1,im
               Diag%rainc(i) = Diag%rainc(i) + frain * rain1(i)
+              Diag%shallow_cnvprcp(i) = Diag%shallow_cnvprcp(i) + frain * rain1(i)
+              Diag%shallow_cnvprcpb(i) = Diag%shallow_cnvprcpb(i) + frain * rain1(i)
             enddo
 ! in shalcnv,  'cnvw' and 'cnvc' are not set to zero
             if (Model%shcnvcw .and. Model%num_p3d == 4 .and. Model%npdf3d == 3) then
@@ -4137,6 +4141,8 @@ subroutine GFS_physics_driver                                 &
 
             do i=1,im
               Diag%rainc(i) = Diag%rainc(i) + frain * rain1(i)
+              Diag%shallow_cnvprcp(i) = Diag%shallow_cnvprcp(i) + frain * rain1(i)
+              Diag%shallow_cnvprcpb(i) = Diag%shallow_cnvprcpb(i) + frain * rain1(i)
             enddo
 ! in  mfshalcnv,  'cnvw' and 'cnvc' are set to zero before computation starts:
             if (Model%shcnvcw .and. Model%num_p3d == 4 .and. Model%npdf3d == 3) then
diff --git a/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 b/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90
index 54e6c10ce..f85b46a1e 100644
--- a/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90
+++ b/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90
@@ -1368,6 +1368,10 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: psmean (:)     => null()   !< surface pressure (kPa)
     real (kind=kind_phys), pointer :: cnvprcp(:)     => null()   !< accumulated convective precipitation (kg/m2)
     real (kind=kind_phys), pointer :: cnvprcpb(:)    => null()   !< accumulated convective precipitation in bucket (kg/m2)
+    real (kind=kind_phys), pointer :: shallow_cnvprcp(:)     => null()   !< accumulated shallow convective precipitation (kg/m2)
+    real (kind=kind_phys), pointer :: shallow_cnvprcpb(:)    => null()   !< accumulated shallow convective precipitation in bucket (kg/m2)
+    real (kind=kind_phys), pointer :: deep_cnvprcp(:)     => null()   !< accumulated deep convective precipitation (kg/m2)
+    real (kind=kind_phys), pointer :: deep_cnvprcpb(:)    => null()   !< accumulated deep convective precipitation in bucket (kg/m2)
     real (kind=kind_phys), pointer :: spfhmin(:)     => null()   !< minimum specific humidity
     real (kind=kind_phys), pointer :: spfhmax(:)     => null()   !< maximum specific humidity
     real (kind=kind_phys), pointer :: u10mmax(:)     => null()   !< maximum u-wind
@@ -5104,6 +5108,10 @@ subroutine diag_create (Diag, IM, Model)
     allocate (Diag%psmean  (IM))
     allocate (Diag%cnvprcp (IM))
     allocate (Diag%cnvprcpb(IM))
+    allocate (Diag%shallow_cnvprcp (IM))
+    allocate (Diag%shallow_cnvprcpb(IM))
+    allocate (Diag%deep_cnvprcp (IM))
+    allocate (Diag%deep_cnvprcpb(IM))
     allocate (Diag%spfhmin (IM))
     allocate (Diag%spfhmax (IM))
     allocate (Diag%u10mmax (IM))
@@ -5474,6 +5482,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     Diag%zmtnblck   = zero
     Diag%totprcpb   = zero
     Diag%cnvprcpb   = zero
+    Diag%shallow_cnvprcpb   = zero
+    Diag%deep_cnvprcpb   = zero
     Diag%toticeb    = zero
     Diag%totsnwb    = zero
     Diag%totgrpb    = zero
@@ -5578,6 +5588,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     if (set_totprcp) then
       Diag%totprcp = zero
       Diag%cnvprcp = zero
+      Diag%shallow_cnvprcp = zero
+      Diag%deep_cnvprcp = zero
       Diag%totice  = zero
       Diag%totsnw  = zero
       Diag%totgrp  = zero