From 932c28816d84493c78b3ba509c21b2013c73d931 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 28 Oct 2022 12:03:28 -0600 Subject: [PATCH 01/15] Add option to convert ocean/wetland to land Resolves ESCOMP/CTSM#1878 --- bld/CLMBuildNamelist.pm | 1 + bld/namelist_files/namelist_defaults_ctsm.xml | 3 ++ .../namelist_definition_ctsm.xml | 7 ++++ src/main/clm_varctl.F90 | 5 +++ src/main/controlMod.F90 | 4 ++ src/main/surfrdMod.F90 | 8 +++- src/main/surfrdUtilsMod.F90 | 42 ++++++++++++++++++- 7 files changed, 67 insertions(+), 3 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index ead6c054f7..797a0acfa3 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2034,6 +2034,7 @@ sub setup_logic_subgrid { my $var = 'run_zero_weight_urban'; add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'convert_ocean_to_land'); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'collapse_urban', 'structure'=>$nl_flags->{'structure'}); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'n_dom_landunits', diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 8ead40701c..82910d7191 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -1754,6 +1754,9 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts .false. +.false. +.true. + .false. .true. diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index a08795dd1f..ad901107d3 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1161,6 +1161,13 @@ Toggle for vancouver specific logic. Toggle for mexico city specific logic. + +If true, any ocean (i.e., wetland) points on the surface dataset are +converted to bare ground (or whatever vegetation is given in that grid +cell - but typically this will be bare ground). + + diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index 9be9af2f73..208baa9be4 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -168,6 +168,11 @@ module clm_varctl ! true => make ALL patches, cols & landunits active (even if weight is 0) logical, public :: all_active = .false. + ! true => any ocean (i.e., "wetland") points on the surface dataset are converted to + ! bare ground (or whatever vegetation is given in that grid cell... but typically this + ! will be bare ground) + logical, public :: convert_ocean_to_land = .false. + logical, public :: collapse_urban = .false. ! true => collapse urban landunits to the dominant urban landunit; default = .false. means "do nothing" i.e. keep all urban landunits as found in the input data integer, public :: n_dom_landunits = -1 ! # of dominant landunits; determines the number of active landunits; default = 0 (set in namelist_defaults_ctsm.xml) means "do nothing" integer, public :: n_dom_pfts = -1 ! # of dominant pfts; determines the number of active pfts; default = 0 (set in namelist_defaults_ctsm.xml) means "do nothing" diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index a07228aa0d..eaa97d8230 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -259,6 +259,8 @@ subroutine control_init(dtime) ! All old cpp-ifdefs are below and have been converted to namelist variables + namelist /clm_inparm/ convert_ocean_to_land + ! Number of dominant pfts and landunits. Enhance ctsm performance by ! reducing the number of active pfts to n_dom_pfts and ! active landunits to n_dom_landunits. @@ -668,6 +670,7 @@ subroutine control_spmd() ! Other subgrid logic call mpi_bcast(run_zero_weight_urban, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast(all_active, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast(convert_ocean_to_land, 1, MPI_LOGICAL, 0, mpicom, ier) ! Number of dominant pfts and landunits. Enhance ctsm performance by ! reducing the number of active pfts to n_dom_pfts and @@ -890,6 +893,7 @@ subroutine control_print () else write(iulog,*) ' land frac data = ',trim(fatmlndfrc) end if + write(iulog,*) ' Convert ocean to land = ', convert_ocean_to_land write(iulog,*) ' Number of ACTIVE PFTS (0 means input pft data NOT collapsed to n_dom_pfts) =', n_dom_pfts write(iulog,*) ' Number of ACTIVE LANDUNITS (0 means input landunit data NOT collapsed to n_dom_landunits) =', n_dom_landunits write(iulog,*) ' Collapse urban landunits; done before collapsing all landunits to n_dom_landunits; .false. means do nothing i.e. keep all the urban landunits, though n_dom_landunits may still remove them =', collapse_urban diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index a8bbf17ec7..d2c37f85f6 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -15,7 +15,7 @@ module surfrdMod use clm_varcon , only : grlnd use clm_varctl , only : iulog use clm_varctl , only : use_cndv, use_crop, use_fates - use surfrdUtilsMod , only : check_sums_equal_1, collapse_crop_types + use surfrdUtilsMod , only : check_sums_equal_1, apply_convert_ocean_to_land, collapse_crop_types use surfrdUtilsMod , only : collapse_to_dominant, collapse_crop_var, collapse_individual_lunits use ncdio_pio , only : file_desc_t, var_desc_t, ncd_pio_openfile, ncd_pio_closefile use ncdio_pio , only : ncd_io, check_var, ncd_inqfdims, check_dim_size, ncd_inqdid, ncd_inqdlen @@ -70,7 +70,7 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) ! o real % abundance PFTs (as a percent of vegetated area) ! ! !USES: - use clm_varctl , only : create_crop_landunit, collapse_urban, & + use clm_varctl , only : create_crop_landunit, convert_ocean_to_land, collapse_urban, & toosmall_soil, toosmall_crop, toosmall_glacier, & toosmall_lake, toosmall_wetland, toosmall_urban, & n_dom_landunits @@ -200,6 +200,10 @@ subroutine surfrd_get_data (begg, endg, ldomain, lfsurdat, actual_numcft) call check_sums_equal_1(wt_lunit, begg, 'wt_lunit', subname) + if (convert_ocean_to_land) then + call apply_convert_ocean_to_land(wt_lunit(begg:endg,:), begg, endg) + end if + ! if collapse_urban = .true. ! collapse urban landunits to the dominant urban landunit diff --git a/src/main/surfrdUtilsMod.F90 b/src/main/surfrdUtilsMod.F90 index 0763d43a16..95303322a4 100644 --- a/src/main/surfrdUtilsMod.F90 +++ b/src/main/surfrdUtilsMod.F90 @@ -20,7 +20,8 @@ module surfrdUtilsMod ! !PUBLIC MEMBER FUNCTIONS: public :: check_sums_equal_1 ! Confirm that sum(arr(n,:)) == 1 for all n public :: renormalize ! Renormalize an array - public :: convert_cft_to_pft ! Conversion of crop CFT to natural veg PFT:w + public :: apply_convert_ocean_to_land ! Apply the conversion of ocean to land points + public :: convert_cft_to_pft ! Conversion of crop CFT to natural veg PFT public :: collapse_crop_types ! Collapse unused crop types into types used in this run public :: collapse_individual_lunits ! Collapse landunits by user-defined thresholds public :: collapse_to_dominant ! Collapse to dominant pfts or landunits @@ -112,6 +113,45 @@ subroutine renormalize(arr, lb, normal) end subroutine renormalize + !----------------------------------------------------------------------- + subroutine apply_convert_ocean_to_land(wt_lunit, begg, endg) + ! + ! !DESCRIPTION: + ! Apply the conversion of ocean points to land, by changing all "wetland" points to + ! natveg; typically this will result in these points becoming bare ground. + ! + ! The motivation for doing this is to avoid the negative runoff that sometimes comes + ! from wetlands. + ! + ! !USES: + use landunit_varcon, only : istsoil, istwet, max_lunit + ! + ! !ARGUMENTS: + integer, intent(in) :: begg ! Beginning grid cell index + integer, intent(in) :: endg ! Ending grid cell index + ! This array is modified in-place: + real(r8), intent(inout) :: wt_lunit(begg:endg, max_lunit) ! Weights of landunits per grid cell + ! + ! !LOCAL VARIABLES: + integer :: g + + character(len=*), parameter :: subname = 'apply_convert_ocean_to_land' + !----------------------------------------------------------------------- + + ! BUG(wjs, 2022-10-27, ESCOMP/CTSM#1886) Ideally we would distinguish between ocean + ! vs. true wetland points on the surface dataset; for now oceans are included in the + ! wetland area on the surface dataset, so we convert all wetlands to land. (Typically + ! there are no true/inland wetlands on the surface dataset, so this is currently okay, + ! but this would become a problem if we started having inland wetlands on the surface + ! dataset again.) + do g = begg, endg + wt_lunit(g,istsoil) = wt_lunit(g,istsoil) + wt_lunit(g,istwet) + wt_lunit(g,istwet) = 0._r8 + end do + + end subroutine apply_convert_ocean_to_land + + !----------------------------------------------------------------------- subroutine convert_cft_to_pft( begg, endg, cftsize, wt_cft ) ! From 5786052bb96df248256a50aa01212c5a6c67ac9e Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 1 Nov 2022 17:33:40 -0600 Subject: [PATCH 02/15] Don't overwrite soil properties just because pctlnd_pft is small It doesn't seem like there's any reason why we need to overwrite soil properties with a default value just because pctlnd_pft is small. Doing so was fine when we turned those grid cells into wetland, but now that we aren't necessarily doing that, it makes sense to keep valid soil data wherever we have it. This could change answers in places where we have valid soil data despite pctlnd_pft being near 0. --- tools/mksurfdata_esmf/src/mksoilcolMod.F90 | 12 ++---------- tools/mksurfdata_esmf/src/mksoiltexMod.F90 | 7 +++---- tools/mksurfdata_esmf/src/mksurfdata.F90 | 4 ++-- 3 files changed, 7 insertions(+), 16 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mksoilcolMod.F90 b/tools/mksurfdata_esmf/src/mksoilcolMod.F90 index 25ba6956b0..6fa9600f3c 100644 --- a/tools/mksurfdata_esmf/src/mksoilcolMod.F90 +++ b/tools/mksurfdata_esmf/src/mksoilcolMod.F90 @@ -29,13 +29,12 @@ module mksoilcolMod contains !================================================================================= - subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pctlnd_pft_o, pioid_o, rc) + subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pioid_o, rc) ! input/output variables character(len=*) , intent(in) :: file_mesh_i ! input mesh file name character(len=*) , intent(in) :: file_data_i ! input data file name type(ESMF_Mesh) , intent(in) :: mesh_o ! model mesho - real(r8) , intent(in) :: pctlnd_pft_o(:) type(file_desc_t) , intent(inout) :: pioid_o integer , intent(out) :: rc @@ -185,13 +184,6 @@ subroutine mksoilcol(file_data_i, file_mesh_i, mesh_o, pctlnd_pft_o, pioid_o, rc end if end do - ! assume medium soil color (15) and loamy texture if pct_lndpft is small - do no = 1,ns_o - if (pctlnd_pft_o(no) < 1.e-6_r8) then - soil_color_o(no) = 15 - end if - end do - ! Write output data if (root_task) write(ndiag, '(a)') trim(subname)//" writing out soil color" call mkfile_output(pioid_o, mesh_o, 'SOIL_COLOR', soil_color_o, rc=rc) @@ -277,7 +269,7 @@ subroutine get_dominant_soilcol(dynamicMaskList, dynamicSrcMaskValue, dynamicDst soil_color_o = maxindex(1) end if - ! If land but no color, set color to 15 (in older dataset generic soil color 4) + ! If no color, set color to 15 (in older dataset generic soil color 4) if (num_soilcolors == 8) then if (soil_color_o == 0) then soil_color_o = 4 diff --git a/tools/mksurfdata_esmf/src/mksoiltexMod.F90 b/tools/mksurfdata_esmf/src/mksoiltexMod.F90 index 295217230f..ae975a9ed9 100644 --- a/tools/mksurfdata_esmf/src/mksoiltexMod.F90 +++ b/tools/mksurfdata_esmf/src/mksoiltexMod.F90 @@ -33,7 +33,7 @@ module mksoiltexMod contains !================================================================================= - subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, pctlnd_pft_o, rc) + subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o, rc) ! ! make %sand, %clay, organic carbon content, coarse fragments, bulk density, ! and pH measured in H2O @@ -44,7 +44,6 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o character(len=*) , intent(in) :: file_lookup_i ! input data file name type(ESMF_Mesh) , intent(in) :: mesh_o ! output mesh type(file_desc_t) , intent(inout) :: pioid_o - real(r8) , intent(in) :: pctlnd_pft_o(:) ! PFT data: % of gridcell for PFTs integer , intent(out) :: rc ! local variables @@ -291,9 +290,9 @@ subroutine mksoiltex(file_mesh_i, file_mapunit_i, file_lookup_i, mesh_o, pioid_o do no = 1,ns_o - if (pctlnd_pft_o(no) < 1.e-6_r8 .or. mapunit_o(no) == 0) then + if (mapunit_o(no) == 0) then - ! Set sand and clay to loam if pctlnd_pft is < 1.e-6 or mapunit is 0 + ! Set sand and clay to loam if mapunit is 0 sand_o(no,:) = 43._r4 clay_o(no,:) = 18._r4 orgc_o(no,:) = 0._r4 diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index 5b216671c3..6a542d36bb 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -458,7 +458,7 @@ program mksurfdata ! ----------------------------------- if (fsurdat /= ' ') then call mksoiltex( mksrf_fsoitex_mesh, file_mapunit_i=mksrf_fsoitex, file_lookup_i=mksrf_fsoitex_lookup, & - mesh_o=mesh_model, pioid_o=pioid, pctlnd_pft_o=pctlnd_pft, rc=rc) + mesh_o=mesh_model, pioid_o=pioid, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mksoiltex') end if @@ -467,7 +467,7 @@ program mksurfdata ! ----------------------------------- if (fsurdat /= ' ') then ! SOIL_COLOR and mxsoil_color is written out in the subroutine - call mksoilcol( mksrf_fsoicol, mksrf_fsoicol_mesh, mesh_model, pctlnd_pft, pioid, rc) + call mksoilcol( mksrf_fsoicol, mksrf_fsoicol_mesh, mesh_model, pioid, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mksoilcol') end if From 27ab891ccf0849cfd83ac58d728bb8e7943b55a5 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 1 Nov 2022 18:13:27 -0600 Subject: [PATCH 03/15] Remove a comment that is no longer true --- tools/mksurfdata_esmf/src/mksurfdata.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index 6a542d36bb..526a999a84 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -1080,7 +1080,6 @@ subroutine normalize_and_check_landuse(ns_o) 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). if (pctlnd_pft(n) < 1.e-6_r8) then if (pctgla(n) < 1.e-6_r8) then From 137bd5d1292514c7ad1b4d7867b039f418c44243 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 7 Nov 2022 17:36:50 -0700 Subject: [PATCH 04/15] Add a note --- bld/namelist_files/namelist_definition_ctsm.xml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index ad901107d3..ffd6ea53eb 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -1165,7 +1165,8 @@ Toggle for mexico city specific logic. group="clm_inparm" value=".false."> If true, any ocean (i.e., wetland) points on the surface dataset are converted to bare ground (or whatever vegetation is given in that grid -cell - but typically this will be bare ground). +cell - but typically this will be bare ground due to lack of vegetation +in grid cells with 100% ocean). Date: Fri, 11 Nov 2022 14:02:14 -0700 Subject: [PATCH 05/15] Add a norm_by_fracs argument to create_routehandle This will allow switching between FRACAREA and DSTAREA normalization; we'll use the latter for PCT areas. --- tools/mksurfdata_esmf/src/mkVICparamsMod.F90 | 3 +- tools/mksurfdata_esmf/src/mkesmfMod.F90 | 42 +++++++++++++++++-- tools/mksurfdata_esmf/src/mkgdpMod.F90 | 3 +- .../src/mkglacierregionMod.F90 | 3 +- tools/mksurfdata_esmf/src/mkglcmecMod.F90 | 6 ++- tools/mksurfdata_esmf/src/mkharvestMod.F90 | 3 +- tools/mksurfdata_esmf/src/mklaiMod.F90 | 3 +- tools/mksurfdata_esmf/src/mklanwatMod.F90 | 9 ++-- tools/mksurfdata_esmf/src/mkpeatMod.F90 | 3 +- tools/mksurfdata_esmf/src/mkpftMod.F90 | 3 +- tools/mksurfdata_esmf/src/mksoildepthMod.F90 | 6 ++- tools/mksurfdata_esmf/src/mksoilfmaxMod.F90 | 3 +- tools/mksurfdata_esmf/src/mkurbanparMod.F90 | 6 ++- tools/mksurfdata_esmf/src/mkvocefMod.F90 | 3 +- 14 files changed, 74 insertions(+), 22 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mkVICparamsMod.F90 b/tools/mksurfdata_esmf/src/mkVICparamsMod.F90 index 7742eb282a..5f69152e2d 100644 --- a/tools/mksurfdata_esmf/src/mkVICparamsMod.F90 +++ b/tools/mksurfdata_esmf/src/mkVICparamsMod.F90 @@ -116,7 +116,8 @@ subroutine mkVICparams(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) do n = 1, ns_o diff --git a/tools/mksurfdata_esmf/src/mkesmfMod.F90 b/tools/mksurfdata_esmf/src/mkesmfMod.F90 index ee31104b80..893b1fec80 100644 --- a/tools/mksurfdata_esmf/src/mkesmfMod.F90 +++ b/tools/mksurfdata_esmf/src/mkesmfMod.F90 @@ -27,11 +27,21 @@ module mkesmfMod contains !=============================================================== - subroutine create_routehandle_r4(mesh_i, mesh_o, routehandle, frac_o, rc) + subroutine create_routehandle_r4(mesh_i, mesh_o, norm_by_fracs, routehandle, frac_o, rc) ! input/output variables type(ESMF_Mesh) , intent(in) :: mesh_i type(ESMF_Mesh) , intent(in) :: mesh_o + ! If norm_by_fracs is .true., then remapping is done using ESMF_NORMTYPE_FRACAREA; + ! otherwise, remapping is done using ESMF_NORMTYPE_DSTAREA. FRACAREA normalization + ! adds a normalization factor of the fraction of the unmasked source grid that + ! overlaps with a destination cell. FRACAREA normalization is appropriate when you + ! want to treat values outside the mask as missing values that shouldn't contribute + ! to the average (this is appropriate for most fields); DSTAREA normalization is + ! appropriate when you want to treat values outside the mask as 0 (this is + ! appropriate for PCT cover fields where we want the final value to be expressed as + ! percent of the entire gridcell area). + logical , intent(in) :: norm_by_fracs type(ESMF_RouteHandle) , intent(inout) :: routehandle real(r4), optional , intent(inout) :: frac_o(:) integer , intent(out) :: rc @@ -40,6 +50,7 @@ subroutine create_routehandle_r4(mesh_i, mesh_o, routehandle, frac_o, rc) integer :: srcMaskValue = 0 ! ignore source points where the mesh mask is 0 integer :: dstMaskValue = -987987 ! don't ingore any destination points integer :: srcTermProcessing_Value = 0 + type(ESMF_NormType_Flag) :: normtype type(ESMF_Field) :: field_i type(ESMF_Field) :: field_o type(ESMF_Field) :: dstfracfield @@ -56,9 +67,15 @@ subroutine create_routehandle_r4(mesh_i, mesh_o, routehandle, frac_o, rc) dstfracfield = ESMF_FieldCreate(mesh_o, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (norm_by_fracs) then + normtype = ESMF_NORMTYPE_FRACAREA + else + normtype = ESMF_NORMTYPE_DSTAREA + end if + ! Create route handle to map field_model to field_data call ESMF_FieldRegridStore(field_i, field_o, routehandle=routehandle, & - regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=ESMF_NORMTYPE_FRACAREA, & + regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=normtype, & srcMaskValues=(/srcMaskValue/), & dstMaskValues=(/dstMaskValue/), & srcTermProcessing=srcTermProcessing_Value, & @@ -83,11 +100,21 @@ subroutine create_routehandle_r4(mesh_i, mesh_o, routehandle, frac_o, rc) end subroutine create_routehandle_r4 !=============================================================== - subroutine create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o, rc) + subroutine create_routehandle_r8(mesh_i, mesh_o, norm_by_fracs, routehandle, frac_o, rc) ! input/output variables type(ESMF_Mesh) , intent(in) :: mesh_i type(ESMF_Mesh) , intent(in) :: mesh_o + ! If norm_by_fracs is .true., then remapping is done using ESMF_NORMTYPE_FRACAREA; + ! otherwise, remapping is done using ESMF_NORMTYPE_DSTAREA. FRACAREA normalization + ! adds a normalization factor of the fraction of the unmasked source grid that + ! overlaps with a destination cell. FRACAREA normalization is appropriate when you + ! want to treat values outside the mask as missing values that shouldn't contribute + ! to the average (this is appropriate for most fields); DSTAREA normalization is + ! appropriate when you want to treat values outside the mask as 0 (this is + ! appropriate for PCT cover fields where we want the final value to be expressed as + ! percent of the entire gridcell area). + logical , intent(in) :: norm_by_fracs type(ESMF_RouteHandle) , intent(inout) :: routehandle real(r8), optional , intent(inout) :: frac_o(:) integer , intent(out) :: rc @@ -96,6 +123,7 @@ subroutine create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o, rc) integer :: srcMaskValue = 0 ! ignore source points where the mesh mask is 0 integer :: dstMaskValue = -987987 ! don't ingore any destination points integer :: srcTermProcessing_Value = 0 + type(ESMF_NormType_Flag) :: normtype type(ESMF_Field) :: field_i type(ESMF_Field) :: field_o type(ESMF_Field) :: dstfracfield @@ -112,9 +140,15 @@ subroutine create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o, rc) dstfracfield = ESMF_FieldCreate(mesh_o, ESMF_TYPEKIND_R8, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + if (norm_by_fracs) then + normtype = ESMF_NORMTYPE_FRACAREA + else + normtype = ESMF_NORMTYPE_DSTAREA + end if + ! Create route handle to map field_model to field_data call ESMF_FieldRegridStore(field_i, field_o, routehandle=routehandle, & - regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=ESMF_NORMTYPE_FRACAREA, & + regridmethod=ESMF_REGRIDMETHOD_CONSERVE, normType=normtype, & srcMaskValues=(/srcMaskValue/), & dstMaskValues=(/dstMaskValue/), & srcTermProcessing=srcTermProcessing_Value, & diff --git a/tools/mksurfdata_esmf/src/mkgdpMod.F90 b/tools/mksurfdata_esmf/src/mkgdpMod.F90 index 801ba5e604..cb37de8bbd 100644 --- a/tools/mksurfdata_esmf/src/mkgdpMod.F90 +++ b/tools/mksurfdata_esmf/src/mkgdpMod.F90 @@ -111,7 +111,8 @@ subroutine mkgdp(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh and get frac_o allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkglacierregionMod.F90 b/tools/mksurfdata_esmf/src/mkglacierregionMod.F90 index 176271a48c..267e46ddfa 100644 --- a/tools/mksurfdata_esmf/src/mkglacierregionMod.F90 +++ b/tools/mksurfdata_esmf/src/mkglacierregionMod.F90 @@ -117,7 +117,8 @@ subroutine mkglacierregion(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) if (ier/=0) call abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkglcmecMod.F90 b/tools/mksurfdata_esmf/src/mkglcmecMod.F90 index 9a20388023..343e772de7 100644 --- a/tools/mksurfdata_esmf/src/mkglcmecMod.F90 +++ b/tools/mksurfdata_esmf/src/mkglcmecMod.F90 @@ -284,7 +284,8 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh and get frac_o allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort(subname//" error in allocating frac_o") - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) @@ -580,7 +581,8 @@ subroutine mkglacier(file_mesh_i, file_data_i, mesh_o, glac_o, rc) ! Create a route handle between the input and output mesh and get frac_o allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkharvestMod.F90 b/tools/mksurfdata_esmf/src/mkharvestMod.F90 index ac5fbf69d0..6e219a2164 100644 --- a/tools/mksurfdata_esmf/src/mkharvestMod.F90 +++ b/tools/mksurfdata_esmf/src/mkharvestMod.F90 @@ -153,7 +153,8 @@ subroutine mkharvest(file_mesh_i, file_data_i, mesh_o, pioid_o, ntime, rc) ! NOTE: this must be done after mask_i is set in mesh_i if (.not. ESMF_RouteHandleIsCreated(routehandle_r8)) then allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle_r8, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle_r8, frac_o=frac_o, rc=rc) call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if diff --git a/tools/mksurfdata_esmf/src/mklaiMod.F90 b/tools/mksurfdata_esmf/src/mklaiMod.F90 index 84069bded4..ba1014c2b9 100644 --- a/tools/mksurfdata_esmf/src/mklaiMod.F90 +++ b/tools/mksurfdata_esmf/src/mklaiMod.F90 @@ -121,7 +121,8 @@ subroutine mklai(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mklanwatMod.F90 b/tools/mksurfdata_esmf/src/mklanwatMod.F90 index 90c8cf2b0a..b6802ea6b6 100644 --- a/tools/mksurfdata_esmf/src/mklanwatMod.F90 +++ b/tools/mksurfdata_esmf/src/mklanwatMod.F90 @@ -116,7 +116,8 @@ subroutine mkpctlak(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, rc) ! Create a route handle between the input and output mesh if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak)) then allocate(frac_o_mklak(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle_mklak, frac_o=frac_o_mklak, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle_mklak, frac_o=frac_o_mklak, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -255,7 +256,8 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) ! Create a route handle between the input and output mesh if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak)) then allocate(frac_o_mklak(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle_mklak, frac_o=frac_o_mklak, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle_mklak, frac_o=frac_o_mklak, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -402,7 +404,8 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkpeatMod.F90 b/tools/mksurfdata_esmf/src/mkpeatMod.F90 index 6dad412427..3acd93b755 100644 --- a/tools/mksurfdata_esmf/src/mkpeatMod.F90 +++ b/tools/mksurfdata_esmf/src/mkpeatMod.F90 @@ -108,7 +108,8 @@ subroutine mkpeat(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh and get frac_o allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkpftMod.F90 b/tools/mksurfdata_esmf/src/mkpftMod.F90 index 0062940513..3ab941c6a0 100644 --- a/tools/mksurfdata_esmf/src/mkpftMod.F90 +++ b/tools/mksurfdata_esmf/src/mkpftMod.F90 @@ -331,7 +331,8 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & if (.not. ESMF_RouteHandleIsCreated(routehandle)) then allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if diff --git a/tools/mksurfdata_esmf/src/mksoildepthMod.F90 b/tools/mksurfdata_esmf/src/mksoildepthMod.F90 index e8f09b20aa..9b7b44f709 100644 --- a/tools/mksurfdata_esmf/src/mksoildepthMod.F90 +++ b/tools/mksurfdata_esmf/src/mksoildepthMod.F90 @@ -87,7 +87,8 @@ subroutine mksoildepth(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Get the landmask from the file and reset the mesh mask based on that allocate(frac_o(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) allocate(frac_i(ns_i), stat=ier) if (ier/=0) call shr_sys_abort() @@ -106,7 +107,8 @@ subroutine mksoildepth(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle between the input and output mesh and get frac_o - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) do no = 1, ns_o diff --git a/tools/mksurfdata_esmf/src/mksoilfmaxMod.F90 b/tools/mksurfdata_esmf/src/mksoilfmaxMod.F90 index 65534ed00a..22cf159849 100644 --- a/tools/mksurfdata_esmf/src/mksoilfmaxMod.F90 +++ b/tools/mksurfdata_esmf/src/mksoilfmaxMod.F90 @@ -102,7 +102,8 @@ subroutine mksoilfmax(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) do n = 1, ns_o diff --git a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 b/tools/mksurfdata_esmf/src/mkurbanparMod.F90 index b25a34a34c..c16cb6682e 100644 --- a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 +++ b/tools/mksurfdata_esmf/src/mkurbanparMod.F90 @@ -184,7 +184,8 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & ! Create a route handle between the input and output mesh if (.not. ESMF_RouteHandleIsCreated(routehandle_mkurban)) then allocate(frac_o_mkurban(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle_mkurban, frac_o=frac_o_mkurban, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle_mkurban, frac_o=frac_o_mkurban, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -938,7 +939,8 @@ subroutine mkurban_topo(file_mesh_i, file_data_i, mesh_o, varname, elev_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o), stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) diff --git a/tools/mksurfdata_esmf/src/mkvocefMod.F90 b/tools/mksurfdata_esmf/src/mkvocefMod.F90 index ae95bb5a14..6aa8ba6ec4 100644 --- a/tools/mksurfdata_esmf/src/mkvocefMod.F90 +++ b/tools/mksurfdata_esmf/src/mkvocefMod.F90 @@ -122,7 +122,8 @@ subroutine mkvocef(file_mesh_i, file_data_i, mesh_o, pioid_o, lat_o, rc) ! Create a route handle between the input and output mesh allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i, mesh_o, routehandle, frac_o=frac_o, rc=rc) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) From 80e89efe96e8fb0b54a6513c45f675b2fb8a8d9f Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Sat, 12 Nov 2022 05:29:51 -0700 Subject: [PATCH 06/15] Map PCT fields using DSTAREA normalization Partially addresses ESCOMP/CTSM#1716 --- tools/mksurfdata_esmf/src/mkglcmecMod.F90 | 46 ++++++++------- tools/mksurfdata_esmf/src/mklanwatMod.F90 | 64 +++++++++++---------- tools/mksurfdata_esmf/src/mkpftMod.F90 | 34 +++++------ tools/mksurfdata_esmf/src/mkurbanparMod.F90 | 28 +++++---- 4 files changed, 94 insertions(+), 78 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mkglcmecMod.F90 b/tools/mksurfdata_esmf/src/mkglcmecMod.F90 index 343e772de7..25fa4a8edc 100644 --- a/tools/mksurfdata_esmf/src/mkglcmecMod.F90 +++ b/tools/mksurfdata_esmf/src/mkglcmecMod.F90 @@ -134,7 +134,7 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) integer , intent(out) :: rc ! local variables: - type(ESMF_RouteHandle) :: routehandle ! nearest neighbor routehandle + type(ESMF_RouteHandle) :: routehandle_nonorm type(ESMF_Mesh) :: mesh_i type(file_desc_t) :: pioid_i type(var_desc_t) :: pio_varid @@ -148,7 +148,7 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) integer :: nlev ! number of levels in input file integer , allocatable :: mask_i(:) real(r8), allocatable :: frac_i(:) - real(r8), allocatable :: frac_o(:) + real(r8), allocatable :: frac_o_nonorm(:) real(r8), allocatable :: area_i(:) real(r8), allocatable :: area_o(:) real(r8), allocatable :: data_pctglc_i(:) @@ -281,11 +281,13 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) call mkpio_iodesc_rawdata(mesh_i, 'PCT_GLC_GIC', pioid_i, pio_varid, pio_vartype, pio_iodesc, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - ! Create a route handle between the input and output mesh and get frac_o - allocate(frac_o(ns_o),stat=ier) - if (ier/=0) call shr_sys_abort(subname//" error in allocating frac_o") - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle, frac_o=frac_o, rc=rc) + ! Create a route handle between the input and output mesh and get frac_o_nonorm + allocate(frac_o_nonorm(ns_o),stat=ier) + if (ier/=0) call shr_sys_abort(subname//" error in allocating frac_o_nonorm") + ! Note that norm_by_fracs is false in the following because this routehandle is + ! used to map fields that are expressed in terms of % of the grid cell. + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_nonorm, frac_o=frac_o_nonorm, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) @@ -308,11 +310,11 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) data_pctglc_i(:) = data_pctglc_gic_i(:) + data_pctglc_icesheet_i(:) ! Map level of data to output grid - call regrid_rawdata(mesh_i, mesh_o, routehandle, data_pctglc_i, data_pctglc_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, data_pctglc_i, data_pctglc_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call regrid_rawdata(mesh_i, mesh_o, routehandle, data_pctglc_gic_i, data_pctglc_gic_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, data_pctglc_gic_i, data_pctglc_gic_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call regrid_rawdata(mesh_i, mesh_o, routehandle, data_pctglc_icesheet_i, data_pctglc_icesheet_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, data_pctglc_icesheet_i, data_pctglc_icesheet_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Compute output variables @@ -470,7 +472,7 @@ subroutine mkglcmec(file_mesh_i, file_data_i, mesh_o, pioid_o, rc) end if ! Deallocate dynamic memory - call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) + call ESMF_RouteHandleDestroy(routehandle_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() @@ -506,14 +508,14 @@ subroutine mkglacier(file_mesh_i, file_data_i, mesh_o, glac_o, rc) integer , intent(out) :: rc ! ! local variables - type(ESMF_RouteHandle) :: routehandle + type(ESMF_RouteHandle) :: routehandle_nonorm type(ESMF_Mesh) :: mesh_i type(file_desc_t) :: pioid integer :: ni,no,k integer :: ns_i, ns_o integer , allocatable :: mask_i(:) real(r8), allocatable :: frac_i(:) - real(r8), allocatable :: frac_o(:) + real(r8), allocatable :: frac_o_nonorm(:) real(r8), allocatable :: area_i(:) real(r8), allocatable :: area_o(:) real(r8), allocatable :: glac_i(:) ! input grid: percent glac @@ -578,18 +580,20 @@ subroutine mkglacier(file_mesh_i, file_data_i, mesh_o, glac_o, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After mkpio_getrawdata in "//trim(subname)) - ! Create a route handle between the input and output mesh and get frac_o - allocate(frac_o(ns_o),stat=ier) + ! Create a route handle between the input and output mesh and get frac_o_nonorm + allocate(frac_o_nonorm(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle, frac_o=frac_o, rc=rc) + ! Note that norm_by_fracs is false in the following because this routehandle is + ! used to map fields that are expressed in terms of % of the grid cell. + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_nonorm, frac_o=frac_o_nonorm, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) ! Area-average percent cover on input grid to output grid (with correction for landmask) ! Note that percent cover is in terms of total grid area. ! Regrid glac_i to glac_o - call regrid_rawdata(mesh_i, mesh_o, routehandle, glac_i, glac_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, glac_i, glac_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return do no = 1,ns_o if (glac_o(no) < 1._r8) then @@ -608,13 +612,13 @@ subroutine mkglacier(file_mesh_i, file_data_i, mesh_o, glac_o, rc) enddo ! Check global areas - call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o, & + call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o_nonorm, & glac_i, glac_o, "pct glacier", percent=.true., ndiag=ndiag, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Deallocate dynamic memory - deallocate (glac_i, frac_o, mask_i) - call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) + deallocate (glac_i, frac_o_nonorm, mask_i) + call ESMF_RouteHandleDestroy(routehandle_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() diff --git a/tools/mksurfdata_esmf/src/mklanwatMod.F90 b/tools/mksurfdata_esmf/src/mklanwatMod.F90 index b6802ea6b6..b4160a830c 100644 --- a/tools/mksurfdata_esmf/src/mklanwatMod.F90 +++ b/tools/mksurfdata_esmf/src/mklanwatMod.F90 @@ -25,8 +25,8 @@ module mklanwatMod public :: mkwetlnd ! make %wetland public :: update_max_array_lake ! Update the maximum lake percent - real(r8), allocatable :: frac_o_mklak(:) - type(ESMF_RouteHandle) :: routehandle_mklak + real(r8), allocatable :: frac_o_mklak_nonorm(:) + type(ESMF_RouteHandle) :: routehandle_mklak_nonorm character(len=*) , parameter :: u_FILE_u = & __FILE__ @@ -114,10 +114,12 @@ subroutine mkpctlak(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle between the input and output mesh - if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak)) then - allocate(frac_o_mklak(ns_o)) - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle_mklak, frac_o=frac_o_mklak, rc=rc) + if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak_nonorm)) then + allocate(frac_o_mklak_nonorm(ns_o)) + ! Note that norm_by_fracs is false in the following because this routehandle is + ! used to map fields that are expressed in terms of % of the grid cell. + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_mklak_nonorm, frac_o=frac_o_mklak_nonorm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -138,14 +140,14 @@ subroutine mkpctlak(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Regrid lake_i to lake_o - call regrid_rawdata(mesh_i, mesh_o, routehandle_mklak, lake_i, lake_o, rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_mklak_nonorm, lake_i, lake_o, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return do no = 1,ns_o if (lake_o(no) < 1.) lake_o(no) = 0. enddo ! Check global areas - call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o_mklak, & + call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o_mklak_nonorm, & lake_i, lake_o, "pct lake", percent=.true., ndiag=ndiag, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -159,8 +161,8 @@ subroutine mkpctlak(file_mesh_i, file_data_i, mesh_o, lake_o, pioid_o, rc) ! Release memory if (mksrf_fdynuse == ' ') then ! ...else we will reuse it - deallocate(frac_o_mklak) - call ESMF_RouteHandleDestroy(routehandle_mklak, nogarbage = .true., rc=rc) + deallocate(frac_o_mklak_nonorm) + call ESMF_RouteHandleDestroy(routehandle_mklak_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() end if call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) @@ -193,9 +195,11 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) ! local variables type(ESMF_Mesh) :: mesh_i + type(ESMF_RouteHandle) :: routehandle type(file_desc_t) :: pioid_i integer , allocatable :: mask_i(:) real(r8), allocatable :: rmask_i(:) + real(r8), allocatable :: frac_o(:) real(r8), allocatable :: lakedepth_i(:) ! iput grid: lake depth (m) real(r8), allocatable :: lakedepth_o(:) ! output grid: lake depth (m) integer :: ni,no,k ! indices @@ -254,13 +258,11 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle between the input and output mesh - if (.not. ESMF_RouteHandleIsCreated(routehandle_mklak)) then - allocate(frac_o_mklak(ns_o)) - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle_mklak, frac_o=frac_o_mklak, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) - end if + allocate(frac_o(ns_o)) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & + routehandle=routehandle, frac_o=frac_o, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) ! ---------------------------------------- ! Create lake parameter (lakedepth) @@ -278,11 +280,11 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) ! regrid lakedepth_i to lakedepth_o - this also returns lakedepth_i to be used in the global sums below allocate (lakedepth_o(ns_o)); lakedepth_o(:) = spval - call regrid_rawdata(mesh_i, mesh_o, routehandle_mklak, lakedepth_i, lakedepth_o, rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle, lakedepth_i, lakedepth_o, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After regrid_rawdata for lakedepth in "//trim(subname)) do no = 1,ns_o - if (frac_o_mklak(no) == 0._r8) then + if (frac_o(no) == 0._r8) then lakedepth_o(no) = 10._r8 end if enddo @@ -300,7 +302,7 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) ! Check global areas for lake depth call output_diagnostics_continuous(mesh_i, mesh_o, & lakedepth_i, lakedepth_o, "lake depth", "m", & - ndiag=ndiag, rc=rc, mask_i=mask_i, frac_o=frac_o_mklak) + ndiag=ndiag, rc=rc, mask_i=mask_i, frac_o=frac_o) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! ---------------------------------------- @@ -310,8 +312,8 @@ subroutine mklakdep(file_mesh_i, file_data_i, mesh_o, pioid_o, fsurdat, rc) call pio_closefile(pioid_i) ! Release memory - deallocate(frac_o_mklak) - call ESMF_RouteHandleDestroy(routehandle_mklak, nogarbage = .true., rc=rc) + deallocate(frac_o) + call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() @@ -339,12 +341,12 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) integer , intent(out) :: rc ! local variables - type(ESMF_RouteHandle) :: routehandle + type(ESMF_RouteHandle) :: routehandle_nonorm type(ESMF_Mesh) :: mesh_i type(file_desc_t) :: pioid_i integer , allocatable :: mask_i(:) real(r8), allocatable :: rmask_i(:) - real(r8), allocatable :: frac_o(:) + real(r8), allocatable :: frac_o_nonorm(:) real(r8), allocatable :: swmp_i(:) ! input grid: percent wetland integer :: ni,no,k ! indices integer :: ns_i,ns_o ! local sizes @@ -403,9 +405,11 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle between the input and output mesh - allocate(frac_o(ns_o)) - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle, frac_o=frac_o, rc=rc) + allocate(frac_o_nonorm(ns_o)) + ! Note that norm_by_fracs is false in the following because this routehandle is + ! used to map fields that are expressed in terms of % of the grid cell. + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_nonorm, frac_o=frac_o_nonorm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) @@ -416,7 +420,7 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! regrid swmp_i to swmp_o - this also returns swmp_i to be used in the global sums below - call regrid_rawdata(mesh_i, mesh_o, routehandle, swmp_i, swmp_o, rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, swmp_i, swmp_o, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After regrid_data for wetland") if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -425,7 +429,7 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) enddo ! Check global areas - call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o, & + call output_diagnostics_area(mesh_i, mesh_o, mask_i, frac_o_nonorm, & swmp_i, swmp_o, "pct wetland", percent=.true., ndiag=ndiag, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -433,7 +437,7 @@ subroutine mkwetlnd(file_mesh_i, file_data_i, mesh_o, swmp_o, rc) call pio_closefile(pioid_i) ! Release memory - call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) + call ESMF_RouteHandleDestroy(routehandle_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() diff --git a/tools/mksurfdata_esmf/src/mkpftMod.F90 b/tools/mksurfdata_esmf/src/mkpftMod.F90 index 3ab941c6a0..cb3899fd44 100644 --- a/tools/mksurfdata_esmf/src/mkpftMod.F90 +++ b/tools/mksurfdata_esmf/src/mkpftMod.F90 @@ -24,8 +24,8 @@ module mkpftMod integer :: m ! index character(len=35) :: veg(0:maxpft) ! vegetation types - real(r8), allocatable :: frac_o(:) - type(ESMF_RouteHandle) :: routehandle + real(r8), allocatable :: frac_o_nonorm(:) + type(ESMF_RouteHandle) :: routehandle_nonorm character(len=*) , parameter :: u_FILE_u = & __FILE__ @@ -250,7 +250,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & real(r8), allocatable :: glob_gpft_o(:) ! output global area pfts integer :: ier, rcode ! error status real(r8) :: relerr = 0.0001_r8 ! max error: sum overlap wts ne 1 - character(len=*), parameter :: subname = 'mkpf' + character(len=*), parameter :: subname = 'mkpft' !----------------------------------------------------------------------- if (root_task) then @@ -326,13 +326,15 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & endif ! ---------------------------------------- - ! Create a route handle between the input and output mesh and get frac_o + ! Create a route handle between the input and output mesh and get frac_o_nonorm ! ---------------------------------------- - if (.not. ESMF_RouteHandleIsCreated(routehandle)) then - allocate(frac_o(ns_o),stat=ier) + if (.not. ESMF_RouteHandleIsCreated(routehandle_nonorm)) then + allocate(frac_o_nonorm(ns_o),stat=ier) if (ier/=0) call shr_sys_abort() - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle, frac_o=frac_o, rc=rc) + ! Note that norm_by_fracs is false in the following because this routehandle is + ! used to map fields that are expressed in terms of % of the grid cell. + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_nonorm, frac_o=frac_o_nonorm, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -340,7 +342,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & ! ---------------------------------------- ! Determine pctlnd_o(:) (in/out argument) ! ---------------------------------------- - pctlnd_o(:) = frac_o(:) * 100._r8 + pctlnd_o(:) = frac_o_nonorm(:) * 100._r8 ! ---------------------------------------- ! Determine pct_nat_pft_o(:,:) @@ -357,7 +359,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & if (chkerr(rc,__LINE__,u_FILE_u)) return ! Regrid to determine pctnatveg_o - call regrid_rawdata(mesh_i, mesh_o, routehandle, pctnatveg_i, pctnatveg_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pctnatveg_i, pctnatveg_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate(pct_nat_pft_i(0:num_natpft,ns_i)) @@ -375,7 +377,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & end do ! Readgrid to determine pct_nat_pft_o - call regrid_rawdata(mesh_i, mesh_o, routehandle, pct_nat_pft_i, pct_nat_pft_o, 0, num_natpft, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pct_nat_pft_i, pct_nat_pft_o, 0, num_natpft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Rescale pct_nat_pft_o @@ -420,7 +422,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & if (ier/=0) call shr_sys_abort() call mkpio_get_rawdata(pioid, 'PCT_CROP', mesh_i, pctcrop_i, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return - call regrid_rawdata(mesh_i, mesh_o, routehandle, pctcrop_i, pctcrop_o, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pctcrop_i, pctcrop_o, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return allocate(pct_cft_i(1:num_cft,ns_i), stat=ier) @@ -464,7 +466,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & end do ! Readgrid pct_cft_i to determine pct_cft_o - call regrid_rawdata(mesh_i, mesh_o, routehandle, pct_cft_i, pct_cft_o, 1, num_cft, rc=rc) + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pct_cft_i, pct_cft_o, 1, num_cft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Rescale pct_nat_pft_o @@ -538,7 +540,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & loc_gpft_o(:) = 0. do no = 1,ns_o do m = 0, numpft_i-1 - loc_gpft_o(m) = loc_gpft_o(m) + pctpft_o(no,m) * area_o(no) * frac_o(no) + loc_gpft_o(m) = loc_gpft_o(m) + pctpft_o(no,m) * area_o(no) * frac_o_nonorm(no) end do end do do m = 0,numpft_i-1 @@ -562,8 +564,8 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & ! Clean up memory if (mksrf_fdynuse == ' ') then ! ...else we will reuse it - deallocate(frac_o) - call ESMF_RouteHandleDestroy(routehandle, nogarbage = .true., rc=rc) + deallocate(frac_o_nonorm) + call ESMF_RouteHandleDestroy(routehandle_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() end if call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) diff --git a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 b/tools/mksurfdata_esmf/src/mkurbanparMod.F90 index c16cb6682e..05ec0cafe5 100644 --- a/tools/mksurfdata_esmf/src/mkurbanparMod.F90 +++ b/tools/mksurfdata_esmf/src/mkurbanparMod.F90 @@ -42,8 +42,8 @@ module mkurbanparMod ! private data members: ! flag to indicate nodata for index variables in output file: integer , parameter :: index_nodata = 0 - real(r8) , allocatable :: frac_o_mkurban(:) - type(ESMF_RouteHandle) :: routehandle_mkurban + real(r8) , allocatable :: frac_o_mkurban_nonorm(:) + type(ESMF_RouteHandle) :: routehandle_mkurban_nonorm character(len=*), parameter :: modname = 'mkurbanparMod' private :: index_nodata @@ -182,10 +182,10 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & if (chkerr(rc,__LINE__,u_FILE_u)) return ! Create a route handle between the input and output mesh - if (.not. ESMF_RouteHandleIsCreated(routehandle_mkurban)) then - allocate(frac_o_mkurban(ns_o)) - call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.true., & - routehandle=routehandle_mkurban, frac_o=frac_o_mkurban, rc=rc) + if (.not. ESMF_RouteHandleIsCreated(routehandle_mkurban_nonorm)) then + allocate(frac_o_mkurban_nonorm(ns_o)) + call create_routehandle_r8(mesh_i=mesh_i, mesh_o=mesh_o, norm_by_fracs=.false., & + routehandle=routehandle_mkurban_nonorm, frac_o=frac_o_mkurban_nonorm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After create routehandle in "//trim(subname)) end if @@ -203,7 +203,9 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & call ESMF_VMLogMemInfo("After mkpio_getrawdata in "//trim(subname)) ! Regrid input data to model resolution - call regrid_rawdata(mesh_i, mesh_o, routehandle_mkurban, data_i, data_o, 1, numurbl, rc) + ! + ! Use a nonorm mapper because we're mapping a field expressed as % of the grid cell area + call regrid_rawdata(mesh_i, mesh_o, routehandle_mkurban_nonorm, data_i, data_o, 1, numurbl, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_VMLogMemInfo("After regrid_data for in "//trim(subname)) @@ -301,7 +303,11 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & if (allocated(data_o)) deallocate(data_o) allocate(data_o(max_regions, ns_o), stat=ier) if (ier/=0) call shr_sys_abort('error allocating data_i(max_regions, ns_o)') - call regrid_rawdata(mesh_i, mesh_o, routehandle_mkurban, data_i, data_o, 1, nregions, rc) + ! This regridding could be done either with or without fracarea normalization, + ! because we just use it to find a dominant value. We use nonorm because we already + ! have a nonorm mapper for the sake of PCTURB and this way we don't need to make a + ! separate mapper with fracarea normalization. + call regrid_rawdata(mesh_i, mesh_o, routehandle_mkurban_nonorm, data_i, data_o, 1, nregions, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Now find dominant region in each output gridcell - this is identical to the maximum index @@ -319,7 +325,7 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & end if ! Output diagnostics - call output_diagnostics_index(mesh_i, mesh_o, mask_i, frac_o_mkurban, & + call output_diagnostics_index(mesh_i, mesh_o, mask_i, frac_o_mkurban_nonorm, & 1, max_regions, region_i, region_o, 'urban region', ndiag, rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() @@ -331,9 +337,9 @@ subroutine mkurban(file_mesh_i, file_data_i, mesh_o, pcturb_o, & ! TODO: determine what to deallocate ! deallocate (urban_classes_gcell_i, urban_classes_gcell_o, region_i) if (mksrf_fdynuse == ' ') then ! ...else we will reuse it - deallocate(frac_o_mkurban) + deallocate(frac_o_mkurban_nonorm) call ESMF_VMLogMemInfo("Before destroy operation in "//trim(subname)) - call ESMF_RouteHandleDestroy(routehandle_mkurban, nogarbage = .true., rc=rc) + call ESMF_RouteHandleDestroy(routehandle_mkurban_nonorm, nogarbage = .true., rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) call shr_sys_abort() end if call ESMF_MeshDestroy(mesh_i, nogarbage = .true., rc=rc) From dc8cb25b930577d97f6559cbca590067fb63cd53 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 15 Nov 2022 13:14:19 -0700 Subject: [PATCH 07/15] Remove / rework some redundant adjustments based on pctlnd_pft Do some rework of adjustments now that we are mapping pctnatveg and pctcrop with DSTAREA rather than FRACAREA normalization - so they end up as percent of the grid cell rather than percent of the land area. The end result won't be identical in all cases, but in practice it should do the same thing as before. - Zero out pctnatveg and pctcrop if their areas are less than 1e-6, rather than the previous logic of zeroing these if pctlnd is less than 1e-6. (This zeroing is now done in mkpft rather than in the main mksurfdata program.) - In deciding whether to set pct_nat_pft or pct_cft to a fixed value, only check for tiny values of pctnatveg / pctcrop, rather than checking for a tiny value of pctlnd (it should now be the case that a tiny value of pctlnd implies tiny values of pctnatveg and pctcrop as well). - Remove the check for no vegetation outside the pft mask. This check no longer seems important and risks doing more harm than good if the threshold for this check differs slightly from the threshold used to zero pctnatveg / pctcrop. --- tools/mksurfdata_esmf/src/mkpftMod.F90 | 16 ++++++---------- tools/mksurfdata_esmf/src/mksurfdata.F90 | 13 ------------- 2 files changed, 6 insertions(+), 23 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mkpftMod.F90 b/tools/mksurfdata_esmf/src/mkpftMod.F90 index cb3899fd44..162b462db2 100644 --- a/tools/mksurfdata_esmf/src/mkpftMod.F90 +++ b/tools/mksurfdata_esmf/src/mkpftMod.F90 @@ -380,16 +380,14 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pct_nat_pft_i, pct_nat_pft_o, 0, num_natpft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Rescale pct_nat_pft_o + ! Rescale pct_nat_pft_o, and set tiny pctnatveg to 0 do no = 1,ns_o - if (pctnatveg_o(no) > 0._r8) then + if (pctnatveg_o(no) >= 1.0e-6_r8) then do m = 0,num_natpft pct_nat_pft_o(m,no) = pct_nat_pft_o(m,no) / (pctnatveg_o(no) * 0.01_r8) end do else - pct_nat_pft_o(0:num_natpft,no) = 0._r8 - end if - if (pctlnd_o(no) < 1.0e-6 .or. pctnatveg_o(no) < 1.0e-6) then + pctnatveg_o(no) = 0._r8 pct_nat_pft_o(0,no) = 100._r8 pct_nat_pft_o(1:num_natpft,no) = 0._r8 end if @@ -469,16 +467,14 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pct_cft_i, pct_cft_o, 1, num_cft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Rescale pct_nat_pft_o + ! Rescale pct_cft_o, and set tiny pctcrop to 0 do no = 1,ns_o - if (pctcrop_o(no) > 0._r8) then + if (pctcrop_o(no) >= 1.0e-6_r8) then do m = 1,num_cft pct_cft_o(m,no) = pct_cft_o(m,no) / (pctcrop_o(no) * 0.01_r8) end do else - pct_cft_o(1:num_cft,no) = 0._r8 - end if - if (pctlnd_o(no) < 1.0e-6 .or. pctcrop_o(no) < 1.0e-6) then + pctcrop_o(no) = 0._r8 pct_cft_o(1,no) = 100._r8 pct_cft_o(2:num_cft,no) = 0._r8 end if diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index 526a999a84..d637df378d 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -393,10 +393,6 @@ program mksurfdata call pctnatpft(n)%set_pct_l2g(0._r8) call pctcft(n)%set_pct_l2g(0._r8) end if - if (pctlnd_pft(n) < 1.e-6_r8) then - call pctnatpft(n)%set_pct_l2g(0._r8) - call pctcft(n)%set_pct_l2g(0._r8) - end if landfrac_pft(n) = pctlnd_pft(n)/100._r8 end do if (fsurdat /= ' ') then @@ -1242,15 +1238,6 @@ subroutine normalize_and_check_landuse(ns_o) end do - ! Make sure that there is no vegetation outside the pft mask - do n = 1,ns_o - if (pctlnd_pft(n) < 1.e-6_r8 .and. (pctnatpft(n)%get_pct_l2g() > 0 .or. pctcft(n)%get_pct_l2g() > 0)) then - write (6,*)'vegetation found outside the pft mask at n=',n - write (6,*)'pctnatveg,pctcrop=', pctnatpft(n)%get_pct_l2g(), pctcft(n)%get_pct_l2g() - call shr_sys_abort() - end if - end do - ! Make sure that sums at the landunit level all add to 100% ! (Note that we don't check pctglcmec here, because it isn't computed at the point ! that this subroutine is called -- but the check of sum(pctglcmec) is done in From d83278cb6b86608fba0b36419db3633f3c29dafd Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 16 Nov 2022 13:01:33 -0700 Subject: [PATCH 08/15] Use LANDFRAC instead of mask for setting pctlnd_pft --- tools/mksurfdata_esmf/src/mkpftMod.F90 | 16 +++++++++-- tools/mksurfdata_esmf/src/mksurfdata.F90 | 35 ++++++++++-------------- 2 files changed, 28 insertions(+), 23 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mkpftMod.F90 b/tools/mksurfdata_esmf/src/mkpftMod.F90 index 162b462db2..54827e7ad3 100644 --- a/tools/mksurfdata_esmf/src/mkpftMod.F90 +++ b/tools/mksurfdata_esmf/src/mkpftMod.F90 @@ -209,7 +209,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & character(len=*) , intent(in) :: file_mesh_i ! input mesh file name character(len=*) , intent(in) :: file_data_i ! input data file name type(ESMF_Mesh) , intent(in) :: mesh_o ! model mesh - real(r8) , intent(inout) :: pctlnd_o(:) ! output grid:%land/gridcell + real(r8) , intent(out) :: pctlnd_o(:) ! output grid:%land/gridcell type(pct_pft_type), intent(inout) :: pctnatpft_o(:) ! natural PFT cover type(pct_pft_type), intent(inout) :: pctcft_o(:) ! crop (CFT) cover @@ -234,6 +234,7 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & real(r8), allocatable :: output_pct_cft_o(:,:) integer , allocatable :: mask_i(:) real(r8), allocatable :: frac_i(:) + real(r8), allocatable :: pctlnd_i(:) ! input land fraction real(r8), allocatable :: pctnatveg_i(:) ! input natural veg percent (% of grid cell) real(r8), allocatable :: pctnatveg_o(:) ! output natural veg percent (% of grid cell) real(r8), allocatable :: pctcrop_i(:) ! input all crop percent (% of grid cell) @@ -340,9 +341,18 @@ subroutine mkpft(file_mesh_i, file_data_i, mesh_o, pctlnd_o, pctnatpft_o, & end if ! ---------------------------------------- - ! Determine pctlnd_o(:) (in/out argument) + ! Determine pctlnd_o(:) ! ---------------------------------------- - pctlnd_o(:) = frac_o_nonorm(:) * 100._r8 + allocate(pctlnd_i(ns_i), stat=ier) + if (ier/=0) call shr_sys_abort('error allocating pctlnd_i') + + call mkpio_get_rawdata(pioid, 'LANDFRAC', mesh_i, pctlnd_i, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + call regrid_rawdata(mesh_i, mesh_o, routehandle_nonorm, pctlnd_i, pctlnd_o, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + ! Convert from fraction to percent: + pctlnd_o(:) = pctlnd_o(:) * 100._r8 ! ---------------------------------------- ! Determine pct_nat_pft_o(:,:) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index d637df378d..a2eecb1af5 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -878,26 +878,21 @@ program mksurfdata if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkpft') call pio_syncfile(pioid) -! ! Consistency check on input land fraction -! ! pctlnd_pft was calculated ABOVE -! ! TODO This error check serves no purpose now that mkpft calls -! ! create_routehandle_r8 only once and, therefore, doesn't update -! ! frac_o and pctlnd_o. Delete? (slevis) -! do n = 1,lsize_o -! if (pctlnd_pft_dyn(n) /= pctlnd_pft(n)) then -! if (root_task) then -! write(ndiag,*) subname,' error: pctlnd_pft for dynamics data = ',& -! pctlnd_pft_dyn(n), ' not equal to pctlnd_pft for surface data = ',& -! pctlnd_pft(n),' at n= ',n -! if ( trim(fname) == ' ' )then -! write(ndiag,*) ' PFT string = ',trim(string) -! else -! write(ndiag,*) ' PFT file = ', fname -! end if -! end if -! call shr_sys_abort() -! end if -! end do + ! Consistency check on input land fraction + ! pctlnd_pft was calculated ABOVE + do n = 1,lsize_o + if (pctlnd_pft_dyn(n) /= pctlnd_pft(n)) then + write(ndiag,*) subname,' error: pctlnd_pft for dynamics data = ',& + pctlnd_pft_dyn(n), ' not equal to pctlnd_pft for surface data = ',& + pctlnd_pft(n),' at n= ',n + if ( trim(fname) == ' ' )then + write(ndiag,*) ' PFT string = ',trim(string) + else + write(ndiag,*) ' PFT file = ', fname + end if + call shr_sys_abort() + end if + end do ! Create harvesting data at model resolution ! Output data is written in mkharvest From e29e26c81d670a9e33300db167ef88fed4e92bb9 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 22 Nov 2022 15:56:05 -0700 Subject: [PATCH 09/15] Rework calculation of wetland area and normalization of landunits This rework makes things work correctly given that the landunit areas coming into normalize_and_check_landuse are in % of grid cell, but the outputs are in % of land area. Addresses much of the remainder of ESCOMP/CTSM#1716. Especially see https://github.com/ESCOMP/CTSM/issues/1716#issuecomment-1322557821 --- tools/mksurfdata_esmf/src/mksurfdata.F90 | 159 ++++++++++++++++++----- 1 file changed, 124 insertions(+), 35 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index a2eecb1af5..eb91be30c7 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -1046,6 +1046,11 @@ 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. ! + ! Coming into this subroutine, landunit areas are expressed as percent of the + ! gridcell and are NOT normalized to sum to 100%. Coming out of this subroutine, + ! landunit areas are expressed as percent of the land portion of the gridcell and + ! ARE normalized to sum to 100%. + ! ! input/output variables: integer, intent(in) :: ns_o @@ -1054,49 +1059,23 @@ subroutine normalize_and_check_landuse(ns_o) 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) :: pct_land ! area considered to be land (% of grid cell) + real(r8) :: frac_land ! area considered to be land (fraction of grid cell) 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 + character(len=32) :: subname = 'normalize_and_check_landuse' ! subroutine name !----------------------------------------------------------------------- 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 + ! ensure that wt is zero (not a very small number such as ! 1e-16) where it really should be zero - 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 - - if (pctlnd_pft(n) < 1.e-6_r8) then - 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) - end if - - ! 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) + 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 if ( pctlak(n) < 0.0_r8 )then write(6,*) subname, ' ERROR: pctlak is negative!' @@ -1124,6 +1103,18 @@ subroutine normalize_and_check_landuse(ns_o) call shr_sys_abort() end if + ! 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) + 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 + 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' @@ -1133,13 +1124,111 @@ subroutine normalize_and_check_landuse(ns_o) call shr_sys_abort() end if - ! Natural vegetated cover is 100% minus the sum of all other landunits + ! Determine the percent of each grid cell considered to be land. (See comments in + ! https://github.com/ESCOMP/CTSM/issues/1716 for details.) + ! + ! Start by using the land fraction field from the PFT raw data set: + pct_land = pctlnd_pft(n) + ! + ! But we don't want to overwrite special landunits or crop with ocean where these + ! special landunits extend beyond the PFT data's land fraction. In essence, this + ! is saying that we'll let special landunit area grow into the natveg area before + ! growing into ocean, but we'll have special landunit area grow into ocean before + ! growing into crop or any other special landunit area. (This check of special + ! landunit area is particularly important for glaciers, where we can have + ! floating ice shelves, so we can have a scenario where pctlnd_pft is 0 but we + ! have non-zero glacier cover and we want the final grid cell to be + ! glacier-covered.) (We could possibly do better by considering the land mask + ! from each special landunit raw dataset, and even better by mapping these + ! various notions of land mask onto each other, but that starts to get messy, and + ! relies on the trustworthiness of each raw dataset's land mask... this + ! formulation seems reasonable enough.) + ! + ! Note that we include pct_crop in the following, but NOT pct_natveg. The + ! assumption behind that is that pct_crop is more reliable and/or more important, + ! and should not be overwritten, whereas pct_natveg can be overwritten with a + ! special landunit. For example, consider a case where the PFT dataset specifies + ! 40% land, with 20% crop and 10% natveg (so, implicitly, 10% special landunits). + ! If the only special landunit is glacier and it has 15% cover, then we want to + ! end up with the glacier overwriting natural veg, so we end up with 20% crop, 5% + ! natveg, 15% glacier and 60% ocean. However, if glacier has 30% cover, then we + ! will assume that some of that glacier extends over the ocean rather than + ! overwriting crop, so we end up with 20% crop, 30% glacier, 50% ocean and 0% + ! natveg. + ! + ! Another reason for excluding pct_natveg from the following is more pragmatic: + ! Typically, we expect the initial sum of all landunit areas to approximately + ! equal pctlnd_pft. This means that, in a coastal grid cell, if we included all + ! landunit areas in the following sum, then in many cases we would take the + ! final pct_land from the landunit sum rather than from pctlnd_pft. But in this + ! scenario where we take pct_land from the landunit sum, it is likely that + ! pct_land will vary from year to year on the landuse timeseries file. This + ! variation from year to year will cause a renormalization of all landunits, + ! leading to changes in the areas of landunits that should really stay fixed + ! from one year to the next. By excluding pct_natveg we give more wiggle room: + ! it will usually be the case that we take the final pct_land from pctlnd_pft, + ! which stays fixed in time, so that the renormalization stays the same from one + ! year to the next. The possible downside is that pct_natveg may end up varying + ! more than is ideal, but this seems better than letting all of the other + ! landunits vary when they should stay fixed. + ! + ! Also, note that this landunit sum agrees with the suma sums elsewhere; this + ! agreement may be important in some cases (so that if we changed the set of + ! landunits included in the sum here, some other changes may be needed below.) + pct_land = max(pct_land, suma) + ! + ! Make sure that we're not ending up with > 100% land area. This can arise if the + ! sum of special landunits exceeds 100% by roundoff; also, due to rounding + ! errors, pctlnd_pft may slightly differ from 100% when it should really be + ! exactly 100%; we want pct_land to be exactly 100% in this case: + if (pct_land > (100._r8 - 1.e-4_r8)) then + pct_land = 100._r8 + end if - new_total_natveg_pct = 100._r8 - suma - ! correct for rounding error: - new_total_natveg_pct = max(new_total_natveg_pct, 0._r8) + if (pct_land < 1.e-6_r8) then + ! If we have essentially 0 land area, set land area to exactly 0 and put all + ! area in wetlands (to simulate ocean). Note that, based on the formulation + ! for pct_land above, this should only arise if the non-natveg landunits + ! already have near-zero area, so the zeroing of these other landunits should + ! only result in changes near the roundoff level. + pct_land = 0._r8 + frac_land = 0._r8 + call pctnatpft(n)%set_pct_l2g(0._r8) + call pctcft(n)%set_pct_l2g(0._r8) + pctlak(n) = 0._r8 + pcturb(n) = 0._r8 + pctgla(n) = 0._r8 + pctwet(n) = 100._r8 + else + ! Fill the rest of the land area with natveg, then renormalize landunits so + ! that they are expressed as percent of the land area rather than percent of + ! the full gridcell area. + + ! First fill the rest of the land area with natveg: + new_total_natveg_pct = pct_land - suma + ! Based on the formulation of pct_land above, pct_land is guaranteed to be at + ! least as large as suma. However, correct it just in case there is rounding + ! error: + new_total_natveg_pct = max(new_total_natveg_pct, 0._r8) + + ! Now renormalize landunit areas so they are expressed as percent of the land + ! area rather than percent of the total gridcell area. Note that we have + ! already corrected pct_land above so that if it was slightly less than 100 + ! it was set to exactly 100, so we can check for any values less than 100 + ! here (rather than having some tolerance): + frac_land = pct_land / 100._r8 + if (pct_land < 100._r8) then + new_total_natveg_pct = new_total_natveg_pct / frac_land + call pctcft(n)%set_pct_l2g(pctcft(n)%get_pct_l2g() / frac_land) + pctlak(n) = pctlak(n) / frac_land + pcturb(n) = pcturb(n) / frac_land + pctgla(n) = pctgla(n) / frac_land + pctwet(n) = pctwet(n) / frac_land + end if - call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct) + ! Finally, set the actual pct_natveg: + call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct) + end if ! Confirm that we have done the rescaling correctly: now the sum of all landunits ! should be 100% within tol_loose From 277db0f9a36892470410e82d0e222c4bcecbaad1 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 1 Dec 2022 16:34:36 -0700 Subject: [PATCH 10/15] Add documentation of the specification of landunit areas --- doc/design/surface_dataset_areas.rst | 39 ++++++++++++++++++++++++ tools/mksurfdata_esmf/src/mksurfdata.F90 | 6 ++-- 2 files changed, 43 insertions(+), 2 deletions(-) create mode 100644 doc/design/surface_dataset_areas.rst diff --git a/doc/design/surface_dataset_areas.rst b/doc/design/surface_dataset_areas.rst new file mode 100644 index 0000000000..ce2c96fac2 --- /dev/null +++ b/doc/design/surface_dataset_areas.rst @@ -0,0 +1,39 @@ +.. sectnum:: + +.. contents:: + +================================== + Overview of this design document +================================== + +This document gives a high-level overview of the specification of sub-grid areas on surface datasets and the raw datasets used to create them. + +See also https://github.com/ESCOMP/CTSM/issues/1716 for related discussion. + +================================================= + Specification of landunit areas in raw datasets +================================================= + +In the high-resolution raw datasets used as input to the mksurfdata tool, landunit areas should be specified as percent of the grid cell, **not** percent of the land area. For example, on the urban raw dataset, if there is a grid cell that is 50% land and 50% ocean, and 60% of the land area is urban, the urban area in that grid cell should be given as 30%. + +One reason for this is that it makes reconciling the different landunit areas more straightforward if different raw datasets disagree about the land mask. In addition, this convention makes it easier to map raw datasets to the model resolution. For example, consider averaging two grid cells into a destination grid cell: one with 2% land area, of which 50% is glacier; and one with 100% land area, of which none is glacier. If we encode these as percent of the land area, we would have 50% glacier and 0% glacier, and then the final glacier cover would be 25%, suggesting that 25% of the land area is glacier, but this is incorrect. If we instead encode these as percent of the total grid cell area, we would have 1% glacier and 0% glacier, and then the final glacier cover would be 0.5%, suggesting that 0.5% of the total grid cell is glacier, which is correct. + +===================================================== + Specification of landunit areas in surface datasets +===================================================== + +In contrast to the raw datasets, landunit areas in surface datasets are specified as percent of the land area, **not** as percent of the total grid cell. This is because we don't know what the model's actual land fraction will be at the time when the surface datasets are created: instead, this land fraction is determined at runtime based on the ocean grid. + +=========================================================================================== + Procedure for converting landunit areas from percent of grid cell to percent of land area +=========================================================================================== + +There are a few important aspects to how we determine final landunit areas in mksurfdata: + +When mapping landunit areas from the raw data resolution to the model resolution, we initially want to maintain areas as percent of the total grid cell area. To achieve this, we set ``norm_by_fracs=.false.`` in the call to ``create_routehandle``, resulting in the use of ``ESMF_NORMTYPE_DSTAREA`` rather than ``ESMF_NORMTYPE_FRACAREA`` as a ``normType`` when computing mapping weights. Using ``FRACAREA`` normalization is appropriate when you want to treat areas outside of the source mask as missing values that do not contribute in any way to the final destination value. Using ``DSTAREA`` normalization, in contrast, essentially treats areas outside of the source mask as zeroes. ``FRACAREA`` normalization is appropriate for many surface dataset fields, but ``DSTAREA`` is appropriate for areas specified as percent of the grid cell. For example, if a source grid cell is entirely ocean, then we want to treat glacier area in that source grid cell as 0%. + +The conversion from percent of the grid cell area to percent of the land area happens in the subroutine ``normalize_and_check_landuse``. An important piece of doing this conversion is to determine an estimate of the land fraction for each model grid cell. This is not straightforward given the disparate land masks used for each raw dataset. We start by using the land fraction from the vegetation (PFT) raw dataset, with the assumption that that is probably the most reliable land mask. However, there are areas where using that land fraction is problematic, particularly where the areas of other landunits extend beyond the PFT's land mask. This is especially an issue for glaciers, where floating ice shelves can extend beyond the land-ocean border. To deal with this problem, if the sum of the areas of special landunits and crops exceeds the land fraction from the PFT data, we instead use that sum as an estimate of the land fraction. Exactly which landunits to include in this sum is a bit arbitrary. Arguably, the natural vegetated landunit should also be included in this sum. However, we ideally want to avoid changes in this estimated land fraction through time in transient datasets, which means that we generally want to use the PFT data's land fraction, only falling back on this landunit sum in exceptional cases. By excluding the natural vegetated area from this sum, we are more likely to use the PFT's land fraction. An implication of this choice is that we are more likely to replace natural vegetated areas with special landunits in cases where there are disagreements between the different raw datasets in coastal grid cells. For more detailed explanation and rationale, see some comments in ``normalize_and_check_landuse``. + +In grid cells where the estimated land fraction is essentially zero, we set the land cover to wetland, as a rough parameterization of ocean. This situation will only arise if the areas of all landunits on the raw datasets are essentially zero for the given grid cell, so we would have no information to choose any particular land cover for the grid cell. (This wetland area may end up being changed to bare ground at runtime, depending on the value of the ``convert_ocean_to_land`` namelist flag.) + +In grid cells where the estimated land fraction is greater than zero, we fill any unclaimed land area with the natural vegetated landunit. We then normalize all landunit areas based on the estimated land fraction so that they now specify areas as percent of the land area rather than as percent of the grid cell. diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index eb91be30c7..09fe291cd5 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -1189,8 +1189,10 @@ subroutine normalize_and_check_landuse(ns_o) ! If we have essentially 0 land area, set land area to exactly 0 and put all ! area in wetlands (to simulate ocean). Note that, based on the formulation ! for pct_land above, this should only arise if the non-natveg landunits - ! already have near-zero area, so the zeroing of these other landunits should - ! only result in changes near the roundoff level. + ! already have near-zero area (and the natveg landunit should also have + ! near-zero area in this case, because its area should be no greater than the + ! land fraction from the PFT raw dataset), so the zeroing of these other + ! landunits should only result in changes near the roundoff level. pct_land = 0._r8 frac_land = 0._r8 call pctnatpft(n)%set_pct_l2g(0._r8) From 4ee2a0012b613334cd204c49bdb119c72327b850 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 1 Dec 2022 16:58:03 -0700 Subject: [PATCH 11/15] Add example scenarios --- doc/design/surface_dataset_areas.rst | 30 ++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/doc/design/surface_dataset_areas.rst b/doc/design/surface_dataset_areas.rst index ce2c96fac2..de12f951da 100644 --- a/doc/design/surface_dataset_areas.rst +++ b/doc/design/surface_dataset_areas.rst @@ -37,3 +37,33 @@ The conversion from percent of the grid cell area to percent of the land area ha In grid cells where the estimated land fraction is essentially zero, we set the land cover to wetland, as a rough parameterization of ocean. This situation will only arise if the areas of all landunits on the raw datasets are essentially zero for the given grid cell, so we would have no information to choose any particular land cover for the grid cell. (This wetland area may end up being changed to bare ground at runtime, depending on the value of the ``convert_ocean_to_land`` namelist flag.) In grid cells where the estimated land fraction is greater than zero, we fill any unclaimed land area with the natural vegetated landunit. We then normalize all landunit areas based on the estimated land fraction so that they now specify areas as percent of the land area rather than as percent of the grid cell. + +=================== + Example scenarios +=================== + +The following example scenarios illustrate the operation of ``normalize_and_check_landuse``; in the following, any landunit not explicitly mentioned has 0% area: + +(a) With pctlnd_pft = 0% and all initial landunit areas 0%: wetland area = 100% + +(b) With pctlnd_pft = 0% and initial glacier area 1%: glacier area = 100% + +(c) With pctlnd_pft > 0% and all initial landunit areas 0%: natural vegetated area = 100% + +(d) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%: crop area = 50%, natural vegetated area = 50% + +(e) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 10%: crop area = 50%, natural vegetated area = 25%, glacier area = 25% + +(f) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 15%: crop area = 50%, natural vegetated area = 12.5%, glacier area = 37.5% + +(g) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 20%: crop area = 50%, glacier area = 50% + +(h) With pctlnd_pft = 40%, initial crop area 20%, natural vegetated area 10%, glacier area 30%: crop area = 40%, glacier area = 60% + +(i) With pctlnd_pft = 40%, initial crop area 0%, natural vegetated area 40%, glacier area 40%: glacier area = 100% + +(j) With pctlnd_pft = 2%, initial natural vegetated area 1%, glacier area 1%: natural vegetated area = 50%, glacier area = 50% + +(k) With pctlnd_pft = 2%, initial natural vegetated area 0%, glacier area 1%: natural vegetated area = 50%, glacier area = 50% + +(l) With pctlnd_pft = 2%, initial natural vegetated area 2%, glacier area 1%: natural vegetated area = 50%, glacier area = 50% From 344c12abe4c081fad3f95bf3eed2504ca8e7a27f Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 2 Dec 2022 05:47:17 -0700 Subject: [PATCH 12/15] Fix minor copy and paste error in error message --- tools/mksurfdata_esmf/src/mksurfdata.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index 09fe291cd5..e3ea832926 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -1008,7 +1008,7 @@ program mksurfdata if (root_task) write(ndiag, '(a,i8)') trim(subname)//" writing PCT_CROP_MAX" call get_pct_l2g_array(pctcft_max, pctcrop) 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 (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_CROP_MAX') 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) From d48239391fed0ed923c705ae5e98ab873b17ce15 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 2 Dec 2022 06:31:42 -0700 Subject: [PATCH 13/15] Write LANDFRAC_MKSURFDATA to surface dataset and landuse timeseries --- tools/mksurfdata_esmf/src/mkfileMod.F90 | 11 +++++++++++ tools/mksurfdata_esmf/src/mksurfdata.F90 | 20 +++++++++++++++++++- 2 files changed, 30 insertions(+), 1 deletion(-) diff --git a/tools/mksurfdata_esmf/src/mkfileMod.F90 b/tools/mksurfdata_esmf/src/mkfileMod.F90 index 6c3046a30b..335a9a508b 100644 --- a/tools/mksurfdata_esmf/src/mkfileMod.F90 +++ b/tools/mksurfdata_esmf/src/mkfileMod.F90 @@ -606,6 +606,17 @@ subroutine mkfile_define_vars(pioid, dynlanduse) call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_PFT', xtype=xtype, & long_name='land fraction from pft dataset (DIFF FROM landfrac USED IN SIMULATION, SHOWN IN HISTORY)', units='unitless') + if (.not. dynlanduse) then + call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_MKSURFDATA', xtype=xtype, & + long_name='land fraction used for renormalization of areas in mksurfdata (DIFF FROM landfrac USED IN SIMULATION)', & + units='unitless') + else + call mkpio_def_spatial_var(pioid=pioid, varname='LANDFRAC_MKSURFDATA', xtype=xtype, & + lev1name='time', & + long_name='land fraction used for renormalization of areas in mksurfdata (DIFF FROM landfrac USED IN SIMULATION)', & + units='unitless') + end if + if (.not. dynlanduse) then call mkpio_def_spatial_var(pioid=pioid, varname='PCT_NATVEG', xtype=xtype, & long_name='total percent natural vegetation landunit', units='unitless') diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index e3ea832926..0b304079bc 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -172,6 +172,9 @@ program mksurfdata real(r8), allocatable :: pctwet_orig(:) ! percent wetland of gridcell before dynamic land use adjustments real(r8), allocatable :: pctgla_orig(:) ! percent glacier of gridcell before dynamic land use adjustments + ! other variables written to file + real(r8), allocatable :: landfrac_mksurfdata(:) ! land fraction used for renormalization of areas + ! pio/esmf variables type(file_desc_t) :: pioid type(var_desc_t) :: pio_varid @@ -595,6 +598,7 @@ program mksurfdata ! Normalize land use and make sure things add up to 100% as well as ! checking that things are as they should be. + allocate(landfrac_mksurfdata(lsize_o)) call normalize_and_check_landuse(lsize_o) ! Write out sum of PFT's @@ -692,6 +696,10 @@ program mksurfdata call mkfile_output(pioid, mesh_model, 'PCT_CFT', pct_cft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for PCT_CFT') end if + + if (root_task) write(ndiag, '(a)') trim(subname)//" writing LANDFRAC_MKSURFDATA" + call mkfile_output(pioid, mesh_model, 'LANDFRAC_MKSURFDATA', landfrac_mksurfdata, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_MKSURFDATA') end if ! ---------------------------------------------------------------------- @@ -779,7 +787,7 @@ program mksurfdata ! landfrac_pft was calculated ABOVE if (root_task) write(ndiag, '(a)') trim(subname)//" writing land fraction calculated in fsurdata calc)" call mkfile_output(pioid, mesh_model, 'LANDFRAC_PFT', landfrac_pft, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output') + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_PFT') ! ----------------------------------------- ! Read in each dynamic pft landuse dataset @@ -993,6 +1001,13 @@ program mksurfdata call pio_syncfile(pioid) end if + if (root_task) write(ndiag, '(a,i8)') trim(subname)//" writing LANDFRAC_MKSURFDATA for year ",year + rcode = pio_inq_varid(pioid, 'LANDFRAC_MKSURFDATA', pio_varid) + call pio_setframe(pioid, pio_varid, int(ntim, kind=Pio_Offset_Kind)) + call mkfile_output(pioid, mesh_model, 'LANDFRAC_MKSURFDATA', landfrac_mksurfdata, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkfile_output for LANDFRAC_MKSURFDATA') + call pio_syncfile(pioid) + if (root_task) then write(ndiag,'(1x,80a1)') ('=',k=1,80) write(ndiag,*) @@ -1232,6 +1247,9 @@ subroutine normalize_and_check_landuse(ns_o) call pctnatpft(n)%set_pct_l2g(new_total_natveg_pct) end if + ! Save landfrac for output to file + landfrac_mksurfdata(n) = frac_land + ! Confirm that we have done the rescaling correctly: now the sum of all landunits ! should be 100% within tol_loose suma = pctlak(n) + pctwet(n) + pctgla(n) + pcturb(n) + pctcft(n)%get_pct_l2g() From 7bb66addd2d071eebbe742550267b3e5c0a7023b Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 5 Dec 2022 07:03:09 -0700 Subject: [PATCH 14/15] Stop hard-coding the settings of pole points --- tools/mksurfdata_esmf/src/mksurfdata.F90 | 37 +----------------------- 1 file changed, 1 insertion(+), 36 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index 0b304079bc..e0f0108d7c 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -389,13 +389,7 @@ program mksurfdata pctlnd_o=pctlnd_pft, pctnatpft_o=pctnatpft, pctcft_o=pctcft, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) call shr_sys_abort('error in calling mkdomain') - ! If have pole points on grid - set south pole to glacier - ! north pole is assumed as non-land do n = 1,lsize_o - if (abs((lat(n) - 90._r8)) < 1.e-6_r8) then - call pctnatpft(n)%set_pct_l2g(0._r8) - call pctcft(n)%set_pct_l2g(0._r8) - end if landfrac_pft(n) = pctlnd_pft(n)/100._r8 end do if (fsurdat /= ' ') then @@ -573,22 +567,9 @@ program mksurfdata end if ! ----------------------------------- - ! Adjust pctlak, pctwet, pcturb and pctgla + ! Save special land unit areas of surface dataset ! ----------------------------------- - do n = 1,lsize_o - - ! If have pole points on grid - set south pole to glacier - ! north pole is assumed as non-land - if (abs((lat(n) - 90._r8)) < 1.e-6_r8) then - pctlak(n) = 0._r8 - pctwet(n) = 0._r8 - pcturb(n) = 0._r8 - pctgla(n) = 100._r8 - end if - - end do - ! Save special land unit areas of surface dataset pctwet_orig(:) = pctwet(:) pctgla_orig(:) = pctgla(:) @@ -929,22 +910,6 @@ program mksurfdata pctwet(:) = pctwet_orig(:) pctgla(:) = pctgla_orig(:) - ! If have pole points on grid - set south pole to glacier - ! north pole is assumed as non-land - ! pctlak, pctwet, pcturb and pctgla were calculated ABOVE - ! pctnatpft and pctcft were calculated ABOVE - do n = 1,lsize_o - if (abs(lat(n) - 90._r8) < 1.e-6_r8) then - pctlak(n) = 0._r8 - pctwet(n) = 0._r8 - pcturb(n) = 0._r8 - pctgla(n) = 100._r8 - call pctnatpft(n)%set_pct_l2g(0._r8) - call pctcft(n)%set_pct_l2g(0._r8) - end if - - end do - ! Normalize land use and make sure things add up to 100% as well as ! checking that things are as they should be. if (root_task) then From 821fdc724199e05949e644f253eda254289e79d3 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Mon, 19 Dec 2022 16:24:22 -0700 Subject: [PATCH 15/15] Tweak a comment --- tools/mksurfdata_esmf/src/mksurfdata.F90 | 30 ++++++++++++++---------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/tools/mksurfdata_esmf/src/mksurfdata.F90 b/tools/mksurfdata_esmf/src/mksurfdata.F90 index e0f0108d7c..9678308800 100644 --- a/tools/mksurfdata_esmf/src/mksurfdata.F90 +++ b/tools/mksurfdata_esmf/src/mksurfdata.F90 @@ -1110,19 +1110,23 @@ subroutine normalize_and_check_landuse(ns_o) ! Start by using the land fraction field from the PFT raw data set: pct_land = pctlnd_pft(n) ! - ! But we don't want to overwrite special landunits or crop with ocean where these - ! special landunits extend beyond the PFT data's land fraction. In essence, this - ! is saying that we'll let special landunit area grow into the natveg area before - ! growing into ocean, but we'll have special landunit area grow into ocean before - ! growing into crop or any other special landunit area. (This check of special - ! landunit area is particularly important for glaciers, where we can have - ! floating ice shelves, so we can have a scenario where pctlnd_pft is 0 but we - ! have non-zero glacier cover and we want the final grid cell to be - ! glacier-covered.) (We could possibly do better by considering the land mask - ! from each special landunit raw dataset, and even better by mapping these - ! various notions of land mask onto each other, but that starts to get messy, and - ! relies on the trustworthiness of each raw dataset's land mask... this - ! formulation seems reasonable enough.) + ! Brief summary of the following: But we don't want to overwrite special + ! landunits or crop with ocean where these special landunits extend beyond the + ! PFT data's land fraction. + ! + ! More details: + ! + ! In essence, this is saying that we'll let special landunit area grow into the + ! natveg area before growing into ocean, but we'll have special landunit area + ! grow into ocean before growing into crop or any other special landunit area. + ! (This check of special landunit area is particularly important for glaciers, + ! where we can have floating ice shelves, so we can have a scenario where + ! pctlnd_pft is 0 but we have non-zero glacier cover and we want the final grid + ! cell to be glacier-covered.) (We could possibly do better by considering the + ! land mask from each special landunit raw dataset, and even better by mapping + ! these various notions of land mask onto each other, but that starts to get + ! messy, and relies on the trustworthiness of each raw dataset's land mask... + ! this formulation seems reasonable enough.) ! ! Note that we include pct_crop in the following, but NOT pct_natveg. The ! assumption behind that is that pct_crop is more reliable and/or more important,