From 70a10382dabd4d31f8bdf7749282842c5b0133ff Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Mon, 21 Mar 2022 12:09:50 -0600 Subject: [PATCH] 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 | 15 ++++ tools/mksurfdata_esmf/src/mkinputMod.F90 | 2 + tools/mksurfdata_esmf/src/mklanwatMod.F90 | 25 ++++++ tools/mksurfdata_esmf/src/mksurfdata.F90 | 96 +++++++++++++++++++-- tools/mksurfdata_esmf/src/mkurbanparMod.F90 | 27 ++++++ 5 files changed, 157 insertions(+), 8 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mkfileMod.F90 b/tools/mksurfdata_esmf/src/mkfileMod.F90 index cdb5192c8f..01d6453000 100644 --- a/tools/mksurfdata_esmf/src/mkfileMod.F90 +++ b/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") diff --git a/tools/mksurfdata_esmf/src/mkinputMod.F90 b/tools/mksurfdata_esmf/src/mkinputMod.F90 index 9afcbfce51..1d5c3a38a9 100644 --- a/tools/mksurfdata_esmf/src/mkinputMod.F90 +++ b/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 diff --git a/tools/mksurfdata_esmf/src/mklanwatMod.F90 b/tools/mksurfdata_esmf/src/mklanwatMod.F90 index f380a0109b..30c53a5ee9 100644 --- a/tools/mksurfdata_esmf/src/mklanwatMod.F90 +++ b/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 diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index d5919c18cb..8295873178 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/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,10 +417,12 @@ 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') @@ -424,6 +430,7 @@ program mksurfdata ! 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 diff --git a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 b/tools/mksurfdata_esmf/src/mkurbanparMod.F90 index 959a75dd7b..0726672642 100644 --- a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 +++ b/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