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/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/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..9a73f3d43 100644
--- a/io/module_wrt_grid_comp.F90
+++ b/io/module_wrt_grid_comp.F90
@@ -3363,6 +3363,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 +3748,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 +3791,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