diff --git a/ccpp/data/GFS_typedefs.F90 b/ccpp/data/GFS_typedefs.F90
index 84302cac4..7007f3e6e 100644
--- a/ccpp/data/GFS_typedefs.F90
+++ b/ccpp/data/GFS_typedefs.F90
@@ -45,7 +45,7 @@ module GFS_typedefs
       integer :: ntrcaer
 
       ! If these are changed to >99, need to adjust formatting string in GFS_diagnostics.F90 (and names in diag_tables)
-      integer, parameter :: naux2dmax = 20 !< maximum number of auxiliary 2d arrays in output (for debugging)
+      integer, parameter :: naux2dmax = 30 !< maximum number of auxiliary 2d arrays in output (for debugging)
       integer, parameter :: naux3dmax = 20 !< maximum number of auxiliary 3d arrays in output (for debugging)
 
       integer, parameter :: dfi_radar_max_intervals = 4 !< Number of radar-derived temperature tendency and/or convection suppression intervals. Do not change.
@@ -390,6 +390,8 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: tsnow_ice(:)       => null()  !< RUC LSM: snow temperature at the bottom of the first snow layer over ice
     real (kind=kind_phys), pointer :: snowfallac_land(:) => null()  !< ruc lsm diagnostics over land
     real (kind=kind_phys), pointer :: snowfallac_ice(:)  => null()  !< ruc lsm diagnostics over ice
+    real (kind=kind_phys), pointer :: acsnow_land(:)     => null()  !< ruc lsm diagnostics over land
+    real (kind=kind_phys), pointer :: acsnow_ice(:)      => null()  !< ruc lsm diagnostics over ice
 
     !  MYNN surface layer
     real (kind=kind_phys), pointer :: ustm (:)         => null()  !u* including drag
@@ -1116,6 +1118,10 @@ module GFS_typedefs
 
 !--- potential temperature definition in surface layer physics
     logical              :: thsfc_loc       !< flag for local vs. standard potential temperature
+!--- flux method in 2-m diagnostics
+    logical              :: diag_flux       !< flag for flux method in 2-m diagnostics
+!--- log method in 2-m diagnostics (for stable conditions)
+    logical              :: diag_log        !< flag for log method in 2-m diagnostics (for stable conditions)
 
 !--- vertical diffusion
     real(kind=kind_phys) :: xkzm_m          !< [in] bkgd_vdif_m  background vertical diffusion for momentum
@@ -1582,15 +1588,21 @@ module GFS_typedefs
                                                !        %upfx0    - clear sky upward lw flux at toa (w/m**2)
 
 ! Input/output - used by physics
-    real (kind=kind_phys), pointer :: srunoff(:)     => null()   !< surface water runoff (from lsm)
-    real (kind=kind_phys), pointer :: evbsa  (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: evcwa  (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: snohfa (:)     => null()   !< noah lsm diagnostics
+    real (kind=kind_phys), pointer :: srunoff(:)     => null()   !< accumulated surface storm runoff (from lsm)
+    real (kind=kind_phys), pointer :: evbsa  (:)     => null()   !< accumulated direct evaporation
+    real (kind=kind_phys), pointer :: evcwa  (:)     => null()   !< accumulated canopy evaporation
+    real (kind=kind_phys), pointer :: snohfa (:)     => null()   !< heat flux for phase change of snow (melting)
+    real (kind=kind_phys), pointer :: transa (:)     => null()   !< accumulated transpiration
+    real (kind=kind_phys), pointer :: sbsnoa (:)     => null()   !< accumulated snow sublimation
+    real (kind=kind_phys), pointer :: snowca (:)     => null()   !< snow cover
+    real (kind=kind_phys), pointer :: sbsno  (:)     => null()   !< instantaneous snow sublimation
+    real (kind=kind_phys), pointer :: evbs(:)        => null()   !< instantaneous direct evaporation
+    real (kind=kind_phys), pointer :: trans  (:)     => null()   !< instantaneous transpiration
+    real (kind=kind_phys), pointer :: evcw(:)        => null()   !< instantaneous canopy evaporation
+    real (kind=kind_phys), pointer :: snowmt_land(:) => null()   !< ruc lsm diagnostics over land
+    real (kind=kind_phys), pointer :: snowmt_ice(:)  => null()   !< ruc lsm diagnostics over ice
+    real (kind=kind_phys), pointer :: soilm  (:)     => null()   !< integrated soil moisture
     real (kind=kind_phys), pointer :: paha   (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: transa (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: sbsnoa (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: snowca (:)     => null()   !< noah lsm diagnostics
-    real (kind=kind_phys), pointer :: soilm  (:)     => null()   !< soil moisture
     real (kind=kind_phys), pointer :: tmpmin (:)     => null()   !< min temperature at 2m height (k)
     real (kind=kind_phys), pointer :: tmpmax (:)     => null()   !< max temperature at 2m height (k)
     real (kind=kind_phys), pointer :: dusfc  (:)     => null()   !< u component of surface stress
@@ -1928,8 +1940,6 @@ module GFS_typedefs
     real (kind=kind_phys), pointer      :: evap_ice(:)        => null()  !<
     real (kind=kind_phys), pointer      :: evap_land(:)       => null()  !<
     real (kind=kind_phys), pointer      :: evap_water(:)      => null()  !<
-    real (kind=kind_phys), pointer      :: evbs(:)            => null()  !<
-    real (kind=kind_phys), pointer      :: evcw(:)            => null()  !<
     real (kind=kind_phys), pointer      :: pah(:)             => null()  !<
     real (kind=kind_phys), pointer      :: ecan(:)            => null()  !<
     real (kind=kind_phys), pointer      :: etran(:)           => null()  !<
@@ -2072,7 +2082,6 @@ module GFS_typedefs
     real (kind=kind_phys), pointer      :: save_tcp(:,:)      => null()  !<
     real (kind=kind_phys), pointer      :: save_u(:,:)        => null()  !<
     real (kind=kind_phys), pointer      :: save_v(:,:)        => null()  !<
-    real (kind=kind_phys), pointer      :: sbsno(:)           => null()  !<
     type (cmpfsw_type),    pointer      :: scmpsw(:)          => null()  !<
     real (kind=kind_phys), pointer      :: sfcalb(:,:)        => null()  !<
     real (kind=kind_phys), pointer      :: sigma(:)           => null()  !<
@@ -2099,7 +2108,6 @@ module GFS_typedefs
     integer                             :: tracers_total                 !<
     integer                             :: tracers_water                 !<
     logical                             :: trans_aero                    !<
-    real (kind=kind_phys), pointer      :: trans(:)           => null()  !<
     real (kind=kind_phys), pointer      :: tseal(:)           => null()  !<
     real (kind=kind_phys), pointer      :: tsfa(:)            => null()  !<
     real (kind=kind_phys), pointer      :: tsfc_water(:)      => null()  !<
@@ -2698,6 +2706,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
        allocate (Sfcprop%tsnow_ice       (IM))
        allocate (Sfcprop%snowfallac_land (IM))
        allocate (Sfcprop%snowfallac_ice  (IM))
+       allocate (Sfcprop%acsnow_land     (IM))
+       allocate (Sfcprop%acsnow_ice      (IM))
        !
        Sfcprop%wetness         = clear_val
        Sfcprop%sh2o            = clear_val
@@ -2714,6 +2724,8 @@ subroutine sfcprop_create (Sfcprop, IM, Model)
        Sfcprop%tsnow_ice       = clear_val
        Sfcprop%snowfallac_land = clear_val
        Sfcprop%snowfallac_ice  = clear_val
+       Sfcprop%acsnow_land     = clear_val
+       Sfcprop%acsnow_ice      = clear_val
        !
        if (Model%rdlai) then
           allocate (Sfcprop%xlaixy (IM))
@@ -3517,6 +3529,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
 !--- potential temperature definition in surface layer physics
     logical              :: thsfc_loc      = .true.          !< flag for local vs. standard potential temperature
+!--- flux method in 2-m diagnostics
+    logical              :: diag_flux      = .false.         !< flag for flux method in 2-m diagnostics
+!--- flux method in 2-m diagnostics (for stable conditions) 
+    logical              :: diag_log       = .false.         !< flag for log method in 2-m diagnostics (for stable conditions)
                                                              !<.true. means use local (gridpoint) surface pressure to define potential temperature
                                                              !<       this is the current GFS physics approach
                                                              !<.false. means use reference pressure of 1000 hPa to define potential temperature
@@ -3711,6 +3727,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                                sfc_z0_type,                                                 &
                           !--- switch beteeen local and standard potential temperature
                                thsfc_loc,                                                   &
+                          !--- switches in 2-m diagnostics
+                               diag_flux, diag_log,                                         &
                           !    vertical diffusion
                                xkzm_m, xkzm_h, xkzm_s, xkzminv, moninq_fac, dspfac,         &
                                bl_upfr, bl_dnfr, rlmx, elmx, sfc_rlm,                       &
@@ -4449,6 +4467,17 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
 
 !--- potential temperature reference in sfc layer
     Model%thsfc_loc        = thsfc_loc
+    if (Model%do_mynnsfclay) then
+!--- flux method in 2-m diagnostics
+      Model%diag_flux        = .true.
+!--- flux method in 2-m diagnostics (for stable conditions)
+      Model%diag_log         = .false.
+    else
+!--- flux method in 2-m diagnostics
+      Model%diag_flux        = diag_flux
+!--- flux method in 2-m diagnostics (for stable conditions)
+      Model%diag_log         = diag_log
+    endif
 
 !--- vertical diffusion
     Model%xkzm_m           = xkzm_m
@@ -5848,6 +5877,8 @@ subroutine control_print(Model)
       print *, ' rbcr              : ', Model%rbcr
       print *, ' do_mynnedmf       : ', Model%do_mynnedmf
       print *, ' do_mynnsfclay     : ', Model%do_mynnsfclay
+      print *, ' diag_flux         : ', Model%diag_flux
+      print *, ' diag_log          : ', Model%diag_log
       print *, ' do_myjsfc         : ', Model%do_myjsfc
       print *, ' do_myjpbl         : ', Model%do_myjpbl
       print *, ' do_ugwp           : ', Model%do_ugwp
@@ -6617,11 +6648,17 @@ subroutine diag_create (Diag, IM, Model)
 !--- In/Out
     allocate (Diag%srunoff (IM))
     allocate (Diag%evbsa   (IM))
-    allocate (Diag%evcwa   (IM))
     allocate (Diag%snohfa  (IM))
+    allocate (Diag%evcwa   (IM))
     allocate (Diag%transa  (IM))
     allocate (Diag%sbsnoa  (IM))
     allocate (Diag%snowca  (IM))
+    allocate (Diag%evbs    (IM))
+    allocate (Diag%evcw    (IM))
+    allocate (Diag%sbsno   (IM))
+    allocate (Diag%trans   (IM))
+    allocate (Diag%snowmt_land (IM))
+    allocate (Diag%snowmt_ice  (IM))
     allocate (Diag%soilm   (IM))
     allocate (Diag%tmpmin  (IM))
     allocate (Diag%tmpmax  (IM))
@@ -6915,8 +6952,14 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
     Diag%evcwa      = zero
     Diag%snohfa     = zero
     Diag%transa     = zero
-    Diag%sbsnoa     = zero
     Diag%snowca     = zero
+    Diag%sbsnoa     = zero
+    Diag%sbsno      = zero
+    Diag%evbs       = zero
+    Diag%evcw       = zero
+    Diag%trans      = zero
+    Diag%snowmt_land= zero
+    Diag%snowmt_ice = zero 
     Diag%soilm      = zero
     Diag%tmpmin     = Model%huge
     Diag%tmpmax     = zero
@@ -7225,8 +7268,6 @@ subroutine interstitial_create (Interstitial, IM, Model)
     allocate (Interstitial%evap_ice        (IM))
     allocate (Interstitial%evap_land       (IM))
     allocate (Interstitial%evap_water      (IM))
-    allocate (Interstitial%evbs            (IM))
-    allocate (Interstitial%evcw            (IM))
     allocate (Interstitial%pah             (IM))
     allocate (Interstitial%ecan            (IM))
     allocate (Interstitial%etran           (IM))
@@ -7319,7 +7360,6 @@ subroutine interstitial_create (Interstitial, IM, Model)
     allocate (Interstitial%save_tcp        (IM,Model%levs))
     allocate (Interstitial%save_u          (IM,Model%levs))
     allocate (Interstitial%save_v          (IM,Model%levs))
-    allocate (Interstitial%sbsno           (IM))
     allocate (Interstitial%scmpsw          (IM))
     allocate (Interstitial%sfcalb          (IM,NF_ALBD))
     allocate (Interstitial%sigma           (IM))
@@ -7339,7 +7379,6 @@ subroutine interstitial_create (Interstitial, IM, Model)
     allocate (Interstitial%tprcp_ice       (IM))
     allocate (Interstitial%tprcp_land      (IM))
     allocate (Interstitial%tprcp_water     (IM))
-    allocate (Interstitial%trans           (IM))
     allocate (Interstitial%tseal           (IM))
     allocate (Interstitial%tsfa            (IM))
     allocate (Interstitial%tsfc_water      (IM))
@@ -7933,8 +7972,6 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     Interstitial%evap_ice        = Model%huge
     Interstitial%evap_land       = Model%huge
     Interstitial%evap_water      = Model%huge
-    Interstitial%evbs            = clear_val
-    Interstitial%evcw            = clear_val
     Interstitial%pah             = clear_val
     Interstitial%ecan            = clear_val
     Interstitial%etran           = clear_val
@@ -8013,7 +8050,6 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     Interstitial%save_tcp        = clear_val
     Interstitial%save_u          = clear_val
     Interstitial%save_v          = clear_val
-    Interstitial%sbsno           = clear_val
     Interstitial%sigma           = clear_val
     Interstitial%sigmaf          = clear_val
     Interstitial%sigmafrac       = clear_val
@@ -8029,7 +8065,6 @@ subroutine interstitial_phys_reset (Interstitial, Model)
     Interstitial%tprcp_ice       = Model%huge
     Interstitial%tprcp_land      = Model%huge
     Interstitial%tprcp_water     = Model%huge
-    Interstitial%trans           = clear_val
     Interstitial%tseal           = clear_val
     Interstitial%tsfc_water      = Model%huge
     Interstitial%tsurf_ice       = Model%huge
diff --git a/ccpp/data/GFS_typedefs.meta b/ccpp/data/GFS_typedefs.meta
index 4cec6d322..4d1f122b7 100644
--- a/ccpp/data/GFS_typedefs.meta
+++ b/ccpp/data/GFS_typedefs.meta
@@ -1496,21 +1496,35 @@
   kind = kind_phys
   active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme)
 [snowfallac_land]
-  standard_name = surface_snow_amount_over_land
-  long_name = run-total snow accumulation on the ground
+  standard_name = surface_snow_amount_vardens_over_land
+  long_name = run-total snow accumulation on the ground with variable snow density over land
   units = kg m-2
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
   active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme)
+[acsnow_land]
+  standard_name = surface_snow_lwe_thickness_amount_over_land
+  long_name = run-total snowfall water equivalent over land
+  units = kg m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
 [snowfallac_ice]
-  standard_name = surface_snow_amount_over_ice
-  long_name = run-total snow accumulation on the ice
+  standard_name = surface_snow_amount_vardens_over_ice
+  long_name = run-total snow accumulation on the ground with variable snow density over ice
   units = kg m-2
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
   active = (control_for_land_surface_scheme == identifier_for_ruc_land_surface_scheme)
+[acsnow_ice]
+  standard_name = surface_snow_lwe_thickness_amount_over_ice
+  long_name = run-total snowfall water equivalent over ice
+  units = kg m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
 [ustm]
   standard_name = surface_friction_velocity_for_momentum
   long_name = friction velocity isolated for momentum only
@@ -4050,6 +4064,18 @@
   units = flag
   dimensions = ()
   type = logical
+[diag_flux]
+  standard_name = flag_for_flux_method_in_2m_diagnostics
+  long_name = flag for flux method in 2-m diagnostics
+  units = flag
+  dimensions = ()
+  type = logical
+[diag_log]
+  standard_name = flag_for_log_method_in_2m_diagnostics
+  long_name = flag for log method in 2-m diagnostics
+  units = flag
+  dimensions = ()
+  type = logical
 [hybedmf]
   standard_name = flag_for_hybrid_edmf_pbl_scheme
   long_name = flag for hybrid edmf pbl scheme (moninedmf)
@@ -6629,6 +6655,34 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
+[sbsno]
+  standard_name = snow_deposition_sublimation_upward_latent_heat_flux
+  long_name = latent heat flux from snow depo/subl
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[evbs]
+  standard_name = soil_upward_latent_heat_flux
+  long_name = soil upward latent heat flux
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[evcw]
+  standard_name = canopy_upward_latent_heat_flux
+  long_name = canopy upward latent heat flux
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[trans]
+  standard_name = transpiration_flux
+  long_name = total plant transpiration rate
+  units = W m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
 [soilm]
   standard_name = soil_moisture_content
   long_name = soil moisture
@@ -6636,6 +6690,20 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
+[snowmt_land]
+  standard_name = surface_snow_melt_over_land
+  long_name = snow melt during timestep over land
+  units = kg m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+[snowmt_ice]
+  standard_name = surface_snow_melt_over_ice
+  long_name = snow melt during timestep over ice
+  units = kg m-2
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
 [tmpmin]
   standard_name = minimum_temperature_at_2m
   long_name = min temperature at 2m height
@@ -8505,20 +8573,6 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
-[evbs]
-  standard_name = soil_upward_latent_heat_flux
-  long_name = soil upward latent heat flux
-  units = W m-2
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
-[evcw]
-  standard_name = canopy_upward_latent_heat_flux
-  long_name = canopy upward latent heat flux
-  units = W m-2
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
 [pah]
   standard_name = total_precipitation_advected_heat
   long_name = precipitation advected heat - total
@@ -9654,13 +9708,6 @@
   dimensions = (horizontal_loop_extent,vertical_layer_dimension)
   type = real
   kind = kind_phys
-[sbsno]
-  standard_name = snow_deposition_sublimation_upward_latent_heat_flux
-  long_name = latent heat flux from snow depo/subl
-  units = W m-2
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
 [scmpsw]
   standard_name = components_of_surface_downward_shortwave_fluxes
   long_name = derived type for special components of surface downward shortwave fluxes
@@ -9861,13 +9908,6 @@
   units = flag
   dimensions = ()
   type = logical
-[trans]
-  standard_name = transpiration_flux
-  long_name = total plant transpiration rate
-  units = W m-2
-  dimensions = (horizontal_loop_extent)
-  type = real
-  kind = kind_phys
 [tseal]
   standard_name = surface_skin_temperature_for_nsst
   long_name = ocean surface skin temperature
diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90
index 6e4b62337..f1e3bfffa 100644
--- a/ccpp/driver/GFS_diagnostics.F90
+++ b/ccpp/driver/GFS_diagnostics.F90
@@ -1044,7 +1044,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'ssrun_acc'
-    ExtDiag(idx)%desc = 'surface storm water runoff - GFS lsm'
+    ExtDiag(idx)%desc = 'Accumulated surface storm water runoff'
     ExtDiag(idx)%unit = 'kg/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     allocate (ExtDiag(idx)%data(nblks))
@@ -1055,7 +1055,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'evbs_ave'
-    ExtDiag(idx)%desc = 'Direct Evaporation from Bare Soil - GFS lsm'
+    ExtDiag(idx)%desc = 'Direct Evaporation from Bare Soil'
     ExtDiag(idx)%unit = 'W/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%time_avg = .TRUE.
@@ -1067,7 +1067,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'evcw_ave'
-    ExtDiag(idx)%desc = 'Canopy water evaporation - GFS lsm'
+    ExtDiag(idx)%desc = 'Canopy water evaporation'
     ExtDiag(idx)%unit = 'W/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%time_avg = .TRUE.
@@ -1079,7 +1079,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'snohf'
-    ExtDiag(idx)%desc = 'Snow Phase Change Heat Flux - GFS lsm'
+    ExtDiag(idx)%desc = 'Snow Phase Change Heat Flux'
     ExtDiag(idx)%unit = 'W/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%time_avg = .TRUE.
@@ -1092,7 +1092,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
      idx = idx + 1
      ExtDiag(idx)%axes = 2
      ExtDiag(idx)%name = 'pah_ave'
-     ExtDiag(idx)%desc = ' Total Precipitation Advected Heat - GFS lsm'
+     ExtDiag(idx)%desc = ' Total Precipitation Advected Heat'
      ExtDiag(idx)%unit = 'W/m**2'
      ExtDiag(idx)%mod_name = 'gfs_phys'
      ExtDiag(idx)%time_avg = .TRUE.
@@ -1105,7 +1105,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'trans_ave'
-    ExtDiag(idx)%desc = 'transpiration - GFS lsm'
+    ExtDiag(idx)%desc = 'transpiration'
     ExtDiag(idx)%unit = 'W/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%time_avg = .TRUE.
@@ -1117,7 +1117,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'sbsno_ave'
-    ExtDiag(idx)%desc = 'Sublimation (evaporation from snow) - GFS lsm'
+    ExtDiag(idx)%desc = 'Sublimation (evaporation from snow)'
     ExtDiag(idx)%unit = 'W/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%time_avg = .TRUE.
@@ -1128,15 +1128,17 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 
     idx = idx + 1
     ExtDiag(idx)%axes = 2
-    ExtDiag(idx)%name = 'snowc_ave'
-    ExtDiag(idx)%desc = 'snow cover - GFS lsm'
+    ExtDiag(idx)%name = 'snowc'
+    !ExtDiag(idx)%name = 'snowc_ave'
+    ExtDiag(idx)%desc = 'snow cover '
     ExtDiag(idx)%unit = '%'
     ExtDiag(idx)%mod_name = 'gfs_phys'
-    ExtDiag(idx)%time_avg = .TRUE.
+    !ExtDiag(idx)%time_avg = .TRUE.
     ExtDiag(idx)%cnvfac = cn_100
     allocate (ExtDiag(idx)%data(nblks))
     do nb = 1,nblks
-      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowca(:)
+      ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%sncovr(:)
+      !ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowca(:)
     enddo
 
     idx = idx + 1
@@ -1279,7 +1281,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'dlwsfc'
-    ExtDiag(idx)%desc = 'time accumulated downward lw flux at surface- GFS physics'
+    ExtDiag(idx)%desc = 'time accumulated downward lw flux at surface'
     ExtDiag(idx)%unit = 'w/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%intpl_method = 'bilinear'
@@ -1291,7 +1293,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'ulwsfc'
-    ExtDiag(idx)%desc = 'time accumulated upward lw flux at surface- GFS physics'
+    ExtDiag(idx)%desc = 'time accumulated upward lw flux at surface'
     ExtDiag(idx)%unit = 'w/m**2'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%intpl_method = 'bilinear'
@@ -2093,7 +2095,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'cpofp'
     ExtDiag(idx)%desc = 'Percent frozen precipitation'
-    ExtDiag(idx)%unit = '%'
+    ExtDiag(idx)%unit = 'fraction'
     ExtDiag(idx)%mod_name = 'gfs_phys'
     ExtDiag(idx)%intpl_method = 'bilinear'
     allocate (ExtDiag(idx)%data(nblks))
@@ -2789,7 +2791,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'canopy'
     ExtDiag(idx)%desc = 'canopy water (cnwat in gfs data)'
-    ExtDiag(idx)%unit = 'XXX'
+    ExtDiag(idx)%unit = 'mm'
     ExtDiag(idx)%mod_name = 'gfs_sfc'
     allocate (ExtDiag(idx)%data(nblks))
     do nb = 1,nblks
@@ -2888,7 +2890,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'hice'
     ExtDiag(idx)%desc = 'sea ice thickness (icetk in gfs_data)'
-    ExtDiag(idx)%unit = 'XXX'
+    ExtDiag(idx)%unit = 'm'
     ExtDiag(idx)%mod_name = 'gfs_sfc'
     allocate (ExtDiag(idx)%data(nblks))
     do nb = 1,nblks
@@ -2940,7 +2942,62 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowd(:)
     enddo
 
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'sbsno'
+    ExtDiag(idx)%desc = 'instantaneous sublimation (evaporation from snow)'
+    ExtDiag(idx)%unit = 'W/m**2'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%sbsno(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'evbs'
+    ExtDiag(idx)%desc = 'instantaneous direct evaporation over land'
+    ExtDiag(idx)%unit = 'W m-2'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%evbs(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'evcw'
+    ExtDiag(idx)%desc = 'instantaneous canopy evaporation'
+    ExtDiag(idx)%unit = 'W m-2'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%evcw(:)
+    enddo
+
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'trans'
+    ExtDiag(idx)%desc = 'instantaneous transpiration'
+    ExtDiag(idx)%unit = 'W m-2'
+    ExtDiag(idx)%mod_name = 'gfs_phys'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%trans(:)
+    enddo
+
     if (Model%lsm == Model%lsm_ruc) then
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'rhofr'
+      ExtDiag(idx)%desc = 'density of frozen precipitation'
+      ExtDiag(idx)%unit = 'kg m-3'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%rhofr(:)
+      enddo
+
       idx = idx + 1
       ExtDiag(idx)%axes = 2
       ExtDiag(idx)%name = 'snowfall_acc_land'
@@ -2952,6 +3009,28 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowfallac_land(:)
       enddo
 
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'acsnow_land'
+      ExtDiag(idx)%desc = 'total accumulated SWE of frozen precipitation over land'
+      ExtDiag(idx)%unit = 'kg m-2'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%acsnow_land(:)
+      enddo
+
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'snowmt_land'
+      ExtDiag(idx)%desc = 'accumulated snow melt over land'
+      ExtDiag(idx)%unit = 'kg m-2'
+      ExtDiag(idx)%mod_name = 'gfs_phys'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowmt_land(:)
+      enddo
+
       idx = idx + 1
       ExtDiag(idx)%axes = 2
       ExtDiag(idx)%name = 'snowfall_acc_ice'
@@ -2962,7 +3041,29 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       do nb = 1,nblks
         ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%snowfallac_ice(:)
       enddo
-    endif
+
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'acsnow_ice'
+      ExtDiag(idx)%desc = 'total accumulated SWE of frozen precipitation over ice'
+      ExtDiag(idx)%unit = 'kg m-2'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%acsnow_ice(:)
+      enddo
+
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'snowmt_ice'
+      ExtDiag(idx)%desc = 'accumulated snow melt over ice'
+      ExtDiag(idx)%unit = 'kg m-2'
+      ExtDiag(idx)%mod_name = 'gfs_phys'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snowmt_ice(:)
+      enddo
+    endif ! RUC lsm
 
     idx = idx + 1
     ExtDiag(idx)%axes = 2
@@ -3074,8 +3175,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'tprcp'
-    ExtDiag(idx)%desc = 'total precipitation'
-    ExtDiag(idx)%unit = 'kg/m**2'
+    ExtDiag(idx)%desc = 'total time-step precipitation'
+    ExtDiag(idx)%unit = 'm'
     ExtDiag(idx)%mod_name = 'gfs_sfc'
     allocate (ExtDiag(idx)%data(nblks))
     do nb = 1,nblks
@@ -3151,6 +3252,18 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%vfrac(:)
     enddo
 
+    idx = idx + 1
+    ExtDiag(idx)%axes = 2
+    ExtDiag(idx)%name = 'wetness'
+    ExtDiag(idx)%desc = 'soil moisture availability in top soil layer'
+    ExtDiag(idx)%unit = 'fraction'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    ExtDiag(idx)%cnvfac = cn_100
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%wetness(:)
+    enddo
+
   if (Model%rdlai) then
     idx = idx + 1
     ExtDiag(idx)%axes = 2
diff --git a/ccpp/physics b/ccpp/physics
index 5d17f5507..4c214a1c1 160000
--- a/ccpp/physics
+++ b/ccpp/physics
@@ -1 +1 @@
-Subproject commit 5d17f550716ac21021b298fbf0e27795e36f5ec8
+Subproject commit 4c214a1c12dd292e5e90f1eba43af306a45116e7
diff --git a/ccpp/suites/suite_FV3_GFS_v15_thompson_mynn.xml b/ccpp/suites/suite_FV3_GFS_v15_thompson_mynn.xml
index 882f4ada0..075b95293 100644
--- a/ccpp/suites/suite_FV3_GFS_v15_thompson_mynn.xml
+++ b/ccpp/suites/suite_FV3_GFS_v15_thompson_mynn.xml
@@ -40,7 +40,7 @@
     </subcycle>
     <!-- Surface iteration loop -->
     <subcycle loop="2">
-      <scheme>sfc_diff</scheme>
+      <scheme>mynnsfc_wrapper</scheme>
       <scheme>GFS_surface_loop_control_part1</scheme>
       <scheme>sfc_nst_pre</scheme>
       <scheme>sfc_nst</scheme>