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)