module subgridWeightsMod #include "shr_assert.h" !----------------------------------------------------------------------- ! !DESCRIPTION: ! Handles modifications, error-checks and diagnostics related to changing subgrid weights ! ! ----- Requirements for subgrid weights that are enforced here ----- ! ! (These requirements are checked in check_weights/weights_okay) ! ! Note: in the following, 'active' refers to a pft, column, landunit or grid cell over ! which computations are performed, and 'inactive' refers to a pft, column or landunit ! where computations are NOT performed (grid cells are always active). ! ! (1) For all columns, landunits and grid cells, the sum of all subgrid weights of its ! children (or grandchildren, etc.) is equal to 1. For example: ! - For all columns, the sum of all patch weights on the column equals 1 ! - For all landunits, the sum of all col weights on the landunit equals 1 ! - For all grid cells, the sum of all patch weights on the grid cell equals 1 ! - etc. ! ! (2) For all ACTIVE columns, landunits and grid cells, the sum of all subgrid weights of ! its ACTIVE children (or grandchildren, etc.) is equal to 1. For example: ! - For all active columns, the sum of all patch weights on the column equals 1 when ! just considering active pfts ! - For all active landunits, the sum of all col weights on the landunit equals 1 when ! just considering active cols ! - For ALL grid cells, the sum of all patch weights on the grid cell equals 1 when ! just considering active pfts -- note that all grid cells are considered active! ! - etc. ! ! (3) For all INACTIVE columns, landunits and grid cells, the sum of all subgrid weights of ! its ACTIVE children, grandchildren, etc. are equal to either 0 or 1. For example: ! - For all inactive columns, the sum of all patch weights on the column equals either 0 ! or 1 when just considering active pfts ! - For all inactive landunits, the sum of all col weights on the landunit equals ! either 0 or 1 when just considering active cols ! - etc. ! ! Another way of stating (2) and (3) is that the sum of weights of all ACTIVE pfts, cols ! or landunits on their parent/grandparent/etc. is always equal to either 0 or 1 -- and ! must be equal to 1 if this parent/grandparent, etc. is itself active. ! ! Note that, together, conditions (1) and (2) imply that any pft, col or landunit whose ! weight on the grid cell is non-zero must be active. In addition, these conditions imply ! that any patch whose weight on the column is non-zero must be active if the column is ! active (and similarly for any patch on an active landunit, and any col on an active ! landunit). ! ! ! ----- Implications of these requirements for computing subgrid averages ----- ! ! The preferred way to average from, say, patch to col is: ! colval(c) = 0 ! do p = pfti(c), pftf(c) ! if (active(p)) colval(c) = colval(c) + pftval(p) * wtcol(p) ! (where wtcol(p) is the weight of the patch on the column) ! If column c is active, then the above conditions guarantee that the pwtcol values ! included in the above sum will sum to 1. If column c is inactive, then the above ! conditions guarantee that the pwtcol values included in the above sum will sum to ! either 1 or 0; if they sum to 0, then colval(c) will remain 0. ! ! Another acceptable method is the following; this method accommodates some unknown ! fraction of pftval's being set to spval, and leaves colval set at spval if there are no ! valid patch values: ! colval(c) = spval ! sumwt(c) = 0 ! do p = pfti(c), pftf(c) ! if (active(p) .and. wtcol(p) /= 0) then ! if (pftval(p) /= spval) then ! if (sumwt(c) == 0) colval(c) = 0 ! colval(c) = colval(c) + pftval(p) * wtcol(p) ! sumwt(c) = sumwt(c) + wtcol(p) ! end if ! end if ! end do ! if (sumwt(c) /= 0) then ! colval(c) = colval(c) / sumwt(c) ! end if ! Note that here we check the condition (active(p) .and. wtcol(p) /= 0). We need to ! include a check for wtcol(p) /= 0 because we don't want to set colval(c) = 0 for zero- ! weight pfts in this line: ! if (sumwt(c) == 0) colval(c) = 0 ! And we include a check for active(p) because we don't want to assume that pftval(p) has ! been set to spval for inactive pfts -- we want to allow for the possibility that ! pftval(p) will be NaN for inactive pfts. ! ! ! !USES: 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, all_active, run_zero_weight_urban, use_fates use clm_varcon , only : nameg, namel, namec, namep use decompMod , only : bounds_type use GridcellType , only : grc use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch use glcBehaviorMod , only : glc_behavior_type ! ! PUBLIC TYPES: implicit none save private ! ! !PUBLIC MEMBER FUNCTIONS: public :: init_subgrid_weights_mod ! initialize stuff in this module public :: compute_higher_order_weights ! given p2c, c2l and l2g weights, compute other weights public :: set_active ! set 'active' flags at pft, column & landunit level public :: check_weights ! check subgrid weights public :: get_landunit_weight ! get the weight of a given landunit on a single grid cell public :: set_landunit_weight ! set the weight of a given landunit on a single grid cell public :: is_gcell_all_ltypeX ! determine whether a grid cell is 100% covered by the given landunit type public :: set_subgrid_diagnostic_fields ! set all subgrid weights diagnostic fields ! ! !REVISION HISTORY: ! Created by Bill Sacks ! ! !PRIVATE TYPES: type subgrid_weights_diagnostics_type ! This type contains diagnostics on subgrid weights, for output to the history file real(r8), pointer :: pct_landunit(:,:) ! % of each landunit on the grid cell [begg:endg, 1:max_lunit] real(r8), pointer :: pct_nat_pft(:,:) ! % of each pft, as % of landunit [begg:endg, natpft_lb:natpft_ub] real(r8), pointer :: pct_cft(:,:) ! % of each crop functional type, as % of landunit [begg:endg, cft_lb:cft_ub] real(r8), pointer :: pct_glc_mec(:,:) ! % of each glacier elevation class, as % of landunit [begg:endg, 1:maxpatch_glcmec] end type subgrid_weights_diagnostics_type type(subgrid_weights_diagnostics_type) :: subgrid_weights_diagnostics ! ! !PRIVATE MEMBER FUNCTIONS: private :: is_active_l ! determine whether the given landunit is active private :: is_active_c ! determine whether the given column is active private :: is_active_p ! determine whether the given patch is active private :: weights_okay ! determine if sum of weights satisfies requirements laid out above private :: set_pct_landunit_diagnostics ! set pct_landunit diagnostic field private :: set_pct_glc_mec_diagnostics ! set pct_glc_mec diagnostic field private :: set_pct_pft_diagnostics ! set pct_nat_pft & pct_cft diagnostic fields character(len=*), parameter, private :: sourcefile = & __FILE__ !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- subroutine init_subgrid_weights_mod(bounds) ! ! !DESCRIPTION: ! Initialize stuff in this module ! ! !USES: use landunit_varcon, only : max_lunit use clm_varpar , only : maxpatch_glcmec, natpft_size, cft_size use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use decompMod , only : BOUNDS_LEVEL_PROC use histFileMod , only : hist_addfld2d ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds ! proc bounds ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'init_subgrid_weights_mod' !----------------------------------------------------------------------- SHR_ASSERT(bounds%level == BOUNDS_LEVEL_PROC, errMsg(sourcefile, __LINE__)) ! ------------------------------------------------------------------------ ! Allocate variables in subgrid_weights_diagnostics ! ------------------------------------------------------------------------ ! Note that, because these variables are output to the history file, it appears that ! their lower bounds need to start at 1 (e.g., 1:natpft_size rather than ! natpft_lb:natpft_ub) allocate(subgrid_weights_diagnostics%pct_landunit(bounds%begg:bounds%endg, 1:max_lunit)) subgrid_weights_diagnostics%pct_landunit(:,:) = nan allocate(subgrid_weights_diagnostics%pct_nat_pft(bounds%begg:bounds%endg, 1:natpft_size)) subgrid_weights_diagnostics%pct_nat_pft(:,:) = nan allocate(subgrid_weights_diagnostics%pct_cft(bounds%begg:bounds%endg, 1:cft_size)) subgrid_weights_diagnostics%pct_cft(:,:) = nan allocate(subgrid_weights_diagnostics%pct_glc_mec(bounds%begg:bounds%endg, 1:maxpatch_glcmec)) subgrid_weights_diagnostics%pct_glc_mec(:,:) = nan ! ------------------------------------------------------------------------ ! Add history fields ! ------------------------------------------------------------------------ call hist_addfld2d (fname='PCT_LANDUNIT', units='%', type2d='ltype', & avgflag='A', long_name='% of each landunit on grid cell', & ptr_lnd=subgrid_weights_diagnostics%pct_landunit) if(.not.use_fates) then call hist_addfld2d (fname='PCT_NAT_PFT', units='%', type2d='natpft', & avgflag='A', long_name='% of each PFT on the natural vegetation (i.e., soil) landunit', & ptr_lnd=subgrid_weights_diagnostics%pct_nat_pft) end if if (cft_size > 0) then call hist_addfld2d (fname='PCT_CFT', units='%', type2d='cft', & avgflag='A', long_name='% of each crop on the crop landunit', & ptr_lnd=subgrid_weights_diagnostics%pct_cft) end if call hist_addfld2d (fname='PCT_GLC_MEC', units='%', type2d='glc_nec', & avgflag='A', long_name='% of each GLC elevation class on the glc_mec landunit', & ptr_lnd=subgrid_weights_diagnostics%pct_glc_mec) end subroutine init_subgrid_weights_mod !----------------------------------------------------------------------- subroutine compute_higher_order_weights(bounds) ! ! !DESCRIPTION: ! Assuming patch%wtcol, col%wtlunit and lun%wtgcell have already been computed, compute ! the "higher-order" weights: patch%wtlunit, patch%wtgcell and col%wtgcell, for all p and c ! ! !USES: ! ! !ARGUMENTS: implicit none type(bounds_type), intent(in) :: bounds ! clump bounds ! ! !LOCAL VARIABLES: integer :: p, c, l ! indices for pft, col & landunit !------------------------------------------------------------------------ do c = bounds%begc, bounds%endc l = col%landunit(c) col%wtgcell(c) = col%wtlunit(c) * lun%wtgcell(l) end do do p = bounds%begp, bounds%endp c = patch%column(p) patch%wtlunit(p) = patch%wtcol(p) * col%wtlunit(c) patch%wtgcell(p) = patch%wtcol(p) * col%wtgcell(c) end do end subroutine compute_higher_order_weights !----------------------------------------------------------------------- subroutine set_active(bounds, glc_behavior) ! ! !DESCRIPTION: ! Set 'active' flags at the pft, column and landunit level ! (note that grid cells are always active) ! ! This should be called whenever any weights change (e.g., patch weights on the column, ! landunit weights on the grid cell, etc.). ! ! Ensures that we don't have any active patch on an inactive column, or an active column on an ! inactive landunit (since these conditions could lead to garbage data) ! ! !USES: ! ! !ARGUMENTS: implicit none type(bounds_type), intent(in) :: bounds ! bounds type(glc_behavior_type), intent(in) :: glc_behavior ! ! !LOCAL VARIABLES: integer :: l,c,p ! loop counters character(len=*), parameter :: subname = 'set_active' !------------------------------------------------------------------------ do l = bounds%begl,bounds%endl lun%active(l) = is_active_l(l, glc_behavior) end do do c = bounds%begc,bounds%endc l = col%landunit(c) col%active(c) = is_active_c(c, glc_behavior) if (col%active(c) .and. .not. lun%active(l)) then write(iulog,*) trim(subname),' ERROR: active column found on inactive landunit', & 'at c = ', c, ', l = ', l call endrun(decomp_index=c, clmlevel=namec, msg=errMsg(sourcefile, __LINE__)) end if end do do p = bounds%begp,bounds%endp c = patch%column(p) patch%active(p) = is_active_p(p) if (patch%active(p) .and. .not. col%active(c)) then write(iulog,*) trim(subname),' ERROR: active patch found on inactive column', & 'at p = ', p, ', c = ', c call endrun(decomp_index=p, clmlevel=namep, msg=errMsg(sourcefile, __LINE__)) end if end do end subroutine set_active !----------------------------------------------------------------------- logical function is_active_l(l, glc_behavior) ! ! !DESCRIPTION: ! Determine whether the given landunit is active ! ! !USES: use landunit_varcon, only : istsoil, istice_mec, isturb_MIN, isturb_MAX ! ! !ARGUMENTS: implicit none integer, intent(in) :: l ! landunit index type(glc_behavior_type), intent(in) :: glc_behavior ! ! !LOCAL VARIABLES: integer :: g ! grid cell index !------------------------------------------------------------------------ if (all_active) then is_active_l = .true. else g =lun%gridcell(l) is_active_l = .false. ! ------------------------------------------------------------------------ ! General conditions under which is_active_l NEEDS to be true in order to satisfy ! the requirements laid out at the top of this module: ! ------------------------------------------------------------------------ if (lun%wtgcell(l) > 0) is_active_l = .true. ! ------------------------------------------------------------------------ ! Conditions under which is_active_p is set to true because we want extra virtual landunits: ! ------------------------------------------------------------------------ if (lun%itype(l) == istice_mec .and. & glc_behavior%has_virtual_columns_grc(g)) then is_active_l = .true. end if if ((lun%itype(l) >= isturb_MIN .and. lun%itype(l) <= isturb_MAX) .and. & run_zero_weight_urban) then is_active_l = .true. end if ! In general, include a virtual natural vegetation landunit. This aids ! initialization of a new landunit; and for runs that are coupled to CISM, this ! provides bare land SMB forcing even if there is no vegetated area. ! ! Also (echoing the similar comment in glcBehaviorMod): 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 vegetated ! 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% bare ! ground 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 vegetated column ! would remain at its cold start initialization values, which would be a Bad ! Thing. Ensuring that all vegetated points within the icemask are active gets ! around this problem - as well as having other benefits, as noted above.) if (lun%itype(l) == istsoil) then is_active_l = .true. end if end if end function is_active_l !----------------------------------------------------------------------- logical function is_active_c(c, glc_behavior) ! ! !DESCRIPTION: ! Determine whether the given column is active ! ! !USES: use landunit_varcon, only : istice_mec, isturb_MIN, isturb_MAX ! ! !ARGUMENTS: implicit none integer, intent(in) :: c ! column index type(glc_behavior_type), intent(in) :: glc_behavior ! ! !LOCAL VARIABLES: integer :: l ! landunit index integer :: g ! grid cell index !------------------------------------------------------------------------ if (all_active) then is_active_c = .true. else l =col%landunit(c) g =col%gridcell(c) is_active_c = .false. ! ------------------------------------------------------------------------ ! General conditions under which is_active_c NEEDS to be true in order to satisfy ! the requirements laid out at the top of this module: ! ------------------------------------------------------------------------ if (lun%active(l) .and. col%wtlunit(c) > 0._r8) is_active_c = .true. ! ------------------------------------------------------------------------ ! Conditions under which is_active_c is set to true because we want extra virtual columns: ! ------------------------------------------------------------------------ if (lun%itype(l) == istice_mec .and. & glc_behavior%has_virtual_columns_grc(g)) then is_active_c = .true. end if ! We don't really need to run over 0-weight urban columns. But because of some ! messiness in the urban code (many loops are over the landunit filter, then drill ! down to columns - so we would need to add 'col%active(c)' conditionals in many ! places) it keeps the code cleaner to run over 0-weight urban columns. This generally ! shouldn't add much computation time, since in most places, all urban columns are ! non-zero weight if the landunit is non-zero weight. if (lun%active(l) .and. (lun%itype(l) >= isturb_MIN .and. lun%itype(l) <= isturb_MAX)) then is_active_c = .true. end if end if end function is_active_c !----------------------------------------------------------------------- logical function is_active_p(p) ! ! !DESCRIPTION: ! Determine whether the given patch is active ! ! !USES: ! ! !ARGUMENTS: implicit none integer, intent(in) :: p ! patch index ! ! !LOCAL VARIABLES: integer :: c ! column index !------------------------------------------------------------------------ if (all_active) then is_active_p = .true. else c =patch%column(p) is_active_p = .false. ! ------------------------------------------------------------------------ ! General conditions under which is_active_p NEEDS to be true in order to satisfy ! the requirements laid out at the top of this module: ! ------------------------------------------------------------------------ if (col%active(c) .and. patch%wtcol(p) > 0._r8) is_active_p = .true. end if end function is_active_p !----------------------------------------------------------------------- function get_landunit_weight(g, ltype) result(weight) ! ! !DESCRIPTION: ! Get the subgrid weight of a given landunit type on a single grid cell ! ! !USES: use clm_varcon, only : ispval ! ! !ARGUMENTS: real(r8) :: weight ! function result integer , intent(in) :: g ! grid cell index integer , intent(in) :: ltype ! landunit type of interest ! ! !LOCAL VARIABLES: integer :: l ! landunit index character(len=*), parameter :: subname = 'get_landunit_weight' !----------------------------------------------------------------------- l = grc%landunit_indices(ltype, g) if (l == ispval) then weight = 0._r8 else weight = lun%wtgcell(l) end if end function get_landunit_weight !----------------------------------------------------------------------- subroutine set_landunit_weight(g, ltype, weight) ! ! !DESCRIPTION: ! Set the subgrid weight of a given landunit type on a single grid cell ! ! !USES: use clm_varcon, only : ispval ! ! !ARGUMENTS: integer , intent(in) :: g ! grid cell index integer , intent(in) :: ltype ! landunit type of interest real(r8), intent(in) :: weight ! new weight of this landunit ! ! !LOCAL VARIABLES: integer :: l ! landunit index character(len=*), parameter :: subname = 'set_landunit_weight' !----------------------------------------------------------------------- l = grc%landunit_indices(ltype, g) if (l /= ispval) then lun%wtgcell(l) = weight else if (weight > 0._r8) then write(iulog,*) subname//' ERROR: Attempt to assign non-zero weight to a non-existent landunit' write(iulog,*) 'g, l, ltype, weight = ', g, l, ltype, weight call endrun(decomp_index=l, clmlevel=namel, msg=errMsg(sourcefile, __LINE__)) end if end subroutine set_landunit_weight !----------------------------------------------------------------------- function is_gcell_all_ltypeX(g, ltype) result(all_ltypeX) ! ! !DESCRIPTION: ! Determine if the given grid cell is 100% covered by the landunit type given by ltype ! ! !USES: ! ! !ARGUMENTS: implicit none logical :: all_ltypeX ! function result integer, intent(in) :: g ! grid cell index integer, intent(in) :: ltype ! landunit type of interest ! ! !LOCAL VARIABLES: real(r8) :: wt_lunit ! subgrid weight of the given landunit real(r8), parameter :: tolerance = 1.e-13_r8 ! tolerance for checking whether landunit's weight is 1 character(len=*), parameter :: subname = 'is_gcell_all_ltypeX' !------------------------------------------------------------------------------ wt_lunit = get_landunit_weight(g, ltype) if (wt_lunit >= (1._r8 - tolerance)) then all_ltypeX = .true. else all_ltypeX = .false. end if end function is_gcell_all_ltypeX !------------------------------------------------------------------------------ subroutine check_weights (bounds, active_only) ! ! !DESCRIPTION: ! Check subgrid weights. ! ! This routine operates in two different modes, depending on the value of active_only. If ! active_only is true, then we check the sum of weights of the ACTIVE children, ! grandchildren, etc. of a given point. If active_only is false, then we check the sum of ! weights of ALL children, grandchildren, etc. of a given point. ! ! Normally this routine will be called twice: once with active_only=false, and once with ! active_only=true. ! ! !USES ! ! !ARGUMENTS implicit none type(bounds_type), intent(in) :: bounds ! bounds logical, intent(in) :: active_only ! true => check sum of weights just of ACTIVE children, grandchildren, etc. ! ! !LOCAL VARIABLES: integer :: g,l,c,p ! loop counters real(r8), allocatable :: sumwtcol(:), sumwtlunit(:), sumwtgcell(:) logical :: error_found ! true if we find an error character(len=*), parameter :: subname = 'check_weights' !------------------------------------------------------------------------------ allocate(sumwtcol(bounds%begc:bounds%endc)) allocate(sumwtlunit(bounds%begl:bounds%endl)) allocate(sumwtgcell(bounds%begg:bounds%endg)) error_found = .false. ! Check patch-level weights sumwtcol(bounds%begc : bounds%endc) = 0._r8 sumwtlunit(bounds%begl : bounds%endl) = 0._r8 sumwtgcell(bounds%begg : bounds%endg) = 0._r8 do p = bounds%begp,bounds%endp c = patch%column(p) l = patch%landunit(p) g = patch%gridcell(p) if ((active_only .and. patch%active(p)) .or. .not. active_only) then sumwtcol(c) = sumwtcol(c) + patch%wtcol(p) sumwtlunit(l) = sumwtlunit(l) + patch%wtlunit(p) sumwtgcell(g) = sumwtgcell(g) + patch%wtgcell(p) end if end do do c = bounds%begc,bounds%endc if (.not. weights_okay(sumwtcol(c), active_only, col%active(c))) then write(iulog,*) trim(subname),' ERROR: at c = ',c,'total PFT weight is ',sumwtcol(c), & 'active_only = ', active_only error_found = .true. end if end do do l = bounds%begl,bounds%endl if (.not. weights_okay(sumwtlunit(l), active_only, lun%active(l))) then write(iulog,*) trim(subname),' ERROR: at l = ',l,'total PFT weight is ',sumwtlunit(l), & 'active_only = ', active_only error_found = .true. end if end do do g = bounds%begg,bounds%endg if (.not. weights_okay(sumwtgcell(g), active_only, i_am_active=.true.)) then write(iulog,*) trim(subname),' ERROR: at g = ',g,'total PFT weight is ',sumwtgcell(g), & 'active_only = ', active_only error_found = .true. end if end do ! Check col-level weights sumwtlunit(bounds%begl : bounds%endl) = 0._r8 sumwtgcell(bounds%begg : bounds%endg) = 0._r8 do c = bounds%begc,bounds%endc l = col%landunit(c) g = col%gridcell(c) if ((active_only .and. col%active(c)) .or. .not. active_only) then sumwtlunit(l) = sumwtlunit(l) + col%wtlunit(c) sumwtgcell(g) = sumwtgcell(g) + col%wtgcell(c) end if end do do l = bounds%begl,bounds%endl if (.not. weights_okay(sumwtlunit(l), active_only, lun%active(l))) then write(iulog,*) trim(subname),' ERROR: at l = ',l,'total col weight is ',sumwtlunit(l), & 'active_only = ', active_only error_found = .true. end if end do do g = bounds%begg,bounds%endg if (.not. weights_okay(sumwtgcell(g), active_only, i_am_active=.true.)) then write(iulog,*) trim(subname),' ERROR: at g = ',g,'total col weight is ',sumwtgcell(g), & 'active_only = ', active_only error_found = .true. end if end do ! Check landunit-level weights sumwtgcell(bounds%begg : bounds%endg) = 0._r8 do l = bounds%begl,bounds%endl g = lun%gridcell(l) if ((active_only .and. lun%active(l)) .or. .not. active_only) then sumwtgcell(g) = sumwtgcell(g) + lun%wtgcell(l) end if end do do g = bounds%begg,bounds%endg if (.not. weights_okay(sumwtgcell(g), active_only, i_am_active=.true.)) then write(iulog,*) trim(subname),' ERROR: at g = ',g,'total lunit weight is ',sumwtgcell(g), & 'active_only = ', active_only error_found = .true. end if end do deallocate(sumwtcol, sumwtlunit, sumwtgcell) if (error_found) then write(iulog,*) ' ' write(iulog,*) 'If you are seeing this message at the beginning of a run with' write(iulog,*) 'use_init_interp = .true. and init_interp_method = "use_finidat_areas",' write(iulog,*) 'and you are seeing weights less than 1, then a likely cause is:' write(iulog,*) 'For the above-mentioned grid cell(s):' write(iulog,*) 'The matching input grid cell had some non-zero-weight subgrid type' write(iulog,*) 'that is not present in memory in the new run.' write(iulog,*) ' ' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Success end subroutine check_weights !----------------------------------------------------------------------- logical function weights_okay(sumwts, active_weights_only, i_am_active) ! ! !DESCRIPTION: ! Determine if sumwts (the sum of weights of children, grandchildren or ! great-grandchildren of a column, landunit or grid cell) satisfies the requirements laid ! out above. ! ! The way this is determined depends on the values of two other variables: ! - active_weights_only: does sumwts just include weights of active children, ! grandchildren or great-grandchilden? (alternative is that it includes weights of ALL ! children, grandchildren or great-grandchildren) ! - i_am_active: true if the column, landunit or grid cell of interest is active ! ! !ARGUMENTS: implicit none real(r8), intent(in) :: sumwts ! sum of weights of children, grandchildren or great-grandchildren logical , intent(in) :: active_weights_only ! true if sumwts just includes active children, etc. logical , intent(in) :: i_am_active ! true if the current point is active ! ! !LOCAL VARIABLES: logical :: weights_equal_1 real(r8), parameter :: tolerance = 1.e-12_r8 ! tolerance for checking whether weights sum to 1 !------------------------------------------------------------------------ weights_equal_1 = (abs(sumwts - 1._r8) <= tolerance) if (active_weights_only) then if (i_am_active) then ! condition (2) above weights_okay = weights_equal_1 else ! condition (3) above weights_okay = (sumwts == 0._r8 .or. weights_equal_1) end if else ! condition (1) above ! (note that i_am_active is irrelevant in this case) weights_okay = weights_equal_1 end if end function weights_okay !----------------------------------------------------------------------- subroutine set_subgrid_diagnostic_fields(bounds) ! ! !DESCRIPTION: ! Set history fields giving diagnostics about subgrid weights ! ! !USES: ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'set_subgrid_diagnostic_fields' !----------------------------------------------------------------------- call set_pct_landunit_diagnostics(bounds) ! Note: (MV, 10-17-14): The following has an use_fates if-block around it since ! the pct_pft_diagnostics referens to patch%itype(p) which is not used by ED ! Note: (SPM, 10-20-15): If this isn't set then debug mode with intel and ! yellowstone will fail when trying to write pct_nat_pft since it contains ! all NaN's. call set_pct_pft_diagnostics(bounds) call set_pct_glc_mec_diagnostics(bounds) end subroutine set_subgrid_diagnostic_fields !----------------------------------------------------------------------- subroutine set_pct_landunit_diagnostics(bounds) ! ! !DESCRIPTION: ! Set pct_landunit diagnostic field: % of each landunit on the grid cell ! ! !USES: ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: g, l ! grid cell & landunit indices integer :: ltype ! landunit type character(len=*), parameter :: subname = 'set_pct_landunit_diagnostics' !----------------------------------------------------------------------- subgrid_weights_diagnostics%pct_landunit(bounds%begg:bounds%endg, :) = 0._r8 do l = bounds%begl, bounds%endl g = lun%gridcell(l) ltype = lun%itype(l) subgrid_weights_diagnostics%pct_landunit(g, ltype) = lun%wtgcell(l) * 100._r8 end do end subroutine set_pct_landunit_diagnostics !----------------------------------------------------------------------- subroutine set_pct_glc_mec_diagnostics(bounds) ! ! !DESCRIPTION: ! Set pct_glc_mec diagnostic field: % of each glc_mec column on the glc_mec landunit ! ! Note: it's safe to call this even if we're not running with glc_mec, but in that ! case it won't do anything. ! ! Note that pct_glc_mec will be 0 for all elevation classes in a grid cell that does ! not have a glc_mec landunit. However, it will still sum to 100% for a grid cell ! that has a 0-weight (i.e., virtual) glc_mec landunit. ! ! !USES: use landunit_varcon, only : istice_mec use column_varcon, only : col_itype_to_icemec_class ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: c,l,g ! indices integer :: icemec_class ! icemec class (1..maxpatch_glcmec) character(len=*), parameter :: subname = 'set_pct_glc_mec_diagnostics' !----------------------------------------------------------------------- subgrid_weights_diagnostics%pct_glc_mec(bounds%begg:bounds%endg, :) = 0._r8 do c = bounds%begc, bounds%endc g = col%gridcell(c) l = col%landunit(c) if (lun%itype(l) == istice_mec) then icemec_class = col_itype_to_icemec_class(col%itype(c)) subgrid_weights_diagnostics%pct_glc_mec(g, icemec_class) = col%wtlunit(c) * 100._r8 end if end do end subroutine set_pct_glc_mec_diagnostics !----------------------------------------------------------------------- subroutine set_pct_pft_diagnostics(bounds) ! ! !DESCRIPTION: ! Set pct_nat_pft & pct_cft diagnostic fields: % of PFTs on their landunit ! ! !USES: use landunit_varcon, only : istsoil, istcrop use clm_varpar, only : natpft_lb, cft_lb ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds ! ! !LOCAL VARIABLES: integer :: p,l,g ! indices integer :: ptype ! patch itype integer :: ptype_1indexing ! patch itype, translated into 1-indexing for the given landunit type character(len=*), parameter :: subname = 'set_pct_pft_diagnostics' !----------------------------------------------------------------------- subgrid_weights_diagnostics%pct_nat_pft(bounds%begg:bounds%endg, :) = 0._r8 ! Note that pct_cft will be 0-size if cft_size is 0 (which can happen if we don't ! have a crop landunit). But it doesn't hurt to have this line setting all elements ! to 0, and doing this always allows us to avoid extra logic which could be a ! maintenance problem. subgrid_weights_diagnostics%pct_cft(bounds%begg:bounds%endg, :) = 0._r8 do p = bounds%begp,bounds%endp g = patch%gridcell(p) l = patch%landunit(p) ptype = patch%itype(p) if (lun%itype(l) == istsoil .and. (.not.use_fates) ) then ptype_1indexing = ptype + (1 - natpft_lb) subgrid_weights_diagnostics%pct_nat_pft(g, ptype_1indexing) = patch%wtlunit(p) * 100._r8 else if (lun%itype(l) == istcrop) then ptype_1indexing = ptype + (1 - cft_lb) subgrid_weights_diagnostics%pct_cft(g, ptype_1indexing) = patch%wtlunit(p) * 100._r8 end if end do end subroutine set_pct_pft_diagnostics end module subgridWeightsMod