From fa89274e2d90ad3273c98a3059184f8d06195528 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Tue, 16 Nov 2021 20:58:44 -0500
Subject: [PATCH 01/18] add PCT_URBAN in landuse.timeseries data

In normalizencheck_landuse, still reduce crop and natveg proportionally when urban exist.
---
 tools/mksurfdata_map/mksurfdata.pl         |   3 +
 tools/mksurfdata_map/src/mkfileMod.F90     |   7 +
 tools/mksurfdata_map/src/mksurfdat.F90     | 185 +++++++++++++--------
 tools/mksurfdata_map/src/mkurbanparMod.F90 |  29 +++-
 4 files changed, 150 insertions(+), 74 deletions(-)

diff --git a/tools/mksurfdata_map/mksurfdata.pl b/tools/mksurfdata_map/mksurfdata.pl
index f9961c7f9e..9d2fe238d3 100755
--- a/tools/mksurfdata_map/mksurfdata.pl
+++ b/tools/mksurfdata_map/mksurfdata.pl
@@ -284,6 +284,9 @@ sub write_transient_timeseries_file {
         my $hrvtypyr = `$scrdir/../../bld/queryDefaultNamelist.pl $queryfilopts $resolhrv -options sim_year='$yr',ssp_rcp=${ssp_rcp}${mkcrop} -var mksrf_fvegtyp -namelist clmexp`;
         chomp( $hrvtypyr );
         printf $fh_landuse_timeseries $dynpft_format, $hrvtypyr, $yr;
+        my $urbanyr = "/glade/scratch/keerzhang/archive/BNU_NoAdjust/05deg_".$yr."_new_NoAdjust.nc";
+        chomp( $urbanyr);
+        printf $fh_landuse_urban_timeseries $dynpft_format, $urbanyr, $yr; # Keer: I don't quiet understand the 'pft_override' option so I made no change to the "landuse_timeseries_override_$desc.txt"
         if ( $yr % 100 == 0 ) {
           print "year: $yr\n";
         }
diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 43bdda4c12..963e386348 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -540,7 +540,14 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
        deallocate(ind1D, ind2D)
 
     else
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN', xtype=xtype, &
+            lev1name='numurbl', lev2name='time', &
+            long_name='percent urban for each density type', units='unitless')
 
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN_MAX', xtype=xtype, &
+            lev1name='numurbl', &
+            long_name='maximum percent urban for each density type', units='unitless')
+         
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
           call ncd_def_spatial_var(ncid=ncid, varname=mkharvest_fieldname(ind1D(j),constant=.false.), xtype=xtype, &
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index bc6ef6f028..d1ae730130 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -29,7 +29,7 @@ program mksurfdat
     use mkharvestMod       , only : mkharvest_numtypes, mkharvest_parse_oride
     use mkharvestMod       , only : harvestDataType
     use mkurbanparCommonMod, only : mkelev
-    use mkurbanparMod      , only : mkurbanInit, mkurban, mkurbanpar, numurbl
+    use mkurbanparMod      , only : mkurbanInit, mkurban, mkurbanpar, numurbl, update_max_array_urban
     use mkutilsMod         , only : normalize_classes_by_gcell
     use mkfileMod          , only : mkfile
     use mkvarpar           , only : nlevsoi, elev_thresh, numstdpft
@@ -80,6 +80,7 @@ program mksurfdat
     character(len=256) :: fdyndat           ! dynamic landuse data file name
     character(len=256) :: fname             ! generic filename
     character(len=256) :: fhrvname          ! generic harvest filename
+    character(len=256) :: furbname          ! generic transient urban land cover filename
     character(len=256) :: string            ! string read in
     integer  :: t1                          ! timer
     real(r8),parameter :: p5  = 0.5_r8      ! constant
@@ -108,6 +109,7 @@ program mksurfdat
     real(r8), allocatable  :: pctlak(:)          ! percent of grid cell that is lake
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
+    real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
     real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: elev(:)            ! glc elevation (m)
@@ -134,6 +136,10 @@ program mksurfdat
     real(r8), allocatable  :: vic_dsmax(:)       ! VIC Dsmax parameter (mm/day)
     real(r8), allocatable  :: vic_ds(:)          ! VIC Ds parameter (unitless)
     real(r8), allocatable  :: lakedepth(:)       ! lake depth (m)
+    real(r8), allocatable  :: pctlak_orig(:)     ! percent lake of gridcell before dynamic land use adjustments            
+    real(r8), allocatable  :: pctwet_orig(:)     ! percent wetland of gridcell before dynamic land use adjustments 
+    real(r8), allocatable  :: pcturb_orig(:)     ! percent urban of gridcell before dynamic land use adjustments 
+    real(r8), allocatable  :: pctgla_orig(:)     ! percent glacier of gridcell before dynamic land use adjustments
 
     real(r8) :: std_elev = -999.99_r8            ! Standard deviation of elevation (m) to use for entire grid
 
@@ -277,7 +283,7 @@ program mksurfdat
     ! ======================================
     ! Optionally specify setting for:
     ! ======================================
-    !    mksrf_fdynuse ----- ASCII text file that lists each year of pft files to use
+    !    mksrf_fdynuse ----- ASCII text file that lists each year of pft and urban files to use
     !    mksrf_gridtype ---- Type of grid (default is 'global')
     !    outnc_double ------ If output should be in double precision
     !    outnc_large_files - If output should be in NetCDF large file format
@@ -442,6 +448,7 @@ program mksurfdat
                pctlak(ns_o)                       , & 
                pctwet(ns_o)                       , & 
                pcturb(ns_o)                       , &
+               pcturb_max(ns_o,numurbl)           , &
                urban_region(ns_o)                 , &
                urbn_classes(ns_o,numurbl)         , &
                urbn_classes_g(ns_o,numurbl)       , &
@@ -459,7 +466,11 @@ program mksurfdat
                vic_dsmax(ns_o)                    , &
                vic_ds(ns_o)                       , &
                lakedepth(ns_o)                    , &
-               glacier_region(ns_o)               )       
+               glacier_region(ns_o)               , &  
+               pctlak_orig(ns_o)                  , &
+               pctwet_orig(ns_o)                  , &
+               pcturb_orig(ns_o)                  , &
+               pctgla_orig(ns_o)                  )       
     landfrac_pft(:)       = spval 
     pctlnd_pft(:)         = spval
     pftdata_mask(:)       = -999
@@ -644,6 +655,12 @@ program mksurfdat
          ndiag=ndiag, zero_out=all_veg, urbn_o=pcturb, urbn_classes_o=urbn_classes, &
          region_o=urban_region)
 
+    ! Save special land unit areas of surface dataset 
+    pctlak_orig(:) = pctlak(:)
+    pctwet_orig(:) = pctwet(:)
+    pcturb_orig(:) = pcturb(:)
+    pctgla_orig(:) = pctgla(:)
+     
     ! Make elevation [elev] from [ftopo, ffrac] dataset
     ! Used only to screen pcturb  
     ! Screen pcturb by elevation threshold from elev dataset
@@ -662,7 +679,6 @@ program mksurfdat
        where (elev .gt. elev_thresh)
          pcturb = 0._r8
        end where
-       deallocate(elev)
     end if
     
     ! Compute topography statistics [topo_stddev, slope] from [ftopostats]
@@ -709,62 +725,6 @@ program mksurfdat
 
     call change_landuse( ldomain, dynpft=.false. )
 
-    do n = 1,ns_o
-
-       ! Truncate all percentage fields on output grid. This is needed to
-       ! insure that wt is zero (not a very small number such as
-       ! 1e-16) where it really should be zero
-       
-       do k = 1,nlevsoi
-          pctsand(n,k) = float(nint(pctsand(n,k)))
-          pctclay(n,k) = float(nint(pctclay(n,k)))
-       end do
-       pctlak(n) = float(nint(pctlak(n)))
-       pctwet(n) = float(nint(pctwet(n)))
-       pctgla(n) = float(nint(pctgla(n)))
-       
-       ! Assume wetland, glacier and/or lake when dataset landmask implies ocean 
-       ! (assume medium soil color (15) and loamy texture).
-       ! Also set pftdata_mask here
-
-       if (pctlnd_pft(n) < 1.e-6_r8) then
-          pftdata_mask(n)  = 0
-          soicol(n)        = 15
-          if (pctgla(n) < 1.e-6_r8) then
-              pctwet(n)    = 100._r8 - pctlak(n)
-              pctgla(n)    = 0._r8
-          else
-              pctwet(n)    = max(100._r8 - pctgla(n) - pctlak(n), 0.0_r8)
-          end if
-          pcturb(n)        = 0._r8
-          call pctnatpft(n)%set_pct_l2g(0._r8)
-          call pctcft(n)%set_pct_l2g(0._r8)
-          pctsand(n,:)     = 43._r8
-          pctclay(n,:)     = 18._r8
-          organic(n,:)   = 0._r8
-       else
-          pftdata_mask(n) = 1
-       end if
-
-       ! Make sure sum of land cover types does not exceed 100. If it does,
-       ! subtract excess from most dominant land cover.
-       
-       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
-       if (suma > 250._r4) then
-          write (6,*) subname, ' error: sum of pctlak, pctwet,', &
-               'pcturb and pctgla is greater than 250%'
-          write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', &
-               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n)
-          call abort()
-       else if (suma > 100._r4) then
-          pctlak(n) = pctlak(n) * 100._r8/suma
-          pctwet(n) = pctwet(n) * 100._r8/suma
-          pcturb(n) = pcturb(n) * 100._r8/suma
-          pctgla(n) = pctgla(n) * 100._r8/suma
-       end if
-       
-    end do
-
     call normalizencheck_landuse(ldomain)
 
     ! Write out sum of PFT's
@@ -1062,14 +1022,11 @@ program mksurfdat
 
     ! Deallocate arrays NOT needed for dynamic-pft section of code
 
-    deallocate ( organic )
     deallocate ( ef1_btr, ef1_fet, ef1_fdt, ef1_shr, ef1_grs, ef1_crp )
     deallocate ( pctglcmec, topoglcmec)
     if ( outnc_3dglc ) deallocate ( pctglc_gic, pctglc_icesheet)
     deallocate ( elevclass )
     deallocate ( fmax )
-    deallocate ( pctsand, pctclay )
-    deallocate ( soicol )
     deallocate ( gdp, fpeat, agfirepkmon )
     deallocate ( soildepth )
     deallocate ( topo_stddev, slope )
@@ -1130,6 +1087,7 @@ program mksurfdat
 
        pctnatpft_max = pctnatpft
        pctcft_max = pctcft
+       pcturb_max = urbn_classes_g
 
        ntim = 0
        do 
@@ -1159,6 +1117,13 @@ program mksurfdat
                 call abort()
              end if
           end if
+          ! Read input urban data
+          read(nfdyn, '(A195,1x,I4)', iostat=ier) furbname, year2
+          if ( year2 /= year ) then
+             write(6,*) subname, ' error: year for urban not equal to year for PFT files'
+             call abort()
+          end if
+      write(6,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname)    
           ntim = ntim + 1
 
           ! Create pctpft data at model resolution
@@ -1187,21 +1152,43 @@ program mksurfdat
              end if
           end do
 
+
+          call mkurban (ldomain, mapfname=map_furban, datfname=furbname, &
+               ndiag=ndiag, zero_out=all_veg, urbn_o=pcturb, urbn_classes_o=urbn_classes, &
+               region_o=urban_region)
+          ! screen pcturb using elevation
+          if ( .not. all_urban .and. .not. all_veg )then
+               where (elev .gt. elev_thresh)
+                    pcturb = 0._r8
+               end where
+          end if
+
+          ! For landunits NOT read each year: reset to their pre-adjustment values in preparation for redoing landunit area normalization
+          pctwet(:) = pctwet_orig(:)
+          pctlak(:) = pctlak_orig(:)
+          pctgla(:) = pctgla_orig(:)
+                                        
           call change_landuse(ldomain, dynpft=.true.)
 
           call normalizencheck_landuse(ldomain)
-	  
+	        call normalize_classes_by_gcell(urbn_classes, pcturb, urbn_classes_g)
+          
           call update_max_array(pctnatpft_max,pctnatpft)
           call update_max_array(pctcft_max,pctcft)
+          call update_max_array_urban(pcturb_max,urbn_classes_g)
 
           ! Output time-varying data for current year
-
+          call check_ret(nf_inq_varid(ncid, 'PCT_NATVEG', varid), subname)
+          call ncd_put_time_slice(ncid, varid, ntim, get_pct_l2g_array(pctnatpft))
+          
           call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname)
           call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctnatpft))
 
           call check_ret(nf_inq_varid(ncid, 'PCT_CROP', varid), subname)
           call ncd_put_time_slice(ncid, varid, ntim, get_pct_l2g_array(pctcft))
-
+          
+          call check_ret(nf_inq_varid(ncid, 'PCT_URBAN', varid), subname)
+          call ncd_put_time_slice(ncid, varid, ntim, urbn_classes_g)
           if (num_cft > 0) then
              call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname)
              call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctcft))
@@ -1241,6 +1228,9 @@ program mksurfdat
        call check_ret(nf_inq_varid(ncid, 'PCT_CROP_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
+       call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
+       call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)
+       
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)
           call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctcft_max)), subname)
@@ -1339,7 +1329,7 @@ subroutine normalizencheck_landuse(ldomain)
 ! Normalize land use and make sure things add up to 100% as well as
 ! checking that things are as they should be.
 !
-! Precondition: pctlak + pctwet + pcturb + pctgla <= 100 (within roundoff)
+! 
 !
 ! !USES:
     use mkpftConstantsMod , only : baregroundindex
@@ -1369,14 +1359,61 @@ subroutine normalizencheck_landuse(ldomain)
     real(r8), parameter :: toosmallPFT = 1.e-10_r8            ! tolerance for PFT's to ignore
     character(len=32) :: subname = 'normalizencheck_landuse'  ! subroutine name
 !-----------------------------------------------------------------------
-
-    ! ------------------------------------------------------------------------
-    ! Normalize vegetated area so that vegetated + special area is 100%
-    ! ------------------------------------------------------------------------
-
     ns_o = ldomain%ns
     do n = 1,ns_o
 
+       ! Truncate all percentage fields on output grid. This is needed to
+       ! insure that wt is zero (not a very small number such as
+       ! 1e-16) where it really should be zero
+ 
+       do k = 1,nlevsoi
+          pctsand(n,k) = float(nint(pctsand(n,k)))
+          pctclay(n,k) = float(nint(pctclay(n,k)))
+       end do
+       pctlak(n) = float(nint(pctlak(n)))
+       pctwet(n) = float(nint(pctwet(n)))
+       pctgla(n) = float(nint(pctgla(n)))
+ 
+       ! Assume wetland, glacier and/or lake when dataset landmask implies ocean 
+       ! (assume medium soil color (15) and loamy texture).
+       ! Also set pftdata_mask here
+ 
+       if (pctlnd_pft(n) < 1.e-6_r8) then
+          pftdata_mask(n)  = 0
+          soicol(n)        = 15
+          if (pctgla(n) < 1.e-6_r8) then
+              pctwet(n)    = 100._r8 - pctlak(n)
+              pctgla(n)    = 0._r8
+          else
+              pctwet(n)    = max(100._r8 - pctgla(n) - pctlak(n), 0.0_r8)
+          end if
+          pcturb(n)        = 0._r8
+          call pctnatpft(n)%set_pct_l2g(0._r8)
+          call pctcft(n)%set_pct_l2g(0._r8)
+          pctsand(n,:)     = 43._r8
+          pctclay(n,:)     = 18._r8
+          organic(n,:)   = 0._r8
+       else
+          pftdata_mask(n) = 1
+       end if
+ 
+       ! Make sure sum of land cover types does not exceed 100. If it does,
+       ! subtract excess from most dominant land cover.
+
+       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
+       if (suma > 250._r4) then
+          write (6,*) subname, ' error: sum of pctlak, pctwet,', &
+               'pcturb and pctgla is greater than 250%'
+          write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', &
+               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n)
+          call abort()
+       else if (suma > 100._r4) then
+          pctlak(n) = pctlak(n) * 100._r8/suma
+          pctwet(n) = pctwet(n) * 100._r8/suma
+          pcturb(n) = pcturb(n) * 100._r8/suma
+          pctgla(n) = pctgla(n) * 100._r8/suma
+       end if
+
        ! Check preconditions
        if ( pctlak(n) < 0.0_r8 )then
           write(6,*) subname, ' ERROR: pctlak is negative!'
@@ -1407,7 +1444,9 @@ subroutine normalizencheck_landuse(ldomain)
                n, pctlak(n), pctwet(n), pcturb(n), pctgla(n)
           call abort()
        end if
-
+       ! ------------------------------------------------------------------------
+       ! Normalize vegetated area so that vegetated + special area is 100%
+       ! ------------------------------------------------------------------------
        ! First normalize vegetated (natural veg + crop) cover so that the total of
        ! (vegetated + (special excluding urban)) is 100%. We'll deal with urban later.
        !
diff --git a/tools/mksurfdata_map/src/mkurbanparMod.F90 b/tools/mksurfdata_map/src/mkurbanparMod.F90
index 07319b1f27..35bbb0ad82 100644
--- a/tools/mksurfdata_map/src/mkurbanparMod.F90
+++ b/tools/mksurfdata_map/src/mkurbanparMod.F90
@@ -23,7 +23,8 @@ module mkurbanparMod
    public :: mkurbanInit
    public :: mkurban
    public :: mkurbanpar
-
+   public :: update_max_array_urban
+   
    ! The following could be private, but because there are associated test routines in a
    ! separate module, it needs to be public
    public :: normalize_urbn_by_tot
@@ -756,4 +757,30 @@ end subroutine lookup_and_check_err
 end subroutine mkurbanpar
 !------------------------------------------------------------------------------
 
+!-----------------------------------------------------------------------
+subroutine update_max_array_urban(pct_urbmax_arr,pct_urban_arr)
+  !
+  ! !DESCRIPTION:
+  ! Update the maximum percent cover of each urban class
+  ! 
+  ! !ARGUMENTS:
+  real(r8)         , intent(inout):: pct_urbmax_arr(:,:)           ! max percent cover of each urban class
+  real(r8)         , intent(in):: pct_urban_arr(:,:)           ! percent cover of each urban class that is used to update the old pct_urbmax_arr
+  !
+  ! !LOCAL VARIABLES:
+  integer :: n,k,ns              ! indices
+
+  character(len=*), parameter :: subname = 'update_max_array_urban'
+  !-----------------------------------------------------------------------
+  ns = size(pct_urban_arr,1)
+  do n = 1, ns
+     do k =1, numurbl
+        if (pct_urban_arr(n,k) > pct_urbmax_arr(n,k)) then 
+           pct_urbmax_arr(n,k) = pct_urban_arr(n,k)
+        end if
+      end do
+  end do
+
+end subroutine update_max_array_urban
+
 end module mkurbanparMod

From 5c69e8ceea71f6b6c6ab6b6728f83ed8eae1d405 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 17 Nov 2021 09:28:53 -0500
Subject: [PATCH 02/18] changed a txt file name

changed a txt file name
---
 tools/mksurfdata_map/mksurfdata.pl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/mksurfdata.pl b/tools/mksurfdata_map/mksurfdata.pl
index 9d2fe238d3..3caa4eb009 100755
--- a/tools/mksurfdata_map/mksurfdata.pl
+++ b/tools/mksurfdata_map/mksurfdata.pl
@@ -286,7 +286,7 @@ sub write_transient_timeseries_file {
         printf $fh_landuse_timeseries $dynpft_format, $hrvtypyr, $yr;
         my $urbanyr = "/glade/scratch/keerzhang/archive/BNU_NoAdjust/05deg_".$yr."_new_NoAdjust.nc";
         chomp( $urbanyr);
-        printf $fh_landuse_urban_timeseries $dynpft_format, $urbanyr, $yr; # Keer: I don't quiet understand the 'pft_override' option so I made no change to the "landuse_timeseries_override_$desc.txt"
+        printf $fh_landuse_timeseries $dynpft_format, $urbanyr, $yr; # Keer: I don't quiet understand the 'pft_override' option so I made no change to the "landuse_timeseries_override_$desc.txt"
         if ( $yr % 100 == 0 ) {
           print "year: $yr\n";
         }

From 71318fa36e1c578a7464095b62623eeaa716d9b9 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 17 Nov 2021 11:36:05 -0500
Subject: [PATCH 03/18] Add PCT_NATVEG for testing purpose

Add PCT_NATVEG for testing purpose
---
 tools/mksurfdata_map/src/mkfileMod.F90 | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 963e386348..bb854a9490 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -547,7 +547,11 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
        call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN_MAX', xtype=xtype, &
             lev1name='numurbl', &
             long_name='maximum percent urban for each density type', units='unitless')
-         
+
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_NATVEG', xtype=xtype, &
+            lev1name='time', &
+            long_name='total percent natural vegetation landunit', units='unitless')
+                          
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
           call ncd_def_spatial_var(ncid=ncid, varname=mkharvest_fieldname(ind1D(j),constant=.false.), xtype=xtype, &

From 88f822a43ce000ab7046931e686a9d4d5a4a2f19 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 17 Nov 2021 20:37:52 -0500
Subject: [PATCH 04/18] Delete variable PCT_NATVEG in landuse.timeseries file

---
 tools/mksurfdata_map/src/mkfileMod.F90     | 3 ---
 tools/mksurfdata_map/src/mksurfdat.F90     | 2 --
 tools/mksurfdata_map/src/mkurbanparMod.F90 | 2 +-
 3 files changed, 1 insertion(+), 6 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index bb854a9490..d77fe838fb 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -548,9 +548,6 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', &
             long_name='maximum percent urban for each density type', units='unitless')
 
-       call ncd_def_spatial_var(ncid=ncid, varname='PCT_NATVEG', xtype=xtype, &
-            lev1name='time', &
-            long_name='total percent natural vegetation landunit', units='unitless')
                           
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index d1ae730130..81a0f731cf 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1178,8 +1178,6 @@ program mksurfdat
           call update_max_array_urban(pcturb_max,urbn_classes_g)
 
           ! Output time-varying data for current year
-          call check_ret(nf_inq_varid(ncid, 'PCT_NATVEG', varid), subname)
-          call ncd_put_time_slice(ncid, varid, ntim, get_pct_l2g_array(pctnatpft))
           
           call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT', varid), subname)
           call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctnatpft))
diff --git a/tools/mksurfdata_map/src/mkurbanparMod.F90 b/tools/mksurfdata_map/src/mkurbanparMod.F90
index 35bbb0ad82..b47adfb61e 100644
--- a/tools/mksurfdata_map/src/mkurbanparMod.F90
+++ b/tools/mksurfdata_map/src/mkurbanparMod.F90
@@ -761,7 +761,7 @@ end subroutine mkurbanpar
 subroutine update_max_array_urban(pct_urbmax_arr,pct_urban_arr)
   !
   ! !DESCRIPTION:
-  ! Update the maximum percent cover of each urban class
+  ! Update the maximum percent cover of each urban class for landuse.timeseries file
   ! 
   ! !ARGUMENTS:
   real(r8)         , intent(inout):: pct_urbmax_arr(:,:)           ! max percent cover of each urban class

From 87b5fbddc8a89a92c73d3b18fd07cc74c9ac72cc Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Mon, 6 Dec 2021 15:49:00 -0500
Subject: [PATCH 05/18] Add hasurban

---
 tools/mksurfdata_map/src/mkfileMod.F90 |  3 +++
 tools/mksurfdata_map/src/mksurfdat.F90 | 14 +++++++++++++-
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index d77fe838fb..746d0b3229 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -548,6 +548,9 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', &
             long_name='maximum percent urban for each density type', units='unitless')
 
+       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=xtype, &
+            lev1name='numurbl', &
+            long_name='whether the grid ever has urban during the considered time period', units='unitless')
                           
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 81a0f731cf..03f64bad17 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -110,6 +110,7 @@ program mksurfdat
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
     real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
+    integer , allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
     real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
     real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: elev(:)            ! glc elevation (m)
@@ -449,6 +450,7 @@ program mksurfdat
                pctwet(ns_o)                       , & 
                pcturb(ns_o)                       , &
                pcturb_max(ns_o,numurbl)           , &
+               hasurban(ns_o,numurbl)             , &
                urban_region(ns_o)                 , &
                urbn_classes(ns_o,numurbl)         , &
                urbn_classes_g(ns_o,numurbl)       , &
@@ -478,6 +480,7 @@ program mksurfdat
     pctlak(:)             = spval
     pctwet(:)             = spval
     pcturb(:)             = spval
+    hasurban(:,:)         = 0
     urban_region(:)       = -999
     urbn_classes(:,:)     = spval
     urbn_classes_g(:,:)   = spval
@@ -1219,13 +1222,22 @@ program mksurfdat
 	  call check_ret(nf_sync(ncid), subname)
 
        end do   ! end of read loop
-
+       do n = 1,ns_o
+          do k =1, numurbl
+             if pcturb_max(n,k) > 1.e-6_r8 then 
+                hasurban(n,k) = 1
+             end if    
+          end do
+       end do        
        call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctnatpft_max)), subname)
 
        call check_ret(nf_inq_varid(ncid, 'PCT_CROP_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
+       call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
+       call check_ret(nf_put_var_double(ncid, varid, hasurban), subname)
+       
        call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)
        

From f506bc92cb9f76e0e2ffa7da39f9374ac0bed337 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Mon, 6 Dec 2021 15:51:44 -0500
Subject: [PATCH 06/18] Fix type (double--int)

---
 tools/mksurfdata_map/src/mksurfdat.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 03f64bad17..33aa782303 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1236,7 +1236,7 @@ program mksurfdat
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
        call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
-       call check_ret(nf_put_var_double(ncid, varid, hasurban), subname)
+       call check_ret(nf_put_var_int(ncid, varid, hasurban), subname)
        
        call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)

From 64abee13ba54736e07f8219789630b7c78424d14 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Mon, 6 Dec 2021 15:53:33 -0500
Subject: [PATCH 07/18] fix syntax error

---
 tools/mksurfdata_map/src/mksurfdat.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 33aa782303..f91cda67f8 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1224,7 +1224,7 @@ program mksurfdat
        end do   ! end of read loop
        do n = 1,ns_o
           do k =1, numurbl
-             if pcturb_max(n,k) > 1.e-6_r8 then 
+             if (pcturb_max(n,k) > 1.e-6_r8) then 
                 hasurban(n,k) = 1
              end if    
           end do

From b3a6e5a164a6aca59853eff8a083d92f8c460a0b Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Tue, 7 Dec 2021 12:04:16 -0500
Subject: [PATCH 08/18] Delete PCT_URBAN_MAX

---
 tools/mksurfdata_map/src/mkfileMod.F90     | 6 +-----
 tools/mksurfdata_map/src/mksurfdat.F90     | 5 +----
 tools/mksurfdata_map/src/mkurbanparMod.F90 | 2 +-
 3 files changed, 3 insertions(+), 10 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 746d0b3229..084848f5d8 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -544,11 +544,7 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', lev2name='time', &
             long_name='percent urban for each density type', units='unitless')
 
-       call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN_MAX', xtype=xtype, &
-            lev1name='numurbl', &
-            long_name='maximum percent urban for each density type', units='unitless')
-
-       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=xtype, &
+       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=nf_int, &
             lev1name='numurbl', &
             long_name='whether the grid ever has urban during the considered time period', units='unitless')
                           
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index f91cda67f8..3d4043de74 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1238,9 +1238,6 @@ program mksurfdat
        call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
        call check_ret(nf_put_var_int(ncid, varid, hasurban), subname)
        
-       call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
-       call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)
-       
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)
           call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctcft_max)), subname)
@@ -1449,7 +1446,7 @@ subroutine normalizencheck_landuse(ldomain)
        suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
        if (suma > (100._r8 + tol_loose)) then
           write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be'
-          write(6,*) '<= 100% before calling this subroutine'
+          write(6,*) '<= 100% before normalizing vegetated area'
           write(6,*) 'n, pctlak, pctwet, pcturb, pctgla = ', &
                n, pctlak(n), pctwet(n), pcturb(n), pctgla(n)
           call abort()
diff --git a/tools/mksurfdata_map/src/mkurbanparMod.F90 b/tools/mksurfdata_map/src/mkurbanparMod.F90
index b47adfb61e..4cbef42b11 100644
--- a/tools/mksurfdata_map/src/mkurbanparMod.F90
+++ b/tools/mksurfdata_map/src/mkurbanparMod.F90
@@ -778,7 +778,7 @@ subroutine update_max_array_urban(pct_urbmax_arr,pct_urban_arr)
         if (pct_urban_arr(n,k) > pct_urbmax_arr(n,k)) then 
            pct_urbmax_arr(n,k) = pct_urban_arr(n,k)
         end if
-      end do
+     end do
   end do
 
 end subroutine update_max_array_urban

From 0dc8ca414b277f4fc76c54ad8e36ad1712b706ed Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Tue, 7 Dec 2021 12:26:07 -0500
Subject: [PATCH 09/18] fixed a wrong indentation

---
 tools/mksurfdata_map/src/mksurfdat.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 3d4043de74..c3732d7751 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1174,7 +1174,7 @@ program mksurfdat
           call change_landuse(ldomain, dynpft=.true.)
 
           call normalizencheck_landuse(ldomain)
-	        call normalize_classes_by_gcell(urbn_classes, pcturb, urbn_classes_g)
+          call normalize_classes_by_gcell(urbn_classes, pcturb, urbn_classes_g)
           
           call update_max_array(pctnatpft_max,pctnatpft)
           call update_max_array(pctcft_max,pctcft)

From f57a3e4bf017021b3478595f4280268bd2dadedf Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 8 Dec 2021 10:04:53 -0500
Subject: [PATCH 10/18] Declare hasurban as a double variable

---
 tools/mksurfdata_map/src/mkfileMod.F90 |  2 +-
 tools/mksurfdata_map/src/mksurfdat.F90 | 10 ++++++----
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 084848f5d8..5cc6bfaf9a 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -544,7 +544,7 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', lev2name='time', &
             long_name='percent urban for each density type', units='unitless')
 
-       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=nf_int, &
+       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=xtype, &
             lev1name='numurbl', &
             long_name='whether the grid ever has urban during the considered time period', units='unitless')
                           
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index c3732d7751..1f124d9bc8 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -110,7 +110,7 @@ program mksurfdat
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
     real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
-    integer , allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
+    real(r8), allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
     real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
     real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: elev(:)            ! glc elevation (m)
@@ -480,7 +480,7 @@ program mksurfdat
     pctlak(:)             = spval
     pctwet(:)             = spval
     pcturb(:)             = spval
-    hasurban(:,:)         = 0
+    hasurban(:,:)         = spval
     urban_region(:)       = -999
     urbn_classes(:,:)     = spval
     urbn_classes_g(:,:)   = spval
@@ -1225,7 +1225,9 @@ program mksurfdat
        do n = 1,ns_o
           do k =1, numurbl
              if (pcturb_max(n,k) > 1.e-6_r8) then 
-                hasurban(n,k) = 1
+                hasurban(n,k) = 1._r8
+             else 
+                hasurban(n,k) = 0._r8             
              end if    
           end do
        end do        
@@ -1236,7 +1238,7 @@ program mksurfdat
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
        call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
-       call check_ret(nf_put_var_int(ncid, varid, hasurban), subname)
+       call check_ret(nf_put_var_double(ncid, varid, hasurban), subname)
        
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)

From da0814061ce7057f04474885c48dd22aceda58cf Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Fri, 10 Dec 2021 12:07:38 -0500
Subject: [PATCH 11/18] Change HASURBAN type as integer

---
 tools/mksurfdata_map/src/mkfileMod.F90 |  2 +-
 tools/mksurfdata_map/src/mksurfdat.F90 | 10 +++++-----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 5cc6bfaf9a..084848f5d8 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -544,7 +544,7 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', lev2name='time', &
             long_name='percent urban for each density type', units='unitless')
 
-       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=xtype, &
+       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=nf_int, &
             lev1name='numurbl', &
             long_name='whether the grid ever has urban during the considered time period', units='unitless')
                           
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 1f124d9bc8..81dc4d9914 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -110,7 +110,7 @@ program mksurfdat
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
     real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
-    real(r8), allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
+    integer , allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
     real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
     real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: elev(:)            ! glc elevation (m)
@@ -480,7 +480,7 @@ program mksurfdat
     pctlak(:)             = spval
     pctwet(:)             = spval
     pcturb(:)             = spval
-    hasurban(:,:)         = spval
+    hasurban(:,:)         = 0
     urban_region(:)       = -999
     urbn_classes(:,:)     = spval
     urbn_classes_g(:,:)   = spval
@@ -1225,9 +1225,9 @@ program mksurfdat
        do n = 1,ns_o
           do k =1, numurbl
              if (pcturb_max(n,k) > 1.e-6_r8) then 
-                hasurban(n,k) = 1._r8
+                hasurban(n,k) = 1
              else 
-                hasurban(n,k) = 0._r8             
+                hasurban(n,k) = 0             
              end if    
           end do
        end do        
@@ -1238,7 +1238,7 @@ program mksurfdat
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
        call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
-       call check_ret(nf_put_var_double(ncid, varid, hasurban), subname)
+       call check_ret(nf_put_var_int(ncid, varid, hasurban), subname)
        
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)

From 5316eb54c5285e8a9a32a383bcaef2fd51ca7596 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Thu, 16 Dec 2021 14:58:42 -0500
Subject: [PATCH 12/18] Add some comments to emphasize hard coding

---
 tools/mksurfdata_map/mksurfdata.pl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tools/mksurfdata_map/mksurfdata.pl b/tools/mksurfdata_map/mksurfdata.pl
index 3caa4eb009..64dbaed738 100755
--- a/tools/mksurfdata_map/mksurfdata.pl
+++ b/tools/mksurfdata_map/mksurfdata.pl
@@ -286,8 +286,8 @@ sub write_transient_timeseries_file {
         printf $fh_landuse_timeseries $dynpft_format, $hrvtypyr, $yr;
         my $urbanyr = "/glade/scratch/keerzhang/archive/BNU_NoAdjust/05deg_".$yr."_new_NoAdjust.nc";
         chomp( $urbanyr);
-        printf $fh_landuse_timeseries $dynpft_format, $urbanyr, $yr; # Keer: I don't quiet understand the 'pft_override' option so I made no change to the "landuse_timeseries_override_$desc.txt"
-        if ( $yr % 100 == 0 ) {
+        printf $fh_landuse_timeseries $dynpft_format, $urbanyr, $yr; # I hard coded this part just to generate a txt file with urban raw data file locations
+        if ( $yr % 100 == 0 ) {                                      # And note that I made no change to the "landuse_timeseries_override_$desc.txt" because I am not sure how to deal with the the 'pft_override' option
           print "year: $yr\n";
         }
       }

From b3d654a187cf93d9bb79ac6468bcf1309d0485bc Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Thu, 16 Dec 2021 15:16:58 -0500
Subject: [PATCH 13/18] Change the pcturb_max threshold for HASURBAN

---
 tools/mksurfdata_map/src/mksurfdat.F90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 81dc4d9914..a753298f44 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1222,9 +1222,9 @@ program mksurfdat
 	  call check_ret(nf_sync(ncid), subname)
 
        end do   ! end of read loop
-       do n = 1,ns_o
+       do n = 1, ns_o
           do k =1, numurbl
-             if (pcturb_max(n,k) > 1.e-6_r8) then 
+             if (pcturb_max(n,k) > 1.e-3_r8) then 
                 hasurban(n,k) = 1
              else 
                 hasurban(n,k) = 0             

From 6f87e1cfde54be72453731c5c7c786a55babf993 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Tue, 21 Dec 2021 09:51:35 -0500
Subject: [PATCH 14/18] Delete HASURBAN

---
 tools/mksurfdata_map/src/mkfileMod.F90 |  4 ++--
 tools/mksurfdata_map/src/mksurfdat.F90 | 17 +++--------------
 2 files changed, 5 insertions(+), 16 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 084848f5d8..2a453b056a 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -544,9 +544,9 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name='numurbl', lev2name='time', &
             long_name='percent urban for each density type', units='unitless')
 
-       call ncd_def_spatial_var(ncid=ncid, varname='HASURBAN', xtype=nf_int, &
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN_MAX', xtype=xtype, &
             lev1name='numurbl', &
-            long_name='whether the grid ever has urban during the considered time period', units='unitless')
+            long_name='maximum percent urban for each density type', units='unitless')
                           
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index a753298f44..ee9dff6831 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -110,7 +110,6 @@ program mksurfdat
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
     real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
-    integer , allocatable  :: hasurban(:,:)      ! whether the urban class should be active in transient urban simulations
     real(r8), allocatable  :: urbn_classes(:,:)  ! percent cover of each urban class, as % of total urban area
     real(r8), allocatable  :: urbn_classes_g(:,:)! percent cover of each urban class, as % of grid cell
     real(r8), allocatable  :: elev(:)            ! glc elevation (m)
@@ -450,7 +449,6 @@ program mksurfdat
                pctwet(ns_o)                       , & 
                pcturb(ns_o)                       , &
                pcturb_max(ns_o,numurbl)           , &
-               hasurban(ns_o,numurbl)             , &
                urban_region(ns_o)                 , &
                urbn_classes(ns_o,numurbl)         , &
                urbn_classes_g(ns_o,numurbl)       , &
@@ -480,7 +478,6 @@ program mksurfdat
     pctlak(:)             = spval
     pctwet(:)             = spval
     pcturb(:)             = spval
-    hasurban(:,:)         = 0
     urban_region(:)       = -999
     urbn_classes(:,:)     = spval
     urbn_classes_g(:,:)   = spval
@@ -1222,23 +1219,15 @@ program mksurfdat
 	  call check_ret(nf_sync(ncid), subname)
 
        end do   ! end of read loop
-       do n = 1, ns_o
-          do k =1, numurbl
-             if (pcturb_max(n,k) > 1.e-3_r8) then 
-                hasurban(n,k) = 1
-             else 
-                hasurban(n,k) = 0             
-             end if    
-          end do
-       end do        
+       
        call check_ret(nf_inq_varid(ncid, 'PCT_NAT_PFT_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctnatpft_max)), subname)
 
        call check_ret(nf_inq_varid(ncid, 'PCT_CROP_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, get_pct_l2g_array(pctcft_max)), subname)
 
-       call check_ret(nf_inq_varid(ncid, 'HASURBAN', varid), subname)
-       call check_ret(nf_put_var_int(ncid, varid, hasurban), subname)
+       call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
+       call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)
        
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)

From b51b3b93a244ef2de942d62e3ea39eb16db60239 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Tue, 21 Dec 2021 09:54:36 -0500
Subject: [PATCH 15/18] Delete pcturb_orig

---
 tools/mksurfdata_map/src/mksurfdat.F90 | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index ee9dff6831..c74e524e50 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -138,7 +138,6 @@ program mksurfdat
     real(r8), allocatable  :: lakedepth(:)       ! lake depth (m)
     real(r8), allocatable  :: pctlak_orig(:)     ! percent lake of gridcell before dynamic land use adjustments            
     real(r8), allocatable  :: pctwet_orig(:)     ! percent wetland of gridcell before dynamic land use adjustments 
-    real(r8), allocatable  :: pcturb_orig(:)     ! percent urban of gridcell before dynamic land use adjustments 
     real(r8), allocatable  :: pctgla_orig(:)     ! percent glacier of gridcell before dynamic land use adjustments
 
     real(r8) :: std_elev = -999.99_r8            ! Standard deviation of elevation (m) to use for entire grid
@@ -469,7 +468,6 @@ program mksurfdat
                glacier_region(ns_o)               , &  
                pctlak_orig(ns_o)                  , &
                pctwet_orig(ns_o)                  , &
-               pcturb_orig(ns_o)                  , &
                pctgla_orig(ns_o)                  )       
     landfrac_pft(:)       = spval 
     pctlnd_pft(:)         = spval
@@ -658,7 +656,6 @@ program mksurfdat
     ! Save special land unit areas of surface dataset 
     pctlak_orig(:) = pctlak(:)
     pctwet_orig(:) = pctwet(:)
-    pcturb_orig(:) = pcturb(:)
     pctgla_orig(:) = pctgla(:)
      
     ! Make elevation [elev] from [ftopo, ffrac] dataset

From 09b35c5fa6e995d546966889ce43bcd95b7747f2 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 22 Dec 2021 11:11:10 -0500
Subject: [PATCH 16/18] Merge #1073 with this PR #1579(Add dynamic lake
 variables)

---
 tools/mksurfdata_map/mksurfdata.pl       |  9 ++++--
 tools/mksurfdata_map/src/mkfileMod.F90   |  8 ++++-
 tools/mksurfdata_map/src/mklanwatMod.F90 | 33 +++++++++++++++++---
 tools/mksurfdata_map/src/mksurfdat.F90   | 38 +++++++++++++++++-------
 4 files changed, 70 insertions(+), 18 deletions(-)

diff --git a/tools/mksurfdata_map/mksurfdata.pl b/tools/mksurfdata_map/mksurfdata.pl
index 64dbaed738..a73f3fdbab 100755
--- a/tools/mksurfdata_map/mksurfdata.pl
+++ b/tools/mksurfdata_map/mksurfdata.pl
@@ -285,9 +285,12 @@ sub write_transient_timeseries_file {
         chomp( $hrvtypyr );
         printf $fh_landuse_timeseries $dynpft_format, $hrvtypyr, $yr;
         my $urbanyr = "/glade/scratch/keerzhang/archive/BNU_NoAdjust/05deg_".$yr."_new_NoAdjust.nc";
-        chomp( $urbanyr);
-        printf $fh_landuse_timeseries $dynpft_format, $urbanyr, $yr; # I hard coded this part just to generate a txt file with urban raw data file locations
-        if ( $yr % 100 == 0 ) {                                      # And note that I made no change to the "landuse_timeseries_override_$desc.txt" because I am not sure how to deal with the the 'pft_override' option
+        chomp( $urbanyr);                                            # !KZ I hard coded this part just to generate a txt file with urban raw data file locations
+        printf $fh_landuse_timeseries $dynpft_format, $urbanyr, $yr; # And note that I made no change to the "landuse_timeseries_override_$desc.txt" because I am not sure how to deal with the the 'pft_override' option
+        my $lakeyr = "/glade/work/ivanderk/inputdata/rawdata/timeseries/mksurf_lake_0.05x0.05_hist_clm5_hydrolakes_2017_c20200305.nc";
+        chomp( $lakeyr);
+        printf $fh_landuse_timeseries $dynpft_format, $lakeyr, $yr;          
+        if ( $yr % 100 == 0 ) {
           print "year: $yr\n";
         }
       }
diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 2a453b056a..69ff8f58d4 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -547,7 +547,13 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
        call ncd_def_spatial_var(ncid=ncid, varname='PCT_URBAN_MAX', xtype=xtype, &
             lev1name='numurbl', &
             long_name='maximum percent urban for each density type', units='unitless')
-                          
+
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_LAKE', xtype=xtype, &
+            lev1name="time", long_name='percent lake', units='unitless') 
+
+       call ncd_def_spatial_var(ncid=ncid, varname='PCT_LAKE_MAX', xtype=xtype, &
+            long_name='maximum percent lake for each density type', units='unitless')
+                                           
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()
           call ncd_def_spatial_var(ncid=ncid, varname=mkharvest_fieldname(ind1D(j),constant=.false.), xtype=xtype, &
diff --git a/tools/mksurfdata_map/src/mklanwatMod.F90 b/tools/mksurfdata_map/src/mklanwatMod.F90
index 49a1485fa7..66159d4bd8 100644
--- a/tools/mksurfdata_map/src/mklanwatMod.F90
+++ b/tools/mksurfdata_map/src/mklanwatMod.F90
@@ -24,10 +24,10 @@ module mklanwatMod
   private
 
 ! !PUBLIC MEMBER FUNCTIONS:
-  public mklakwat           ! make % lake
-  public mkwetlnd           ! make % wetland
-  public mklakparams        ! make lake parameters
-
+  public mklakwat                ! make % lake
+  public mkwetlnd                ! make % wetland
+  public mklakparams             ! make lake parameters
+  public update_max_array_lake   ! Update the maximum lake percent
 !EOP
 !===============================================================
 contains
@@ -499,5 +499,30 @@ subroutine mklakparams(ldomain, mapfname, datfname, ndiag, &
   call shr_sys_flush(6)
 
 end subroutine mklakparams
+!------------------------------------------------------------------------------
+
+!-----------------------------------------------------------------------
+subroutine update_max_array_lake(pct_lakmax_arr,pct_lake_arr)
+  !
+  ! !DESCRIPTION:
+  ! Update the maximum lake percent for landuse.timeseries file
+  ! 
+  ! !ARGUMENTS:
+  real(r8)         , intent(inout):: pct_lakmax_arr(:)       ! max lake percent
+  real(r8)         , intent(in):: pct_lake_arr(:)            ! lake percent that is used to update the old pct_lakmax_arr
+  !
+  ! !LOCAL VARIABLES:
+  integer :: n,ns              ! indices
+
+  character(len=*), parameter :: subname = 'update_max_array_lake'
+  !-----------------------------------------------------------------------
+  ns = size(pct_lake_arr,1)
+  do n = 1, ns
+     if (pct_lake_arr(n) > pct_lakmax_arr(n)) then 
+        pct_lakmax_arr(n) = pct_lake_arr(n)
+     end if
+  end do
+
+end subroutine update_max_array_lake
 
 end module mklanwatMod
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index c74e524e50..2e9ffaab64 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -22,7 +22,7 @@ program mksurfdat
                                     soil_color, mksoilcol, mkorganic, &
                                     soil_fmax, mkfmax
     use mkvocefMod         , only : mkvocef
-    use mklanwatMod        , only : mklakwat, mkwetlnd, mklakparams
+    use mklanwatMod        , only : mklakwat, mkwetlnd, mklakparams, update_max_array_lake
     use mkglacierregionMod , only : mkglacierregion
     use mkglcmecMod        , only : nglcec, mkglcmec, mkglcmecInit, mkglacier
     use mkharvestMod       , only : mkharvest, mkharvest_init, mkharvest_fieldname
@@ -81,6 +81,7 @@ program mksurfdat
     character(len=256) :: fname             ! generic filename
     character(len=256) :: fhrvname          ! generic harvest filename
     character(len=256) :: furbname          ! generic transient urban land cover filename
+    character(len=256) :: flakname          ! generic lake filename
     character(len=256) :: string            ! string read in
     integer  :: t1                          ! timer
     real(r8),parameter :: p5  = 0.5_r8      ! constant
@@ -107,6 +108,7 @@ program mksurfdat
     real(r8), allocatable  :: elevclass(:)       ! glacier_mec elevation classes
     integer,  allocatable  :: glacier_region(:)  ! glacier region ID
     real(r8), allocatable  :: pctlak(:)          ! percent of grid cell that is lake
+    real(r8), allocatable  :: pctlak_max(:)      ! maximum percent of grid cell that is lake  
     real(r8), allocatable  :: pctwet(:)          ! percent of grid cell that is wetland  
     real(r8), allocatable  :: pcturb(:)          ! percent of grid cell that is urbanized (total across all urban classes)
     real(r8), allocatable  :: pcturb_max(:,:)    ! maximum percent cover of each urban class, as % of grid cell
@@ -135,8 +137,7 @@ program mksurfdat
     real(r8), allocatable  :: vic_ws(:)          ! VIC Ws parameter (unitless)
     real(r8), allocatable  :: vic_dsmax(:)       ! VIC Dsmax parameter (mm/day)
     real(r8), allocatable  :: vic_ds(:)          ! VIC Ds parameter (unitless)
-    real(r8), allocatable  :: lakedepth(:)       ! lake depth (m)
-    real(r8), allocatable  :: pctlak_orig(:)     ! percent lake of gridcell before dynamic land use adjustments            
+    real(r8), allocatable  :: lakedepth(:)       ! lake depth (m)            
     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
 
@@ -282,7 +283,7 @@ program mksurfdat
     ! ======================================
     ! Optionally specify setting for:
     ! ======================================
-    !    mksrf_fdynuse ----- ASCII text file that lists each year of pft and urban files to use
+    !    mksrf_fdynuse ----- ASCII text file that lists each year of pft, urban, and lake files to use
     !    mksrf_gridtype ---- Type of grid (default is 'global')
     !    outnc_double ------ If output should be in double precision
     !    outnc_large_files - If output should be in NetCDF large file format
@@ -444,7 +445,8 @@ program mksurfdat
                pctcft(ns_o)                       , &
                pctcft_max(ns_o)                   , &
                pctgla(ns_o)                       , & 
-               pctlak(ns_o)                       , & 
+               pctlak(ns_o)                       , &
+               pctlak_max(ns_o)                   , & 
                pctwet(ns_o)                       , & 
                pcturb(ns_o)                       , &
                pcturb_max(ns_o,numurbl)           , &
@@ -466,7 +468,6 @@ program mksurfdat
                vic_ds(ns_o)                       , &
                lakedepth(ns_o)                    , &
                glacier_region(ns_o)               , &  
-               pctlak_orig(ns_o)                  , &
                pctwet_orig(ns_o)                  , &
                pctgla_orig(ns_o)                  )       
     landfrac_pft(:)       = spval 
@@ -654,7 +655,6 @@ program mksurfdat
          region_o=urban_region)
 
     ! Save special land unit areas of surface dataset 
-    pctlak_orig(:) = pctlak(:)
     pctwet_orig(:) = pctwet(:)
     pctgla_orig(:) = pctgla(:)
      
@@ -1085,6 +1085,7 @@ program mksurfdat
        pctnatpft_max = pctnatpft
        pctcft_max = pctcft
        pcturb_max = urbn_classes_g
+       pctlak_max = pctlak
 
        ntim = 0
        do 
@@ -1120,7 +1121,14 @@ program mksurfdat
              write(6,*) subname, ' error: year for urban not equal to year for PFT files'
              call abort()
           end if
-      write(6,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname)    
+      write(6,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname) 
+          read(nfdyn, '(A195,1x,I4)', iostat=ier) flakname, year2
+          if ( year2 /= year ) then
+             write(6,*) subname, ' error: year for lake not equal to year for PFT files'
+             call abort()
+          end if
+      write(6,*)'input lake dynamic dataset for year ', year2, ' is : ', trim(flakname)    
+           
           ntim = ntim + 1
 
           ! Create pctpft data at model resolution
@@ -1149,7 +1157,10 @@ program mksurfdat
              end if
           end do
 
-
+          ! Create pctlak data at model resolution (use original mapping file from lake data)
+          call mklakwat (ldomain, mapfname=map_flakwat, datfname=flakname, &
+               ndiag=ndiag, zero_out=all_urban.or.all_veg, lake_o=pctlak)
+               
           call mkurban (ldomain, mapfname=map_furban, datfname=furbname, &
                ndiag=ndiag, zero_out=all_veg, urbn_o=pcturb, urbn_classes_o=urbn_classes, &
                region_o=urban_region)
@@ -1162,7 +1173,6 @@ program mksurfdat
 
           ! For landunits NOT read each year: reset to their pre-adjustment values in preparation for redoing landunit area normalization
           pctwet(:) = pctwet_orig(:)
-          pctlak(:) = pctlak_orig(:)
           pctgla(:) = pctgla_orig(:)
                                         
           call change_landuse(ldomain, dynpft=.true.)
@@ -1173,6 +1183,7 @@ program mksurfdat
           call update_max_array(pctnatpft_max,pctnatpft)
           call update_max_array(pctcft_max,pctcft)
           call update_max_array_urban(pcturb_max,urbn_classes_g)
+          call update_max_array_lake(pctlak_max,pctlak)
 
           ! Output time-varying data for current year
           
@@ -1184,6 +1195,10 @@ program mksurfdat
           
           call check_ret(nf_inq_varid(ncid, 'PCT_URBAN', varid), subname)
           call ncd_put_time_slice(ncid, varid, ntim, urbn_classes_g)
+          
+          call check_ret(nf_inq_varid(ncid, 'PCT_LAKE', varid), subname)
+          call ncd_put_time_slice(ncid, varid, ntim, pctlak)
+                    
           if (num_cft > 0) then
              call check_ret(nf_inq_varid(ncid, 'PCT_CFT', varid), subname)
              call ncd_put_time_slice(ncid, varid, ntim, get_pct_p2l_array(pctcft))
@@ -1226,6 +1241,9 @@ program mksurfdat
        call check_ret(nf_inq_varid(ncid, 'PCT_URBAN_MAX', varid), subname)
        call check_ret(nf_put_var_double(ncid, varid, pcturb_max), subname)
        
+       call check_ret(nf_inq_varid(ncid, 'PCT_LAKE_MAX', varid), subname)
+       call check_ret(nf_put_var_double(ncid, varid, pctlak_max), subname)       
+       
        if (num_cft > 0) then
           call check_ret(nf_inq_varid(ncid, 'PCT_CFT_MAX', varid), subname)
           call check_ret(nf_put_var_double(ncid, varid, get_pct_p2l_array(pctcft_max)), subname)

From fc7119da410c9f106f2239f1a4d435cff2470b54 Mon Sep 17 00:00:00 2001
From: keerzhang1 <86242752+keerzhang1@users.noreply.github.com>
Date: Wed, 22 Dec 2021 14:25:17 -0500
Subject: [PATCH 17/18] Correct copy and paste error

---
 tools/mksurfdata_map/src/mkfileMod.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tools/mksurfdata_map/src/mkfileMod.F90 b/tools/mksurfdata_map/src/mkfileMod.F90
index 69ff8f58d4..abba061b69 100644
--- a/tools/mksurfdata_map/src/mkfileMod.F90
+++ b/tools/mksurfdata_map/src/mkfileMod.F90
@@ -552,7 +552,7 @@ subroutine mkfile(domain, fname, harvdata, dynlanduse)
             lev1name="time", long_name='percent lake', units='unitless') 
 
        call ncd_def_spatial_var(ncid=ncid, varname='PCT_LAKE_MAX', xtype=xtype, &
-            long_name='maximum percent lake for each density type', units='unitless')
+            long_name='maximum percent lake', units='unitless')
                                            
        call harvdata%getFieldsIdx( ind1D, ind2D )
        do j = 1, harvdata%num1Dfields()

From 4b3ef5a77bf514e66bd6c7a5b72947e1e009ca75 Mon Sep 17 00:00:00 2001
From: Keith Oleson <oleson@ucar.edu>
Date: Tue, 28 Dec 2021 14:50:38 -0700
Subject: [PATCH 18/18] Incorporate changes from PR1564 into PR1579

---
 tools/mksurfdata_map/src/mkpftUtilsMod.F90 |  40 ----
 tools/mksurfdata_map/src/mksurfdat.F90     | 203 ++++++---------------
 2 files changed, 52 insertions(+), 191 deletions(-)

diff --git a/tools/mksurfdata_map/src/mkpftUtilsMod.F90 b/tools/mksurfdata_map/src/mkpftUtilsMod.F90
index 4a9ea12f97..0faede43a1 100644
--- a/tools/mksurfdata_map/src/mkpftUtilsMod.F90
+++ b/tools/mksurfdata_map/src/mkpftUtilsMod.F90
@@ -24,7 +24,6 @@ module mkpftUtilsMod
   ! !PUBLIC MEMBER FUNCTIONS:
   !
   public :: convert_from_p2g      ! Convert a p2g array into pct_pft_type objects
-  public :: adjust_total_veg_area ! Adjust the total vegetated area (natural veg & crop) to a new specified total
 
   !
   ! !PRIVATE MEMBER FUNCTIONS:
@@ -213,45 +212,6 @@ function get_default_cft() result(default_cft)
 
   end function get_default_cft
 
-
-  !-----------------------------------------------------------------------
-  subroutine adjust_total_veg_area(new_total_pct, pctnatpft, pctcft)
-    !
-    ! !DESCRIPTION:
-    ! Adjust the total vegetated area on the grid cell (natural veg & crop) to a new
-    ! specified total.
-    !
-    ! If the old areas are 0%, then all the new area goes into pctnatpft.
-    !
-    ! !USES:
-    use mkpctPftTypeMod, only : pct_pft_type
-    !
-    ! !ARGUMENTS:
-    real(r8), intent(in) :: new_total_pct ! new total % of natural veg + crop landunits
-    class(pct_pft_type), intent(inout) :: pctnatpft ! natural veg cover information
-    class(pct_pft_type), intent(inout) :: pctcft    ! crop cover information
-    !
-    ! !LOCAL VARIABLES:
-    real(r8) :: natpft_l2g ! grid cell % cover of nat. veg.
-    real(r8) :: cft_l2g    ! grid cell % cover of crop
-    real(r8) :: old_total  ! old total % cover of natural veg + crop landunits
-    
-    character(len=*), parameter :: subname = 'adjust_total_veg_area'
-    !-----------------------------------------------------------------------
-    
-    natpft_l2g = pctnatpft%get_pct_l2g()
-    cft_l2g = pctcft%get_pct_l2g()
-    old_total = natpft_l2g + cft_l2g
-    if (old_total > 0._r8) then
-       call pctnatpft%set_pct_l2g(natpft_l2g * new_total_pct / old_total)
-       call pctcft%set_pct_l2g(cft_l2g * new_total_pct / old_total)
-    else
-       call pctnatpft%set_pct_l2g(new_total_pct)
-    end if
-
-  end subroutine adjust_total_veg_area
-
-
 end module mkpftUtilsMod
 
 
diff --git a/tools/mksurfdata_map/src/mksurfdat.F90 b/tools/mksurfdata_map/src/mksurfdat.F90
index 32390a5188..ba51a4480c 100644
--- a/tools/mksurfdata_map/src/mksurfdat.F90
+++ b/tools/mksurfdata_map/src/mksurfdat.F90
@@ -1120,13 +1120,14 @@ program mksurfdat
              write(6,*) subname, ' error: year for urban not equal to year for PFT files'
              call abort()
           end if
-      write(6,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname) 
+          write(6,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname)
+          ! Read input lake data
           read(nfdyn, '(A195,1x,I4)', iostat=ier) flakname, year2
           if ( year2 /= year ) then
              write(6,*) subname, ' error: year for lake not equal to year for PFT files'
              call abort()
           end if
-      write(6,*)'input lake dynamic dataset for year ', year2, ' is : ', trim(flakname)    
+          write(6,*)'input lake dynamic dataset for year ', year2, ' is : ', trim(flakname)
            
           ntim = ntim + 1
 
@@ -1170,7 +1171,8 @@ program mksurfdat
                end where
           end if
 
-          ! For landunits NOT read each year: reset to their pre-adjustment values in preparation for redoing landunit area normalization
+          ! For landunits NOT read each year: reset to their pre-adjustment values in 
+          ! preparation for redoing landunit area normalization
           pctwet(:) = pctwet_orig(:)
           pctgla(:) = pctgla_orig(:)
                                         
@@ -1344,8 +1346,6 @@ subroutine normalizencheck_landuse(ldomain)
 ! 
 !
 ! !USES:
-    use mkpftConstantsMod , only : baregroundindex
-    use mkpftUtilsMod     , only : adjust_total_veg_area
     implicit none
 ! !ARGUMENTS:
     type(domain_type)   :: ldomain
@@ -1360,13 +1360,7 @@ subroutine normalizencheck_landuse(ldomain)
     integer  :: nsmall                      ! number of small PFT values for a single check
     integer  :: nsmall_tot                  ! total number of small PFT values in all grid cells
     real(r8) :: suma                        ! sum for error check
-    real(r8) :: suma2                       ! another sum for error check
-    real(r8) :: new_total_veg_pct           ! new % veg (% of grid cell, total of natural veg & crop)
-    real(r8) :: bare_pct_p2g                ! % of bare soil, as % of grid cell
-    real(r8) :: bare_urb_diff               ! difference between bare soil and urban %
-    real(r8) :: pcturb_excess               ! excess urban % not accounted for by bare soil
-    real(r8) :: sum8, sum8a                 ! sum for error check
-    real(r4) :: sum4a                       ! sum for error check
+    real(r8) :: new_total_natveg_pct        ! new % veg (% of grid cell, natural veg)
     real(r8), parameter :: tol_loose = 1.e-4_r8               ! tolerance for some 'loose' error checks
     real(r8), parameter :: toosmallPFT = 1.e-10_r8            ! tolerance for PFT's to ignore
     character(len=32) :: subname = 'normalizencheck_landuse'  ! subroutine name
@@ -1409,21 +1403,17 @@ subroutine normalizencheck_landuse(ldomain)
           pftdata_mask(n) = 1
        end if
  
-       ! Make sure sum of land cover types does not exceed 100. If it does,
-       ! subtract excess from most dominant land cover.
+       ! Make sure sum of all land cover types except natural vegetation
+       ! does not exceed 100. If it does, subtract excess from these land cover
+       ! types proportionally.
 
-       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
-       if (suma > 250._r4) then
-          write (6,*) subname, ' error: sum of pctlak, pctwet,', &
-               'pcturb and pctgla is greater than 250%'
-          write (6,*)'n,pctlak,pctwet,pcturb,pctgla= ', &
-               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n)
-          call abort()
-       else if (suma > 100._r4) then
+       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
+       if (suma > 100._r4) then
           pctlak(n) = pctlak(n) * 100._r8/suma
           pctwet(n) = pctwet(n) * 100._r8/suma
           pcturb(n) = pcturb(n) * 100._r8/suma
           pctgla(n) = pctgla(n) * 100._r8/suma
+          call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma)
        end if
 
        ! Check preconditions
@@ -1447,81 +1437,33 @@ subroutine normalizencheck_landuse(ldomain)
           write(6,*) 'n, pctgla = ', n, pctgla(n)
           call abort()
        end if
-
-       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
-       if (suma > (100._r8 + tol_loose)) then
-          write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be'
-          write(6,*) '<= 100% before normalizing vegetated area'
-          write(6,*) 'n, pctlak, pctwet, pcturb, pctgla = ', &
-               n, pctlak(n), pctwet(n), pcturb(n), pctgla(n)
+       if ( pctcft(n)%get_pct_l2g() < 0.0_r8 )then
+          write(6,*) subname, ' ERROR: pctcrop is negative!'
+          write(6,*) 'n, pctcrop = ', n, pctcft(n)%get_pct_l2g()
           call abort()
        end if
-       ! ------------------------------------------------------------------------
-       ! Normalize vegetated area so that vegetated + special area is 100%
-       ! ------------------------------------------------------------------------
-       ! First normalize vegetated (natural veg + crop) cover so that the total of
-       ! (vegetated + (special excluding urban)) is 100%. We'll deal with urban later.
-       !
-       ! Note that, in practice, the total area of natural veg + crop is typically 100%
-       ! going into this routine. However, the following code does NOT rely on this, and
-       ! will work properly regardless of the initial area of natural veg + crop (even if
-       ! that initial area is 0%).
-       
-       suma = pctlak(n)+pctwet(n)+pctgla(n)
-       new_total_veg_pct = 100._r8 - suma
-       ! correct for rounding error:
-       new_total_veg_pct = max(new_total_veg_pct, 0._r8)
 
-       call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n))
-
-       ! Make sure we did the above rescaling correctly
-
-       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
-       if (abs(suma - 100._r8) > tol_loose) then
-          write(6,*) subname, ' ERROR in rescaling veg based on (special excluding urban'
-          write(6,*) 'suma = ', suma
+       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
+       if (suma > (100._r8 + tol_loose)) then
+          write(6,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla + pctcrop must be'
+          write(6,*) '<= 100% before normalizing natural vegetation area'
+          write(6,*) 'n, pctlak, pctwet, pcturb, pctgla, pctcrop = ', &
+               n, pctlak(n), pctwet(n), pcturb(n), pctgla(n), pctcft(n)%get_pct_l2g()
           call abort()
        end if
 
-       ! Now decrease the vegetated area to account for urban area. Urban needs to be
-       ! handled specially because we replace bare soil preferentially with urban, rather
-       ! than rescaling all PFTs equally.
+       ! Natural vegetated cover is 100% minus the sum of all other landunits
 
-       if (pcturb(n) > 0._r8) then
-          
-          ! Replace bare soil preferentially with urban
-          bare_pct_p2g = pctnatpft(n)%get_one_pct_p2g(baregroundindex)
-          bare_urb_diff = bare_pct_p2g - pcturb(n)
-          bare_pct_p2g = max(0._r8, bare_urb_diff)
-          call pctnatpft(n)%set_one_pct_p2g(baregroundindex, bare_pct_p2g)
-          pcturb_excess = abs(min(0._r8,bare_urb_diff))
-          
-          ! For any urban not accounted for by bare soil, replace other PFTs
-          ! proportionally
-          if (pcturb_excess > 0._r8) then
-             ! Note that, in this case, we will have already reduced bare ground to 0%
-
-             new_total_veg_pct = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g() - pcturb_excess
-             if (new_total_veg_pct < 0._r8) then
-                if (abs(new_total_veg_pct) < tol_loose) then
-                   ! only slightly less than 0; correct it
-                   new_total_veg_pct = 0._r8
-                else
-                   write(6,*) subname, ' ERROR: trying to replace veg with urban,'
-                   write(6,*) 'but pcturb_excess exceeds current vegetation percent'
-                   call abort()
-                end if
-             end if
-
-             call adjust_total_veg_area(new_total_veg_pct, pctnatpft=pctnatpft(n), pctcft=pctcft(n))
-          end if
-
-       end if ! pcturb(n) > 0
+       new_total_natveg_pct = 100._r8 - suma
+       ! correct for rounding error:
+       new_total_natveg_pct = max(new_total_natveg_pct, 0._r8)
 
+       call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct)
+       
        ! Confirm that we have done the rescaling correctly: now the sum of all landunits
-       ! should be 100%
-       suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)
-       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
+       ! should be 100% within tol_loose
+       suma = pctlak(n)+pctwet(n)+pctgla(n)+pcturb(n)+pctcft(n)%get_pct_l2g()
+       suma = suma + pctnatpft(n)%get_pct_l2g()
        if (abs(suma - 100._r8) > tol_loose) then
           write(6,*) subname, ' ERROR: landunits do not sum to 100%'
           write(6,*) 'n, suma, pctlak, pctwet, pctgla, pcturb, pctnatveg, pctcrop = '
@@ -1558,30 +1500,34 @@ subroutine normalizencheck_landuse(ldomain)
           call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() * 100._r8/suma)
        end if
        
-       ! Roundoff error fix
-       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
-       suma2 = pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
+       ! This roundoff error fix is needed to handle the situation where new_total_natveg_pct
+       ! ends up near 0 but not exactly 0 due to roundoff issues. In this situation, we set the
+       ! natveg landunit area to exactly 0 and put the remainder into some other landunit. Since
+       ! the remainder is very small, it doesn't really matter which other landunit we add it to,
+       ! so we just pick some landunit that already has at least 1% cover.
+       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
        if ( (suma < 100._r8 .and. suma > (100._r8 - 1.e-6_r8)) .or. &
-            (suma2 > 0.0_r8 .and. suma2 <  1.e-6_r8) ) then
-          write (6,*) 'Special land units near 100%, but not quite for n,suma =',n,suma
-          write (6,*) 'Adjusting special land units to 100%'
-          if (pctlak(n) >= 25._r8) then
-             pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n))
-          else if (pctwet(n) >= 25._r8) then
-             pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n))
-          else if (pcturb(n) >= 25._r8) then
-             pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n))
-          else if (pctgla(n) >= 25._r8) then
-             pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n))
+            (pctnatpft(n)%get_pct_l2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_l2g() <  1.e-6_r8) ) then
+          write (6,*) 'Special plus crop land units near 100%, but not quite for n,suma =',n,suma
+          write (6,*) 'Adjusting special plus crop land units to 100%'
+          if (pctlak(n) >= 1.0_r8) then
+             pctlak(n) = 100._r8 - (pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g())
+          else if (pctwet(n) >= 1.0_r8) then
+             pctwet(n) = 100._r8 - (pctlak(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g())
+          else if (pcturb(n) >= 1.0_r8) then
+             pcturb(n) = 100._r8 - (pctlak(n) + pctwet(n) + pctgla(n) + pctcft(n)%get_pct_l2g())
+          else if (pctgla(n) >= 1.0_r8) then
+             pctgla(n) = 100._r8 - (pctlak(n) + pctwet(n) + pcturb(n) + pctcft(n)%get_pct_l2g())
+          else if (pctcft(n)%get_pct_l2g() >= 1.0_r8) then
+             call pctcft(n)%set_pct_l2g(100._r8 - (pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)))
           else
-             write (6,*) subname, 'Error: sum of special land units nearly 100% but none is >= 25% at ', &
+             write (6,*) subname, 'Error: sum of special plus crop land units nearly 100% but none is >= 1.0% at ', &
                   'n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),pctnatveg(n),pctcrop(n),suma = ', &
                   n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
                   pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(),suma
              call abort()
           end if
           call pctnatpft(n)%set_pct_l2g(0._r8)
-          call pctcft(n)%set_pct_l2g(0._r8)
        end if
        if ( any(pctnatpft(n)%get_pct_p2g() > 0.0_r8 .and. pctnatpft(n)%get_pct_p2g() < toosmallPFT ) .or. &
             any(pctcft(n)%get_pct_p2g()    > 0.0_r8 .and. pctcft(n)%get_pct_p2g()    < toosmallPFT )) then
@@ -1593,14 +1539,8 @@ subroutine normalizencheck_landuse(ldomain)
           call abort()
        end if
        
-       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n)
-       if (suma < 100._r8-epsilon(suma) .and. suma > (100._r8 - 4._r8*epsilon(suma))) then
-          write (6,*) subname, 'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', &
-               n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
-               pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g()
-          call abort()
-       end if
-       suma = suma + pctnatpft(n)%get_pct_l2g() + pctcft(n)%get_pct_l2g()
+       suma = pctlak(n) + pctwet(n) + pcturb(n) + pctgla(n) + pctcft(n)%get_pct_l2g()
+       suma = suma + pctnatpft(n)%get_pct_l2g()
        if ( abs(suma-100._r8) > 1.e-10_r8) then
           write (6,*) subname, ' error: sum of pctlak, pctwet,', &
                'pcturb, pctgla, pctnatveg and pctcrop is NOT equal to 100'
@@ -1612,45 +1552,6 @@ subroutine normalizencheck_landuse(ldomain)
        
     end do
 
-    ! Check that when pctnatveg+pctcrop identically zero, sum of special landunits is identically 100%
-
-    if ( .not. outnc_double )then
-       do n = 1,ns_o
-          sum8  =         real(pctlak(n),r4)
-          sum8  = sum8  + real(pctwet(n),r4)
-          sum8  = sum8  + real(pcturb(n),r4)
-          sum8  = sum8  + real(pctgla(n),r4)
-          sum4a =         real(pctnatpft(n)%get_pct_l2g(),r4)
-          sum4a = sum4a + real(pctcft(n)%get_pct_l2g(),r4)
-          if ( sum4a==0.0_r4 .and. sum8 < 100._r4-2._r4*epsilon(sum4a) )then
-             write (6,*) subname, ' error: sum of pctlak, pctwet,', &
-                  'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8
-             write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop= ', &
-                  n,pctlak(n),pctwet(n),pcturb(n),pctgla(n), &
-                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g()
-             call abort()
-          end if
-       end do
-    else
-       do n = 1,ns_o
-          sum8  =         pctlak(n)
-          sum8  = sum8  + pctwet(n)
-          sum8  = sum8  + pcturb(n)
-          sum8  = sum8  + pctgla(n)
-          sum8a =         pctnatpft(n)%get_pct_l2g()
-          sum8a = sum8a + pctcft(n)%get_pct_l2g()
-          if ( sum8a==0._r8 .and. sum8 < (100._r8-4._r8*epsilon(sum8)) )then
-             write (6,*) subname, ' error: sum of pctlak, pctwet,', &
-                  'pcturb, pctgla is < 100% when pctnatveg+pctcrop==0 sum = ', sum8
-             write (6,*) 'Total error, error/epsilon = ',100._r8-sum8, ((100._r8-sum8)/epsilon(sum8))
-             write (6,*)'n,pctlak,pctwet,pcturb,pctgla,pctnatveg,pctcrop,epsilon= ', &
-                  n,pctlak(n),pctwet(n),pcturb(n),pctgla(n),&
-                  pctnatpft(n)%get_pct_l2g(),pctcft(n)%get_pct_l2g(), epsilon(sum8)
-             call abort()
-          end if
-       end do
-    end if
-
     ! Make sure that there is no vegetation outside the pft mask
     do n = 1,ns_o
        if (pftdata_mask(n) == 0 .and. (pctnatpft(n)%get_pct_l2g() > 0 .or. pctcft(n)%get_pct_l2g() > 0)) then