diff --git a/tools/mksurfdata_esmf/src/mkfileMod.F90 b/tools/mksurfdata_esmf/src/mkfileMod.F90
index 6c3046a30b..335a9a508b 100644
--- a/tools/mksurfdata_esmf/src/mkfileMod.F90
+++ b/tools/mksurfdata_esmf/src/mkfileMod.F90
@@ -606,6 +606,17 @@ subroutine mkfile_define_vars(pioid, dynlanduse)
     call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_PFT', xtype=xtype, &
          long_name='land fraction from pft dataset (DIFF FROM landfrac USED IN SIMULATION, SHOWN IN HISTORY)', units='unitless')
 
+    if (.not. dynlanduse) then
+       call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_MKSURFDATA', xtype=xtype, &
+            long_name='land fraction used for renormalization of areas in mksurfdata (DIFF FROM landfrac USED IN SIMULATION)', &
+            units='unitless')
+    else
+       call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_MKSURFDATA', xtype=xtype, &
+            lev1name='time', &
+            long_name='land fraction used for renormalization of areas in mksurfdata (DIFF FROM landfrac USED IN SIMULATION)', &
+            units='unitless')
+    end if
+
     if (.not. dynlanduse) then
        call mkpio_def_spatial_var(pioid=pioid, varname='PCT_NATVEG', xtype=xtype, &
             long_name='total percent natural vegetation landunit', units='unitless')
diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90
index e3ea832926..0b304079bc 100644
--- a/tools/mksurfdata_esmf/src/mksurfdata.F90
+++ b/tools/mksurfdata_esmf/src/mksurfdata.F90
@@ -172,6 +172,9 @@ program mksurfdata
   real(r8), allocatable           :: pctwet_orig(:)          ! percent wetland of gridcell before dynamic land use adjustments
   real(r8), allocatable           :: pctgla_orig(:)          ! percent glacier of gridcell before dynamic land use adjustments
 
+  ! other variables written to file
+  real(r8), allocatable           :: landfrac_mksurfdata(:)  ! land fraction used for renormalization of areas
+
   ! pio/esmf variables
   type(file_desc_t)               :: pioid
   type(var_desc_t)                :: pio_varid
@@ -595,6 +598,7 @@ program mksurfdata
 
   ! Normalize land use and make sure things add up to 100% as well as
   ! checking that things are as they should be.
+  allocate(landfrac_mksurfdata(lsize_o))
   call normalize_and_check_landuse(lsize_o)
 
   ! Write out sum of PFT's
@@ -692,6 +696,10 @@ program mksurfdata
         call mkfile_output(pioid, mesh_model, 'PCT_CFT', pct_cft, rc=rc)
         if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_CFT')
      end if
+
+     if (root_task)  write(ndiag, '(a)') trim(subname)//" writing LANDFRAC_MKSURFDATA"
+     call mkfile_output(pioid, mesh_model, 'LANDFRAC_MKSURFDATA', landfrac_mksurfdata, rc=rc)
+     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_MKSURFDATA')
   end if
 
   ! ----------------------------------------------------------------------
@@ -779,7 +787,7 @@ program mksurfdata
      ! landfrac_pft was calculated ABOVE
      if (root_task)  write(ndiag, '(a)') trim(subname)//" writing land fraction calculated in fsurdata calc)"
      call mkfile_output(pioid,  mesh_model, 'LANDFRAC_PFT', landfrac_pft, rc=rc)
-     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output')
+     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_PFT')
 
      ! -----------------------------------------
      ! Read in each dynamic pft landuse dataset
@@ -993,6 +1001,13 @@ program mksurfdata
            call pio_syncfile(pioid)
         end if
 
+        if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing LANDFRAC_MKSURFDATA for year ",year
+        rcode = pio_inq_varid(pioid, 'LANDFRAC_MKSURFDATA', pio_varid)
+        call pio_setframe(pioid, pio_varid, int(ntim, kind=Pio_Offset_Kind))
+        call mkfile_output(pioid, mesh_model, 'LANDFRAC_MKSURFDATA', landfrac_mksurfdata, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_MKSURFDATA')
+        call pio_syncfile(pioid)
+
         if (root_task) then
            write(ndiag,'(1x,80a1)') ('=',k=1,80)
            write(ndiag,*)
@@ -1232,6 +1247,9 @@ subroutine normalize_and_check_landuse(ns_o)
             call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct)
          end if
 
+         ! Save landfrac for output to file
+         landfrac_mksurfdata(n) = frac_land
+
          ! Confirm that we have done the rescaling correctly: now the sum of all landunits
          ! should be 100% within tol_loose
          suma = pctlak(n) + pctwet(n) + pctgla(n) + pcturb(n) + pctcft(n)%get_pct_l2g()