diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 820975655d..210e4bf243 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -3012,7 +3012,7 @@ snow melt of Brock et al. (2006) group="clm_initinterp_inparm" valid_values="" > If FALSE (which is the default): If an output type cannot be found in the input for initInterp, code aborts -If TRUE: If an output type cannot be found in the input, fill with closest natural veg column +If TRUE: If a non-urban output type cannot be found in the input, fill with closest natural veg column (using bare soil for patch-level variables) NOTE: Natural vegetation and crop landunits always behave as if this were true. e.g., if @@ -3021,6 +3021,14 @@ always fill with the closest natural veg patch / column, regardless of the value flag. So interpolation from non-crop to crop cases can be done without setting this flag. + +If FALSE (which is the default): If an urban output type cannot be found in the input for initInterp, +code aborts +If TRUE: If an urban output type cannot be found in the input, fill with closest urban high density +(HD) landunit + + diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 94c6dd4572..2f5738f4dc 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -173,6 +173,16 @@ + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/include_user_mods new file mode 100644 index 0000000000..fe0e18cf88 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/include_user_mods @@ -0,0 +1 @@ +../default diff --git a/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/user_nl_clm new file mode 100644 index 0000000000..499b6026ea --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/f09_FillMissingW_Urban/user_nl_clm @@ -0,0 +1,6 @@ +! NOTE: Using an initial file that does NOT have TBD on it and 5.4 landuse timeseries dataset that has TBD on it +fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_esmf/ctsm5.4.0/surfdata_0.9x1.25_hist_1850_78pfts_c250428.nc' +flanduse_timeseries = '$DIN_LOC_ROOT/lnd/clm2/surfdata_esmf/ctsm5.4.0/landuse.timeseries_0.9x1.25_hist_1850-2023_78pfts_c250428.nc' +finidat = '$DIN_LOC_ROOT/lnd/clm2/initdata_esmf/ctsm5.4/ctsm53041_54surfdata_snowTherm_100_pSASU.clm2.r.0161-01-01-00000.nc' +init_interp_fill_missing_urban_with_HD = .true. +use_init_interp = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/run_self_tests/shell_commands b/cime_config/testdefs/testmods_dirs/clm/run_self_tests/shell_commands index e9a2f934d5..7762f69e36 100755 --- a/cime_config/testdefs/testmods_dirs/clm/run_self_tests/shell_commands +++ b/cime_config/testdefs/testmods_dirs/clm/run_self_tests/shell_commands @@ -5,11 +5,11 @@ ./xmlchange ROF_NCPL='$ATM_NCPL' # Turn MEGAN off to run faster -./xmlchange CLM_BLDNML_OPTS='--no-megan' -append +./xmlchange CLM_BLDNML_OPTS='--no-megan' --append # Use fast structure and NWP configuration for speed ./xmlchange CLM_STRUCTURE="fast" ./xmlchange CLM_CONFIGURATION="nwp" # Turn cpl history off -./xmlchagne HIST_OPTION="never" \ No newline at end of file +./xmlchange HIST_OPTION="never" \ No newline at end of file diff --git a/src/init_interp/initInterp.F90 b/src/init_interp/initInterp.F90 index c8f80e16ae..3ccdcd9b58 100644 --- a/src/init_interp/initInterp.F90 +++ b/src/init_interp/initInterp.F90 @@ -75,6 +75,9 @@ module initInterpMod ! patch-level variables) logical :: init_interp_fill_missing_with_natveg + ! If true, fill missing urban landunit type with closest urban high density (HD) landunit + logical :: init_interp_fill_missing_urban_with_HD + character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -106,11 +109,13 @@ subroutine initInterp_readnl(NLFilename) !----------------------------------------------------------------------- namelist /clm_initinterp_inparm/ & - init_interp_method, init_interp_fill_missing_with_natveg + init_interp_method, init_interp_fill_missing_with_natveg, & + init_interp_fill_missing_urban_with_HD ! Initialize options to default values, in case they are not specified in the namelist init_interp_method = ' ' init_interp_fill_missing_with_natveg = .false. + init_interp_fill_missing_urban_with_HD = .false. if (masterproc) then unitn = getavu() @@ -130,6 +135,7 @@ subroutine initInterp_readnl(NLFilename) call shr_mpi_bcast (init_interp_method, mpicom) call shr_mpi_bcast (init_interp_fill_missing_with_natveg, mpicom) + call shr_mpi_bcast (init_interp_fill_missing_urban_with_HD, mpicom) if (masterproc) then write(iulog,*) ' ' @@ -287,12 +293,36 @@ subroutine initInterp (filei, fileo, bounds, glc_behavior) status = pio_get_att(ncidi, pio_global, & 'icol_vegetated_or_bare_soil', & subgrid_special_indices%icol_vegetated_or_bare_soil) + status = pio_get_att(ncidi, pio_global, & + 'icol_urban_roof', & + subgrid_special_indices%icol_urban_roof) + status = pio_get_att(ncidi, pio_global, & + 'icol_urban_sunwall', & + subgrid_special_indices%icol_urban_sunwall) + status = pio_get_att(ncidi, pio_global, & + 'icol_urban_shadewall', & + subgrid_special_indices%icol_urban_shadewall) + status = pio_get_att(ncidi, pio_global, & + 'icol_urban_impervious_road', & + subgrid_special_indices%icol_urban_impervious_road) + status = pio_get_att(ncidi, pio_global, & + 'icol_urban_pervious_road', & + subgrid_special_indices%icol_urban_pervious_road) status = pio_get_att(ncidi, pio_global, & 'ilun_vegetated_or_bare_soil', & subgrid_special_indices%ilun_vegetated_or_bare_soil) status = pio_get_att(ncidi, pio_global, & 'ilun_crop', & subgrid_special_indices%ilun_crop) + status = pio_get_att(ncidi, pio_global, & + 'ilun_urban_tbd', & + subgrid_special_indices%ilun_urban_TBD) + status = pio_get_att(ncidi, pio_global, & + 'ilun_urban_hd', & + subgrid_special_indices%ilun_urban_HD) + status = pio_get_att(ncidi, pio_global, & + 'ilun_urban_md', & + subgrid_special_indices%ilun_urban_MD) ! BACKWARDS_COMPATIBILITY(wjs, 2021-04-16) ilun_landice_multiple_elevation_classes has ! been renamed to ilun_landice. For now we need to handle both possibilities for the @@ -321,10 +351,26 @@ subroutine initInterp (filei, fileo, bounds, glc_behavior) subgrid_special_indices%ipft_not_vegetated write(iulog,*)'icol_vegetated_or_bare_soil = ' , & subgrid_special_indices%icol_vegetated_or_bare_soil + write(iulog,*)'icol_urban_roof = ' , & + subgrid_special_indices%icol_urban_roof + write(iulog,*)'icol_urban_sunwall = ' , & + subgrid_special_indices%icol_urban_sunwall + write(iulog,*)'icol_urban_shadewall = ' , & + subgrid_special_indices%icol_urban_shadewall + write(iulog,*)'icol_urban_impervious_road = ' , & + subgrid_special_indices%icol_urban_impervious_road + write(iulog,*)'icol_urban_pervious_road = ' , & + subgrid_special_indices%icol_urban_pervious_road write(iulog,*)'ilun_vegetated_or_bare_soil = ' , & subgrid_special_indices%ilun_vegetated_or_bare_soil write(iulog,*)'ilun_crop = ' , & subgrid_special_indices%ilun_crop + write(iulog,*)'ilun_urban_tbd = ' , & + subgrid_special_indices%ilun_urban_TBD + write(iulog,*)'ilun_urban_hd = ' , & + subgrid_special_indices%ilun_urban_HD + write(iulog,*)'ilun_urban_md = ' , & + subgrid_special_indices%ilun_urban_MD write(iulog,*)'ilun_landice = ' , & subgrid_special_indices%ilun_landice write(iulog,*)'create_glacier_mec_landunits = ', & @@ -839,6 +885,7 @@ subroutine findMinDist( dimname, begi, endi, bego, endo, ncidi, ncido, & glc_behavior=glc_behavior, & glc_elevclasses_same=glc_elevclasses_same, & fill_missing_with_natveg=init_interp_fill_missing_with_natveg, & + fill_missing_urban_with_HD=init_interp_fill_missing_urban_with_HD, & mindist_index=minindx) case (interp_method_finidat_areas) if (masterproc) then diff --git a/src/init_interp/initInterpMindist.F90 b/src/init_interp/initInterpMindist.F90 index f6853b1cd3..9fdc9f81dd 100644 --- a/src/init_interp/initInterpMindist.F90 +++ b/src/init_interp/initInterpMindist.F90 @@ -32,11 +32,20 @@ module initInterpMindist type, public :: subgrid_special_indices_type integer :: ipft_not_vegetated integer :: icol_vegetated_or_bare_soil + integer :: icol_urban_roof + integer :: icol_urban_sunwall + integer :: icol_urban_shadewall + integer :: icol_urban_impervious_road + integer :: icol_urban_pervious_road integer :: ilun_vegetated_or_bare_soil integer :: ilun_crop integer :: ilun_landice + integer :: ilun_urban_TBD + integer :: ilun_urban_HD + integer :: ilun_urban_MD contains procedure :: is_vegetated_landunit ! returns true if the given landunit type is natural veg or crop + procedure :: is_urban_landunit ! returns true if the given landunit type is urban end type subgrid_special_indices_type type, public :: subgrid_type @@ -58,8 +67,10 @@ module initInterpMindist private :: set_glc_must_be_same_type private :: set_ice_adjustable_type private :: do_fill_missing_with_natveg + private :: do_fill_missing_urban_with_HD private :: is_sametype private :: is_baresoil + private :: is_urban_HD character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -147,7 +158,7 @@ end subroutine destroy_subgrid_type subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgrido, & subgrid_special_indices, glc_behavior, glc_elevclasses_same, & - fill_missing_with_natveg, mindist_index) + fill_missing_with_natveg, fill_missing_urban_with_HD, mindist_index) ! -------------------------------------------------------------------- ! arguments @@ -165,7 +176,7 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr logical , intent(in) :: glc_elevclasses_same ! If false: if an output type cannot be found in the input, code aborts - ! If true: if an output type cannot be found in the input, fill with closest natural + ! If true: if a non-urban output type cannot be found in the input, fill with closest natural ! veg column (using bare soil for patch-level variables) ! ! NOTE: always treated as true for natural veg and crop landunits/columns/patches in @@ -173,6 +184,11 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr ! use the closest natural veg column, regardless of the value of this flag. logical , intent(in) :: fill_missing_with_natveg + + ! If false: if an urban output type cannot be found in the input, code aborts + ! If true: if an urban output type cannot be found in the input, fill with closest urban HD + logical , intent(in) :: fill_missing_urban_with_HD + integer , intent(out) :: mindist_index(bego:endo) ! ! local variables @@ -187,6 +203,8 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr ! considered the same type. This is only valid for glc points, and is only valid ! for subgrid name = 'pft' or 'column'. logical :: glc_must_be_same_type_o(bego:endo) + + character(len=*), parameter :: subname = 'set_mindist' ! -------------------------------------------------------------------- if (associated(subgridi%topoglc) .and. associated(subgrido%topoglc)) then @@ -221,7 +239,8 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr subgridi = subgridi, subgrido = subgrido, & subgrid_special_indices = subgrid_special_indices, & glc_must_be_same_type = glc_must_be_same_type_o(no), & - veg_patch_just_considers_ptype = .true.)) then + veg_patch_just_considers_ptype = .true., & + do_fill_missing_urban_with_HD = .false.)) then dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * & 0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni)) @@ -260,7 +279,11 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr end if end do - ! If output type is not contained in input dataset, then use closest bare soil, + ! Note that do_fill_missing_with_natveg below will return .false. for pfts and columnns associated + ! with urban landunits so that the fill missing with bare soil will be implemented only for + ! non-urban types (pfts, columns, landunits, gridcells). + + ! If non-urban output type is not contained in input dataset, then use closest bare soil, ! if this point is one for which we fill missing with natveg. if ( distmin == spval .and. & do_fill_missing_with_natveg( & @@ -279,6 +302,50 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr end if end if end do + + ! If urban output type is not contained in input dataset, then use closest urban HD, + ! if this point is one for which we fill missing urban with urban HD. + else if (distmin == spval & + .and. do_fill_missing_urban_with_HD( & + fill_missing_urban_with_HD, no, subgrido, subgrid_special_indices)) then + do ni = begi, endi + if (activei(ni)) then + ! We need to call is_sametype for pfts and columns here to make sure that each + ! urban input pft and column type matches the output pft and column type. We don't + ! want to call it for landunits because they intentionally won't be the same type + ! (since we are filling missing urban landunits with HD) + if (subgrido%name .eq. 'landunit') then + if ( is_urban_HD(ni, subgridi, subgrid_special_indices)) then + dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re + dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * & + 0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni)) + dist = dx*dx + dy*dy + if ( dist < distmin )then + distmin = dist + nmin = ni + end if + end if + else + if (is_sametype(ni = ni, no = no, & + subgridi = subgridi, subgrido = subgrido, & + subgrid_special_indices = subgrid_special_indices, & + glc_must_be_same_type = glc_must_be_same_type_o(no), & + veg_patch_just_considers_ptype = .false., & + do_fill_missing_urban_with_HD = .true.)) then + if ( is_urban_HD(ni, subgridi, subgrid_special_indices)) then + dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re + dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * & + 0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni)) + dist = dx*dx + dy*dy + if ( dist < distmin )then + distmin = dist + nmin = ni + end if + end if + end if + end if + end if + end do end if ! Error conditions @@ -287,13 +354,29 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr &Cannot find any input points matching output point:' call subgrido%print_point(no, iulog) write(iulog,*) ' ' - write(iulog,*) 'Consider rerunning with the following in user_nl_clm:' + write(iulog,*) 'If this is an urban type' + write(iulog,*) '(ltype = ', subgrid_special_indices%ilun_urban_TBD, & + ',', subgrid_special_indices%ilun_urban_HD, & + ', or', subgrid_special_indices%ilun_urban_MD, ')' + write(iulog,*) 'then consider rerunning with the following in user_nl_clm:' + write(iulog,*) 'init_interp_fill_missing_urban_with_HD = .true.' + write(iulog,*) 'However, note that this will fill all urban missing types in the output' + write(iulog,*) 'with the closest urban high density (HD) type in the input' + write(iulog,*) 'So, you should consider whether that is what you want.' + write(iulog,*) ' ' + write(iulog,*) 'If this is a non-urban type' + write(iulog,*) '(ltype \= ',subgrid_special_indices%ilun_urban_TBD, & + ',', subgrid_special_indices%ilun_urban_HD, & + ', or', subgrid_special_indices%ilun_urban_MD, ')' + write(iulog,*) 'consider rerunning with the following in user_nl_clm:' write(iulog,*) 'init_interp_fill_missing_with_natveg = .true.' - write(iulog,*) 'However, note that this will fill all missing types in the output' + write(iulog,*) 'However, note that this will fill all non-urban missing types in the output' write(iulog,*) 'with the closest natural veg column in the input' write(iulog,*) '(using bare soil for patch-level variables).' write(iulog,*) 'So, you should consider whether that is what you want.' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(iulog,*) errMsg(sourcefile, __LINE__) + call endrun(msg=subname// & + ' ERROR: Cannot find any input points matching output point') end if mindist_index(no) = nmin @@ -378,7 +461,8 @@ subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, subgridi = subgridi, subgrido = subgrido, & subgrid_special_indices = subgrid_special_indices, & glc_must_be_same_type = glc_must_be_same_type_o(no), & - veg_patch_just_considers_ptype = .false.) + veg_patch_just_considers_ptype = .false., & + do_fill_missing_urban_with_HD = .false.) if (ni_sametype) then if (found) then write(iulog,*) subname// & @@ -555,7 +639,7 @@ function do_fill_missing_with_natveg(fill_missing_with_natveg, & no, subgrido, subgrid_special_indices) ! ! !DESCRIPTION: - ! Returns true if the given output point, if missing, should be filled with the + ! Returns true if the given non-urban output point, if missing, should be filled with the ! closest natural veg point. ! ! !ARGUMENTS: @@ -576,8 +660,8 @@ function do_fill_missing_with_natveg(fill_missing_with_natveg, & if (subgrido%name == 'gridcell') then ! It makes no sense to try to fill missing with natveg for gridcell-level values do_fill_missing_with_natveg = .false. - else if (fill_missing_with_natveg) then - ! User has asked for all missing points to be filled with natveg + else if (fill_missing_with_natveg .and. .not. subgrid_special_indices%is_urban_landunit(subgrido%ltype(no))) then + ! User has asked for all non-urban missing points to be filled with natveg do_fill_missing_with_natveg = .true. else if (subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then ! Even if user hasn't asked for it, we fill missing vegetated points (natural veg @@ -591,11 +675,46 @@ function do_fill_missing_with_natveg(fill_missing_with_natveg, & end function do_fill_missing_with_natveg + !----------------------------------------------------------------------- + function do_fill_missing_urban_with_HD(fill_missing_urban_with_HD, & + no, subgrido, subgrid_special_indices) + ! + ! !DESCRIPTION: + ! Returns true if the given urban output point, if missing, should be filled with the + ! closest urban HD point. + ! + ! !ARGUMENTS: + logical :: do_fill_missing_urban_with_HD ! function result + + ! whether we should fill ALL missing points with urban HD + logical, intent(in) :: fill_missing_urban_with_HD + + integer , intent(in) :: no + type(subgrid_type), intent(in) :: subgrido + type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'do_fill_missing_urban_with_HD' + !----------------------------------------------------------------------- + + if (subgrido%name == 'gridcell') then + ! It makes no sense to try to fill missing with urban HD for gridcell-level values + do_fill_missing_urban_with_HD = .false. + else if (fill_missing_urban_with_HD) then + ! User has asked for all missing urban points to be filled with urban HD + do_fill_missing_urban_with_HD = .true. + else + do_fill_missing_urban_with_HD = .false. + end if + + end function do_fill_missing_urban_with_HD !======================================================================= logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indices, & - glc_must_be_same_type, veg_patch_just_considers_ptype) + glc_must_be_same_type, veg_patch_just_considers_ptype, & + do_fill_missing_urban_with_HD) ! -------------------------------------------------------------------- ! arguments @@ -620,6 +739,12 @@ logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indice ! If false, then they need to have the same column and landunit types, too (as is the ! general case). logical, intent(in) :: veg_patch_just_considers_ptype + + ! If True, we allow for landunits to be different when checking if pft and column are + ! the same type, to allow for HD fill of missing urban output points. + logical, intent(in) :: do_fill_missing_urban_with_HD + + ! For urban columns/patches ! -------------------------------------------------------------------- is_sametype = .false. @@ -644,6 +769,10 @@ logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indice subgridi%ptype(ni) == subgrido%ptype(no)) then is_sametype = .true. end if + else if (subgridi%ptype(ni) == subgrido%ptype(no) .and. & + subgridi%ctype(ni) == subgrido%ctype(no) .and. & + do_fill_missing_urban_with_HD) then + is_sametype = .true. else if (subgridi%ptype(ni) == subgrido%ptype(no) .and. & subgridi%ctype(ni) == subgrido%ctype(no) .and. & subgridi%ltype(ni) == subgrido%ltype(no)) then @@ -654,6 +783,9 @@ logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indice subgridi%ltype(ni) == subgrid_special_indices%ilun_landice .and. & subgrido%ltype(no) == subgrid_special_indices%ilun_landice ) then is_sametype = .true. + else if (subgridi%ctype(ni) == subgrido%ctype(no) .and. & + do_fill_missing_urban_with_HD) then + is_sametype = .true. else if (subgridi%ctype(ni) == subgrido%ctype(no) .and. & subgridi%ltype(ni) == subgrido%ltype(no)) then is_sametype = .true. @@ -712,6 +844,31 @@ logical function is_baresoil (n, subgrid, subgrid_special_indices) end function is_baresoil + !----------------------------------------------------------------------- + logical function is_urban_HD (n, subgrid, subgrid_special_indices) + + ! -------------------------------------------------------------------- + ! arguments + integer , intent(in) :: n + type(subgrid_type), intent(in) :: subgrid + type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices + ! -------------------------------------------------------------------- + + is_urban_HD = .false. + + if (subgrid%name == 'pft' .or. subgrid%name == 'column' .or. subgrid%name == 'landunit') then + if (subgrid%ltype(n) == subgrid_special_indices%ilun_urban_HD) then + is_urban_HD = .true. + end if + else + if (masterproc) then + write(iulog,*)'ERROR interpinic: is_urban_HD subgrid type ',subgrid%name,' not supported' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + end function is_urban_HD + !----------------------------------------------------------------------- function is_vegetated_landunit(this, ltype) ! @@ -739,5 +896,30 @@ function is_vegetated_landunit(this, ltype) end function is_vegetated_landunit + function is_urban_landunit(this, ltype) + ! + ! !DESCRIPTION: + ! Returns true if the given landunit type is urban + ! + ! !USES: + ! + ! !ARGUMENTS: + logical :: is_urban_landunit ! function result + class(subgrid_special_indices_type), intent(in) :: this + integer, intent(in) :: ltype ! landunit type of interest + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'is_urban_landunit' + !----------------------------------------------------------------------- + + if (ltype == this%ilun_urban_TBD .or. ltype == this%ilun_urban_HD & + .or. ltype == this%ilun_urban_MD) then + is_urban_landunit = .true. + else + is_urban_landunit = .false. + end if + + end function is_urban_landunit end module initInterpMindist diff --git a/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf b/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf index 04f09cb55d..7d277566af 100644 --- a/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf +++ b/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf @@ -19,12 +19,20 @@ module initInterpMindistTestUtils subgrid_special_indices_type( & ipft_not_vegetated = 0, & icol_vegetated_or_bare_soil = 10, & + icol_urban_roof = 71, & + icol_urban_sunwall = 72, & + icol_urban_shadewall = 73, & + icol_urban_impervious_road = 74, & + icol_urban_pervious_road = 75, & ilun_vegetated_or_bare_soil = 3, & ilun_crop = 4, & - ilun_landice = 5) + ilun_landice = 5, & + ilun_urban_TBD = 7, & + ilun_urban_HD = 8, & + ilun_urban_MD = 9) - ! value we can use for a special landunit; note that this just needs to differ from - ! ilun_vegetated_or_bare_soil and from ilun_crop + ! value we can use for a special landunit; note that this needs to differ from + ! ilun_vegetated_or_bare_soil, ilun_crop, ilun_urban_TBD, ilun_urban_HD, ilun_urban_MD integer, parameter, public :: ilun_special = 6 contains diff --git a/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf b/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf index 06ce20d7de..7a2c51456d 100644 --- a/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf +++ b/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf @@ -10,6 +10,7 @@ module test_set_mindist use clm_varcon , only: spval use unittestSimpleSubgridSetupsMod use unittestSubgridMod + use unittestUtils, only : endrun_msg use glcBehaviorMod, only: glc_behavior_type implicit none @@ -41,7 +42,7 @@ contains end subroutine tearDown subroutine wrap_set_mindist(subgridi, subgrido, mindist_index, activei, activeo, & - glc_behavior, glc_elevclasses_same, fill_missing_with_natveg) + glc_behavior, glc_elevclasses_same, fill_missing_with_natveg, fill_missing_urban_with_HD) ! Wrap the call to set_mindist. ! ! If activei / activeo are not provided, they are assumed to be .true. for all points. @@ -52,6 +53,7 @@ contains ! If glc_elevclasses_same is not present, it is assumed to be true. ! ! If fill_missing_with_natveg is not provided, it is assumed to be false + ! If fill_missing_urban_with_HD is not provided, it is assumed to be false ! Arguments: type(subgrid_type), intent(in) :: subgridi @@ -62,6 +64,7 @@ contains type(glc_behavior_type), intent(in), optional :: glc_behavior logical, intent(in), optional :: glc_elevclasses_same logical, intent(in), optional :: fill_missing_with_natveg + logical, intent(in), optional :: fill_missing_urban_with_HD ! Local variables: integer :: npts_i, npts_o @@ -71,6 +74,7 @@ contains type(glc_behavior_type) :: l_glc_behavior logical :: l_glc_elevclasses_same logical :: l_fill_missing_with_natveg + logical :: l_fill_missing_urban_with_HD !----------------------------------------------------------------------- @@ -115,12 +119,19 @@ contains l_fill_missing_with_natveg = .false. end if + if (present(fill_missing_urban_with_HD)) then + l_fill_missing_urban_with_HD = fill_missing_urban_with_HD + else + l_fill_missing_urban_with_HD = .false. + end if + call set_mindist(begi = 1, endi = npts_i, bego = bego, endo = endo, & activei = l_activei, activeo = l_activeo, subgridi = subgridi, subgrido = subgrido, & subgrid_special_indices = subgrid_special_indices, & glc_behavior = l_glc_behavior, & glc_elevclasses_same = l_glc_elevclasses_same, & fill_missing_with_natveg = l_fill_missing_with_natveg, & + fill_missing_urban_with_HD = l_fill_missing_urban_with_HD, & mindist_index = mindist_index) end subroutine wrap_set_mindist @@ -724,6 +735,186 @@ contains end associate end subroutine newveg_usesBaresoil + @Test + subroutine TBDurban_usesHDurban(this) + ! If there's a new urban TBD type, this should take inputs from the closest + ! HD type, if fill_missing_urban_with_HD = .true and fill_missing_with_natveg = .false. + ! + class(TestSetMindist), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: i + integer :: mindist_index(1) + + associate( & + icol_urban_roof => subgrid_special_indices%icol_urban_roof, & + icol_urban_sunwall => subgrid_special_indices%icol_urban_sunwall, & + icol_urban_shadewall => subgrid_special_indices%icol_urban_shadewall, & + icol_urban_impervious_road => subgrid_special_indices%icol_urban_impervious_road, & + icol_urban_pervious_road => subgrid_special_indices%icol_urban_pervious_road, & + ilun_urban_TBD => subgrid_special_indices%ilun_urban_TBD, & + ilun_urban_HD => subgrid_special_indices%ilun_urban_HD, & + ilun_urban_MD => subgrid_special_indices%ilun_urban_MD & + ) + + call setup_landunit_ncols(ltype=ilun_urban_TBD, & + ctypes=[icol_urban_roof,icol_urban_sunwall,icol_urban_shadewall, & + icol_urban_impervious_road,icol_urban_pervious_road], & + cweights=[0.6_r8,0.1_r8,0.1_r8,0.1_r8,0.1_r8], & + ptype=0) + + call create_subgrid_info( & + subgrid_info = subgrido, & + npts = 1, & + beg = 1, & + name = 'landunit', & + ltype = [ilun_urban_TBD], & + lat = [my_lat], & + lon = [my_lon]) + + ! Input points differ in landunit type + call create_subgrid_info( & + subgrid_info = subgridi, & + npts = 2, & + name = 'landunit', & + ltype = [ilun_urban_MD, ilun_urban_HD], & + lat = [(my_lat, i=1,2)], & + lon = [(my_lon, i=1,2)]) + + call wrap_set_mindist(subgridi, subgrido, mindist_index, & + fill_missing_urban_with_HD = .true., & + fill_missing_with_natveg = .false.) + + ! Note that the mindist_index should return the second index of the + ! ltype array (2), not the actual value of ilun_urban_HD + @assertEqual(2, mindist_index(1)) + + end associate + end subroutine TBDurban_usesHDurban + + @Test + subroutine TBDurban_usesHDurban_aborts(this) + ! If there's a new urban TBD type, this should take inputs from the closest + ! HD type. This test will abort correctly if fill_missing_urban_with_HD = .false. + ! + class(TestSetMindist), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: i + integer :: mindist_index(1) + character(len=:), allocatable :: expected_msg + + associate( & + icol_urban_roof => subgrid_special_indices%icol_urban_roof, & + icol_urban_sunwall => subgrid_special_indices%icol_urban_sunwall, & + icol_urban_shadewall => subgrid_special_indices%icol_urban_shadewall, & + icol_urban_impervious_road => subgrid_special_indices%icol_urban_impervious_road, & + icol_urban_pervious_road => subgrid_special_indices%icol_urban_pervious_road, & + ilun_urban_TBD => subgrid_special_indices%ilun_urban_TBD, & + ilun_urban_HD => subgrid_special_indices%ilun_urban_HD, & + ilun_urban_MD => subgrid_special_indices%ilun_urban_MD & + ) + + call setup_landunit_ncols(ltype=ilun_urban_TBD, & + ctypes=[icol_urban_roof,icol_urban_sunwall,icol_urban_shadewall, & + icol_urban_impervious_road,icol_urban_pervious_road], & + cweights=[0.6_r8,0.1_r8,0.1_r8,0.1_r8,0.1_r8], & + ptype=0) + + call create_subgrid_info( & + subgrid_info = subgrido, & + npts = 1, & + beg = 1, & + name = 'landunit', & + ltype = [ilun_urban_TBD], & + lat = [my_lat], & + lon = [my_lon]) + + ! Input points differ in landunit type + call create_subgrid_info( & + subgrid_info = subgridi, & + npts = 2, & + name = 'landunit', & + ltype = [ilun_urban_MD, ilun_urban_HD], & + lat = [(my_lat, i=1,2)], & + lon = [(my_lon, i=1,2)]) + + call wrap_set_mindist(subgridi, subgrido, mindist_index, & + fill_missing_urban_with_HD = .false.) + + expected_msg = endrun_msg( & + 'set_mindist ERROR: Cannot find any input points matching output point') + @assertExceptionRaised(expected_msg) + + end associate + end subroutine TBDurban_usesHDurban_aborts + + @Test + subroutine urbanlandunits_NotFilled_with_natveg_aborts(this) + ! With fill_missing_urban_with_HD = .false. and fill_missing_with_natveg = .true., + ! urban landunit should not be filled with natveg, and an error in set_mindist will be + ! thrown, and this test should pass. + ! + class(TestSetMindist), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: i + integer :: mindist_index(1) + character(len=:), allocatable :: expected_msg + + associate( & + ipft_bare => subgrid_special_indices%ipft_not_vegetated, & + icol_urban_roof => subgrid_special_indices%icol_urban_roof, & + icol_urban_sunwall => subgrid_special_indices%icol_urban_sunwall, & + icol_urban_shadewall => subgrid_special_indices%icol_urban_shadewall, & + icol_urban_impervious_road => subgrid_special_indices%icol_urban_impervious_road, & + icol_urban_pervious_road => subgrid_special_indices%icol_urban_pervious_road, & + icol_natveg => subgrid_special_indices%icol_vegetated_or_bare_soil, & + ilun_natveg => subgrid_special_indices%ilun_vegetated_or_bare_soil, & + ilun_urban_TBD => subgrid_special_indices%ilun_urban_TBD & + ) + + call setup_landunit_ncols(ltype=ilun_urban_TBD, & + ctypes=[icol_urban_roof,icol_urban_sunwall,icol_urban_shadewall, & + icol_urban_impervious_road,icol_urban_pervious_road], & + cweights=[0.6_r8,0.1_r8,0.1_r8,0.1_r8,0.1_r8], & + ptype=0) + + call create_subgrid_info( & + subgrid_info = subgrido, & + npts = 1, & + beg = 1, & + name = 'pft', & + ptype = [0], & + ctype = [icol_urban_roof], & + ltype = [ilun_urban_TBD], & + lat = [my_lat], & + lon = [my_lon]) + + call create_subgrid_info( & + subgrid_info = subgridi, & + npts = 1, & + name = 'pft', & + ptype = [ipft_bare], & + ctype = [icol_natveg], & + ltype = [ilun_natveg], & + lat = [my_lat], & + lon = [my_lon]) + + call wrap_set_mindist(subgridi, subgrido, mindist_index, & + fill_missing_urban_with_HD = .false., & + fill_missing_with_natveg = .true.) + + expected_msg = endrun_msg( & + 'set_mindist ERROR: Cannot find any input points matching output point') + @assertExceptionRaised(expected_msg) + + end associate + end subroutine urbanlandunits_NotFilled_with_natveg_aborts + @Test subroutine baresoil_ignoresSpecialLandunits(this) ! This test ensures that, when finding a match for a bare soil patch, we ignore