diff --git a/ccpp/driver/GFS_diagnostics.F90 b/ccpp/driver/GFS_diagnostics.F90
index 4c0e7cdf7..01d49593e 100644
--- a/ccpp/driver/GFS_diagnostics.F90
+++ b/ccpp/driver/GFS_diagnostics.F90
@@ -1,11 +1,11 @@
 module GFS_diagnostics
 
 !-----------------------------------------------------------------------
-!    GFS_diagnostics_mod defines a data type and contains the routine 
+!    GFS_diagnostics_mod defines a data type and contains the routine
 !    to populate said type with diagnostics from the GFS physics for
 !    use by the modeling system for output
 !-----------------------------------------------------------------------
- 
+
   use machine,            only: kind_phys
 
   !--- GFS_typedefs ---
@@ -51,7 +51,7 @@ module GFS_diagnostics
   CONTAINS
 
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    
+
  ! Helper function for GFS_externaldiag_populate to handle the massive dtend(:,:,dtidx(:,:)) array
     subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit)
     implicit none
@@ -62,7 +62,7 @@ subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit)
     integer, intent(inout) :: idx
     real(kind=kind_phys), pointer :: dtend(:,:,:) ! Assumption: dtend is null iff all(dtidx <= 1)
     character(len=*), intent(in), optional :: desc, unit
-    
+
     integer :: idtend, nb
 
     idtend = Model%dtidx(itrac,iprocess)
@@ -88,17 +88,17 @@ subroutine add_dtend(Model,ExtDiag,IntDiag,idx,nblks,itrac,iprocess,desc,unit)
        enddo
     endif
   end subroutine add_dtend
-  
-!-------------------------------------------------------------------------      
+
+!-------------------------------------------------------------------------
 !--- GFS_externaldiag_populate ---
-!-------------------------------------------------------------------------      
-!    creates and populates a data type with GFS physics diagnostic 
+!-------------------------------------------------------------------------
+!    creates and populates a data type with GFS physics diagnostic
 !    variables which is then handed off to the IPD for use by the model
-!    infrastructure layer to output as needed.  The data type includes 
-!    names, units, conversion factors, etc.  There is no copying of data, 
-!    but instead pointers are associated to the internal representation 
+!    infrastructure layer to output as needed.  The data type includes
+!    names, units, conversion factors, etc.  There is no copying of data,
+!    but instead pointers are associated to the internal representation
 !    of each individual physics diagnostic.
-!-------------------------------------------------------------------------      
+!-------------------------------------------------------------------------
   subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop, Coupling,  &
                                         Grid, Tbd, Cldprop, Radtend, IntDiag, Init_parm)
 !---------------------------------------------------------------------------------------------!
@@ -158,7 +158,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     ExtDiag(:)%name = ''
     ExtDiag(:)%intpl_method = 'nearest_stod'
 
-    idx = 0 
+    idx = 0
 
     idx = idx + 1
     ExtDiag(idx)%axes = 2
@@ -949,7 +949,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 !
 !--- accumulated diagnostics ---
     do num = 1,NFXR
-      write (xtra,'(I2.2)') num 
+      write (xtra,'(I2.2)') num
       idx = idx + 1
       ExtDiag(idx)%axes = 2
       ExtDiag(idx)%name = 'fluxr_'//trim(xtra)
@@ -965,7 +965,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 !--- the next two appear to be appear to be coupling fields in gloopr
 !--- each has four elements
 !rab    do num = 1,4
-!rab      write (xtra,'(I1)') num 
+!rab      write (xtra,'(I1)') num
 !rab      idx = idx + 1
 !rab      ExtDiag(idx)%axes = 2
 !rab      ExtDiag(idx)%name = 'dswcmp_'//trim(xtra)
@@ -978,7 +978,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 !rab    enddo
 !rab
 !rab    do num = 1,4
-!rab      write (xtra,'(I1)') num 
+!rab      write (xtra,'(I1)') num
 !rab      idx = idx + 1
 !rab      ExtDiag(idx)%axes = 2
 !rab      ExtDiag(idx)%name = 'uswcmp_'//trim(xtra)
@@ -1103,7 +1103,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     do nb = 1,nblks
       ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%snohfa(:)
     enddo
-    
+
     if (Model%lsm == Model%lsm_noahmp) then
      idx = idx + 1
      ExtDiag(idx)%axes = 2
@@ -1383,7 +1383,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%tedir(:)
     enddo
 
-    if (Model%lsm == Model%lsm_noahmp) then   
+    if (Model%lsm == Model%lsm_noahmp) then
      idx = idx + 1
      ExtDiag(idx)%axes = 2
      ExtDiag(idx)%name = 'wa_acc'
@@ -2197,7 +2197,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%smcref2(:)
     enddo
 
-    if (Model%lsm == Model%lsm_noahmp) then   
+    if (Model%lsm == Model%lsm_noahmp) then
      idx = idx + 1
      ExtDiag(idx)%axes = 2
      ExtDiag(idx)%name = 'pahi'
@@ -2468,7 +2468,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_pbl(:,:)
       enddo
     endif
-    
+
     if (Model%do_spp) then
       idx = idx + 1
       ExtDiag(idx)%axes = 3
@@ -2481,7 +2481,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_sfc(:,:)
       enddo
     endif
-    
+
     if (Model%do_spp) then
       idx = idx + 1
       ExtDiag(idx)%axes = 3
@@ -2494,7 +2494,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_mp(:,:)
       enddo
     endif
-    
+
     if (Model%do_spp) then
       idx = idx + 1
       ExtDiag(idx)%axes = 3
@@ -2507,7 +2507,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var3 => Coupling(nb)%spp_wts_gwd(:,:)
       enddo
     endif
-    
+
     if (Model%do_spp) then
       idx = idx + 1
       ExtDiag(idx)%axes = 3
@@ -2677,7 +2677,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       do nb = 1,nblks
         ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%use_lake_model(:)
       enddo
-    
+
       if(Model%iopt_lake==Model%iopt_lake_clm) then
 
         ! Populate the 3D arrays separately since the code is complicated:
@@ -2704,7 +2704,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         do nb = 1,nblks
            ExtDiag(idx)%data(nb)%int2 => Sfcprop(nb)%lake_cannot_freeze(:)
         enddo
-        
+
         idx = idx + 1
         ExtDiag(idx)%axes = 2
         ExtDiag(idx)%name = 'lake_t2m'
@@ -2812,7 +2812,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         do nb = 1,nblks
           ExtDiag(idx)%data(nb)%var2 => Sfcprop(nb)%lake_ht(:)
         enddo
-        
+
       endif
 
   endif
@@ -2909,8 +2909,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     do nb = 1,nblks
       ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%du3dt_pbl(:,:)
     enddo
-!    
-! dv3dt_pbl     
+!
+! dv3dt_pbl
     idx = idx + 1
     ExtDiag(idx)%axes = 3
     ExtDiag(idx)%name = 'dv3dt_pbl_ugwp'
@@ -2921,8 +2921,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     do nb = 1,nblks
       ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%dv3dt_pbl(:,:)
     enddo
-!    
-! dt3dt_pbl     
+!
+! dt3dt_pbl
     idx = idx + 1
     ExtDiag(idx)%axes = 3
     ExtDiag(idx)%name = 'dt3dt_pbl_ugwp'
@@ -2934,8 +2934,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%dt3dt_pbl(:,:)
     enddo
 !
-! uav_ugwp 
-! 
+! uav_ugwp
+!
     idx = idx + 1
     ExtDiag(idx)%axes = 3
     ExtDiag(idx)%name = 'uav_ugwp'
@@ -2947,8 +2947,8 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%uav_ugwp(:,:)
     enddo
 !
-! tav_ugwp 
-! 
+! tav_ugwp
+!
     idx = idx + 1
     ExtDiag(idx)%axes = 3
     ExtDiag(idx)%name = 'tav_ugwp'
@@ -2982,7 +2982,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
       ExtDiag(idx)%data(nb)%var3 => IntDiag(nb)%du3dt_ngw(:,:)
     enddo
 !
-!    
+!
     idx = idx + 1
     ExtDiag(idx)%axes = 3
     ExtDiag(idx)%name = 'du3dt_mtb'
@@ -3454,7 +3454,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
           endif
         enddo
       enddo
-      
+
       if_qdiag3d: if(Model%qdiag3d) then
 
         idx = idx + 1
@@ -3499,7 +3499,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
 
 !rab
 !rab    do num = 1,5+Mdl_parms%pl_coeff
-!rab      write (xtra,'(I1)') num 
+!rab      write (xtra,'(I1)') num
 !rab      idx = idx + 1
 !rab      ExtDiag(idx)%axes = 3
 !rab      ExtDiag(idx)%name = 'dtend_'//trim(xtra)
@@ -3877,7 +3877,7 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
     idx = idx + 1
     ExtDiag(idx)%axes = 2
     ExtDiag(idx)%name = 'scolor'
-    ExtDiag(idx)%desc = 'soil color in integer 1-20' 
+    ExtDiag(idx)%desc = 'soil color in integer 1-20'
     ExtDiag(idx)%unit = 'number'
     ExtDiag(idx)%mod_name = 'gfs_sfc'
     allocate (ExtDiag(idx)%data(nblks))
@@ -4203,6 +4203,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%sh2o(:,num)
       enddo
     enddo
+    idx = idx + 1
+    ExtDiag(idx)%axes = 3
+    ExtDiag(idx)%name = 'soill'
+    ExtDiag(idx)%desc = 'liquid soil moisture'
+    ExtDiag(idx)%unit = 'm**3/m**3'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%sh2o(:,:)
+    enddo
   else
     do num = 1,Model%lsoil_lsm
       write (xtra,'(i1)') num
@@ -4225,6 +4235,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%slc(:,num)
       enddo
     enddo
+    idx = idx + 1
+    ExtDiag(idx)%axes = 3
+    ExtDiag(idx)%name = 'soill'
+    ExtDiag(idx)%desc = 'liquid soil moisture'
+    ExtDiag(idx)%unit = 'm**3/m**3'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%slc(:,:)
+    enddo
   endif
 
   if (Model%lsm == Model%lsm_ruc) then
@@ -4241,6 +4261,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%smois(:,num)
       enddo
     enddo
+    idx = idx + 1
+    ExtDiag(idx)%axes = 3
+    ExtDiag(idx)%name = 'soilw'
+    ExtDiag(idx)%desc = 'volumetric soil moisture'
+    ExtDiag(idx)%unit = 'fraction'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%smois(:,:)
+    enddo
   else
     do num = 1,Model%lsoil_lsm
       write (xtra,'(i1)') num
@@ -4255,6 +4285,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
         ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%smc(:,num)
       enddo
     enddo
+    idx = idx + 1
+    ExtDiag(idx)%axes = 3
+    ExtDiag(idx)%name = 'soilw'
+    ExtDiag(idx)%desc = 'volumetric soil moisture'
+    ExtDiag(idx)%unit = 'fraction'
+    ExtDiag(idx)%mod_name = 'gfs_sfc'
+    allocate (ExtDiag(idx)%data(nblks))
+    do nb = 1,nblks
+      ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%smc(:,:)
+    enddo
   endif
 
     if (Model%lsm == Model%lsm_ruc) then
@@ -4271,6 +4311,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
           ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%tslb(:,num)
         enddo
       enddo
+      idx = idx + 1
+      ExtDiag(idx)%axes = 3
+      ExtDiag(idx)%name = 'soilt'
+      ExtDiag(idx)%desc = 'soil temperature'
+      ExtDiag(idx)%unit = 'K'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%tslb(:,:)
+      enddo
     else
       do num = 1,Model%lsoil_lsm
         write (xtra,'(i1)') num
@@ -4285,6 +4335,16 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
           ExtDiag(idx)%data(nb)%var2 => sfcprop(nb)%stc(:,num)
         enddo
       enddo
+      idx = idx + 1
+      ExtDiag(idx)%axes = 3
+      ExtDiag(idx)%name = 'soilt'
+      ExtDiag(idx)%desc = 'soil temperature'
+      ExtDiag(idx)%unit = 'K'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var3 => sfcprop(nb)%stc(:,:)
+      enddo
     endif
 
 !--------------------------nsst variables
@@ -5396,7 +5456,7 @@ subroutine clm_lake_externaldiag_populate(ExtDiag, Model, Sfcprop, idx, cn_one,
     character(:), allocatable :: fullname
 
     integer :: nk, idx0, iblk
-    
+
     do iblk=1,nblks
       call link_all_levels(Sfcprop(iblk)%lake_snow_z3d, 'lake_snow_z3d', 'lake snow level depth', 'm')
     enddo
@@ -5524,6 +5584,6 @@ function soil_layer_depth(lsm, lsm_ruc, lsm_noah, layer) result(layer_depth)
      !
   end function soil_layer_depth
 
-!-------------------------------------------------------------------------      
+!-------------------------------------------------------------------------
 
 end module GFS_diagnostics
diff --git a/ccpp/framework b/ccpp/framework
index 219f2e9c8..f0b9a18b0 160000
--- a/ccpp/framework
+++ b/ccpp/framework
@@ -1 +1 @@
-Subproject commit 219f2e9c88b7b774becac2bd1453696e105af1c4
+Subproject commit f0b9a18b005d950cb9b0038fbc827b6b37500f43
diff --git a/ccpp/physics b/ccpp/physics
index 6b0599ef6..b2a733836 160000
--- a/ccpp/physics
+++ b/ccpp/physics
@@ -1 +1 @@
-Subproject commit 6b0599ef6c7ea076336b1f073899c5b97e37d584
+Subproject commit b2a733836e4a335a817120baec915fd21ea64f76
diff --git a/cpl/module_block_data.F90 b/cpl/module_block_data.F90
index 7bd0a71b3..b5dc6cc80 100644
--- a/cpl/module_block_data.F90
+++ b/cpl/module_block_data.F90
@@ -365,7 +365,7 @@ subroutine block_array_copy_2d_r8_to_2d_r8(destin_ptr, source_arr, block, block_
         jb = block%index(block_index)%jj(ix)
         i = ib - block%isc + 1
         j = jb - block%jsc + 1
-        destin_ptr(i,j) = factor * source_arr(ib,jb)
+        destin_ptr(i,j) = factor * source_arr(i,j)
       enddo
       localrc = ESMF_SUCCESS
     end if
@@ -440,7 +440,7 @@ subroutine block_array_copy_3d_r8_to_3d_r8(destin_ptr, source_arr, block, block_
           jb = block%index(block_index)%jj(ix)
           i = ib - block%isc + 1
           j = jb - block%jsc + 1
-          destin_ptr(i,j,k) = factor * source_arr(ib,jb,k)
+          destin_ptr(i,j,k) = factor * source_arr(i,j,k)
         enddo
       enddo
       localrc = ESMF_SUCCESS
@@ -523,7 +523,7 @@ subroutine block_array_copy_3dslice_r8_to_3d_r8(destin_ptr, source_arr, slice, b
             jb = block%index(block_index)%jj(ix)
             i = ib - block%isc + 1
             j = jb - block%jsc + 1
-            destin_ptr(i,j,k) = factor * source_arr(ib,jb,k,slice)
+            destin_ptr(i,j,k) = factor * source_arr(i,j,k,slice)
           enddo
         enddo
         localrc = ESMF_SUCCESS
@@ -1004,7 +1004,7 @@ subroutine block_array_copy_2d_r4_to_2d_r8(destin_ptr, source_arr, block, block_
         jb = block%index(block_index)%jj(ix)
         i = ib - block%isc + 1
         j = jb - block%jsc + 1
-        destin_ptr(i,j) = factor * source_arr(ib,jb)
+        destin_ptr(i,j) = factor * source_arr(i,j)
       enddo
       localrc = ESMF_SUCCESS
     end if
@@ -1079,7 +1079,7 @@ subroutine block_array_copy_3d_r4_to_3d_r8(destin_ptr, source_arr, block, block_
           jb = block%index(block_index)%jj(ix)
           i = ib - block%isc + 1
           j = jb - block%jsc + 1
-          destin_ptr(i,j,k) = factor * source_arr(ib,jb,k)
+          destin_ptr(i,j,k) = factor * source_arr(i,j,k)
         enddo
       enddo
       localrc = ESMF_SUCCESS
@@ -1162,7 +1162,7 @@ subroutine block_array_copy_3dslice_r4_to_3d_r8(destin_ptr, source_arr, slice, b
             jb = block%index(block_index)%jj(ix)
             i = ib - block%isc + 1
             j = jb - block%jsc + 1
-            destin_ptr(i,j,k) = factor * source_arr(ib,jb,k,slice)
+            destin_ptr(i,j,k) = factor * source_arr(i,j,k,slice)
           enddo
         enddo
         localrc = ESMF_SUCCESS
diff --git a/io/fv3atm_history_io.F90 b/io/fv3atm_history_io.F90
index 7c73fe296..d7ee4e808 100644
--- a/io/fv3atm_history_io.F90
+++ b/io/fv3atm_history_io.F90
@@ -49,9 +49,10 @@ module fv3atm_history_io_mod
   type history_type
     integer :: tot_diag_idx = 0
 
-    integer :: isco=0,ieco=0,jsco=0,jeco=0,levo=0,num_axes_phys=0
-    integer :: fhzero=0, ncld=0, nsoil=0, imp_physics=0, landsfcmdl=0
+    integer :: isco=0,ieco=0,jsco=0,jeco=0,num_axes_phys=0
+    integer :: fhzero=0, ncld=0, nsoil=0, nsoil_lsm=0, imp_physics=0, landsfcmdl=0
     real(4) :: dtp=0
+    integer,dimension(:),        pointer         :: levo => null()
     integer,dimension(:),        pointer         :: nstt => null()
     integer,dimension(:),        pointer         :: nstt_vctbl => null()
     integer,dimension(:),        pointer         :: all_axes => null()
@@ -182,11 +183,11 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat,
     hist%ieco   = Atm_block%iec
     hist%jsco   = Atm_block%jsc
     hist%jeco   = Atm_block%jec
-    hist%levo   = model%levs
     hist%fhzero = nint(Model%fhzero)
     !   hist%ncld   = Model%ncld
     hist%ncld   = Model%imp_physics
     hist%nsoil  = Model%lsoil
+    hist%nsoil_lsm  = Model%lsoil_lsm
     hist%dtp    = Model%dtp
     hist%imp_physics  = Model%imp_physics
     hist%landsfcmdl  = Model%lsm
@@ -208,7 +209,9 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat,
       call mpp_error(fatal, 'fv3atm_io::fv3atm_diag_register - need to increase parameter DIAG_SIZE')
     endif
 
+    allocate(hist%levo(hist%tot_diag_idx))
     allocate(hist%nstt(hist%tot_diag_idx), hist%nstt_vctbl(hist%tot_diag_idx))
+    hist%levo          = 0
     hist%nstt          = 0
     hist%nstt_vctbl    = 0
     nrgst_bl      = 0
@@ -224,6 +227,7 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat,
            trim(Diag(idx)%unit), missing_value=real(missing_value))
       if(Diag(idx)%id > 0) then
         if (Diag(idx)%axes == 2) then
+          hist%levo(idx) = 1
           if( index(trim(Diag(idx)%intpl_method),'bilinear') > 0 ) then
             nrgst_bl = nrgst_bl + 1
             hist%nstt(idx) = nrgst_bl
@@ -239,17 +243,18 @@ subroutine history_type_register(hist, Diag, Time, Atm_block, Model, xlon, xlat,
             endif
           endif
         else if (diag(idx)%axes == 3) then
+          hist%levo(idx) = size(Diag(idx)%data(1)%var3, dim=2)
           if( index(trim(diag(idx)%intpl_method),'bilinear') > 0 ) then
             hist%nstt(idx) = nrgst_bl + 1
-            nrgst_bl  = nrgst_bl + hist%levo
+            nrgst_bl  = nrgst_bl + hist%levo(idx)
           else if (trim(diag(idx)%intpl_method) == 'nearest_stod' ) then
             hist%nstt(idx) = nrgst_nb + 1
-            nrgst_nb  = nrgst_nb + hist%levo
+            nrgst_nb  = nrgst_nb + hist%levo(idx)
           endif
           if(trim(diag(idx)%intpl_method) == 'vector_bilinear') then
             if(diag(idx)%name(1:1) == 'v' .or. diag(idx)%name(1:1) == 'V') then
               hist%nstt_vctbl(idx) = nrgst_vctbl + 1
-              nrgst_vctbl = nrgst_vctbl + hist%levo
+              nrgst_vctbl = nrgst_vctbl + hist%levo(idx)
               !             print *,'in phy_setup, vector_bilinear, name=', trim(diag(idx)%name),' nstt_vctbl=', hist%nstt_vctbl(idx), 'idx=',idx
             endif
           endif
@@ -289,14 +294,14 @@ subroutine history_type_output(hist, time, diag, atm_block, nx, ny, levs, ntcw,
     real(kind=kind_phys),      intent(in) :: time_radsw
     real(kind=kind_phys),      intent(in) :: time_radlw
     !--- local variables
-    integer :: i, j, k, idx, nb, ix, ii, jj
+    integer :: i, j, k, idx, nb, ix, ii, jj, levo_3d
     character(len=2) :: xtra
 #ifdef CCPP_32BIT
     real, dimension(nx,ny)      :: var2
-    real, dimension(nx,ny,levs) :: var3
+    real, dimension(:,:,:), allocatable :: var3
 #else
     real(kind=kind_phys), dimension(nx,ny)      :: var2
-    real(kind=kind_phys), dimension(nx,ny,levs) :: var3
+    real(kind=kind_phys), dimension(:,:,:), allocatable :: var3
 #endif
     real(kind=kind_phys) :: rtime_int, rtime_intfull, lcnvfac
     real(kind=kind_phys) :: rtime_radsw, rtime_radlw
@@ -444,31 +449,37 @@ subroutine history_type_output(hist, time, diag, atm_block, nx, ny, levs, ntcw,
               enddo
             endif if_mask
           endif int_or_real
-          !           used=send_data(Diag(idx)%id, var2, Time)
-          !           print *,'in phys, after store_data, idx=',idx,' var=', trim(Diag(idx)%name)
+
           call hist%store_data(Diag(idx)%id, var2, Time, idx, Diag(idx)%intpl_method, Diag(idx)%name)
-          !           if(trim(Diag(idx)%name) == 'totprcp_ave' ) print *,'in gfs_io, totprcp=',Diag(idx)%data(1)%var2(1:3), &
-          !             ' lcnvfac=', lcnvfac
+
         elseif (Diag(idx)%axes == 3) then
           !---
           !--- skipping other 3D variables with the following else statement
           !---
-          !         if(mpp_pe()==mpp_root_pe())print *,'in,fv3atm_io. 3D fields, idx=',idx,'varname=',trim(diag(idx)%name), &
-          !             'lcnvfac=',lcnvfac, 'hist%levo=',hist%levo,'nx=',nx,'ny=',ny
-          do k=1, hist%levo
+
+          levo_3d = hist%levo(idx)
+          allocate(var3(nx,ny,levo_3d))
+
+          do k=1, levo_3d
             do j = 1, ny
               jj = j + Atm_block%jsc -1
               do i = 1, nx
                 ii = i + Atm_block%isc -1
                 nb = Atm_block%blkno(ii,jj)
                 ix = Atm_block%ixp(ii,jj)
-                !         if(mpp_pe()==mpp_root_pe())print *,'in,fv3atm_io,sze(Diag(idx)%data(nb)%var3)=',  &
-                !             size(Diag(idx)%data(nb)%var3,1),size(Diag(idx)%data(nb)%var3,2)
-                var3(i,j,k) = Diag(idx)%data(nb)%var3(ix,hist%levo-k+1)*lcnvfac
+                ! flip only 3d variables with vertical dimension == levs (atm model levels)
+                if (levo_3d == levs) then
+                  var3(i,j,k) = Diag(idx)%data(nb)%var3(ix,levo_3d-k+1)*lcnvfac
+                else
+                  var3(i,j,k) = Diag(idx)%data(nb)%var3(ix,        k  )*lcnvfac
+                endif
               enddo
             enddo
           enddo
+
           call hist%store_data3D(Diag(idx)%id, var3, Time, idx, Diag(idx)%intpl_method, Diag(idx)%name)
+          deallocate(var3)
+
         endif if_2d
       endif has_id
     end do history_loop
@@ -574,7 +585,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
 #ifdef CCPP_32BIT
     real, intent(in)                    :: work(:,:,:)
 #else
-    real(kind=kind_phys), intent(in)    :: work(hist%ieco-hist%isco+1,hist%jeco-hist%jsco+1,hist%levo)
+    real(kind=kind_phys), intent(in)    :: work(hist%ieco-hist%isco+1,hist%jeco-hist%jsco+1,hist%levo(idx))
 #endif
     type(time_type), intent(in)         :: Time
     character(*), intent(in)            :: intpl_method
@@ -584,11 +595,12 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
     integer k,j,i,nv,i1,j1
     logical used
     !
+    write(0,*) ' history_type_store_data3D kinds ', kind_phys, kind(work), lbound(work), ubound(work), size(work)
     if( id > 0 ) then
       if( hist%use_wrtgridcomp_output ) then
         if( trim(intpl_method) == 'bilinear') then
           !$omp parallel do default(shared) private(i,j,k)
-          do k= 1,hist%levo
+          do k= 1,hist%levo(idx)
             do j= hist%jsco,hist%jeco
               do i= hist%isco,hist%ieco
                 hist%buffer_phys_bl(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k)
@@ -597,7 +609,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
           enddo
         else if(trim(intpl_method) == 'nearest_stod') then
           !$omp parallel do default(shared) private(i,j,k)
-          do k= 1,hist%levo
+          do k= 1,hist%levo(idx)
             do j= hist%jsco,hist%jeco
               do i= hist%isco,hist%ieco
                 hist%buffer_phys_nb(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k)
@@ -607,7 +619,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
         else if(trim(intpl_method) == 'vector_bilinear') then
           !first save the data
           !$omp parallel do default(shared) private(i,j,k)
-          do k= 1,hist%levo
+          do k= 1,hist%levo(idx)
             do j= hist%jsco,hist%jeco
               do i= hist%isco,hist%ieco
                 hist%buffer_phys_bl(i,j,hist%nstt(idx)+k-1) = work(i-hist%isco+1,j-hist%jsco+1,k)
@@ -615,9 +627,9 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
             enddo
           enddo
           if( fldname(1:1) == 'u' .or. fldname(1:1) == 'U') then
-            if(.not.associated(hist%uwork3d)) allocate(hist%uwork3d(hist%isco:hist%ieco,hist%jsco:hist%jeco,hist%levo))
+            if(.not.associated(hist%uwork3d)) allocate(hist%uwork3d(hist%isco:hist%ieco,hist%jsco:hist%jeco,hist%levo(idx)))
             !$omp parallel do default(shared) private(i,j,k)
-            do k= 1, hist%levo
+            do k= 1, hist%levo(idx)
               do j= hist%jsco,hist%jeco
                 do i= hist%isco,hist%ieco
                   hist%uwork3d(i,j,k) = work(i-hist%isco+1,j-hist%jsco+1,k)
@@ -642,7 +654,7 @@ subroutine history_type_store_data3D(hist, id, work, Time, idx, intpl_method, fl
                 enddo
               enddo
               !$omp parallel do default(shared) private(i,j,k,nv,i1,j1)
-              do k= 1, hist%levo
+              do k= 1, hist%levo(idx)
                 nv = hist%nstt_vctbl(idx)+k-1
                 do j= hist%jsco,hist%jeco
                   j1 = j-hist%jsco+1
@@ -895,6 +907,36 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q
         deallocate(axis_data)
       enddo
     endif
+
+    ! add zsoil axis
+    call ESMF_AttributeAdd(fcst_grid, convention="NetCDF", purpose="FV3",  &
+         attrList=(/"zsoil               ", &
+                    "zsoil:long_name     ", &
+                    "zsoil:units         ", &
+                    "zsoil:cartesian_axis", &
+                    "zsoil:positive      "/), rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+    call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", &
+         name="zsoil", valueList=(/ (i, i=1,hist%nsoil_lsm) /), rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+    call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", &
+         name="zsoil:long_name", value="soil level", rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+    call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", &
+         name="zsoil:units", value="level", rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+    call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", &
+         name="zsoil:cartesian_axis", value="Z", rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
+    call ESMF_AttributeSet(fcst_grid, convention="NetCDF", purpose="FV3", &
+         name="zsoil:positive", value="down", rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+
     !   print *,'in setup fieldbundle_phys, hist%num_axes_phys=',hist%num_axes_phys,'hist%tot_diag_idx=',hist%tot_diag_idx, &
     !       'nbdlphys=',nbdlphys
     !
@@ -913,7 +955,7 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q
 
             !add origin field
             call hist%add_field_to_phybundle(trim(output_name),trim(Diag(idx)%desc),trim(Diag(idx)%unit), "time: point",         &
-                 axes(1:Diag(idx)%axes), fcst_grid, hist%nstt(idx), phys_bundle(ibdl), outputfile(ibdl),  &
+                 axes(1:Diag(idx)%axes), fcst_grid, hist%nstt(idx), hist%levo(idx), phys_bundle(ibdl), outputfile(ibdl),  &
                  bdl_intplmethod(ibdl), rcd=rc)
             !           if( mpp_pe() == mpp_root_pe()) print *,'phys, add field,',trim(Diag(idx)%name),'idx=',idx,'ibdl=',ibdl
             !
@@ -923,7 +965,7 @@ subroutine history_type_bundle_setup(hist, Diag, axes, phys_bundle, fcst_grid, q
                 output_name = 'wind'//trim(output_name)//'vector'
                 outputfile1 = 'none'
                 call hist%add_field_to_phybundle(trim(output_name),trim(Diag(idx)%desc),trim(Diag(idx)%unit), "time: point",       &
-                     axes(1:Diag(idx)%axes), fcst_grid, hist%nstt_vctbl(idx),phys_bundle(ibdl), outputfile1, &
+                     axes(1:Diag(idx)%axes), fcst_grid, hist%nstt_vctbl(idx), hist%levo(idx), phys_bundle(ibdl), outputfile1, &
                      bdl_intplmethod(ibdl),l2dvector=l2dvector,  rcd=rc)
                 !               if( mpp_pe() == mpp_root_pe()) print *,'in phys, add vector field,',trim(Diag(idx)%name),' idx=',idx,' ibdl=',ibdl
               endif
@@ -951,7 +993,7 @@ end subroutine history_type_bundle_setup
   !! It sets attributes for and logs information about a single ESMF field. Do not call this subroutine directly.
   !! Call fv_phys_bundle_setup instead.
   subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cell_methods, axes,phys_grid, &
-       kstt,phys_bundle,output_file,intpl_method,range,l2dvector,rcd)
+       kstt,levo,phys_bundle,output_file,intpl_method,range,l2dvector,rcd)
     !
     use esmf
     !
@@ -961,7 +1003,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel
     character(*), intent(in)             :: output_file, intpl_method
     integer, intent(in)                  :: axes(:)
     type(esmf_grid), intent(in)          :: phys_grid
-    integer, intent(in)                  :: kstt
+    integer, intent(in)                  :: kstt, levo
     type(esmf_fieldbundle),intent(inout) :: phys_bundle
     real, intent(in), optional           :: range(2)
     logical, intent(in), optional        :: l2dvector
@@ -1027,7 +1069,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel
              line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
 
       else if(size(axes) == 3) then
-        temp_r3d => hist%buffer_phys_nb(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+hist%levo-1)
+        temp_r3d => hist%buffer_phys_nb(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+levo-1)
         field = ESMF_FieldCreate(phys_grid, temp_r3d, datacopyflag=copyflag, &
              name=var_name, indexFlag=ESMF_INDEX_DELOCAL, rc=rc)
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,       &
@@ -1043,7 +1085,7 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,       &
              line=__LINE__, file=__FILE__)) call ESMF_Finalize(endflag=ESMF_END_ABORT)
       else if(size(axes) == 3) then
-        temp_r3d => hist%buffer_phys_bl(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+hist%levo-1)
+        temp_r3d => hist%buffer_phys_bl(hist%isco:hist%ieco,hist%jsco:hist%jeco,kstt:kstt+levo-1)
         field = ESMF_FieldCreate(phys_grid, temp_r3d, datacopyflag=copyflag, &
              name=var_name, indexFlag=ESMF_INDEX_DELOCAL, rc=rc)
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU,       &
@@ -1128,8 +1170,14 @@ subroutine history_type_add_field_to_phybundle(hist,var_name,long_name,units,cel
           call ESMF_AttributeAdd(field, convention="NetCDF", purpose="FV3", &
                attrList=(/"ESMF:ungridded_dim_labels"/), rc=rc)
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
-          call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
-               name="ESMF:ungridded_dim_labels", valueList=(/trim(hist%axis_name(idx))/), rc=rc)
+
+          if (levo == hist%nsoil_lsm) then
+            call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
+                 name="ESMF:ungridded_dim_labels", valueList=(/"zsoil"/), rc=rc)
+          else
+            call ESMF_AttributeSet(field, convention="NetCDF", purpose="FV3", &
+                 name="ESMF:ungridded_dim_labels", valueList=(/trim(hist%axis_name(idx))/), rc=rc)
+          endif
           if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
         endif
       enddo
diff --git a/io/fv3atm_restart_io.F90 b/io/fv3atm_restart_io.F90
index d32be0586..a567f992f 100644
--- a/io/fv3atm_restart_io.F90
+++ b/io/fv3atm_restart_io.F90
@@ -14,7 +14,8 @@ module fv3atm_restart_io_mod
                                 register_axis, register_restart_field, &
                                 register_variable_attribute, register_field, &
                                 read_restart, write_restart, write_data,     &
-                                get_global_io_domain_indices, get_dimension_size
+                                get_global_io_domain_indices, get_dimension_size, &
+                                global_att_exists, get_global_attribute
   use mpp_domains_mod,    only: domain2d
   use fv3atm_common_io,   only: create_2d_field_and_add_to_bundle, &
        create_3d_field_and_add_to_bundle, copy_from_gfs_data, axis_type
@@ -515,6 +516,7 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
     !--- directory of the input files
     character(5)  :: indir='INPUT'
     character(37) :: infile
+    character(2)  :: file_ver
     !--- fms2_io file open logic
     logical :: amiopen
     logical :: override_frac_grid
@@ -644,8 +646,19 @@ subroutine sfc_prop_restart_read (Sfcprop, Atm_block, Model, fv_domain, warm_sta
     amiopen=open_file(Sfc_restart, trim(infile), "read", domain=fv_domain, is_restart=.true., dont_add_res_to_filename=.true.)
     if( .not.amiopen ) call mpp_error(FATAL, 'Error opening file'//trim(infile))
 
+    if (global_att_exists(Sfc_restart, "file_version")) then
+      call get_global_attribute(Sfc_restart, "file_version", file_ver)
+      if (file_ver == "V2") then
+        sfc%is_v2_file=.true.
+      endif
+    endif
+
     if(sfc%allocate_arrays(Model, Atm_block, .true., warm_start)) then
-      call sfc%fill_2d_names(Model, warm_start)
+      if (sfc%is_v2_file) then
+        call sfc%fill_2d_names_v2(Model, warm_start)
+      else
+        call sfc%fill_2d_names(Model, warm_start)
+      endif
       call sfc%register_axes(Model, Sfc_restart, .true., warm_start)
 
       ! Tell CLM Lake to allocate data, and register its axes and fields
diff --git a/io/fv3atm_sfc_io.F90 b/io/fv3atm_sfc_io.F90
index c0bfcf6d9..95957682a 100644
--- a/io/fv3atm_sfc_io.F90
+++ b/io/fv3atm_sfc_io.F90
@@ -21,7 +21,8 @@ module fv3atm_sfc_io
   private
 
   public :: Sfc_io_data_type
-  public :: Sfc_io_fill_2d_names, Sfc_io_fill_3d_names, Sfc_io_allocate_arrays, &
+  public :: Sfc_io_fill_2d_names, Sfc_io_fill_2d_names_v2, &
+       Sfc_io_fill_3d_names, Sfc_io_allocate_arrays, &
        Sfc_io_register_axes, Sfc_io_write_axes, Sfc_io_register_2d_fields, &
        Sfc_io_register_3d_fields, Sfc_io_copy_to_grid, Sfc_io_copy_from_grid, &
        Sfc_io_apply_safeguards, Sfc_io_transfer, Sfc_io_final
@@ -49,6 +50,8 @@ module fv3atm_sfc_io
     ! The lsoil flag is only meaningful when reading:;
     logical, public :: is_lsoil = .false.
 
+    logical, public :: is_v2_file = .false.
+
     ! SYNONYMS: Some nvar variables had two names in fv3atm_io.F90. They have
     ! only one name here. The "_s" is redundant because this file only has
     ! surface restart variables.
@@ -79,6 +82,7 @@ module fv3atm_sfc_io
     procedure, public :: bundle_2d_fields => Sfc_io_bundle_2d_fields
     procedure, public :: bundle_3d_fields => Sfc_io_bundle_3d_fields
     procedure, public :: fill_2d_names => Sfc_io_fill_2d_names
+    procedure, public :: fill_2d_names_v2 => Sfc_io_fill_2d_names_v2
     procedure, public :: fill_3d_names => Sfc_io_fill_3d_names
     procedure, public :: init_fields => Sfc_io_init_fields
     procedure, public :: transfer => Sfc_io_transfer
@@ -558,6 +562,174 @@ subroutine Sfc_io_fill_2d_names(sfc,Model,warm_start)
     endif
   end subroutine Sfc_io_fill_2d_names
 
+  !>@ Fills the name2d array with all surface 2D field names. Updates nvar2m if needed.
+  !!  This routine is for v2 coldstart files.
+  subroutine Sfc_io_fill_2d_names_v2(sfc,Model,warm_start)
+    implicit none
+    class(Sfc_io_data_type)           :: sfc
+    type(GFS_control_type),    intent(in) :: Model
+    logical, intent(in) :: warm_start
+    integer :: nt
+
+    !--- names of the 2D variables to save
+    nt=0
+    nt=nt+1 ; sfc%name2(nt) = 'slmsk'
+    nt=nt+1 ; sfc%name2(nt) = 'tsea'    ! tsfc
+    nt=nt+1 ; sfc%name2(nt) = 'sheleg'  ! weasd in file. Optional for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'tg3'
+    nt=nt+1 ; sfc%name2(nt) = 'zorl'    ! Optional for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'alvsf'
+    nt=nt+1 ; sfc%name2(nt) = 'alvwf'
+    nt=nt+1 ; sfc%name2(nt) = 'alnsf'
+    nt=nt+1 ; sfc%name2(nt) = 'alnwf'
+    nt=nt+1 ; sfc%name2(nt) = 'facsf'
+    nt=nt+1 ; sfc%name2(nt) = 'facwf'
+    nt=nt+1 ; sfc%name2(nt) = 'vfrac'
+    nt=nt+1 ; sfc%name2(nt) = 'canopy'
+    nt=nt+1 ; sfc%name2(nt) = 'f10m'
+    nt=nt+1 ; sfc%name2(nt) = 't2m'
+    nt=nt+1 ; sfc%name2(nt) = 'q2m'
+    nt=nt+1 ; sfc%name2(nt) = 'vtype'
+    nt=nt+1 ; sfc%name2(nt) = 'stype'
+    nt=nt+1 ; sfc%name2(nt) = 'uustar'
+    nt=nt+1 ; sfc%name2(nt) = 'ffmm'
+    nt=nt+1 ; sfc%name2(nt) = 'ffhh'
+    nt=nt+1 ; sfc%name2(nt) = 'hice'
+    nt=nt+1 ; sfc%name2(nt) = 'fice'
+    nt=nt+1 ; sfc%name2(nt) = 'tisfc'
+    nt=nt+1 ; sfc%name2(nt) = 'tprcp'
+    nt=nt+1 ; sfc%name2(nt) = 'srflag'
+    nt=nt+1 ; sfc%name2(nt) = 'snwdph'  ! snowd in file. Optional for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'shdmin'
+    nt=nt+1 ; sfc%name2(nt) = 'shdmax'
+    nt=nt+1 ; sfc%name2(nt) = 'slope'
+    nt=nt+1 ; sfc%name2(nt) = 'snoalb'
+    !--- variables below here are optional, unless indicated.
+    nt=nt+1 ; sfc%name2(nt) = 'scolor'
+    nt=nt+1 ; sfc%name2(nt) = 'sncovr'
+    nt=nt+1 ; sfc%name2(nt) = 'snodl' ! snowd on land portion of a cell.
+                                      ! Mandatory for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'weasdl'! weasd on land portion of a cell.
+                                      ! Mandatory for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'tsfc'  ! tsfc composite
+    nt=nt+1 ; sfc%name2(nt) = 'tsfcl' ! temp on land portion of a cell
+    nt=nt+1 ; sfc%name2(nt) = 'zorlw' ! zorl on water portion of a cell
+                                      ! Mandatory for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'zorll' ! zorl on land portion of a cell
+    nt=nt+1 ; sfc%name2(nt) = 'zorli' ! zorl on ice portion of a cell
+                                      ! Mandatory for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'albdirvis_lnd'
+    nt=nt+1 ; sfc%name2(nt) = 'albdirnir_lnd'
+    nt=nt+1 ; sfc%name2(nt) = 'albdifvis_lnd'
+    nt=nt+1 ; sfc%name2(nt) = 'albdifnir_lnd'
+    nt=nt+1 ; sfc%name2(nt) = 'emis_lnd'
+    nt=nt+1 ; sfc%name2(nt) = 'emis_ice'
+    nt=nt+1 ; sfc%name2(nt) = 'sncovr_ice'
+    nt=nt+1 ; sfc%name2(nt) = 'snodi' ! snowd on ice portion of a cell.
+                                      ! Mandatory for cold starts.
+    nt=nt+1 ; sfc%name2(nt) = 'weasdi'! weasd on ice portion of a cell
+                                      ! Mandatory for cold starts.
+
+    if (Model%use_cice_alb .or. Model%lsm == Model%lsm_ruc) then
+      nt=nt+1 ; sfc%name2(nt) = 'albdirvis_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'albdifvis_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'albdirnir_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'albdifnir_ice'
+    endif
+
+    if(Model%cplwav) then
+      nt=nt+1 ; sfc%name2(nt) = 'zorlwav' !zorl from wave component
+      sfc%nvar2m = nt
+    endif
+
+    if (Model%nstf_name(1) > 0) then
+      !--- NSSTM inputs only needed when (nstf_name(1) > 0) .and. (nstf_name(2)) == 0)
+      nt=nt+1 ; sfc%name2(nt) = 'tref'
+      nt=nt+1 ; sfc%name2(nt) = 'z_c'
+      nt=nt+1 ; sfc%name2(nt) = 'c_0'
+      nt=nt+1 ; sfc%name2(nt) = 'c_d'
+      nt=nt+1 ; sfc%name2(nt) = 'w_0'
+      nt=nt+1 ; sfc%name2(nt) = 'w_d'
+      nt=nt+1 ; sfc%name2(nt) = 'xt'
+      nt=nt+1 ; sfc%name2(nt) = 'xs'
+      nt=nt+1 ; sfc%name2(nt) = 'xu'
+      nt=nt+1 ; sfc%name2(nt) = 'xv'
+      nt=nt+1 ; sfc%name2(nt) = 'xz'
+      nt=nt+1 ; sfc%name2(nt) = 'zm'
+      nt=nt+1 ; sfc%name2(nt) = 'xtts'
+      nt=nt+1 ; sfc%name2(nt) = 'xzts'
+      nt=nt+1 ; sfc%name2(nt) = 'd_conv'
+      nt=nt+1 ; sfc%name2(nt) = 'ifd'
+      nt=nt+1 ; sfc%name2(nt) = 'dt_cool'
+      nt=nt+1 ; sfc%name2(nt) = 'qrain'
+    endif
+    !
+    ! Only needed when Noah MP LSM is used - 29 2D
+    !
+    if (Model%lsm == Model%lsm_noahmp) then
+      nt=nt+1 ; sfc%name2(nt) = 'snowxy'
+      nt=nt+1 ; sfc%name2(nt) = 'tvxy'
+      nt=nt+1 ; sfc%name2(nt) = 'tgxy'
+      nt=nt+1 ; sfc%name2(nt) = 'canicexy'
+      nt=nt+1 ; sfc%name2(nt) = 'canliqxy'
+      nt=nt+1 ; sfc%name2(nt) = 'eahxy'
+      nt=nt+1 ; sfc%name2(nt) = 'tahxy'
+      nt=nt+1 ; sfc%name2(nt) = 'cmxy'
+      nt=nt+1 ; sfc%name2(nt) = 'chxy'
+      nt=nt+1 ; sfc%name2(nt) = 'fwetxy'
+      nt=nt+1 ; sfc%name2(nt) = 'sneqvoxy'
+      nt=nt+1 ; sfc%name2(nt) = 'alboldxy'
+      nt=nt+1 ; sfc%name2(nt) = 'qsnowxy'
+      nt=nt+1 ; sfc%name2(nt) = 'wslakexy'
+      nt=nt+1 ; sfc%name2(nt) = 'zwtxy'
+      nt=nt+1 ; sfc%name2(nt) = 'waxy'
+      nt=nt+1 ; sfc%name2(nt) = 'wtxy'
+      nt=nt+1 ; sfc%name2(nt) = 'lfmassxy'
+      nt=nt+1 ; sfc%name2(nt) = 'rtmassxy'
+      nt=nt+1 ; sfc%name2(nt) = 'stmassxy'
+      nt=nt+1 ; sfc%name2(nt) = 'woodxy'
+      nt=nt+1 ; sfc%name2(nt) = 'stblcpxy'
+      nt=nt+1 ; sfc%name2(nt) = 'fastcpxy'
+      nt=nt+1 ; sfc%name2(nt) = 'xsaixy'
+      nt=nt+1 ; sfc%name2(nt) = 'xlaixy'
+      nt=nt+1 ; sfc%name2(nt) = 'taussxy'
+      nt=nt+1 ; sfc%name2(nt) = 'smcwtdxy'
+      nt=nt+1 ; sfc%name2(nt) = 'deeprechxy'
+      nt=nt+1 ; sfc%name2(nt) = 'rechxy'
+    else if (Model%lsm == Model%lsm_ruc .and. warm_start) then
+      nt=nt+1 ; sfc%name2(nt) = 'wetness'
+      nt=nt+1 ; sfc%name2(nt) = 'clw_surf_land'
+      nt=nt+1 ; sfc%name2(nt) = 'clw_surf_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'qwv_surf_land'
+      nt=nt+1 ; sfc%name2(nt) = 'qwv_surf_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'tsnow_land'
+      nt=nt+1 ; sfc%name2(nt) = 'tsnow_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'snowfall_acc_land'
+      nt=nt+1 ; sfc%name2(nt) = 'snowfall_acc_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'sfalb_lnd'
+      nt=nt+1 ; sfc%name2(nt) = 'sfalb_lnd_bck'
+      nt=nt+1 ; sfc%name2(nt) = 'sfalb_ice'
+      if (Model%rdlai) then
+        nt=nt+1 ; sfc%name2(nt) = 'lai'
+      endif
+    else if (Model%lsm == Model%lsm_ruc .and. Model%rdlai) then
+      nt=nt+1 ; sfc%name2(nt) = 'lai'
+    endif
+
+    if (Model%lkm > 0 .and. Model%iopt_lake==Model%iopt_lake_flake) then
+      nt=nt+1 ; sfc%name2(nt) = 'T_snow'
+      nt=nt+1 ; sfc%name2(nt) = 'T_ice'
+      nt=nt+1 ; sfc%name2(nt) = 'h_ML'
+      nt=nt+1 ; sfc%name2(nt) = 't_ML'
+      nt=nt+1 ; sfc%name2(nt) = 't_mnw'
+      nt=nt+1 ; sfc%name2(nt) = 'h_talb'
+      nt=nt+1 ; sfc%name2(nt) = 't_talb'
+      nt=nt+1 ; sfc%name2(nt) = 't_bot1'
+      nt=nt+1 ; sfc%name2(nt) = 't_bot2'
+      nt=nt+1 ; sfc%name2(nt) = 'c_t'
+    endif
+  end subroutine Sfc_io_fill_2d_names_v2
+
   !>@ Registers 2D fields with FMS for reading or writing non-quilt restart files
   subroutine Sfc_io_register_2d_fields(sfc,Model,Sfc_restart,reading,warm_start)
     implicit none
@@ -590,33 +762,62 @@ subroutine Sfc_io_register_2d_fields(sfc,Model,Sfc_restart,reading,warm_start)
     endif
 
     !--- register the 2D fields
-    do num = 1,sfc%nvar2m
-      var2_p => sfc%var2(:,:,num)
-      if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'tsfcl' .or. trim(sfc%name2(num)) == 'zorll' &
-           .or. trim(sfc%name2(num)) == 'zorli' .or. trim(sfc%name2(num)) == 'zorlwav' &
-           .or. trim(sfc%name2(num)) == 'snodl' .or. trim(sfc%name2(num)) == 'weasdl'  &
-           .or. trim(sfc%name2(num)) == 'snodi' .or. trim(sfc%name2(num)) == 'weasdi'  &
-           .or. trim(sfc%name2(num)) == 'tsfc'  .or. trim(sfc%name2(num)) ==  'zorlw'  &
+    if (sfc%is_v2_file) then
+      do num = 1,sfc%nvar2m
+        var2_p => sfc%var2(:,:,num)
+        if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'zorll' &
+           .or. trim(sfc%name2(num)) == 'zorl' .or. trim(sfc%name2(num)) == 'zorlwav' &
+           .or. trim(sfc%name2(num)) == 'snwdph' .or. trim(sfc%name2(num)) == 'sheleg'  &
+           .or. trim(sfc%name2(num)) == 'tsfc' &
            .or. trim(sfc%name2(num)) == 'albdirvis_lnd' .or. trim(sfc%name2(num)) == 'albdirnir_lnd' &
            .or. trim(sfc%name2(num)) == 'albdifvis_lnd' .or. trim(sfc%name2(num)) == 'albdifnir_lnd' &
            .or. trim(sfc%name2(num)) == 'albdirvis_ice' .or. trim(sfc%name2(num)) == 'albdirnir_ice' &
            .or. trim(sfc%name2(num)) == 'albdifvis_ice' .or. trim(sfc%name2(num)) == 'albdifnir_ice' &
            .or. trim(sfc%name2(num)) == 'emis_lnd'      .or. trim(sfc%name2(num)) == 'emis_ice'      &
            .or. trim(sfc%name2(num)) == 'sncovr_ice'    .or. trim(sfc%name2(num)) == 'scolor') then
-        if(reading .and. sfc%is_lsoil) then
-          call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.)
+          if(reading .and. sfc%is_lsoil) then
+            call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.)
+          else
+            call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,&
+               & chunksizes=chunksizes2d, is_optional=.true.)
+          end if
         else
-          call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,&
+          if(reading .and. sfc%is_lsoil) then
+            call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/))
+          else
+            call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d)
+          end if
+        endif
+      enddo
+    else
+      do num = 1,sfc%nvar2m
+        var2_p => sfc%var2(:,:,num)
+        if (trim(sfc%name2(num)) == 'sncovr' .or. trim(sfc%name2(num)) == 'tsfcl' .or. trim(sfc%name2(num)) == 'zorll' &
+             .or. trim(sfc%name2(num)) == 'zorli' .or. trim(sfc%name2(num)) == 'zorlwav' &
+             .or. trim(sfc%name2(num)) == 'snodl' .or. trim(sfc%name2(num)) == 'weasdl'  &
+             .or. trim(sfc%name2(num)) == 'snodi' .or. trim(sfc%name2(num)) == 'weasdi'  &
+             .or. trim(sfc%name2(num)) == 'tsfc'  .or. trim(sfc%name2(num)) ==  'zorlw'  &
+             .or. trim(sfc%name2(num)) == 'albdirvis_lnd' .or. trim(sfc%name2(num)) == 'albdirnir_lnd' &
+             .or. trim(sfc%name2(num)) == 'albdifvis_lnd' .or. trim(sfc%name2(num)) == 'albdifnir_lnd' &
+             .or. trim(sfc%name2(num)) == 'albdirvis_ice' .or. trim(sfc%name2(num)) == 'albdirnir_ice' &
+             .or. trim(sfc%name2(num)) == 'albdifvis_ice' .or. trim(sfc%name2(num)) == 'albdifnir_ice' &
+             .or. trim(sfc%name2(num)) == 'emis_lnd'      .or. trim(sfc%name2(num)) == 'emis_ice'      &
+             .or. trim(sfc%name2(num)) == 'sncovr_ice'    .or. trim(sfc%name2(num)) == 'scolor') then
+          if(reading .and. sfc%is_lsoil) then
+            call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=(/'lat','lon'/), is_optional=.true.)
+          else
+            call register_restart_field(Sfc_restart, sfc%name2(num), var2_p, dimensions=time2d,&
                & chunksizes=chunksizes2d, is_optional=.true.)
-        end if
-      else
-        if(reading .and. sfc%is_lsoil) then
-          call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/))
+          end if
         else
-          call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d)
-        end if
-      endif
-    enddo
+          if(reading .and. sfc%is_lsoil) then
+            call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=(/'lat','lon'/))
+          else
+            call register_restart_field(Sfc_restart,sfc%name2(num),var2_p, dimensions=time2d, chunksizes=chunksizes2d)
+          end if
+        endif
+      enddo
+    endif
 
     if (Model%nstf_name(1) > 0) then
       mand = .false.
@@ -702,7 +903,11 @@ subroutine Sfc_io_register_3d_fields(sfc,Model,Sfc_restart,reading,warm_start)
 
     !--- register the 3D fields
     var3_p => sfc%var3ice(:,:,:)
-    call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.true.)
+    if (sfc%is_v2_file) then
+      call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.false.)
+    else
+      call register_restart_field(Sfc_restart, sfc%name3(0), var3_p, dimensions=xyz1_time, chunksizes=chunksizes3d, is_optional=.true.)
+    endif
 
     if(reading) then
       do num = 1,sfc%nvar3
@@ -1309,6 +1514,157 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
+    if (sfc%is_v2_file) then
+
+    if (sfc%var2(i,j,27) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snowd')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          if (Sfcprop(nb)%fice(ix) > zero) then
+            Sfcprop(nb)%snowd(ix)  = Sfcprop(nb)%snodi(ix)
+          elseif (Sfcprop(nb)%landfrac(ix) > zero) then
+            Sfcprop(nb)%snowd(ix)  = Sfcprop(nb)%snodl(ix)
+          else
+            Sfcprop(nb)%snowd(ix)  = zero
+          endif
+        enddo
+      enddo
+    endif
+
+    if (sfc%var2(i,j,3) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing weasd')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          if (Sfcprop(nb)%fice(ix) > zero) then
+            Sfcprop(nb)%weasd(ix) = Sfcprop(nb)%weasdi(ix)
+          elseif (Sfcprop(nb)%landfrac(ix) > zero) then
+            Sfcprop(nb)%weasd(ix) = Sfcprop(nb)%weasdl(ix)
+          else
+            Sfcprop(nb)%weasd(ix) = zero
+          endif
+        enddo
+      enddo
+    endif
+
+! Needed for first time step in radiation before Noah/NoahMP sets it from look up table.
+! Just use a nominal value.
+    if (sfc%var2(i,j,39) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorll')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          if (Sfcprop(nb)%landfrac(ix) > zero) then
+            Sfcprop(nb)%zorll(ix) = 25.0
+          endif
+        enddo
+      enddo
+    endif
+
+    if (sfc%var2(i,j,5) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorl')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          if (Sfcprop(nb)%fice(ix) > zero) then
+            Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorli(ix)
+          elseif (Sfcprop(nb)%landfrac(ix) > zero) then
+            Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorll(ix)
+          else
+            Sfcprop(nb)%zorl(ix) = Sfcprop(nb)%zorlw(ix)
+          endif
+        enddo
+      enddo
+    endif
+
+    if (sfc%var2(i,j,46) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing emis_ice')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          Sfcprop(nb)%emis_ice(ix) = 0.96
+        enddo
+      enddo
+    endif
+
+    if (sfc%var2(i,j,47) < -9990.0_kind_phys .and. Model%lsm /= Model%lsm_ruc) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing sncovr_ice')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          !         Sfcprop(nb)%sncovr_ice(ix) = Sfcprop(nb)%sncovr(ix)
+          Sfcprop(nb)%sncovr_ice(ix) = zero
+        enddo
+      enddo
+    endif
+
+    if (Model%use_cice_alb) then
+      if (sfc%var2(i,j,50) < -9990.0_kind_phys) then
+        !$omp parallel do default(shared) private(nb, ix)
+        do nb = 1, Atm_block%nblks
+          do ix = 1, Atm_block%blksz(nb)
+            if (Sfcprop(nb)%oceanfrac(ix) > zero .and. &
+                 Sfcprop(nb)%fice(ix) >= Model%min_seaice) then
+              Sfcprop(nb)%albdirvis_ice(ix) = 0.6_kind_phys
+              Sfcprop(nb)%albdifvis_ice(ix) = 0.6_kind_phys
+              Sfcprop(nb)%albdirnir_ice(ix) = 0.6_kind_phys
+              Sfcprop(nb)%albdifnir_ice(ix) = 0.6_kind_phys
+            endif
+          enddo
+        enddo
+      endif
+
+    endif
+
+    ! Fill in composite tsfc for coldstart runs - must happen after tsfcl is computed
+    compute_tsfc_for_coldstart: if (sfc%var2(i,j,36) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing composite tsfc')
+      if(Model%frac_grid) then ! 3-way composite
+        !$omp parallel do default(shared) private(nb, ix, tem, tem1)
+        do nb = 1, Atm_block%nblks
+          do ix = 1, Atm_block%blksz(nb)
+            Sfcprop(nb)%tsfco(ix) = max(con_tice, Sfcprop(nb)%tsfco(ix)) ! this may break restart reproducibility
+            tem1 = one - Sfcprop(nb)%landfrac(ix)
+            tem  = tem1 * Sfcprop(nb)%fice(ix) ! tem = ice fraction wrt whole cell
+            Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix) * Sfcprop(nb)%landfrac(ix) &
+                 + Sfcprop(nb)%tisfc(ix) * tem                      &
+                 + Sfcprop(nb)%tsfco(ix) * (tem1-tem)
+          enddo
+        enddo
+      else
+        !$omp parallel do default(shared) private(nb, ix, tem)
+        do nb = 1, Atm_block%nblks
+          do ix = 1, Atm_block%blksz(nb)
+            if (Sfcprop(nb)%slmsk(ix) == 1) then
+              Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfcl(ix)
+              if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc land ',nb,ix,Sfcprop(nb)%tsfcl(ix)
+            elseif(Sfcprop(nb)%fice(ix) > 0.0)then
+              Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tisfc(ix)
+              if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc ice  ',nb,ix,Sfcprop(nb)%tisfc(ix)
+            else
+              Sfcprop(nb)%tsfc(ix) = Sfcprop(nb)%tsfco(ix)
+              if (Sfcprop(nb)%tsfc(ix) < -99 .or. Sfcprop(nb)%tsfc(ix) > 999.) print*,'bad tsfc water ',nb,ix,Sfcprop(nb)%tsfco(ix)
+            endif
+          enddo
+        enddo
+      endif
+    endif compute_tsfc_for_coldstart
+
+    if (sfc%var2(i,j,sfc%nvar2m) < -9990.0_kind_phys) then
+      if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing zorlwav')
+      !$omp parallel do default(shared) private(nb, ix)
+      do nb = 1, Atm_block%nblks
+        do ix = 1, Atm_block%blksz(nb)
+          Sfcprop(nb)%zorlwav(ix) = Sfcprop(nb)%zorl(ix) !--- compute zorlwav from existing variables
+        enddo
+      enddo
+    endif
+
+
+    else ! old verion of coldstart file
+
+
     if (sfc%var2(i,j,34) < -9990.0_kind_phys) then
       if (Model%me == Model%master ) call mpp_error(NOTE, 'gfs_driver::surface_props_input - computing snodl')
       !$omp parallel do default(shared) private(nb, ix, tem)
@@ -1503,6 +1859,8 @@ subroutine Sfc_io_apply_safeguards(sfc, Model, Atm_block, Sfcprop)
       enddo
     endif
 
+    endif ! check on which version of the surface file.
+
   end subroutine Sfc_io_apply_safeguards
 
   !>@ destructor for Sfc_io_data_type
diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90
index 082b38839..d9d8ff970 100644
--- a/io/module_write_netcdf.F90
+++ b/io/module_write_netcdf.F90
@@ -57,7 +57,7 @@ subroutine write_netcdf(wrtfb, filename, &
 !
 !** local vars
     integer :: i,j,t, istart,iend,jstart,jend
-    integer :: im, jm, lm
+    integer :: im, jm, lm, lsoil
 
     integer, dimension(:), allocatable              :: fldlev
 
@@ -95,9 +95,10 @@ subroutine write_netcdf(wrtfb, filename, &
     integer :: ncerr,ierr
     integer :: ncid
     integer :: oldMode
-    integer :: im_dimid, jm_dimid, tile_dimid, pfull_dimid, phalf_dimid, time_dimid, ch_dimid
+    integer :: dim_len
+    integer :: im_dimid, jm_dimid, tile_dimid, pfull_dimid, phalf_dimid, time_dimid, ch_dimid, lsoil_dimid
     integer :: im_varid, jm_varid, tile_varid, lon_varid, lat_varid, timeiso_varid
-    integer, dimension(:), allocatable :: dimids_2d, dimids_3d, dimids, chunksizes
+    integer, dimension(:), allocatable :: dimids_2d, dimids_3d, dimids_soil, dimids, chunksizes
     integer, dimension(:), allocatable :: varids
     integer :: xtype
     integer :: quant_mode
@@ -204,17 +205,24 @@ subroutine write_netcdf(wrtfb, filename, &
     end do
 
     lm = maxval(fldlev(:))
+    call get_dimlen_if_exists(ncid, "zsoil", wrtgrid, lsoil, rc)
+    if (lsoil > 0 .and. (.not. any(fldlev(:) == lsoil))) then
+       lsoil = 0
+    end if
+    if (lsoil > 0 .and. (.not. any(fldlev(:) > 1 .and. fldlev(:) /= lsoil))) then
+       lm = 1
+    end if
 
     ! for serial output allocate 'global' arrays
     if (.not. par) then
-       allocate(array_r4(im,jm))
+       ! allocate(array_r4(im,jm))
        allocate(array_r8(im,jm))
-       allocate(array_r4_3d(im,jm,lm))
+       ! allocate(array_r4_3d(im,jm,lm))
        allocate(array_r8_3d(im,jm,lm))
        if (is_cubed_sphere) then
-          allocate(array_r4_cube(im,jm,tileCount))
+          ! allocate(array_r4_cube(im,jm,tileCount))
           allocate(array_r8_cube(im,jm,tileCount))
-          allocate(array_r4_3d_cube(im,jm,lm,tileCount))
+          ! allocate(array_r4_3d_cube(im,jm,lm,tileCount))
           allocate(array_r8_3d_cube(im,jm,lm,tileCount))
        end if
     end if
@@ -243,6 +251,9 @@ subroutine write_netcdf(wrtfb, filename, &
          call add_dim(ncid, "pfull", pfull_dimid, wrtgrid, mype, rc)
          call add_dim(ncid, "phalf", phalf_dimid, wrtgrid, mype, rc)
        end if
+       if (lsoil > 1) then
+         call add_dim(ncid, "zsoil", lsoil_dimid, wrtgrid, mype, rc)
+       end if
        if (is_cubed_sphere) then
           ncerr = nf90_def_dim(ncid, "tile", tileCount, tile_dimid); NC_ERR_STOP(ncerr)
        end if
@@ -318,14 +329,16 @@ subroutine write_netcdf(wrtfb, filename, &
        ! define variables (fields)
        if (is_cubed_sphere) then
           allocate(dimids_2d(4))
-          allocate(dimids_3d(5))
-          dimids_2d =             [im_dimid,jm_dimid,            tile_dimid,time_dimid]
-          if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid,tile_dimid,time_dimid]
+          allocate(dimids_3d(5), dimids_soil(5))
+          dimids_2d =                  [im_dimid,jm_dimid,            tile_dimid,time_dimid]
+          if (lm > 1) dimids_3d =      [im_dimid,jm_dimid,pfull_dimid,tile_dimid,time_dimid]
+          if (lsoil > 1) dimids_soil = [im_dimid,jm_dimid,lsoil_dimid,tile_dimid,time_dimid]
        else
           allocate(dimids_2d(3))
-          allocate(dimids_3d(4))
-          dimids_2d =             [im_dimid,jm_dimid,                       time_dimid]
-          if (lm > 1) dimids_3d = [im_dimid,jm_dimid,pfull_dimid,           time_dimid]
+          allocate(dimids_3d(4), dimids_soil(4))
+          dimids_2d =                  [im_dimid,jm_dimid,                       time_dimid]
+          if (lm > 1) dimids_3d =      [im_dimid,jm_dimid,pfull_dimid,           time_dimid]
+          if (lsoil > 1) dimids_soil = [im_dimid,jm_dimid,lsoil_dimid,           time_dimid]
        end if
 
        do i=1, fieldCount
@@ -336,7 +349,11 @@ subroutine write_netcdf(wrtfb, filename, &
          if (rank == 2) then
            dimids = dimids_2d
          else if (rank == 3) then
-           dimids = dimids_3d
+           if (fldlev(i) == lsoil) then
+             dimids = dimids_soil
+           else if (fldlev(i) == lm) then
+             dimids = dimids_3d
+           endif
          else
            if (mype==0) write(0,*)'Unsupported rank ', rank
            call ESMF_Finalize(endflag=ESMF_END_ABORT)
@@ -368,7 +385,7 @@ subroutine write_netcdf(wrtfb, filename, &
                if (is_cubed_sphere) then
                   chunksizes = [im, jm, lm, tileCount, 1]
                else
-                  chunksizes = [ichunk3d(grid_id), jchunk3d(grid_id), kchunk3d(grid_id),            1]
+                  chunksizes = [ichunk3d(grid_id), jchunk3d(grid_id), fldlev(i), 1]
                end if
                ncerr = nf90_def_var_chunking(ncid, varids(i), NF90_CHUNKED, chunksizes) ; NC_ERR_STOP(ncerr)
             end if
@@ -610,17 +627,21 @@ subroutine write_netcdf(wrtfb, filename, &
             else
                if (is_cubed_sphere) then
                   call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc)
+                  allocate(array_r4_cube(im,jm,tileCount))
                   do t=1,tileCount
                      call ESMF_ArrayGather(array, array_r4_cube(:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc)
                   end do
                   if (do_io) then
                      ncerr = nf90_put_var(ncid, varids(i), values=array_r4_cube, start=start_idx); NC_ERR_STOP(ncerr)
                   end if
+                  deallocate(array_r4_cube)
                else
+                  allocate(array_r4(im,jm))
                   call ESMF_FieldGather(fcstField(i), array_r4, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc)
                   if (do_io) then
                      ncerr = nf90_put_var(ncid, varids(i), values=array_r4, start=start_idx); NC_ERR_STOP(ncerr)
                   end if
+                  deallocate(array_r4)
                end if
             end if
          else if (typekind == ESMF_TYPEKIND_R8) then
@@ -663,17 +684,21 @@ subroutine write_netcdf(wrtfb, filename, &
             else
                if (is_cubed_sphere) then
                   call ESMF_FieldGet(fcstField(i), array=array, rc=rc); ESMF_ERR_RETURN(rc)
+                  allocate(array_r4_3d_cube(im,jm,fldlev(i),tileCount))
                   do t=1,tileCount
                      call ESMF_ArrayGather(array, array_r4_3d_cube(:,:,:,t), rootPet=0, tile=t, rc=rc); ESMF_ERR_RETURN(rc)
                   end do
                   if (mype==0) then
                      ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d_cube, start=start_idx); NC_ERR_STOP(ncerr)
                   end if
+                  deallocate(array_r4_3d_cube)
                else
+                  allocate(array_r4_3d(im,jm,fldlev(i)))
                   call ESMF_FieldGather(fcstField(i), array_r4_3d, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc)
                   if (mype==0) then
                      ncerr = nf90_put_var(ncid, varids(i), values=array_r4_3d, start=start_idx); NC_ERR_STOP(ncerr)
                   end if
+                  deallocate(array_r4_3d)
                end if
             end if
          else if (typekind == ESMF_TYPEKIND_R8) then
@@ -708,14 +733,14 @@ subroutine write_netcdf(wrtfb, filename, &
     end do ! end fieldCount
 
     if (.not. par) then
-       deallocate(array_r4)
+       ! deallocate(array_r4)
        deallocate(array_r8)
-       deallocate(array_r4_3d)
+       ! deallocate(array_r4_3d)
        deallocate(array_r8_3d)
        if (is_cubed_sphere) then
-          deallocate(array_r4_cube)
+          ! deallocate(array_r4_cube)
           deallocate(array_r8_cube)
-          deallocate(array_r4_3d_cube)
+          ! deallocate(array_r4_3d_cube)
           deallocate(array_r8_3d_cube)
        end if
     end if
@@ -723,6 +748,7 @@ subroutine write_netcdf(wrtfb, filename, &
     if (do_io) then
        deallocate(dimids_2d)
        deallocate(dimids_3d)
+       deallocate(dimids_soil)
     end if
 
     deallocate(fcstField)
@@ -883,6 +909,22 @@ subroutine get_grid_attr(grid, prefix, ncid, varid, rc)
 
   end subroutine get_grid_attr
 
+!----------------------------------------------------------------------------------------
+  subroutine get_dimlen_if_exists(ncid, dim_name, grid, dim_len, rc)
+
+    integer, intent(in)             :: ncid
+    character(len=*), intent(in)    :: dim_name
+    type(ESMF_Grid), intent(in)     :: grid
+    integer, intent(out)            :: dim_len
+    integer, intent(out)            :: rc
+
+    dim_len = 0
+
+    call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
+                           name=dim_name, itemCount=dim_len, rc=rc); ESMF_ERR_RETURN(rc)
+
+  end subroutine get_dimlen_if_exists
+
   !> Add a dimension.
   !>
   !> @param[in] ncid NetCDF file ID.
@@ -907,6 +949,7 @@ subroutine add_dim(ncid, dim_name, dimid, grid, mype, rc)
     integer :: ncerr
     type(ESMF_TypeKind_Flag)   :: typekind
 
+    real(ESMF_KIND_I4), allocatable  :: valueListI4(:)
     real(ESMF_KIND_R4), allocatable  :: valueListR4(:)
     real(ESMF_KIND_R8), allocatable  :: valueListR8(:)
 !
@@ -944,6 +987,15 @@ subroutine add_dim(ncid, dim_name, dimid, grid, mype, rc)
        ncerr = nf90_put_var(ncid, dim_varid, values=valueListR4); NC_ERR_STOP(ncerr)
        ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
        deallocate(valueListR4)
+     else if (typekind==ESMF_TYPEKIND_I4) then
+       ncerr = nf90_def_var(ncid, dim_name, NF90_INT4, dimids=[dimid], varid=dim_varid); NC_ERR_STOP(ncerr)
+       allocate(valueListI4(n))
+       call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
+                              name=trim(dim_name), valueList=valueListI4, rc=rc); ESMF_ERR_RETURN(rc)
+       ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
+       ncerr = nf90_put_var(ncid, dim_varid, values=valueListI4); NC_ERR_STOP(ncerr)
+       ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
+       deallocate(valueListI4)
      else
         if (mype==0) write(0,*)'Error in module_write_netcdf.F90(add_dim) unknown typekind for ',trim(dim_name)
         call ESMF_Finalize(endflag=ESMF_END_ABORT)
diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90
index b7e93e28f..de0cedb6f 100644
--- a/io/module_wrt_grid_comp.F90
+++ b/io/module_wrt_grid_comp.F90
@@ -28,7 +28,9 @@ module module_wrt_grid_comp
 !
       use mpi
       use esmf
+      use fms_mod, only : uppercase
       use fms
+      use mpp_mod, only : mpp_init, mpp_error
 
       use write_internal_state
       use module_fv3_io_def,   only : num_pes_fcst,                             &
@@ -3363,6 +3365,7 @@ subroutine ioCompRun(comp, importState, exportState, clock, rc)
 
     integer                          :: localPet, petCount, i, j, k, ind
     type(ESMF_Grid)                  :: grid
+    real(ESMF_KIND_I4), allocatable  :: valueListi4(:)
     real(ESMF_KIND_R4), allocatable  :: valueListr4(:)
     real(ESMF_KIND_R8), allocatable  :: valueListr8(:)
     integer                          :: valueCount, fieldCount, udimCount
@@ -3747,6 +3750,12 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc)
         call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
                               name=trim(dimLabel), valueList=valueListr8, rc=rc)
 
+        if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
+      else if ( typekind == ESMF_TYPEKIND_I4) then
+        allocate(valueListi4(valueCount))
+        call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
+                              name=trim(dimLabel), valueList=valueListi4, rc=rc)
+
         if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return
       else
         write(0,*) 'in write_out_ungridded_dim_atts: ERROR unknown typekind'
@@ -3784,6 +3793,17 @@ subroutine write_out_ungridded_dim_atts(dimLabel, rc)
         ncerr = nf90_put_var(ncid, varid, values=valueListr8)
         if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return
         deallocate(valueListr8)
+      else if(typekind == ESMF_TYPEKIND_I4) then
+        ncerr = nf90_def_var(ncid, trim(dimLabel), NF90_INT4, &
+                             dimids=(/dimid/), varid=varid)
+        if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return
+
+        ncerr = nf90_enddef(ncid=ncid)
+        if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return
+
+        ncerr = nf90_put_var(ncid, varid, values=valueListi4)
+        if (ESMF_LogFoundNetCDFError(ncerr, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__, rcToReturn=rc)) return
+        deallocate(valueListi4)
       endif
       ! add attributes to this vertical variable
       call ESMF_AttributeGet(grid, convention="NetCDF", purpose="FV3", &
diff --git a/io/post_fv3.F90 b/io/post_fv3.F90
index 97962bdd9..30eae6ba4 100644
--- a/io/post_fv3.F90
+++ b/io/post_fv3.F90
@@ -3562,6 +3562,49 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            ! soilt
+            if(trim(fieldname)=='soilt') then
+              !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,stc,arrayr43d,sm,sice,fillvalue,spval)
+              do l=1,nsoil
+                do j=jsta,jend
+                  do i=ista, iend
+                    stc(i,j,l) = arrayr43d(i,j,l)
+                    if( abs(arrayr43d(i,j,l)-fillValue) < small) stc(i,j,l) = spval
+                    !mask open water areas, combine with sea ice tmp
+                    if (sm(i,j) /= 0.0 .and. sice(i,j) ==0.) stc(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! soilw
+            if(trim(fieldname)=='soilw') then
+              !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,smc,arrayr43d,sm,fillvalue,spval)
+              do l=1,nsoil
+                do j=jsta,jend
+                  do i=ista, iend
+                    smc(i,j,l) = arrayr43d(i,j,l)
+                    if( abs(arrayr43d(i,j,l)-fillValue) < small) smc(i,j,l) = spval
+                    if (sm(i,j) /= 0.0) smc(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
+            ! soill
+            if(trim(fieldname)=='soill') then
+              !$omp parallel do default(none) private(i,j,l) shared(nsoil,jsta,jend,ista,iend,sh2o,arrayr43d,sm,fillvalue,spval)
+              do l=1,nsoil
+                do j=jsta,jend
+                  do i=ista, iend
+                    sh2o(i,j,l) = arrayr43d(i,j,l)
+                    if( abs(arrayr43d(i,j,l)-fillValue) < small) sh2o(i,j,l) = spval
+                    if (sm(i,j) /= 0.0) sh2o(i,j,l) = spval
+                  enddo
+                enddo
+              enddo
+            endif
+
             ! model level ozone mixing ratio
 #ifdef MULTI_GASES
             if(trim(fieldname)=='spo3') then
diff --git a/upp b/upp
index 78f369b01..945cb2cef 160000
--- a/upp
+++ b/upp
@@ -1 +1 @@
-Subproject commit 78f369b011ec41dca87a1260298d0228f2c4f38f
+Subproject commit 945cb2cef5e8bd5949afd4f0fc35c4fb6e95a1bf