Manually introduce the changes found in #1586

- The code builds and a non-transient case completes:
mpiexec_mpt -p "%g:" -np 288 ./src/mksurfdata < surfdata_1.9x2.5_hist_78pfts_CMIP6_1850_c220316.namelist
- Transient case gives error:
mpiexec_mpt -p "%g:" -np 288 ./src/mksurfdata < surfdata_1.9x2.5_hist_78pfts_CMIP6_1850-1899_c220316.namelist
mksrfdata error: year for urban not equal to year for PFT files

...but I don't, yet, see a mistake in my implementation:
from mksurfdata_map/src/mksurfdat.F90 lines 1117-1130
to mksurfdata_esmf/src/mksurfdat.F90 lines 883-908

tools/mksurfdata_esmf/src/mkfileMod.F90

@@ -526,6 +526,21 @@ subroutine mkfile_define_vars(pioid, dynlanduse)
 
     else
 
+       call mkpio_def_spatial_var(pioid=pioid, varname='PCT_URBAN', xtype=xtype, &
+            lev1name = 'numurbl', lev2name='time', &
+            long_name = "percent urban for each density type", units = "unitless")
+
+       call mkpio_def_spatial_var(pioid=pioid, varname='PCT_URBAN_MAX', xtype=xtype, &
+            lev1name = 'numurbl', &
+            long_name = "maximum percent urban for each density type", units = "unitless")
+
+       call mkpio_def_spatial_var(pioid=pioid, varname='PCT_LAKE', xtype=xtype, &
+            lev1name = 'time', &
+            long_name = "percent lake", units = "unitless")
+
+       call mkpio_def_spatial_var(pioid=pioid, varname='PCT_LAKE_MAX', xtype=xtype, &
+            long_name = "maximum percent lake", units = "unitless")
+
        call mkpio_def_spatial_var(pioid=pioid, varname='HARVEST_VH1', xtype=xtype, &
             lev1name='time', &
             long_name = "harvest from primary forest", units = "gC/m2/yr")

tools/mksurfdata_esmf/src/mkinputMod.F90

@@ -31,6 +31,8 @@ module mkinputMod
   character(CL) , public    :: fsurlog  ! output surface log file name
   character(CL) , public    :: fdyndat  ! dynamic landuse data file name
   character(CL) , public    :: fhrvname ! generic harvest filename
+  character(CL) , public    :: furbname ! generic transient urban land cover filename
+  character(CL) , public    :: flakname ! generic lake filename
 
   character(CX) , public    :: mksrf_fgrid_mesh             = ' ' ! land grid file name to use
   integer       , public    :: mksrf_fgrid_mesh_nx          = -999

tools/mksurfdata_esmf/src/mklanwatMod.F90

@@ -22,6 +22,7 @@ module mklanwatMod
 
   public :: mklakwat    ! make %lake and lake dpeth
   public :: mkwetlnd    ! make %wetland
+  public :: update_max_array_lake  ! Update the maximum lake percent
 
   character(len=*) , parameter :: u_FILE_u = &
        __FILE__
@@ -328,4 +329,28 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc)
 
   end subroutine mkwetlnd
 
+!===============================================================
+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

tools/mksurfdata_esmf/src/mksurfdata.F90

@@ -60,7 +60,7 @@ program mksurfdata
   ! ======================================
   ! 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, 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
@@ -106,8 +106,8 @@ program mksurfdata
   use mksoilfmaxMod      , only : mksoilfmax
   use mksoildepthMod     , only : mksoildepth
   use mksoilcolMod       , only : mksoilcol
-  use mkurbanparMod      , only : mkurbanInit, mkurban, mkurbanpar, mkurban_topo, numurbl
-  use mklanwatMod        , only : mklakwat, mkwetlnd
+  use mkurbanparMod      , only : mkurbanInit, mkurban, mkurbanpar, mkurban_topo, numurbl, update_max_array_urban
+  use mklanwatMod        , only : mklakwat, mkwetlnd, update_max_array_lake
   use mkorganicMod       , only : mkorganic
   use mkutilsMod         , only : normalize_classes_by_gcell, chkerr
   use mkfileMod          , only : mkfile_define_dims, mkfile_define_atts, mkfile_define_vars
@@ -161,13 +161,17 @@ program mksurfdata
 
   ! inland water data, glacier data and urban data
   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           :: pctgla(:)               ! percent of grid cell that is glacier
   integer , allocatable           :: urban_region(:)         ! urban region ID
   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           :: urban_classes(:,:)      ! percent cover of each urban class, as % of total urban area
   real(r8), allocatable           :: urban_classes_g(:,:)    ! percent cover of each urban class, as % of grid cell
   real(r8), allocatable           :: elev(:)                 ! glc elevation (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
 
   ! pio/esmf variables
   type(file_desc_t)               :: pioid
@@ -413,17 +417,20 @@ program mksurfdata
   ! LAKEDEPTH is written out in the subroutine
   ! Need to keep pctlak and pctwet external for use below
   allocate ( pctlak(lsize_o)) ; pctlak(:) = spval
+  allocate ( pctlak_max(lsize_o)) ; pctlak_max(:) = spval
   call mklakwat(mksrf_flakwat_mesh, mksrf_flakwat, mesh_model, pctlak, pioid, fsurdat, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mklatwat')
 
   allocate ( pctwet(lsize_o)) ; pctwet(:) = spval
+  allocate ( pctwet_orig(lsize_o)) ; pctwet_orig(:) = spval
   call mkwetlnd(mksrf_fwetlnd_mesh, mksrf_fwetlnd, mesh_model, pctwet, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkwetlnd')
 
   ! -----------------------------------
   ! Make glacier fraction [pctgla] from [fglacier] dataset
   ! -----------------------------------
   allocate (pctgla(lsize_o)) ; pctgla(:) = spval
+  allocate (pctgla_orig(lsize_o)) ; pctgla_orig(:) = spval
   call mkglacier (mksrf_fglacier_mesh, mksrf_fglacier, mesh_model, glac_o=pctgla, rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkglacier')
 
@@ -500,6 +507,7 @@ program mksurfdata
   ! Make urban fraction [pcturb] from [furban] dataset and
   ! -----------------------------------
   allocate (pcturb(lsize_o))                 ; pcturb(:)            = spval
+  allocate (pcturb_max(lsize_o, numurbl))    ; pcturb_max(:,:)      = spval
   allocate (urban_classes(lsize_o,numurbl))  ; urban_classes(:,:)   = spval
   allocate (urban_region(lsize_o))           ; urban_region(:)      = -999
   call mkurban(mksrf_furban_mesh, mksrf_furban, mesh_model, pcturb, urban_classes, urban_region, rc=rc)
@@ -526,7 +534,6 @@ program mksurfdata
   where (elev > elev_thresh)
      pcturb = 0._r8
   end where
-  deallocate(elev)
 
   ! -----------------------------------
   ! Compute topography statistics [topo_stddev, slope] from [ftopostats]
@@ -611,6 +618,10 @@ program mksurfdata
      end if
   end do
 
+  ! Save special land unit areas of surface dataset 
+  pctwet_orig(:) = pctwet(:)
+  pctgla_orig(:) = pctgla(:)
+
   ! -----------------------------------
   ! Perform other normalizations
   ! -----------------------------------
@@ -835,6 +846,8 @@ program mksurfdata
 
      pctnatpft_max = pctnatpft
      pctcft_max = pctcft
+     pcturb_max = urban_classes_g
+     pctlak_max = pctlak
 
      end_of_fdynloop = .false.
      ntim = 0
@@ -867,6 +880,33 @@ program mksurfdata
            end if
            call shr_sys_abort()
         end if
+        ! Read input urban data
+        if (root_task) then
+           read(nfdyn, '(A195,1x,I4)', iostat=ier) furbname, year2
+           write(ndiag,*)'input urban dynamic dataset for year ', year2, ' is : ', trim(furbname)
+        end if
+        call mpi_bcast (furbname, len(furbname), MPI_CHARACTER, 0, mpicom, ier)
+        call mpi_bcast (year2, 1, MPI_INTEGER, 0, mpicom, ier)
+        if ( year2 /= year ) then
+           if (root_task) then
+              write(ndiag,*) subname, ' error: year for urban not equal to year for PFT files'
+           end if
+           call shr_sys_abort()
+        end if
+        ! Read input lake data
+        if (root_task) then
+           read(nfdyn, '(A195,1x,I4)', iostat=ier) flakname, year2
+           write(ndiag,*)'input lake dynamic dataset for year ', year2, ' is : ', trim(flakname)
+        end if
+        call mpi_bcast (flakname, len(flakname), MPI_CHARACTER, 0, mpicom, ier)
+        call mpi_bcast (year2, 1, MPI_INTEGER, 0, mpicom, ier)
+        if ( year2 /= year ) then
+           if (root_task) then
+              write(ndiag,*) subname, ' error: year for lake not equal to year for PFT files'
+           end if
+           call shr_sys_abort()
+        end if
+
         ntim = ntim + 1
         if (root_task) then
            write(ndiag,'(a,i8)')subname//' ntime = ',ntim
@@ -911,6 +951,23 @@ program mksurfdata
         if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkharvest')
         call pio_syncfile(pioid)
 
+        ! Create pctlak data at model resolution (use original mapping file from lake data)
+        call mklakwat(mksrf_flakwat_mesh, flakname, mesh_model, pctlak, pioid, fsurdat, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mklakwat')
+        call pio_syncfile(pioid)
+
+        call mkurban(mksrf_furban_mesh, furbname, mesh_model, pcturb, urban_classes, urban_region, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkurban')
+        call pio_syncfile(pioid)
+        ! screen pcturb using elevation
+        where (elev > elev_thresh)
+           pcturb = 0._r8
+        end where
+
+        ! 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(:)
+
         ! Do landuse changes such as for the poles, etc.
         ! If have pole points on grid - set south pole to glacier
         ! north pole is assumed as non-land
@@ -935,10 +992,13 @@ program mksurfdata
            write(ndiag,'(a)')' calling normalize_and_check_landuse'
         end if
         call normalize_and_check_landuse(lsize_o)
+        call normalize_classes_by_gcell(urban_classes, pcturb, urban_classes_g)
 
         ! Given an array of pct_pft_type variables, update all the max_p2l variables.
         call update_max_array(pctnatpft_max, pctnatpft)
         call update_max_array(pctcft_max, pctcft)
+        call update_max_array_urban(pcturb_max,urban_classes_g)
+        call update_max_array_lake(pctlak_max,pctlak)
 
         if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_NAT_PFT for year ",year
         rcode = pio_inq_varid(pioid, 'PCT_NAT_PFT', pio_varid)
@@ -956,6 +1016,20 @@ program mksurfdata
         if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_CROP')
         call pio_syncfile(pioid)
 
+        if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_URBAN for year ",year
+        rcode = pio_inq_varid(pioid, 'PCT_URBAN', pio_varid)
+        call pio_setframe(pioid, pio_varid, int(ntim, kind=Pio_Offset_Kind))
+        call mkfile_output(pioid, mesh_model, 'PCT_URBAN', urban_classes_g, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_URBAN')
+        call pio_syncfile(pioid)
+
+        if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_LAKE for year ",year
+        rcode = pio_inq_varid(pioid, 'PCT_LAKE', pio_varid)
+        call pio_setframe(pioid, pio_varid, int(ntim, kind=Pio_Offset_Kind))
+        call mkfile_output(pioid, mesh_model, 'PCT_LAKE', pctlak, rc=rc)
+        if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_LAKE')
+        call pio_syncfile(pioid)
+
         if (num_cft > 0) then
            if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_CFT for year ",year
            rcode = pio_inq_varid(pioid, 'PCT_CFT', pio_varid)
@@ -983,6 +1057,14 @@ program mksurfdata
      call mkfile_output(pioid, mesh_model, 'PCT_CROP_MAX', pctcrop, rc=rc)
      if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_CROP')
 
+     if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_URBAN_MAX"
+     call mkfile_output(pioid, mesh_model, 'PCT_URBAN_MAX', pcturb_max, rc=rc)
+     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_URBAN_MAX')
+
+     if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_LAKE_MAX"
+     call mkfile_output(pioid, mesh_model, 'PCT_LAKE_MAX', pctlak_max, rc=rc)
+     if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_LAKE_MAX')
+
      if (num_cft > 0) then
         if (root_task)  write(ndiag, '(a,i8)') trim(subname)//" writing PCT_CFT_MAX"
         call get_pct_p2l_array(pctcft_max, ndim1=lsize_o, ndim2=num_cft, pct_p2l=pct_cft)
@@ -1011,8 +1093,6 @@ subroutine normalize_and_check_landuse(ns_o)
       ! 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)
-      !
       use mkpftConstantsMod , only : baregroundindex
       use mkpftUtilsMod     , only : adjust_total_veg_area
 
@@ -1068,12 +1148,12 @@ subroutine normalize_and_check_landuse(ns_o)
          if (suma > (100._r8 + tol_loose)) then
             if (root_task) then
                write(ndiag,*) subname, ' ERROR: pctlak + pctwet + pcturb + pctgla must be'
-               write(ndiag,*) '<= 100% before calling this subroutine'
+               write(ndiag,*) '<= 100% before normalizing vegetated area'
                write(ndiag,*) 'n, pctlak, pctwet, pcturb, pctgla = ', &
                     n, pctlak(n), pctwet(n), pcturb(n), pctgla(n)
             else
                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)
             end if

tools/mksurfdata_esmf/src/mkurbanparMod.F90

@@ -24,6 +24,7 @@ module mkurbanparMod
   public :: mkurbanInit
   public :: mkurban
   public :: mkurbanpar
+  public :: update_max_array_urban
   public :: normalize_urbn_by_tot
   public :: mkurban_topo                   ! Get elevation to reduce urban for high elevation areas
   public :: mkurban_pct_diagnostics        ! print diagnostics related to pct urban
@@ -954,4 +955,30 @@ subroutine mkurban_topo(file_mesh_i, file_data_i, mesh_o, varname, elev_o, rc)
 
   end subroutine mkurban_topo
 
+  !===============================================================
+  subroutine update_max_array_urban(pct_urbmax_arr,pct_urban_arr)
+    !
+    ! !DESCRIPTION:
+    ! 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
+    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