module glcBehaviorMod !----------------------------------------------------------------------- ! !DESCRIPTION: ! Determines a number of aspects of the behavior of glacier_mec classes in each grid ! cell. ! ! !USES: #include "shr_assert.h" use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : iulog use landunit_varcon, only : istice_mec use clm_instur , only : wt_lunit, wt_glc_mec use decompMod , only : bounds_type use filterColMod , only : filter_col_type use ColumnType , only : col use PatchType , only : patch ! !PUBLIC TYPES: implicit none private save type, public :: glc_behavior_type private ! ------------------------------------------------------------------------ ! Public data ! ------------------------------------------------------------------------ ! If has_virtual_columns_grc(g) is true, then grid cell g has virtual columns for ! all possible glc_mec columns. ! ! For the sake of coupling with CISM, this should only be needed within the icemask, ! where we need virtual columns for the sake of coupling with CISM. This is needed in ! order to (1) provide SMB in all elevation classes, in case it is being used with ! 1-way coupling (or to force a later TG run); (2) even with two-way coupling, ! provide SMB in the elevation classes above and below existing elevation classes, ! for the sake of vertical interpolation; (3) provide place-holder columns (which are ! already spun-up) for dynamic landunits; (4) ensure that all glacier columns are ! given spun-up initial conditions by init_interp. ! ! More details on (4) (echoing the similar comment in subgridWeightsMod): We need all ! glacier and vegetated points to be active in the icemask region for the sake of ! init_interp - since we only interpolate onto active points, and we don't know which ! points will have non-zero area until after initialization (as long as we can't send ! information from glc to clm in initialization). (If we had an inactive glacier ! point in the icemask region, according to the weights on the surface dataset, and ! ran init_interp, this point would keep its cold start initialization values. Then, ! in the first time step of the run loop, it's possible that this point would become ! active because, according to glc, there is actually > 0% glacier in that grid ! cell. We don't do any state / flux adjustments in the first time step after ! init_interp due to glacier area changes, so this glacier column would remain at its ! cold start initialization values, which would be a Bad Thing. Ensuring that all ! glacier points within the icemask are active gets around this problem - as well as ! having other benefits, as noted above.) ! ! However, by making this part of the user-modifiable "glc behavior", we make it easy ! for the user to add virtual columns, if this is desired for diagnostic ! purposes. One important reason why this may be desired is to produce coupler ! history forcings to force a later TG run, with SMB forcings outside the original ! CISM area. (Also, we cannot use icemask for all purposes, because it isn't known at ! initialization.) logical, allocatable, public :: has_virtual_columns_grc(:) ! If allow_multiple_columns_grc(g) is true, then grid cell g may have multiple ! glacier columns, for the different elevation classes. If ! allow_multiple_columns_grc(g) is false, then grid cell g is guaranteed to have at ! most one glacier column. 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. logical, allocatable, public :: melt_replaced_by_ice_grc(:) ! If ice_runoff_melted_grc(g) is true, then ice runoff generated by the ! CLM physics over glacier columns in gridcell g is melted (generating a negative ! sensible heat flux) and runs off as liquid. If it is false, then ice runoff is ! sent to the river model as ice (a crude parameterization of iceberg calving). logical, allocatable, public :: ice_runoff_melted_grc(:) ! If rain_to_snow_runs_off_grc(g) is true, then any rain that would be converted to ! snow through precipitation repartitioning/downscaling instead runs off immediately. ! ! Note that, unlike the other options here, this applies to ALL landunits in the ! given region, so that the physics are consistent between different landunit types. logical, allocatable, public :: rain_to_snow_runs_off_grc(:) ! ------------------------------------------------------------------------ ! Private data ! ------------------------------------------------------------------------ ! If collapse_to_atm_topo_grc(g) is true, then grid cell g has at most one glc_mec ! column, whose topographic height exactly matches the atmosphere's topographic ! height for that grid cell (so that there is no adjustment of atmospheric ! forcings). ! ! Note that has_virtual_columns_grc(g) is guaranteed to be false if ! collapse_to_atm_topo_grc(g) is true. logical, allocatable :: collapse_to_atm_topo_grc(:) contains ! ------------------------------------------------------------------------ ! Public routines ! ------------------------------------------------------------------------ procedure, public :: Init ! version of Init meant for production use procedure, public :: InitFromInputs ! version of Init meant for unit testing (and called by other code in this class) procedure, public :: InitSetDirectly ! version of Init meant for unit testing ! get number of subgrid units in glc_mec landunit on one grid cell procedure, public :: get_num_glc_mec_subgrid ! returns true if memory should be allocated for the given glc_mec column, and its ! weight on the landunit procedure, public :: glc_mec_col_exists ! returns true if glc_mec columns on the given grid cell have dynamic type (type ! potentially changing at runtime) procedure, public :: cols_have_dynamic_type ! Sets a column-level logical array to true for any ice_mec column that needs ! downscaling, false for any ice_mec column that does not need downscaling procedure, public :: icemec_cols_need_downscaling ! Sets a column-level logical array to true for any ice_mec column that has ! dynamic type, false for any ice_mec column that does not have dynamic type procedure, public :: cols_have_dynamic_type_array ! Sets a patch-level logical array to true for any ice_mec column that has ! dynamic type, false for any ice_mec column that does not have dynamic type procedure, public :: patches_have_dynamic_type_array ! update the column class types of any glc_mec columns that need to be updated procedure, public :: update_glc_classes ! ------------------------------------------------------------------------ ! Public routines, for unit tests only ! ------------------------------------------------------------------------ ! get the value of collapse_to_atm_topo at a given grid cell procedure, public :: get_collapse_to_atm_topo ! ------------------------------------------------------------------------ ! Private routines ! ------------------------------------------------------------------------ procedure, private :: InitAllocate ! reads GLACIER_REGION field from surface dataset procedure, private, nopass :: read_surface_dataset ! reads local namelist items procedure, private, nopass :: read_namelist ! returns a column-level filter of ice_mec columns with the collapse_to_atm_topo ! behavior procedure, private :: collapse_to_atm_topo_icemec_filterc ! update class of glc_mec columns in regions where these are collapsed to a single ! column, given a filter procedure, private :: update_collapsed_columns_classes end type glc_behavior_type ! !PRIVATE MEMBER DATA: ! Longest name allowed for glacier_region_behavior, glacier_region_melt_behavior, ! glacier_region_ice_runoff_behavior and glacier_region_rain_to_snow_behavior integer, parameter :: max_behavior_name_len = 32 ! Smallest and largest allowed values for a glacier region ID integer, parameter :: min_glacier_region_id = 0 integer, parameter :: max_glacier_region_id = 10 ! Value indicating that a namelist item has not been set character(len=*), parameter :: behavior_str_unset = 'UNSET' character(len=*), parameter, private :: sourcefile = & __FILE__ contains !----------------------------------------------------------------------- subroutine Init(this, begg, endg, NLFilename) ! ! !DESCRIPTION: ! Initialize a glc_behavior_type object. ! ! This version of Init is the one intended for production code use. It reads the ! information it needs from the surface dataset and namelist. ! ! !USES: ! ! !ARGUMENTS: class(glc_behavior_type), intent(inout) :: this integer, intent(in) :: begg ! beginning gridcell index integer, intent(in) :: endg ! ending gridcell index character(len=*), intent(in) :: NLFilename ! Namelist filename ! ! !LOCAL VARIABLES: integer, allocatable :: glacier_region_map(:) character(len=max_behavior_name_len) :: glacier_region_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len) :: glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len) :: glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len) :: glacier_region_rain_to_snow_behavior(min_glacier_region_id:max_glacier_region_id) character(len=*), parameter :: subname = 'Init' !----------------------------------------------------------------------- allocate(glacier_region_map(begg:endg)) call this%read_surface_dataset(begg, endg, glacier_region_map(begg:endg)) call this%read_namelist(NLFilename, glacier_region_behavior, & glacier_region_melt_behavior, glacier_region_ice_runoff_behavior, & glacier_region_rain_to_snow_behavior) call this%InitFromInputs(begg, endg, & glacier_region_map(begg:endg), glacier_region_behavior, & glacier_region_melt_behavior, glacier_region_ice_runoff_behavior, & glacier_region_rain_to_snow_behavior) end subroutine Init !----------------------------------------------------------------------- subroutine InitFromInputs(this, begg, endg, & glacier_region_map, glacier_region_behavior_str, glacier_region_melt_behavior_str, & glacier_region_ice_runoff_behavior_str, glacier_region_rain_to_snow_behavior_str) ! ! !DESCRIPTION: ! Initialize a glc_behavior_type object given a map of glacier region IDs and an ! array of behavior specifications for each of these IDs. ! ! This version should generally only be called directly by tests, but it is also used ! by the main production Init method. ! ! !USES: ! ! !ARGUMENTS: class(glc_behavior_type), intent(inout) :: this integer, intent(in) :: begg ! beginning gridcell index integer, intent(in) :: endg ! ending gridcell index ! map of glacier region IDs integer, intent(in) :: glacier_region_map(begg:) ! string giving behavior for each glacier region ID ! allowed values are: ! - 'multiple': grid cells can potentially have multiple glacier elevation classes, ! but no virtual columns ! - 'virtual': grid cells have virtual columns: values are computed for every glacier ! elevation class, even those with 0 area ! - 'single_at_atm_topo': glacier landunits in these grid cells have a single column, ! whose elevation matches the atmosphere's topographic height (so that there is no ! adjustment due to downscaling) character(len=*), intent(in) :: glacier_region_behavior_str(min_glacier_region_id:) ! string giving treatment of ice melt for each glacier region ID ! allowed values are: ! - 'replaced_by_ice' ! - 'remains_in_place' character(len=*), intent(in) :: glacier_region_melt_behavior_str(min_glacier_region_id:) ! string giving treatment of ice runoff for each glacier region ID ! allowed values are: ! - 'remains_ice' ! - 'melted' character(len=*), intent(in) :: glacier_region_ice_runoff_behavior_str(min_glacier_region_id:) ! string giving treatment of rain-to-snow conversion for each glacier region ID ! allowed values are: ! - 'converted_to_snow' ! - 'runs_off' character(len=*), intent(in) :: glacier_region_rain_to_snow_behavior_str(min_glacier_region_id:) ! ! !LOCAL VARIABLES: ! whether each glacier region ID is present in the glacier_region_map logical :: glacier_region_present(min_glacier_region_id:max_glacier_region_id) ! integer codes corresponding to glacier_region_behavior_str integer :: glacier_region_behavior(min_glacier_region_id:max_glacier_region_id) ! integer codes corresponding to glacier_region_melt_behavior_str integer :: glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id) ! integer codes corresponding to glacier_region_ice_runoff_behavior_str integer :: glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id) ! integer codes corresponding to glacier_region_rain_to_snow_behavior_str integer :: glacier_region_rain_to_snow_behavior(min_glacier_region_id:max_glacier_region_id) integer :: g integer :: my_id integer :: my_behavior integer :: my_melt_behavior integer :: my_ice_runoff_behavior integer :: my_rain_to_snow_behavior ! possible glacier_region_behavior codes integer, parameter :: BEHAVIOR_MULTIPLE = 1 integer, parameter :: BEHAVIOR_VIRTUAL = 2 integer, parameter :: BEHAVIOR_SINGLE_AT_ATM_TOPO = 3 ! possible glacier_region_melt_behavior codes integer, parameter :: MELT_BEHAVIOR_REPLACED_BY_ICE = 1 integer, parameter :: MELT_BEHAVIOR_REMAINS_IN_PLACE = 2 ! possible glacier_region_ice_runoff_behavior codes integer, parameter :: ICE_RUNOFF_BEHAVIOR_REMAINS_ICE = 1 integer, parameter :: ICE_RUNOFF_BEHAVIOR_MELTED = 2 ! possible glacier_region_rain_to_snow_behavior codes integer, parameter :: RAIN_TO_SNOW_BEHAVIOR_CONVERTED_TO_SNOW = 1 integer, parameter :: RAIN_TO_SNOW_BEHAVIOR_RUNS_OFF = 2 ! value indicating that a behavior code has not been set (for glacier_region_behavior, ! glacier_region_melt_behavior, glacier_region_ice_runoff_behavior, or ! glacier_region_rain_to_snow_behavior) integer, parameter :: BEHAVIOR_UNSET = -1 character(len=*), parameter :: subname = 'InitFromInputs' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(glacier_region_map) == (/endg/)), errMsg(sourcefile, __LINE__)) call check_glacier_region_map call determine_region_presence call translate_glacier_region_behavior call translate_glacier_region_melt_behavior call translate_glacier_region_ice_runoff_behavior call translate_glacier_region_rain_to_snow_behavior call this%InitAllocate(begg, endg) do g = begg, endg my_id = glacier_region_map(g) my_behavior = glacier_region_behavior(my_id) my_melt_behavior = glacier_region_melt_behavior(my_id) my_ice_runoff_behavior = glacier_region_ice_runoff_behavior(my_id) my_rain_to_snow_behavior = glacier_region_rain_to_snow_behavior(my_id) ! This should only happen due to a programming error, not due to a user input error SHR_ASSERT(my_behavior /= BEHAVIOR_UNSET, errMsg(sourcefile, __LINE__)) SHR_ASSERT(my_melt_behavior /= BEHAVIOR_UNSET, errMsg(sourcefile, __LINE__)) SHR_ASSERT(my_ice_runoff_behavior /= BEHAVIOR_UNSET, errMsg(sourcefile, __LINE__)) SHR_ASSERT(my_rain_to_snow_behavior /= BEHAVIOR_UNSET, errMsg(sourcefile, __LINE__)) if (my_behavior == BEHAVIOR_VIRTUAL) then this%has_virtual_columns_grc(g) = .true. else this%has_virtual_columns_grc(g) = .false. end if if (my_melt_behavior == MELT_BEHAVIOR_REMAINS_IN_PLACE) then this%melt_replaced_by_ice_grc(g) = .false. else this%melt_replaced_by_ice_grc(g) = .true. end if if (my_ice_runoff_behavior == ICE_RUNOFF_BEHAVIOR_MELTED) then this%ice_runoff_melted_grc(g) = .true. else this%ice_runoff_melted_grc(g) = .false. end if if (my_rain_to_snow_behavior == RAIN_TO_SNOW_BEHAVIOR_RUNS_OFF) then this%rain_to_snow_runs_off_grc(g) = .true. else this%rain_to_snow_runs_off_grc(g) = .false. end if ! For now, allow_multiple_columns_grc is simply the opposite of ! collapse_to_atm_topo_grc. However, we maintain the separate ! allow_multiple_columns_grc so that the public interface can stay the same if we ! differentiate between the two in the future - e.g., allowing for the possibility ! of a behavior where we have at most one glacier column, but not forced to the ! atmosphere's elevation. if (my_behavior == BEHAVIOR_SINGLE_AT_ATM_TOPO) then this%collapse_to_atm_topo_grc(g) = .true. this%allow_multiple_columns_grc(g) = .false. else this%collapse_to_atm_topo_grc(g) = .false. this%allow_multiple_columns_grc(g) = .true. end if end do contains subroutine check_glacier_region_map if (minval(glacier_region_map) < min_glacier_region_id) then write(iulog,*) subname//' ERROR: Expect GLACIER_REGION to be >= ', min_glacier_region_id write(iulog,*) 'minval = ', minval(glacier_region_map) call endrun(msg=' ERROR: GLACIER_REGION smaller than expected'// & errMsg(sourcefile, __LINE__)) end if if (maxval(glacier_region_map) > max_glacier_region_id) then write(iulog,*) subname//' ERROR: Max GLACIER_REGION is ', & maxval(glacier_region_map) write(iulog,*) 'but max_glacier_region_id is only ', max_glacier_region_id write(iulog,*) 'Try increasing max_glacier_region_id in ', sourcefile call endrun(msg=' ERROR: GLACIER_REGION larger than expected'// & errMsg(sourcefile, __LINE__)) end if end subroutine check_glacier_region_map subroutine determine_region_presence integer :: g integer :: my_id glacier_region_present(:) = .false. do g = begg, endg my_id = glacier_region_map(g) glacier_region_present(my_id) = .true. end do end subroutine determine_region_presence subroutine translate_glacier_region_behavior integer :: i do i = min_glacier_region_id, max_glacier_region_id glacier_region_behavior(i) = BEHAVIOR_UNSET if (glacier_region_present(i)) then SHR_ASSERT_ALL((ubound(glacier_region_behavior_str) >= (/i/)), errMsg(sourcefile, __LINE__)) select case (glacier_region_behavior_str(i)) case ('multiple') glacier_region_behavior(i) = BEHAVIOR_MULTIPLE case ('virtual') glacier_region_behavior(i) = BEHAVIOR_VIRTUAL case ('single_at_atm_topo') glacier_region_behavior(i) = BEHAVIOR_SINGLE_AT_ATM_TOPO case (behavior_str_unset) write(iulog,*) ' ERROR: glacier_region_behavior not specified for ID ', i write(iulog,*) 'You probably need to extend the glacier_region_behavior namelist array' call endrun(msg=' ERROR: glacier_region_behavior not specified for ID '// & errMsg(sourcefile, __LINE__)) case default write(iulog,*) ' ERROR: Unknown glacier_region_behavior for ID ', i write(iulog,*) glacier_region_behavior_str(i) write(iulog,*) 'Allowable values are: multiple, virtual, single_at_atm_topo' call endrun(msg=' ERROR: Unknown glacier_region_behavior'// & errMsg(sourcefile, __LINE__)) end select end if end do end subroutine translate_glacier_region_behavior subroutine translate_glacier_region_melt_behavior integer :: i do i = min_glacier_region_id, max_glacier_region_id glacier_region_melt_behavior(i) = BEHAVIOR_UNSET if (glacier_region_present(i)) then SHR_ASSERT_ALL((ubound(glacier_region_melt_behavior_str) >= (/i/)), errMsg(sourcefile, __LINE__)) select case (glacier_region_melt_behavior_str(i)) case ('replaced_by_ice') glacier_region_melt_behavior(i) = MELT_BEHAVIOR_REPLACED_BY_ICE case ('remains_in_place') glacier_region_melt_behavior(i) = MELT_BEHAVIOR_REMAINS_IN_PLACE case (behavior_str_unset) write(iulog,*) ' ERROR: glacier_region_melt_behavior not specified for ID ', i write(iulog,*) 'You probably need to extend the glacier_region_melt_behavior namelist array' call endrun(msg=' ERROR: glacier_region_melt_behavior not specified for ID '// & errMsg(sourcefile, __LINE__)) case default write(iulog,*) ' ERROR: Unknown glacier_region_melt_behavior for ID ', i write(iulog,*) glacier_region_melt_behavior_str(i) write(iulog,*) 'Allowable values are: replaced_by_ice, remains_in_place' call endrun(msg=' ERROR: Unknown glacier_region_melt_behavior'// & errMsg(sourcefile, __LINE__)) end select end if end do end subroutine translate_glacier_region_melt_behavior subroutine translate_glacier_region_ice_runoff_behavior integer :: i do i = min_glacier_region_id, max_glacier_region_id glacier_region_ice_runoff_behavior(i) = BEHAVIOR_UNSET if (glacier_region_present(i)) then SHR_ASSERT_ALL((ubound(glacier_region_ice_runoff_behavior_str) >= (/i/)), errMsg(sourcefile, __LINE__)) select case (glacier_region_ice_runoff_behavior_str(i)) case ('remains_ice') glacier_region_ice_runoff_behavior(i) = ICE_RUNOFF_BEHAVIOR_REMAINS_ICE case('melted') glacier_region_ice_runoff_behavior(i) = ICE_RUNOFF_BEHAVIOR_MELTED case (behavior_str_unset) write(iulog,*) ' ERROR: glacier_region_ice_runoff_behavior not specified for ID ', i write(iulog,*) 'You probably need to extend the glacier_region_ice_runoff_behavior namelist array' call endrun(msg=' ERROR: glacier_region_ice_runoff_behavior not specified for ID '// & errMsg(sourcefile, __LINE__)) case default write(iulog,*) ' ERROR: Unknown glacier_region_ice_runoff_behavior for ID ', i write(iulog,*) glacier_region_ice_runoff_behavior_str(i) write(iulog,*) 'Allowable values are: remains_ice, melted' call endrun(msg=' ERROR: Unknown glacier_region_ice_runoff_behavior'// & errMsg(sourcefile, __LINE__)) end select end if end do end subroutine translate_glacier_region_ice_runoff_behavior subroutine translate_glacier_region_rain_to_snow_behavior integer :: i do i = min_glacier_region_id, max_glacier_region_id glacier_region_rain_to_snow_behavior(i) = BEHAVIOR_UNSET if (glacier_region_present(i)) then SHR_ASSERT_ALL((ubound(glacier_region_rain_to_snow_behavior_str) >= [i]), errMsg(sourcefile, __LINE__)) select case (glacier_region_rain_to_snow_behavior_str(i)) case ('converted_to_snow') glacier_region_rain_to_snow_behavior(i) = RAIN_TO_SNOW_BEHAVIOR_CONVERTED_TO_SNOW case ('runs_off') glacier_region_rain_to_snow_behavior(i) = RAIN_TO_SNOW_BEHAVIOR_RUNS_OFF case (behavior_str_unset) write(iulog,*) ' ERROR: glacier_region_rain_to_snow_behavior not specified for ID ', i write(iulog,*) 'You probably need to extend the glacier_region_rain_to_snow_behavior namelist array' call endrun(msg=' ERROR: glacier_region_rain_to_snow_behavior not specified for ID '// & errMsg(sourcefile, __LINE__)) case default write(iulog,*) ' ERROR: Unknown glacier_region_rain_to_snow_behavior for ID ', i write(iulog,*) glacier_region_rain_to_snow_behavior_str(i) write(iulog,*) 'Allowable values are: converted_to_snow, runs_off' call endrun(msg=' ERROR: Unknown glacier_region_rain_to_snow_behavior'// & errMsg(sourcefile, __LINE__)) end select end if end do end subroutine translate_glacier_region_rain_to_snow_behavior end subroutine InitFromInputs !----------------------------------------------------------------------- subroutine InitSetDirectly(this, begg, endg, & has_virtual_columns, collapse_to_atm_topo, rain_to_snow_runs_off) ! ! !DESCRIPTION: ! Initialize a glc_behavior_type object by directly setting provided fields ! ! This version is meant for testing ! ! !USES: ! ! !ARGUMENTS: class(glc_behavior_type), intent(inout) :: this integer, intent(in) :: begg ! beginning gridcell index integer, intent(in) :: endg ! ending gridcell index logical, intent(in), optional :: has_virtual_columns(begg:) logical, intent(in), optional :: collapse_to_atm_topo(begg:) logical, intent(in), optional :: rain_to_snow_runs_off(begg:) ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'InitForTesting' !----------------------------------------------------------------------- if (present(has_virtual_columns)) then SHR_ASSERT_ALL((ubound(has_virtual_columns) == (/endg/)), errMsg(sourcefile, __LINE__)) end if if (present(collapse_to_atm_topo)) then SHR_ASSERT_ALL((ubound(collapse_to_atm_topo) == (/endg/)), errMsg(sourcefile, __LINE__)) end if if (present(rain_to_snow_runs_off)) then SHR_ASSERT_ALL((ubound(rain_to_snow_runs_off) == [endg]), errMsg(sourcefile, __LINE__)) end if call this%InitAllocate(begg, endg) if (present(has_virtual_columns)) then this%has_virtual_columns_grc(:) = has_virtual_columns(:) end if if (present(collapse_to_atm_topo)) then this%collapse_to_atm_topo_grc(:) = collapse_to_atm_topo(:) end if if (present(rain_to_snow_runs_off)) then this%rain_to_snow_runs_off_grc(:) = rain_to_snow_runs_off(:) end if end subroutine InitSetDirectly !----------------------------------------------------------------------- subroutine InitAllocate(this, begg, endg) ! ! !DESCRIPTION: ! Allocate variables in this object ! ! !ARGUMENTS: class(glc_behavior_type), intent(inout) :: this integer, intent(in) :: begg ! beginning gridcell index integer, intent(in) :: endg ! ending gridcell index ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'InitAllocate' !----------------------------------------------------------------------- allocate(this%has_virtual_columns_grc (begg:endg)); this%has_virtual_columns_grc (:) = .false. allocate(this%allow_multiple_columns_grc(begg:endg)); this%allow_multiple_columns_grc(:) = .false. allocate(this%melt_replaced_by_ice_grc(begg:endg)); this%melt_replaced_by_ice_grc(:) = .false. allocate(this%collapse_to_atm_topo_grc(begg:endg)); this%collapse_to_atm_topo_grc(:) = .false. allocate(this%ice_runoff_melted_grc(begg:endg)); this%ice_runoff_melted_grc(:) = .false. allocate(this%rain_to_snow_runs_off_grc(begg:endg)); this%rain_to_snow_runs_off_grc(:) = .false. end subroutine InitAllocate !----------------------------------------------------------------------- subroutine read_surface_dataset(begg, endg, glacier_region_map) ! ! !DESCRIPTION: ! Reads GLACIER_REGION field from surface dataset, returns it in glacier_region_map ! ! !USES: use clm_varctl , only : fsurdat use fileutils , only : getfil use ncdio_pio , only : file_desc_t, ncd_io, ncd_pio_openfile, ncd_pio_closefile use spmdMod , only : masterproc use clm_varcon , only : grlnd ! ! !ARGUMENTS: integer, intent(in) :: begg ! beginning grid cell index integer, intent(in) :: endg ! ending grid cell index integer, intent(out) :: glacier_region_map(begg:) ! ! !LOCAL VARIABLES: integer, pointer :: glacier_region_map_ptr(:) ! pointer version needed for ncd_io interface character(len=256) :: locfn ! local filename type(file_desc_t) :: ncid ! netcdf id logical :: readvar character(len=*), parameter :: subname = 'read_surface_dataset' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(glacier_region_map) == (/endg/)), errMsg(sourcefile, __LINE__)) if (masterproc) then write(iulog,*) 'Attempting to read GLACIER_REGION...' end if call getfil(fsurdat, locfn, 0) call ncd_pio_openfile(ncid, locfn, 0) allocate(glacier_region_map_ptr(begg:endg)) call ncd_io(ncid=ncid, varname='GLACIER_REGION', flag='read', & data=glacier_region_map_ptr, dim1name=grlnd, readvar=readvar) if (.not. readvar) then call endrun(msg=' ERROR: GLACIER_REGION NOT on surfdata file'//errMsg(sourcefile, __LINE__)) end if call ncd_pio_closefile(ncid) glacier_region_map(begg:endg) = glacier_region_map_ptr(begg:endg) deallocate(glacier_region_map_ptr) end subroutine read_surface_dataset !----------------------------------------------------------------------- subroutine read_namelist(NLFilename, glacier_region_behavior, & glacier_region_melt_behavior, glacier_region_ice_runoff_behavior, & glacier_region_rain_to_snow_behavior) ! ! !DESCRIPTION: ! Read local namelist items ! ! !USES: use fileutils , only : getavu, relavu, opnfil use shr_nl_mod , only : shr_nl_find_group_name use clm_nlUtilsMod , only : find_nlgroup_name use spmdMod , only : masterproc, mpicom use shr_mpi_mod , only : shr_mpi_bcast ! ! !ARGUMENTS: character(len=*), intent(in) :: NLFilename ! Namelist filename character(len=max_behavior_name_len), intent(out) :: & glacier_region_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len), intent(out) :: & glacier_region_melt_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len), intent(out) :: & glacier_region_ice_runoff_behavior(min_glacier_region_id:max_glacier_region_id) character(len=max_behavior_name_len), intent(out) :: & glacier_region_rain_to_snow_behavior(min_glacier_region_id:max_glacier_region_id) ! ! !LOCAL VARIABLES: integer :: unitn ! unit for namelist file integer :: nml_error ! namelist i/o error flag character(len=*), parameter :: subname = 'read_namelist' !----------------------------------------------------------------------- namelist /clm_glacier_behavior/ & glacier_region_behavior, glacier_region_melt_behavior, & glacier_region_ice_runoff_behavior, glacier_region_rain_to_snow_behavior ! Initialize options to default values glacier_region_behavior(:) = behavior_str_unset glacier_region_melt_behavior(:) = behavior_str_unset glacier_region_ice_runoff_behavior(:) = behavior_str_unset glacier_region_rain_to_snow_behavior(:) = behavior_str_unset if (masterproc) then unitn = getavu() call opnfil(NLFilename, unitn, 'F') call shr_nl_find_group_name(unitn, 'clm_glacier_behavior', status=nml_error) if (nml_error == 0) then read(unitn, nml=clm_glacier_behavior, iostat=nml_error) if (nml_error /= 0) then call endrun(msg='ERROR reading clm_glacier_behavior namelist'// & errMsg(sourcefile, __LINE__)) end if else call endrun(msg='ERROR finding clm_glacier_behavior namelist'// & errMsg(sourcefile, __LINE__)) end if call relavu( unitn ) endif call shr_mpi_bcast(glacier_region_behavior, mpicom) call shr_mpi_bcast(glacier_region_melt_behavior, mpicom) call shr_mpi_bcast(glacier_region_ice_runoff_behavior, mpicom) call shr_mpi_bcast(glacier_region_rain_to_snow_behavior, mpicom) if (masterproc) then write(iulog,*) ' ' write(iulog,*) 'clm_glacier_behavior settings:' write(iulog,nml=clm_glacier_behavior) write(iulog,*) ' ' end if end subroutine read_namelist !----------------------------------------------------------------------- subroutine get_num_glc_mec_subgrid(this, gi, atm_topo, npatches, ncols, nlunits) ! ! !DESCRIPTION: ! Get number of subgrid units in glc_mec landunit on one grid cell ! ! !USES: use clm_varpar , only : maxpatch_glcmec ! ! !ARGUMENTS: class(glc_behavior_type), intent(in) :: this integer , intent(in) :: gi ! grid cell index real(r8), intent(in) :: atm_topo ! atmosphere's topographic height for this grid cell (m) integer , intent(out) :: npatches ! number of glacier_mec patches in this grid cell integer , intent(out) :: ncols ! number of glacier_mec columns in this grid cell integer , intent(out) :: nlunits ! number of glacier_mec landunits in this grid cell ! ! !LOCAL VARIABLES: integer :: m ! loop index logical :: col_exists real(r8) :: col_wt_lunit character(len=*), parameter :: subname = 'get_num_glc_mec_subgrid' !----------------------------------------------------------------------- ncols = 0 do m = 1, maxpatch_glcmec call this%glc_mec_col_exists(gi = gi, elev_class = m, atm_topo = atm_topo, & exists = col_exists, col_wt_lunit = col_wt_lunit) if (col_exists) then ncols = ncols + 1 end if end do if (this%collapse_to_atm_topo_grc(gi) .and. & wt_lunit(gi, istice_mec) > 0.0_r8) then ! For grid cells with the collapse_to_atm_topo behavior, with a non-zero weight ! ice_mec landunit, we expect exactly one column SHR_ASSERT(ncols == 1, errMsg(sourcefile, __LINE__)) end if if (ncols > 0) then npatches = ncols nlunits = 1 else npatches = 0 nlunits = 0 end if end subroutine get_num_glc_mec_subgrid !----------------------------------------------------------------------- subroutine glc_mec_col_exists(this, gi, elev_class, atm_topo, exists, col_wt_lunit) ! ! !DESCRIPTION: ! For the given glc_mec column, with elevation class index elev_class, in grid cell ! gi: sets exists to true if memory should be allocated for this column, and sets ! col_wt_lunit to the column's weight on the icemec landunit. ! ! If exists is false, then col_wt_lunit is arbitrary and should be ignored. ! ! !USES: use glc_elevclass_mod, only : glc_get_elevation_class, GLC_ELEVCLASS_ERR_NONE use glc_elevclass_mod, only : GLC_ELEVCLASS_ERR_TOO_LOW, GLC_ELEVCLASS_ERR_TOO_HIGH use glc_elevclass_mod, only : glc_errcode_to_string ! ! !ARGUMENTS: class(glc_behavior_type), intent(in) :: this integer, intent(in) :: gi ! grid cell index integer, intent(in) :: elev_class ! elevation class index real(r8), intent(in) :: atm_topo ! atmosphere's topographic height for this grid cell (m) logical, intent(out) :: exists ! whether memory should be allocated for this column real(r8), intent(out) :: col_wt_lunit ! column's weight on the icemec landunit ! ! !LOCAL VARIABLES: integer :: atm_elev_class ! elevation class corresponding to atmosphere topographic height integer :: err_code character(len=*), parameter :: subname = 'glc_mec_col_exists' !----------------------------------------------------------------------- ! Set default outputs exists = .false. col_wt_lunit = wt_glc_mec(gi, elev_class) if (this%collapse_to_atm_topo_grc(gi)) then if (wt_lunit(gi, istice_mec) > 0.0_r8) then call glc_get_elevation_class(atm_topo, atm_elev_class, err_code) if ( err_code == GLC_ELEVCLASS_ERR_NONE .or. & err_code == GLC_ELEVCLASS_ERR_TOO_LOW .or. & err_code == GLC_ELEVCLASS_ERR_TOO_HIGH) then ! These are all acceptable "errors" - it is even okay for these purposes if ! the elevation is lower than the lower bound of elevation class 1, or ! higher than the upper bound of the top elevation class. ! Do nothing else write(iulog,*) subname, ': ERROR getting elevation class for topo = ', atm_topo write(iulog,*) glc_errcode_to_string(err_code) call endrun(msg=subname//': ERROR getting elevation class') end if if (elev_class == atm_elev_class) then exists = .true. col_wt_lunit = 1._r8 else exists = .false. col_wt_lunit = 0._r8 end if end if else ! collapse_to_atm_topo_grc .false. if (this%has_virtual_columns_grc(gi)) then exists = .true. else if (wt_lunit(gi, istice_mec) > 0.0_r8 .and. & wt_glc_mec(gi, elev_class) > 0.0_r8) then ! If the landunit has non-zero weight on the grid cell, and this column has ! non-zero weight on the landunit... exists = .true. end if end if end subroutine glc_mec_col_exists !----------------------------------------------------------------------- function cols_have_dynamic_type(this, gi) ! ! !DESCRIPTION: ! Returns true if glc_mec columns on the given grid cell have dynamic type (i.e., ! type potentially changing at runtime) ! ! !USES: ! ! !ARGUMENTS: logical :: cols_have_dynamic_type ! function result class(glc_behavior_type), intent(in) :: this integer, intent(in) :: gi ! grid cell index ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'cols_have_dynamic_type' !----------------------------------------------------------------------- if (this%collapse_to_atm_topo_grc(gi)) then cols_have_dynamic_type = .true. else cols_have_dynamic_type = .false. end if end function cols_have_dynamic_type !----------------------------------------------------------------------- subroutine icemec_cols_need_downscaling(this, bounds, num_icemecc, filter_icemecc, & needs_downscaling_col) ! ! !DESCRIPTION: ! Sets needs_downscaling_col to true for any ice_mec column that needs downscaling, ! false for any ice_mec column that does not need downscaling. ! ! Outside of filter_icemecc, leaves needs_downscaling_col untouched. ! ! !USES: ! ! !ARGUMENTS: class(glc_behavior_type) , intent(in) :: this type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_icemecc ! number of points in filter_icemecc integer , intent(in) :: filter_icemecc(:) ! col filter for ice_mec logical , intent(inout) :: needs_downscaling_col( bounds%begc: ) ! ! !LOCAL VARIABLES: integer :: fc integer :: c integer :: g character(len=*), parameter :: subname = 'icemec_cols_need_downscaling' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(needs_downscaling_col) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) do fc = 1, num_icemecc c = filter_icemecc(fc) g = col%gridcell(c) if (this%collapse_to_atm_topo_grc(g)) then needs_downscaling_col(c) = .false. else needs_downscaling_col(c) = .true. end if end do end subroutine icemec_cols_need_downscaling !----------------------------------------------------------------------- subroutine cols_have_dynamic_type_array(this, begc, endc, has_dynamic_type_col) ! ! !DESCRIPTION: ! Sets a column-level logical array to true for any ice_mec column that has ! dynamic type, false for any ice_mec column that does not have dynamic type. ! ! The value is undefined for non-ice_mec columns. ! ! !ARGUMENTS: class(glc_behavior_type) , intent(in) :: this integer, intent(in) :: begc ! beginning column index integer, intent(in) :: endc ! ending column index logical, intent(inout) :: has_dynamic_type_col( begc: ) ! ! !LOCAL VARIABLES: integer :: c integer :: g character(len=*), parameter :: subname = 'cols_have_dynamic_type_array' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(has_dynamic_type_col) == (/endc/)), errMsg(sourcefile, __LINE__)) do c = begc, endc g = col%gridcell(c) ! Users shouldn't rely on the values set for non-ice_mec columns, but it's simpler ! just to set this for all column types. if (this%collapse_to_atm_topo_grc(g)) then has_dynamic_type_col(c) = .true. else has_dynamic_type_col(c) = .false. end if end do end subroutine cols_have_dynamic_type_array !----------------------------------------------------------------------- subroutine patches_have_dynamic_type_array(this, begp, endp, has_dynamic_type_patch) ! ! !DESCRIPTION: ! Sets a patch-level logical array to true for any ice_mec patch that has ! dynamic type, false for any ice_mec patch that does not have dynamic type. ! ! The value is undefined for non-ice_mec patches. ! ! !ARGUMENTS: class(glc_behavior_type) , intent(in) :: this integer, intent(in) :: begp ! beginning patch index integer, intent(in) :: endp ! ending patch index logical, intent(inout) :: has_dynamic_type_patch( begp: ) ! ! !LOCAL VARIABLES: integer :: p integer :: g character(len=*), parameter :: subname = 'patches_have_dynamic_type_array' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(has_dynamic_type_patch) == (/endp/)), errMsg(sourcefile, __LINE__)) do p = begp, endp g = patch%gridcell(p) ! Users shouldn't rely on the values set for non-ice_mec patches, but it's simpler ! just to set this for all patch types. if (this%collapse_to_atm_topo_grc(g)) then has_dynamic_type_patch(p) = .true. else has_dynamic_type_patch(p) = .false. end if end do end subroutine patches_have_dynamic_type_array !----------------------------------------------------------------------- subroutine update_glc_classes(this, bounds, topo_col) ! ! !DESCRIPTION: ! Update the column class types of any glc_mec columns that need to be updated. ! ! Assumes that topo_col has already been set appropriately. ! ! !USES: ! ! !ARGUMENTS: class(glc_behavior_type), intent(in) :: this type(bounds_type), intent(in) :: bounds real(r8), intent(in) :: topo_col( bounds%begc: ) ! ! !LOCAL VARIABLES: type(filter_col_type) :: collapse_filterc character(len=*), parameter :: subname = 'update_glc_classes' !----------------------------------------------------------------------- collapse_filterc = this%collapse_to_atm_topo_icemec_filterc(bounds) call this%update_collapsed_columns_classes(bounds, collapse_filterc, topo_col) end subroutine update_glc_classes !----------------------------------------------------------------------- subroutine update_collapsed_columns_classes(this, bounds, collapse_filterc, topo_col) ! ! !DESCRIPTION: ! Update class of glc_mec columns in regions where these are collapsed to a single ! column, given a filter. ! ! Assumes that topo_col has already been updated appropriately for these columns. ! ! !USES: use glc_elevclass_mod, only : glc_get_elevation_class, GLC_ELEVCLASS_ERR_NONE use glc_elevclass_mod, only : GLC_ELEVCLASS_ERR_TOO_LOW, GLC_ELEVCLASS_ERR_TOO_HIGH use glc_elevclass_mod, only : glc_errcode_to_string use column_varcon , only : icemec_class_to_col_itype ! ! !ARGUMENTS: class(glc_behavior_type), intent(in) :: this type(bounds_type), intent(in) :: bounds type(filter_col_type), intent(in) :: collapse_filterc real(r8), intent(in) :: topo_col( bounds%begc: ) ! ! !LOCAL VARIABLES: integer :: fc ! filter index integer :: c ! column index integer :: elev_class ! elevation class of the single column on the ice_mec landunit integer :: err_code integer :: new_coltype character(len=*), parameter :: subname = 'update_collapsed_columns_classes' !----------------------------------------------------------------------- SHR_ASSERT_ALL((ubound(topo_col) == (/bounds%endc/)), errMsg(sourcefile, __LINE__)) do fc = 1, collapse_filterc%num c = collapse_filterc%indices(fc) call glc_get_elevation_class(topo_col(c), elev_class, err_code) if ( err_code == GLC_ELEVCLASS_ERR_NONE .or. & err_code == GLC_ELEVCLASS_ERR_TOO_LOW .or. & err_code == GLC_ELEVCLASS_ERR_TOO_HIGH) then ! These are all acceptable "errors" - it is even okay for these purposes if ! the elevation is lower than the lower bound of elevation class 1, or ! higher than the upper bound of the top elevation class. ! Do nothing else write(iulog,*) subname, ': ERROR getting elevation class for topo = ', & topo_col(c) write(iulog,*) glc_errcode_to_string(err_code) call endrun(msg=subname//': ERROR getting elevation class') end if new_coltype = icemec_class_to_col_itype(elev_class) if (new_coltype /= col%itype(c)) then call col%update_itype(c = c, itype = new_coltype) end if end do end subroutine update_collapsed_columns_classes !----------------------------------------------------------------------- function collapse_to_atm_topo_icemec_filterc(this, bounds) result(filter) ! ! !DESCRIPTION: ! Returns a column-level filter of ice_mec columns with the collapse_to_atm_topo behavior ! ! !USES: use filterColMod, only : filter_col_type, col_filter_from_grcflags_ltypes ! ! !ARGUMENTS: class(glc_behavior_type), intent(in) :: this type(filter_col_type) :: filter ! function result type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'collapse_to_atm_topo_icemec_filterc' !----------------------------------------------------------------------- ! Currently this creates the filter on the fly, recreating it every time this ! function is called. Since this is a static filter, we could just compute it once ! and save it, returning the already-computed filter when this function is called. ! However, the problem with that is the need to have a different filter for each ! clump (and potentially another filter for calls from outside a clump loop). This ! will become easier to handle if we rework CLM's threading so that there is a ! separate instance of each object for each clump: in that case, we'll have multiple ! instances of glc_behavior_type, each corresponding to one clump, each with its own ! filter. filter = col_filter_from_grcflags_ltypes( & bounds = bounds, & grcflags = this%collapse_to_atm_topo_grc(bounds%begg:bounds%endg), & ltypes = [istice_mec], & include_inactive = .true.) end function collapse_to_atm_topo_icemec_filterc !----------------------------------------------------------------------- function get_collapse_to_atm_topo(this, gi) result(collapse_to_atm_topo) ! ! !DESCRIPTION: ! Get the value of collapse_to_atm_topo at a given grid cell ! ! This function just exists to support unit testing, and should not be called from ! production code. ! ! !USES: ! ! !ARGUMENTS: logical :: collapse_to_atm_topo ! function result class(glc_behavior_type), intent(in) :: this integer, intent(in) :: gi ! grid cell index ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'get_collapse_to_atm_topo' !----------------------------------------------------------------------- collapse_to_atm_topo = this%collapse_to_atm_topo_grc(gi) end function get_collapse_to_atm_topo end module glcBehaviorMod