From a14758bede89c544c45a63915fe03ed2139bbbd6 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 22 Aug 2024 17:05:38 -0600 Subject: [PATCH] Change logic to account for new allowed / disallowed behavior combos - Weaken the consistency checks between icemask and glc_behavior: We no longer require has_virtual_columns and melt_replaced_by_ice - we now only require that we do NOT have (melt_replaced_by_ice_grc .and. .not. has_virtual_columns_grc). - Prevent users from setting the combination of glacier_region_melt_behavior = "replaced_by_ice" with glacier_region_ice_runoff_behavior = "melted". (While there is nothing fundamentally wrong with this combination, it can result in problematic, non-physical fluxes - particularly, a large positive sensible heat flux during glacial melt in regions where the icesheet is not fully dynamic and two-way-coupled; see https://github.com/ESCOMP/ctsm/issues/423 for details.) - Only update glacier areas and topo values where the glacier region behavior is 'virtual', because that's the only region where we are guaranteed to have all of the elevation classes we need in order to remain in sync. (Note that, for conservation purposes, it's important that we update areas in all regions where we're fully-two-way-coupled to the icesheet and we're computing SMB; this requirement is checked in check_glc2lnd_icemask.) This change is needed now that we no longer require grid cells within the ice mask to have the 'virtual' behavior. - Ensure that glc_dyn_runoff_routing is 0 wherever we're not computing SMB. This change isn't strictly necessary with the current code, because it appears that glc_dyn_runoff_routing is only used within the do_smb filter. However, this change makes the code more robust to future changes. This change is needed now that we no longer require grid cells within the ice mask to have the melt_replaced_by_ice behavior. Also fixes / adds unit tests covering these behavior changes --- src/main/glc2lndMod.F90 | 206 ++++++++++-------- src/main/glcBehaviorMod.F90 | 31 ++- src/main/lnd2glcMod.F90 | 20 +- .../test/glcBehavior_test/test_glcBehavior.pf | 2 +- src/main/test/topo_test/test_topo.pf | 69 ++++-- 5 files changed, 214 insertions(+), 114 deletions(-) diff --git a/src/main/glc2lndMod.F90 b/src/main/glc2lndMod.F90 index 2d0dbb5791..c2e6290300 100644 --- a/src/main/glc2lndMod.F90 +++ b/src/main/glc2lndMod.F90 @@ -35,7 +35,7 @@ module glc2lndMod ! Public data ! ------------------------------------------------------------------------ - ! Where we should do runoff routing that is appropriate for having a dynamic icesheet underneath. + ! Where we should do SMB-related runoff routing that is appropriate for having a dynamic icesheet underneath. real(r8), pointer :: glc_dyn_runoff_routing_grc (:) => null() ! ------------------------------------------------------------------------ @@ -383,32 +383,32 @@ subroutine check_glc2lnd_icemask(this, bounds) if (this%icemask_grc(g) > 0._r8) then - ! Ensure that icemask is a subset of has_virtual_columns. This is needed because - ! we allocated memory based on has_virtual_columns, so it is a problem if the - ! ice sheet tries to expand beyond the area defined by has_virtual_columns. - if (.not. this%glc_behavior%has_virtual_columns_grc(g)) then - write(iulog,'(a)') subname//' ERROR: icemask must be a subset of has_virtual_columns.' - write(iulog,'(a)') 'Ensure that the glacier_region_behavior namelist item is set correctly.' - write(iulog,'(a)') '(It should specify "virtual" for the region corresponding to the GLC domain.)' - write(iulog,'(a)') 'If glacier_region_behavior is set correctly, then you can fix this problem' - write(iulog,'(a)') 'by modifying GLACIER_REGION on the surface dataset.' - write(iulog,'(a)') '(Expand the region that corresponds to the GLC domain' - write(iulog,'(a)') '- i.e., the region specified as "virtual" in glacier_region_behavior.)' - call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg=errMsg(sourcefile, __LINE__)) - end if - - ! Ensure that icemask is a subset of melt_replaced_by_ice. This is needed - ! because we only compute SMB in the region given by melt_replaced_by_ice - ! (according to the logic for building the do_smb filter), and we need SMB - ! everywhere inside the icemask. - if (.not. this%glc_behavior%melt_replaced_by_ice_grc(g)) then - write(iulog,'(a)') subname//' ERROR: icemask must be a subset of melt_replaced_by_ice.' - write(iulog,'(a)') 'Ensure that the glacier_region_melt_behavior namelist item is set correctly.' - write(iulog,'(a)') '(It should specify "replaced_by_ice" for the region corresponding to the GLC domain.)' - write(iulog,'(a)') 'If glacier_region_behavior is set correctly, then you can fix this problem' - write(iulog,'(a)') 'by modifying GLACIER_REGION on the surface dataset.' - write(iulog,'(a)') '(Expand the region that corresponds to the GLC domain' - write(iulog,'(a)') '- i.e., the region specified as "replaced_by_ice" in glacier_region_melt_behavior.)' + ! Ensure that, within the icemask, there are no points that have (non-virtual + ! and compute-SMB). This is important for two reasons: + ! + ! (1) To ensure that, in grid cells where we're producing SMB, we have SMB for + ! all elevation classes, so that the downscaling / vertical interpolation + ! can be done correctly. + ! + ! (2) To avoid conservation issues, we want to ensure that, in grid cells where + ! we're producing SMB and are dynamically coupled to the ice sheet (if 2-way + ! coupling is enabled), glacier areas are remaining in-sync with glc. (Note + ! that has_virtual_columns_grc dictates where we're able to keep glacier + ! areas in sync with glc.) (In principle, I think this one could check + ! icemask_coupled_fluxes rather than icemask; we check icemask because we + ! needed to check icemask for the other reason anyway; this is okay because + ! icemask_coupled_fluxes is a subset of icemask.) + if (this%glc_behavior%melt_replaced_by_ice_grc(g) .and. & + .not. this%glc_behavior%has_virtual_columns_grc(g)) then + write(iulog,'(a)') subname//' ERROR: Within the icemask, there cannot be any points that have' + write(iulog,'(a)') '(non-virtual and compute-SMB).' + write(iulog,'(a)') 'Ensure that GLACIER_REGION on the surface dataset and the namelist items,' + write(iulog,'(a)') 'glacier_region_behavior and glacier_region_melt_behavior are all set correctly:' + write(iulog,'(a)') 'Typically, the region encompassing the active GLC domain should specify' + write(iulog,'(a)') 'glacier_region_behavior="virtual" and glacier_region_melt_behavior="replaced_by_ice".' + write(iulog,'(a)') '(But it is also okay for part of the GLC domain to have' + write(iulog,'(a)') 'glacier_region_melt_behavior="remains_in_place"; this part of the domain can have' + write(iulog,'(a)') 'any setting for glacier_region_behavior.)' call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg=errMsg(sourcefile, __LINE__)) end if @@ -437,10 +437,12 @@ subroutine check_glc2lnd_icemask_coupled_fluxes(this, bounds) do g = bounds%begg, bounds%endg - ! Ensure that icemask_coupled_fluxes is a subset of icemask. Although there - ! currently is no code in CLM that depends on this relationship, it seems helpful - ! to ensure that this intuitive relationship holds, so that code developed in the - ! future can rely on it. + ! Ensure that icemask_coupled_fluxes is a subset of icemask. This is helpful to + ! ensure that the consistency checks that are done on glc behaviors within the + ! icemask (in check_glc2lnd_icemask) also apply within the icemask_coupled_fluxes + ! region. Other than that convenience, there currently is no code in CLM that + ! depends on this relationship, but it seems helpful to ensure that this intuitive + ! relationship holds, so that code developed in the future can rely on it. if (this%icemask_coupled_fluxes_grc(g) > 0._r8 .and. this%icemask_grc(g) == 0._r8) then write(iulog,*) subname//' ERROR: icemask_coupled_fluxes must be a subset of icemask.' call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg=errMsg(sourcefile, __LINE__)) @@ -477,70 +479,73 @@ subroutine update_glc2lnd_dyn_runoff_routing(this, bounds) ! where CISM is running in diagnostic-only mode and therefore is not sending a calving flux - ! we have glc_dyn_runoff_routing = 0, and the snowcap flux goes to the runoff model. ! This is needed to conserve water correctly in the absence of a calving flux. + ! + ! In places where we are not computing SMB, we also have glc_dyn_runoff_routing = 0. + ! Currently glc_dyn_runoff_routing is only used where we're computing SMB, but if it + ! were ever used elsewhere, it seems best to have it set to 0 there: this seems + ! consistent with the fact that we zero out the SMB flux sent to GLC in that region. + ! (However, it's possible that, once we start actually using glc_dyn_runoff_routing + ! for some purpose outside the do_smb filter, we'll discover that this logic should + ! be changed.) do g = bounds%begg, bounds%endg - ! Set glc_dyn_runoff_routing_grc(g) to a value in the range [0,1]. - ! - ! This value gives the grid cell fraction that is deemed to be coupled to the - ! dynamic ice sheet model. For this fraction of the grid cell, snowcap fluxes are - ! sent to the ice sheet model. The remainder of the grid cell sends snowcap fluxes - ! to the runoff model. - ! - ! Note: The coupler (in prep_glc_mod.F90) assumes that the fraction coupled to the - ! dynamic ice sheet model is min(lfrac, Sg_icemask_l), where lfrac is the - ! "frac" component of fraction_lx, and Sg_icemask_l is obtained by mapping - ! Sg_icemask_g from the glc to the land grid. Here, ldomain%frac is - ! equivalent to lfrac, and this%icemask_grc is equivalent to Sg_icemask_l. - ! However, here we use icemask_coupled_fluxes_grc, so that we route all snow - ! capping to runoff in areas where the ice sheet is not generating calving - ! fluxes. In addition, here we need to divide by lfrac, because the coupler - ! multiplies by it later (and, for example, if lfrac = 0.1 and - ! icemask_coupled_fluxes = 1, we want all snow capping to go to the ice - ! sheet model, not to the runoff model). - ! - ! Note: In regions where CLM overlaps the CISM domain, this%icemask_grc(g) typically - ! is nearly equal to ldomain%frac(g). So an alternative would be to simply set - ! glc_dyn_runoff_routing_grc(g) = icemask_grc(g). - ! The reason to cap glc_dyn_runoff_routing at lfrac is to avoid sending the - ! ice sheet model a greater mass of water (in the form of snowcap fluxes) - ! than is allowed to fall on a CLM grid cell that is part ocean. - - ! TODO(wjs, 2017-05-08) Ideally, we wouldn't have this duplication in logic - ! between the coupler and CLM. The best solution would be to have the coupler - ! itself do the partitioning of the snow capping flux between the ice sheet model - ! and the runoff model. A next-best solution would be to have the coupler send a - ! field to CLM telling it what fraction of snow capping should go to the runoff - ! model in each grid cell. - - if (ldomain%frac(g) == 0._r8) then - ! Avoid divide by 0; note that, in this case, the amount going to runoff isn't - ! important for system-wide conservation, so we could really choose anything we - ! want. - this%glc_dyn_runoff_routing_grc(g) = this%icemask_coupled_fluxes_grc(g) - else - this%glc_dyn_runoff_routing_grc(g) = & - min(ldomain%frac(g), this%icemask_coupled_fluxes_grc(g)) / & - ldomain%frac(g) - end if - - if (this%glc_dyn_runoff_routing_grc(g) > 0.0_r8) then - - ! Ensure that glc_dyn_runoff_routing is a subset of melt_replaced_by_ice. This - ! is needed because glacial melt is only sent to the runoff stream in the region - ! given by melt_replaced_by_ice (because the latter is used to create the do_smb - ! filter, and the do_smb filter controls where glacial melt is computed). - if (.not. this%glc_behavior%melt_replaced_by_ice_grc(g)) then - write(iulog,'(a)') subname//' ERROR: icemask_coupled_fluxes must be a subset of melt_replaced_by_ice.' - write(iulog,'(a)') 'Ensure that the glacier_region_melt_behavior namelist item is set correctly.' - write(iulog,'(a)') '(It should specify "replaced_by_ice" for the region corresponding to the GLC domain.)' - write(iulog,'(a)') 'If glacier_region_behavior is set correctly, then you can fix this problem' - write(iulog,'(a)') 'by modifying GLACIER_REGION on the surface dataset.' - write(iulog,'(a)') '(Expand the region that corresponds to the GLC domain' - write(iulog,'(a)') '- i.e., the region specified as "replaced_by_ice" in glacier_region_melt_behavior.)' - call endrun(subgrid_index=g, subgrid_level=subgrid_level_gridcell, msg=errMsg(sourcefile, __LINE__)) + if (this%glc_behavior%melt_replaced_by_ice_grc(g)) then + ! As noted in the comments at the top of this routine, we only set + ! glc_dyn_runoff_routing where we are computing SMB + + ! Set glc_dyn_runoff_routing_grc(g) to a value in the range [0,1]. + ! + ! This value gives the grid cell fraction that is deemed to be coupled to the + ! dynamic ice sheet model. For this fraction of the grid cell, snowcap fluxes are + ! sent to the ice sheet model. The remainder of the grid cell sends snowcap fluxes + ! to the runoff model. + ! + ! Note: The coupler (in prep_glc_mod.F90) assumes that the fraction coupled to the + ! dynamic ice sheet model is min(lfrac, Sg_icemask_l), where lfrac is the + ! "frac" component of fraction_lx, and Sg_icemask_l is obtained by mapping + ! Sg_icemask_g from the glc to the land grid. Here, ldomain%frac is + ! equivalent to lfrac, and this%icemask_grc is equivalent to Sg_icemask_l. + ! However, here we use icemask_coupled_fluxes_grc, so that we route all snow + ! capping to runoff in areas where the ice sheet is not generating calving + ! fluxes. In addition, here we need to divide by lfrac, because the coupler + ! multiplies by it later (and, for example, if lfrac = 0.1 and + ! icemask_coupled_fluxes = 1, we want all snow capping to go to the ice + ! sheet model, not to the runoff model). + ! + ! Note: In regions where CLM overlaps the CISM domain, this%icemask_grc(g) typically + ! is nearly equal to ldomain%frac(g). So an alternative would be to simply set + ! glc_dyn_runoff_routing_grc(g) = icemask_grc(g). + ! The reason to cap glc_dyn_runoff_routing at lfrac is to avoid sending the + ! ice sheet model a greater mass of water (in the form of snowcap fluxes) + ! than is allowed to fall on a CLM grid cell that is part ocean. + + ! TODO(wjs, 2017-05-08) Ideally, we wouldn't have this duplication in logic + ! between the coupler and CLM. The best solution would be to have the coupler + ! itself do the partitioning of the snow capping flux between the ice sheet model + ! and the runoff model. A next-best solution would be to have the coupler send a + ! field to CLM telling it what fraction of snow capping should go to the runoff + ! model in each grid cell. + + if (ldomain%frac(g) == 0._r8) then + ! Avoid divide by 0; note that, in this case, the amount going to runoff isn't + ! important for system-wide conservation, so we could really choose anything we + ! want. + this%glc_dyn_runoff_routing_grc(g) = this%icemask_coupled_fluxes_grc(g) + else + this%glc_dyn_runoff_routing_grc(g) = & + min(ldomain%frac(g), this%icemask_coupled_fluxes_grc(g)) / & + ldomain%frac(g) end if + + else ! .not. this%glc_behavior%melt_replaced_by_ice_grc(g) + ! As noted in the comments at the top of this routine, we set + ! glc_dyn_runoff_routing to 0 where we are not computing SMB. (This assumes that + ! gridcells where we compute SMB are the same as gridcells for which + ! melt_replaced_by_ice is true.) + this%glc_dyn_runoff_routing_grc(g) = 0._r8 end if + end do end subroutine update_glc2lnd_dyn_runoff_routing @@ -578,8 +583,15 @@ subroutine update_glc2lnd_fracs(this, bounds) if (glc_do_dynglacier) then do g = bounds%begg, bounds%endg - ! Values from GLC are only valid within the icemask, so we only update CLM's areas there - if (this%icemask_grc(g) > 0._r8) then + ! Values from GLC are only valid within the icemask, so we only update CLM's + ! areas there. Also, we only update areas where the glacier region behavior is + ! 'virtual', because that's the only region where we are guaranteed to have all + ! of the elevation classes we need in order to remain in sync. (Note that, for + ! conservation purposes, it's important that we update areas in all regions + ! where we're fully-two-way-coupled to the icesheet and we're computing SMB; + ! this requirement is checked in check_glc2lnd_icemask.) (This conditional + ! should be kept consistent with the conditional in update_glc2lnd_topo.) + if (this%icemask_grc(g) > 0._r8 .and. this%glc_behavior%has_virtual_columns_grc(g)) then ! Set total ice landunit area area_ice = sum(this%frac_grc(g, 1:maxpatch_glc)) @@ -623,7 +635,7 @@ subroutine update_glc2lnd_fracs(this, bounds) msg="at least one glc column has non-zero area from cpl but has no slot in memory") end if ! error end if ! area_ice > 0 - end if ! this%icemask_grc(g) > 0 + end if ! this%icemask_grc(g) > 0 .and. this%glc_behavior%has_virtual_columns_grc(g) end do ! g end if ! glc_do_dynglacier @@ -667,8 +679,12 @@ subroutine update_glc2lnd_topo(this, bounds, topo_col, needs_downscaling_col) l = col%landunit(c) g = col%gridcell(c) - ! Values from GLC are only valid within the icemask, so we only update CLM's topo values there - if (this%icemask_grc(g) > 0._r8) then + ! Values from GLC are only valid within the icemask, so we only update CLM's + ! topo values there. Also, consistently with the conditional in + ! update_glc2lnd_fracs, we only update topo values where the glacier region + ! behavior is 'virtual': it could be problematic to update topo values in a + ! grid cell where we're not updating areas. + if (this%icemask_grc(g) > 0._r8 .and. this%glc_behavior%has_virtual_columns_grc(g)) then if (lun%itype(l) == istice) then ice_class = col_itype_to_ice_class(col%itype(c)) else diff --git a/src/main/glcBehaviorMod.F90 b/src/main/glcBehaviorMod.F90 index d2719a6533..315aa88c67 100644 --- a/src/main/glcBehaviorMod.F90 +++ b/src/main/glcBehaviorMod.F90 @@ -72,9 +72,10 @@ module glcBehaviorMod logical, allocatable, public :: allow_multiple_columns_grc(:) ! If melt_replaced_by_ice_grc(g) is true, then any glacier ice melt in gridcell g - ! runs off and is replaced by ice. Note that SMB cannot be computed in gridcell g if - ! melt_replaced_by_ice_grc(g) is false, since we can't compute a sensible negative - ! smb in that case. + ! runs off and is replaced by ice. This flag is also used to determine where we + ! compute SMB: We compute SMB in all grid cells for which melt_replaced_by_ice_grc is + ! true. (SMB cannot be computed in gridcells where melt_replaced_by_ice_grc is false, + ! since we can't compute a sensible negative SMB in that case.) logical, allocatable, public :: melt_replaced_by_ice_grc(:) ! If ice_runoff_melted_grc(g) is true, then ice runoff generated by the @@ -310,6 +311,7 @@ subroutine InitFromInputs(this, begg, endg, & call translate_glacier_region_behavior call translate_glacier_region_melt_behavior call translate_glacier_region_ice_runoff_behavior + call check_behaviors call this%InitAllocate(begg, endg) @@ -484,7 +486,28 @@ subroutine translate_glacier_region_ice_runoff_behavior end do end subroutine translate_glacier_region_ice_runoff_behavior - end subroutine InitFromInputs + subroutine check_behaviors + ! Check the various behaviors for validity / consistency + integer :: i + + do i = min_glacier_region_id, max_glacier_region_id + if (glacier_region_melt_behavior(i) == MELT_BEHAVIOR_REPLACED_BY_ICE .and. & + glacier_region_ice_runoff_behavior(i) == ICE_RUNOFF_BEHAVIOR_MELTED) then + write(iulog,*) ' ERROR: Bad glacier region behavior combination for ID ', i + write(iulog,*) 'You cannot combine glacier_region_melt_behavior = "replaced_by_ice"' + write(iulog,*) 'with glacier_region_ice_runoff_behavior = "melted".' + write(iulog,*) 'While there is nothing fundamentally wrong with this combination,' + write(iulog,*) 'it can result in problematic, non-physical fluxes (particularly,' + write(iulog,*) 'a large positive sensible heat flux during glacial melt in regions' + write(iulog,*) 'where the icesheet is not fully dynamic and two-way-coupled;' + write(iulog,*) 'see https://github.com/ESCOMP/ctsm/issues/423 for details).' + call endrun(msg=' ERROR: Bad glacier region behavior combination '// & + errMsg(sourcefile, __LINE__)) + end if + end do + end subroutine check_behaviors + + end subroutine InitFromInputs !----------------------------------------------------------------------- diff --git a/src/main/lnd2glcMod.F90 b/src/main/lnd2glcMod.F90 index 27fa7639d7..26359ff261 100644 --- a/src/main/lnd2glcMod.F90 +++ b/src/main/lnd2glcMod.F90 @@ -169,10 +169,28 @@ subroutine update_lnd2glc(this, bounds, num_do_smb_c, filter_do_smb_c, & character(len=*), parameter :: subname = 'update_lnd2glc' !------------------------------------------------------------------------------ - ! Initialize to reasonable defaults + ! Initialize to reasonable defaults. These values will be sent for gridcells / + ! columns outside the do_smb filter. + ! NOTE(wjs, 2018-07-03) qice should be 0 outside the do_smb filter to ensure conservation this%qice_grc(bounds%begg : bounds%endg, :) = 0._r8 + + ! NOTE(wjs, 2018-07-03) tsrf can be anything outside the do_smb filter; 0 deg C seems + ! as reasonable as anything (based on input from Bill Lipscomb and Gunter Leguy) this%tsrf_grc(bounds%begg : bounds%endg, :) = tfrz + + ! NOTE(wjs, 2018-07-03) The topo values outside the do_smb filter could matter for + ! gridcells where we compute SMB for some but not all elevation classes (possibly + ! because we haven't even allocated memory for some elevation classes - i.e., if we're + ! not using the 'virtual' behavior in that gridcell). In glc2lndMod, we ensure that + ! this cannot occur for gridcells within the icemask (i.e., within the icemask, we + ! ensure that there are no points that have (non-virtual and compute-SMB)), so this + ! isn't a conservation issue, but it could still be important, e.g., for generating + ! appropriate forcings for a later dlnd-driven T compset. I'm not sure what is "right" + ! here. We've historically used 0 for this, and maybe that's as good as anything, + ! because it may lead to the 0 SMB values being ignored for the sake of vertical + ! interpolation, but I'm not sure about this. But maybe it would be better to use + ! topo at the center of each elevation class? this%topo_grc(bounds%begg : bounds%endg, :) = 0._r8 ! Fill the lnd->glc data on the clm grid diff --git a/src/main/test/glcBehavior_test/test_glcBehavior.pf b/src/main/test/glcBehavior_test/test_glcBehavior.pf index ff104458b1..37271d2f98 100644 --- a/src/main/test/glcBehavior_test/test_glcBehavior.pf +++ b/src/main/test/glcBehavior_test/test_glcBehavior.pf @@ -170,7 +170,7 @@ contains call glc_behavior%InitFromInputs(bounds%begg, bounds%endg, & glacier_region_map = [0], & glacier_region_behavior_str = ['single_at_atm_topo'], & - glacier_region_melt_behavior_str = ['replaced_by_ice'], & + glacier_region_melt_behavior_str = ['remains_in_place'], & glacier_region_ice_runoff_behavior_str = ['melted']) @assertTrue(glc_behavior%ice_runoff_melted_grc(bounds%begg)) diff --git a/src/main/test/topo_test/test_topo.pf b/src/main/test/topo_test/test_topo.pf index 196ce34763..82c3b2cb90 100644 --- a/src/main/test/topo_test/test_topo.pf +++ b/src/main/test/topo_test/test_topo.pf @@ -251,13 +251,33 @@ contains expected_filter = col_filter_empty(bounds) call this%topo%Init(bounds) - ! Need icemask 0, because we can't have single_at_atm_topo inside the icemask - call this%do_UpdateTopo(glc_behavior, icemask_grc = grc_array(0._r8)) + call this%do_UpdateTopo(glc_behavior, icemask_grc = grc_array(1._r8)) filter = this%topo%DownscaleFilterc(bounds) @assertTrue(filter == expected_filter) end subroutine downscaleFilter_afterUpdate_doesNotContain_singleAtAtmTopo + @Test + subroutine downscaleFilter_afterUpdate_contains_vegInsideIcemaskAndVirtual(this) + ! We expect the downscaleFilter to contain vegetated points if they are both (1) + ! inside the icemask, and (2) in the 'virtual' region - because topo is updated in + ! that region. + class(TestTopo), intent(inout) :: this + type(glc_behavior_type) :: glc_behavior + type(filter_col_type) :: filter + type(filter_col_type) :: expected_filter + + call setup_single_veg_patch(pft_type = 1) + glc_behavior = create_glc_behavior_all_virtual() + expected_filter = col_filter_from_index_array(bounds, [bounds%begc]) + + call this%topo%Init(bounds) + call this%do_UpdateTopo(glc_behavior, icemask_grc = grc_array(1._r8)) + filter = this%topo%DownscaleFilterc(bounds) + + @assertTrue(filter == expected_filter) + end subroutine downscaleFilter_afterUpdate_contains_vegInsideIcemaskAndVirtual + @Test subroutine downscaleFilter_afterUpdate_doesNotContain_vegOutsideIcemask(this) class(TestTopo), intent(inout) :: this @@ -266,8 +286,6 @@ contains type(filter_col_type) :: expected_filter call setup_single_veg_patch(pft_type = 1) - ! Use 'virtual' behavior, to make sure that we're not accidentally trying to - ! downscale vegetation over virtual columns. glc_behavior = create_glc_behavior_all_virtual() expected_filter = col_filter_empty(bounds) @@ -279,7 +297,10 @@ contains end subroutine downscaleFilter_afterUpdate_doesNotContain_vegOutsideIcemask @Test - subroutine downscaleFilter_afterUpdate_contains_vegInsideIcemask(this) + subroutine downscaleFilter_afterUpdate_doesNotContain_vegNonVirtual(this) + ! Since topo is only updated in the 'virtual' region, we expect the downscale filter + ! to NOT include vegetated points outside the 'virtual' region, because topo + ! shouldn't be updated for those vegetated points. class(TestTopo), intent(inout) :: this type(glc_behavior_type) :: glc_behavior type(filter_col_type) :: filter @@ -287,14 +308,36 @@ contains call setup_single_veg_patch(pft_type = 1) glc_behavior = create_glc_behavior_all_multiple() - expected_filter = col_filter_from_index_array(bounds, [bounds%begc]) + expected_filter = col_filter_empty(bounds) call this%topo%Init(bounds) call this%do_UpdateTopo(glc_behavior, icemask_grc = grc_array(1._r8)) filter = this%topo%DownscaleFilterc(bounds) @assertTrue(filter == expected_filter) - end subroutine downscaleFilter_afterUpdate_contains_vegInsideIcemask + end subroutine downscaleFilter_afterUpdate_doesNotContain_vegNonVirtual + + @Test + subroutine topo_changes_for_glcmecInsideIcemaskAndVirtual(this) + class(TestTopo), intent(inout) :: this + type(glc_behavior_type) :: glc_behavior + real(r8), parameter :: topo_orig = 7._r8 + real(r8), parameter :: atm_topo = 23._r8 + + ! our column should get set to this: + real(r8), parameter :: glc_topo = 27._r8 + + call setup_single_ice_column(elev_class = 1) + glc_behavior = create_glc_behavior_all_virtual() + topo_glc_mec(:,:) = topo_orig + + call this%topo%Init(bounds) + call this%do_UpdateTopo(glc_behavior, icemask_grc = grc_array(1._r8), & + atm_topo_grc = grc_array(atm_topo), & + glc_topo = glc_topo) + + @assertEqual(glc_topo, this%topo%topo_col(bounds%begc)) + end subroutine topo_changes_for_glcmecInsideIcemaskAndVirtual @Test subroutine topo_noChange_for_glcmecOutsideIcemask(this) @@ -319,17 +362,17 @@ contains end subroutine topo_noChange_for_glcmecOutsideIcemask @Test - subroutine topo_changes_for_glcmecInsideIcemask(this) + subroutine topo_noChange_for_glcmecNonVirtual(this) class(TestTopo), intent(inout) :: this type(glc_behavior_type) :: glc_behavior real(r8), parameter :: topo_orig = 7._r8 - real(r8), parameter :: atm_topo = 23._r8 - ! our column should get set to this: + ! our column should NOT get set to either of these: + real(r8), parameter :: atm_topo = 23._r8 real(r8), parameter :: glc_topo = 27._r8 call setup_single_ice_column(elev_class = 1) - glc_behavior = create_glc_behavior_all_virtual() + glc_behavior = create_glc_behavior_all_multiple() topo_glc_mec(:,:) = topo_orig call this%topo%Init(bounds) @@ -337,8 +380,8 @@ contains atm_topo_grc = grc_array(atm_topo), & glc_topo = glc_topo) - @assertEqual(glc_topo, this%topo%topo_col(bounds%begc)) - end subroutine topo_changes_for_glcmecInsideIcemask + @assertEqual(topo_orig, this%topo%topo_col(bounds%begc)) + end subroutine topo_noChange_for_glcmecNonVirtual @Test subroutine topo_changes_for_singleAtAtmTopo(this)