module EDCanopyStructureMod ! ===================================================================================== ! Code to determine whether the canopy is closed, and which plants are either in the ! understorey or overstorey. This is obviosuly far too complicated for it's own good ! ===================================================================================== use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : itrue, ifalse use FatesConstantsMod , only : tinyr8 use FatesConstantsMod , only : nearzero use FatesConstantsMod , only : rsnbl_math_prec use FatesGlobals , only : fates_log use EDPftvarcon , only : EDPftvarcon_inst use FatesAllometryMod , only : carea_allom use EDCohortDynamicsMod , only : copy_cohort, terminate_cohorts, fuse_cohorts use EDCohortDynamicsMod , only : InitPRTObject use EDCohortDynamicsMod , only : InitPRTBoundaryConditions use EDCohortDynamicsMod , only : SendCohortToLitter use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : nclmax use EDTypesMod , only : nlevleaf use EDtypesMod , only : AREA use FatesGlobals , only : endrun => fates_endrun use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : numpft use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use EDTypesMod , only : maxCohortsPerPatch use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : fnrt_organ use PRTGenericMod, only : sapw_organ use PRTGenericMod, only : store_organ use PRTGenericMod, only : repro_organ use PRTGenericMod, only : struct_organ use PRTGenericMod, only : SetState ! CIME Globals use shr_log_mod , only : errMsg => shr_log_errMsg implicit none private public :: canopy_structure public :: canopy_spread public :: calc_areaindex public :: canopy_summarization public :: update_hlm_dynamics logical, parameter :: debug=.false. character(len=*), parameter, private :: sourcefile = & __FILE__ real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation ! will attempt to reduce errors ! below this level real(r8), parameter :: area_check_precision = 1.0E-7_r8 ! Area conservation checks must ! be within this absolute tolerance real(r8), parameter :: area_check_rel_precision = 1.0E-4_r8 ! Area conservation checks must ! be within this relative tolerance real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm ! can be roughly considered the same right? ! 10/30/09: Created by Rosie Fisher ! 2017/2018: Modifications and updates by Ryan Knox ! ============================================================================ contains ! ============================================================================ subroutine canopy_structure( currentSite , bc_in ) ! ! !DESCRIPTION: ! create cohort instance ! ! This routine allocates the 'canopy_layer' attribute to each cohort ! All top leaves in the same canopy layer get the same light resources. ! The first canopy layer is the 'canopy' or 'overstorey'. The second is the 'understorey'. ! More than two layers is not permitted at the moment ! Seeds germinating into the 3rd or higher layers are automatically removed. ! ! ------Perfect Plasticity----- ! The idea of these canopy layers derives originally from Purves et al. 2009 ! Their concept is that, given enoughplasticity in canopy position, size, shape and depth ! all of the gound area will be filled perfectly by leaves, and additional leaves will have ! to exist in the understorey. ! Purves et al. use the concept of 'Z*' to assume that the height required to attain a place in the ! canopy is spatially uniform. In this implementation, described in Fisher et al. (2010, New Phyt) we ! extent that concept to assume that position in the canopy has some random element, and that BOTH height ! and chance combine to determine whether trees get into the canopy. ! Thus, when the canopy is closed and there is excess area, some of it must be demoted ! If we demote -all- the trees less than a given height, there is a massive advantage in being the cohort that is ! the biggest when the canopy is closed. ! In this implementation, the amount demoted, ('weight') is a function of the height weighted by the competitive exclusion ! parameter (ED_val_comp_excln). ! Complexity in this routine results from a few things. ! Firstly, the complication of the demotion amount sometimes being larger than the cohort area (for a very small, short cohort) ! Second, occasionaly, disturbance (specifically fire) can cause the canopy layer to become less than closed, ! without changing the area of the patch. If this happens, then some of the plants in the lower layer need to be 'promoted' so ! all of the routine has to happen in both the downwards and upwards directions. ! ! The order of events here is therefore: ! (The entire subroutine has a single outer 'patch' loop. ! Section 1: figure out the total area, and whether there are >1 canopy layers at all. ! ! Sorts out cohorts into canopy and understorey layers... ! ! !USES: use EDParamsMod, only : ED_val_comp_excln use EDTypesMod , only : min_patch_area use FatesInterfaceMod, only : bc_in_type ! ! !ARGUMENTS type(ed_site_type) , intent(inout), target :: currentSite type(bc_in_type), intent(in) :: bc_in ! ! !LOCAL VARIABLES: type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort integer :: i_lyr ! current layer index integer :: z ! Current number of canopy layers. (1= canopy, 2 = understorey) integer :: ipft real(r8) :: arealayer(nclmax+2) ! Amount of plant area currently in each canopy layer integer :: patch_area_counter ! count iterations used to solve canopy areas logical :: area_not_balanced ! logical controlling if the patch layer areas ! have successfully been redistributed integer :: return_code ! math checks on variables will return>0 if problems exist ! We only iterate because of possible imprecisions generated by the cohort ! termination process. These should be super small, so at the most ! try to re-balance 3 times. If that doesn't give layer areas ! within tolerance of canopy area, there is something wrong integer, parameter :: max_patch_iterations = 10 !---------------------------------------------------------------------- currentPatch => currentSite%oldest_patch ! ! zero site-level demotion / promotion tracking info currentSite%demotion_rate(:) = 0._r8 currentSite%promotion_rate(:) = 0._r8 currentSite%demotion_carbonflux = 0._r8 currentSite%promotion_carbonflux = 0._r8 ! ! Section 1: Check total canopy area. ! do while (associated(currentPatch)) ! Patch loop ! ------------------------------------------------------------------------------ ! Perform numerical checks on some cohort and patch structures ! ------------------------------------------------------------------------------ ! canopy layer has a special bounds check currentCohort => currentPatch%tallest do while (associated(currentCohort)) if( currentCohort%canopy_layer < 1 .or. currentCohort%canopy_layer > nclmax+1 ) then write(fates_log(),*) 'lat:',currentSite%lat write(fates_log(),*) 'lon:',currentSite%lon write(fates_log(),*) 'BOGUS CANOPY LAYER: ',currentCohort%canopy_layer call endrun(msg=errMsg(sourcefile, __LINE__)) end if currentCohort => currentCohort%shorter enddo ! Does any layer have excess area in it? Keep going until it does not... patch_area_counter = 0 area_not_balanced = .true. do while(area_not_balanced) ! --------------------------------------------------------------------------- ! Demotion Phase: Identify upper layers that are too full, and demote them to ! the layers below. ! --------------------------------------------------------------------------- ! Its possible that before we even enter this scheme ! some cohort numbers are very low. Terminate them. call terminate_cohorts(currentSite, currentPatch, 1, 12) ! Calculate how many layers we have in this canopy ! This also checks the understory to see if its crown ! area is large enough to warrant a temporary sub-understory layer z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.) do i_lyr = 1,z ! Loop around the currently occupied canopy layers. call DemoteFromLayer(currentSite, currentPatch, i_lyr) end do ! After demotions, we may then again have cohorts that are very very ! very sparse, remove them call terminate_cohorts(currentSite, currentPatch, 1,13) call fuse_cohorts(currentSite, currentPatch, bc_in) ! Remove cohorts for various other reasons call terminate_cohorts(currentSite, currentPatch, 2,13) ! --------------------------------------------------------------------------------------- ! Promotion Phase: Identify if any upper-layers are underful and layers below them ! have cohorts that can be split and promoted to the layer above. ! --------------------------------------------------------------------------------------- ! Re-calculate Number of layers without the false substory z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.) ! We only promote if we have at least two layers if (z>1) then do i_lyr=1,z-1 call PromoteIntoLayer(currentSite, currentPatch, i_lyr) end do ! Remove cohorts that are incredibly sparse call terminate_cohorts(currentSite, currentPatch, 1,14) call fuse_cohorts(currentSite, currentPatch, bc_in) ! Remove cohorts for various other reasons call terminate_cohorts(currentSite, currentPatch, 2,14) end if ! --------------------------------------------------------------------------------------- ! Check on Layer Area (if the layer differences are not small ! Continue trying to demote/promote. Its possible on the first pass through, ! that cohort fusion has nudged the areas a little bit. ! --------------------------------------------------------------------------------------- z = NumPotentialCanopyLayers(currentPatch,currentSite%spread,include_substory=.false.) area_not_balanced = .false. do i_lyr = 1,z call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer(i_lyr)) if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. & ((arealayer(i_lyr)-currentPatch%area) > area_check_precision ) ) then area_not_balanced = .true. endif enddo ! --------------------------------------------------------------------------------------- ! Gracefully exit if too many iterations have gone by ! --------------------------------------------------------------------------------------- patch_area_counter = patch_area_counter + 1 if(patch_area_counter > max_patch_iterations .and. area_not_balanced) then write(fates_log(),*) 'PATCH AREA CHECK NOT CLOSING' write(fates_log(),*) 'patch area:',currentpatch%area do i_lyr = 1,z write(fates_log(),*) 'layer: ',i_lyr,' area: ',arealayer(i_lyr) write(fates_log(),*) 'rel error: ',(arealayer(i_lyr)-currentPatch%area)/currentPatch%area write(fates_log(),*) 'abs error: ',arealayer(i_lyr)-currentPatch%area enddo write(fates_log(),*) 'lat:',currentSite%lat write(fates_log(),*) 'lon:',currentSite%lon write(fates_log(),*) 'spread:',currentSite%spread currentCohort => currentPatch%tallest do while (associated(currentCohort)) write(fates_log(),*) 'coh ilayer:',currentCohort%canopy_layer write(fates_log(),*) 'coh dbh:',currentCohort%dbh write(fates_log(),*) 'coh pft:',currentCohort%pft write(fates_log(),*) 'coh n:',currentCohort%n write(fates_log(),*) 'coh carea:',currentCohort%c_area ipft=currentCohort%pft write(fates_log(),*) 'maxh:',EDPftvarcon_inst%allom_dbh_maxheight(ipft) write(fates_log(),*) 'lmode: ',EDPftvarcon_inst%allom_lmode(ipft) write(fates_log(),*) 'd2bl2: ',EDPftvarcon_inst%allom_d2bl2(ipft) write(fates_log(),*) 'd2bl_ediff: ',EDPftvarcon_inst%allom_blca_expnt_diff(ipft) write(fates_log(),*) 'd2ca_min: ',EDPftvarcon_inst%allom_d2ca_coefficient_min(ipft) write(fates_log(),*) 'd2ca_max: ',EDPftvarcon_inst%allom_d2ca_coefficient_max(ipft) currentCohort => currentCohort%shorter enddo call endrun(msg=errMsg(sourcefile, __LINE__)) end if enddo ! do while(area_not_balanced) ! Set current canopy layer occupancy indicator. currentPatch%NCL_p = min(nclmax,z) ! ------------------------------------------------------------------------------------------- ! if we are using "strict PPA", then calculate a z_star value as ! the height of the smallest tree in the canopy ! loop from top to bottom and locate the shortest cohort in level 1 whose shorter ! neighbor is in level 2 set zstar as the ehight of that shortest level 1 cohort ! ------------------------------------------------------------------------------------------- if ( ED_val_comp_excln .lt. 0.0_r8) then currentPatch%zstar = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer .eq. 2)then if (associated(currentCohort%taller)) then if (currentCohort%taller%canopy_layer .eq. 1 ) then currentPatch%zstar = currentCohort%taller%hite endif endif endif currentCohort => currentCohort%shorter enddo endif currentPatch => currentPatch%younger enddo !patch return end subroutine canopy_structure ! ============================================================================================== subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) use EDParamsMod, only : ED_val_comp_excln use SFParamsMod, only : SF_val_CWD_frac ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), intent(inout), target :: currentPatch integer, intent(in) :: i_lyr ! Current canopy layer of interest ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: currentCohort type(ed_cohort_type), pointer :: copyc type(ed_cohort_type), pointer :: nextc ! The next cohort in line integer :: i_cwd ! Index for CWD pool real(r8) :: cc_loss ! cohort crown area loss in demotion (m2) real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] real(r8) :: scale_factor ! for prob. exclusion - scales weight to a fraction real(r8) :: scale_factor_min ! "" minimum before exeedance of 1 real(r8) :: scale_factor_res ! "" applied to residual areas real(r8) :: area_res ! residual area to demote after weakest cohort hits max real(r8) :: newarea real(r8) :: demote_area real(r8) :: sumweights real(r8) :: sumequal ! for rank-ordered same-size cohorts ! this tallies their excluded area real(r8) :: arealayer ! the area of the current canopy layer logical :: tied_size_with_neighbors real(r8) :: total_crownarea_of_tied_cohorts ! First, determine how much total canopy area we have in this layer call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer) demote_area = arealayer - currentPatch%area if ( demote_area > area_target_precision ) then ! Is this layer currently over-occupied? ! In that case, we need to work out which cohorts to demote. ! We go in order from shortest to tallest for ranked demotion sumweights = 0.0_r8 currentCohort => currentPatch%shortest do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) if(debug) then if(currentCohort%c_area<0._r8)then write(fates_log(),*) 'negative c_area stage 1d: ',currentCohort%dbh,i_lyr,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if if( currentCohort%canopy_layer == i_lyr)then if (ED_val_comp_excln .ge. 0.0_r8 ) then ! ---------------------------------------------------------- ! Stochastic method. ! Weight cohort demotion by inverse size to a constant power. ! In this hypothesis, it is assumed that even the tallest ! cohorts have a chance (although smaller) of being forced ! to the understory. ! ---------------------------------------------------------- currentCohort%excl_weight = 1._r8 / (currentCohort%hite**ED_val_comp_excln) sumweights = sumweights + currentCohort%excl_weight else ! ----------------------------------------------------------- ! Rank ordered deterministic method ! ----------------------------------------------------------- ! If there are cohorts that have the exact same height (which is possible, really) ! we don't want to unilaterally promote/demote one before the others. ! So we <>mote them as a unit ! now we need to go through and figure out how many equal-size cohorts there are. ! then we need to go through, add up the collective crown areas of all equal-sized ! and equal-canopy-layer cohorts, ! and then demote from each as if they were a single group total_crownarea_of_tied_cohorts = currentCohort%c_area tied_size_with_neighbors = .false. nextc => currentCohort%taller do while (associated(nextc)) if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then tied_size_with_neighbors = .true. total_crownarea_of_tied_cohorts = & total_crownarea_of_tied_cohorts + nextc%c_area end if else exit endif nextc => nextc%taller end do if ( tied_size_with_neighbors ) then currentCohort%excl_weight = & max(0.0_r8,min(currentCohort%c_area, & (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) )) sumequal = currentCohort%excl_weight nextc => currentCohort%taller do while (associated(nextc)) if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then ! now we know the total crown area of all equal-sized, ! equal-canopy-layer cohorts nextc%excl_weight = & max(0.0_r8,min(nextc%c_area, & (nextc%c_area/total_crownarea_of_tied_cohorts) * & (demote_area - sumweights) )) sumequal = sumequal + nextc%excl_weight end if else exit endif nextc => nextc%taller end do ! Update the current cohort pointer to the last similar cohort ! Its ok if this is not in the right layer if(associated(nextc))then currentCohort => nextc%shorter else currentCohort => currentPatch%tallest end if sumweights = sumweights + sumequal else currentCohort%excl_weight = & max(min(currentCohort%c_area, demote_area - sumweights ), 0._r8) sumweights = sumweights + currentCohort%excl_weight end if endif endif currentCohort => currentCohort%taller enddo ! If this is probabalistic demotion, we need to do a round of normalization. ! And then a few rounds where we pre-calculate the demotion areas ! and adjust things if the demoted area wants to be greater than ! what is available. The math is too hard to explain here, see ! the tech note section on promotion/demotion. if (ED_val_comp_excln .ge. 0.0_r8 ) then scale_factor_min = 1.e10_r8 scale_factor = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then currentCohort%excl_weight = currentCohort%excl_weight/sumweights if( 1._r8/currentCohort%excl_weight < scale_factor_min ) & scale_factor_min = 1._r8/currentCohort%excl_weight scale_factor = scale_factor + currentCohort%excl_weight * currentCohort%c_area endif currentCohort => currentCohort%shorter enddo ! This is the factor by which we need to multiply ! the demotion probabilities, so the sum result equals ! the total amount to demote scale_factor = demote_area/scale_factor if(scale_factor <= scale_factor_min) then ! Trivial case, all of the demotion fractions are less than 1. currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then currentCohort%excl_weight = currentCohort%c_area * currentCohort%excl_weight * scale_factor if(debug) then if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%excl_weight < 0._r8) ) then write(fates_log(),*) 'exclusion area too big (1)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'dbh: ',currentCohort%dbh write(fates_log(),*) 'n: ',currentCohort%n write(fates_log(),*) 'spread: ',currentSite%spread write(fates_log(),*) 'pft: ',currentCohort%pft write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if endif currentCohort => currentCohort%shorter enddo else ! Non-trivial case, at least 1 cohort's demotion ! rate would exceed its area, given the trivial scale factor area_res = 0._r8 scale_factor_res = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then area_res = area_res + & currentCohort%c_area * currentCohort%excl_weight * & scale_factor_min scale_factor_res = scale_factor_res + & currentCohort%c_area * & (1._r8 - (currentCohort%excl_weight * scale_factor_min)) endif currentCohort => currentCohort%shorter enddo area_res = demote_area - area_res scale_factor_res = area_res / scale_factor_res currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then currentCohort%excl_weight = currentCohort%c_area * & (currentCohort%excl_weight * scale_factor_min + & (1._r8 - (currentCohort%excl_weight*scale_factor_min) ) * scale_factor_res) if(debug)then if((currentCohort%excl_weight > & (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%excl_weight < 0._r8) ) then write(fates_log(),*) 'exclusion area error (2)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'currentCohort%excl_weight: ', & currentCohort%excl_weight write(fates_log(),*) 'excess: ', & currentCohort%excl_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if endif currentCohort => currentCohort%shorter enddo end if end if ! perform a check and see if the demotions meet the demand sumweights = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then sumweights = sumweights + currentCohort%excl_weight end if currentCohort => currentCohort%shorter end do if (abs(sumweights - demote_area) > area_check_precision ) then write(fates_log(),*) 'demotions dont add up' write(fates_log(),*) 'sum demotions: ',sumweights write(fates_log(),*) 'area needed to be demoted: ',demote_area write(fates_log(),*) 'excess: ',sumweights - demote_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Weights have been calculated. Now move them to the lower layer currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr )then cc_loss = currentCohort%excl_weight leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) if ( (cc_loss-currentCohort%c_area) > -nearzero .and. & (cc_loss-currentCohort%c_area) < area_target_precision ) then ! If the whole cohort is being demoted, just change its ! layer index currentCohort%canopy_layer = i_lyr+1 ! keep track of number and biomass of demoted cohort currentSite%demotion_rate(currentCohort%size_class) = & currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + & (leaf_c + store_c + fnrt_c + sapw_c + struct_c) * currentCohort%n elseif( (cc_loss < currentCohort%c_area) .and. & (cc_loss > area_target_precision) ) then ! If only part of the cohort is demoted ! then it must be split (little more complicated) ! Make a copy of the current cohort. The copy and the original ! conserve total number density of the original. The copy ! remains in the upper-story. The original is the one ! demoted to the understory allocate(copyc) ! Initialize the PARTEH object and point to the ! correct boundary condition fields copyc%prt => null() call InitPRTObject(copyc%prt) call InitPRTBoundaryConditions(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(currentSite,copyc) endif call copy_cohort(currentCohort, copyc) newarea = currentCohort%c_area - cc_loss copyc%n = currentCohort%n*newarea/currentCohort%c_area currentCohort%n = currentCohort%n - copyc%n copyc%canopy_layer = i_lyr !the taller cohort is the copy ! Demote the current cohort to the understory. currentCohort%canopy_layer = i_lyr + 1 ! keep track of number and biomass of demoted cohort currentSite%demotion_rate(currentCohort%size_class) = & currentSite%demotion_rate(currentCohort%size_class) + currentCohort%n currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + & (leaf_c + store_c + fnrt_c + sapw_c + struct_c) * currentCohort%n call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,copyc%c_area) call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) !----------- Insert copy into linked list ------------------------! copyc%shorter => currentCohort if(associated(currentCohort%taller))then copyc%taller => currentCohort%taller currentCohort%taller%shorter => copyc else currentPatch%tallest => copyc copyc%taller => null() endif currentCohort%taller => copyc elseif(cc_loss > currentCohort%c_area)then write(fates_log(),*) 'more area than the cohort has is being demoted' write(fates_log(),*) 'loss:',cc_loss write(fates_log(),*) 'existing area:',currentCohort%c_area write(fates_log(),*) 'excess: ',cc_loss - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! kill the ones which go into canopy layers that are not allowed if(currentCohort%canopy_layer>nclmax )then ! put the litter from the terminated cohorts ! straight into the fragmenting pools call SendCohortToLitter(currentSite,currentPatch, & currentCohort,currentCohort%n) currentCohort%n = 0.0_r8 currentCohort%c_area = 0.0_r8 currentCohort%canopy_layer = i_lyr end if call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) endif !canopy layer = i_ly currentCohort => currentCohort%shorter enddo !currentCohort ! Update the area calculations of the current layer ! And the layer below that may or may not had recieved ! Demotions call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer) if ( (abs(arealayer - currentPatch%area)/arealayer > area_check_rel_precision ) .or. & (abs(arealayer - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'demotion did not trim area within tolerance' write(fates_log(),*) 'arealayer:',arealayer write(fates_log(),*) 'patch%area:',currentPatch%area write(fates_log(),*) 'ilayer: ',i_lyr write(fates_log(),*) 'bias:',arealayer - currentPatch%area write(fates_log(),*) 'rel bias:',(arealayer - currentPatch%area)/arealayer write(fates_log(),*) 'demote_area:',demote_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if return end subroutine DemoteFromLayer ! ============================================================================================== subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! ------------------------------------------------------------------------------------------- ! Check whether the intended 'full' layers are actually filling all the space. ! If not, promote some fraction of cohorts upwards. ! THIS SECTION MIGHT BE TRIGGERED BY A FIRE OR MORTALITY EVENT, FOLLOWED BY A PATCH FUSION, ! SO THE TOP LAYER IS NO LONGER FULL. ! ------------------------------------------------------------------------------------------- use EDParamsMod, only : ED_val_comp_excln ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), intent(inout), target :: currentPatch integer, intent(in) :: i_lyr ! Current canopy layer of interest ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: currentCohort type(ed_cohort_type), pointer :: copyc type(ed_cohort_type), pointer :: nextc ! the next cohort, or used for looping ! cohorts against the current real(r8) :: scale_factor ! for prob. exclusion - scales weight to a fraction real(r8) :: scale_factor_min ! "" minimum before exeedance of 1 real(r8) :: scale_factor_res ! "" applied to residual areas real(r8) :: area_res ! residual area to demote after weakest cohort hits max real(r8) :: promote_area real(r8) :: newarea real(r8) :: sumweights real(r8) :: sumequal ! for tied cohorts, the sum of weights in ! their group real(r8) :: cc_gain ! cohort crown area gain in promotion (m2) real(r8) :: arealayer_current ! area (m2) of the current canopy layer real(r8) :: arealayer_below ! area (m2) of the layer below the current layer real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] logical :: tied_size_with_neighbors real(r8) :: total_crownarea_of_tied_cohorts call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr+1,arealayer_below) ! how much do we need to gain? promote_area = currentPatch%area - arealayer_current if( promote_area > area_target_precision ) then if(arealayer_below <= promote_area ) then ! --------------------------------------------------------------------------- ! Promote all cohorts from layer below if that whole layer has area smaller ! than the tolerance on the gains needed into current layer ! --------------------------------------------------------------------------- currentCohort => currentPatch%tallest do while (associated(currentCohort)) !look at the cohorts in the canopy layer below... if(currentCohort%canopy_layer == i_lyr+1)then leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) currentCohort%canopy_layer = i_lyr call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(currentCohort%size_class) = & currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + & (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * currentCohort%n endif currentCohort => currentCohort%shorter enddo else ! --------------------------------------------------------------------------- ! This is the non-trivial case where the lower layer can accomodate ! more than what is necessary. ! --------------------------------------------------------------------------- ! figure out with what weighting we need to promote cohorts. ! This is the opposite of the demotion weighting... sumweights = 0.0_r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) if(currentCohort%canopy_layer == i_lyr+1)then !look at the cohorts in the canopy layer below... if (ED_val_comp_excln .ge. 0.0_r8 ) then ! ------------------------------------------------------------------ ! Stochastic case, as above (in demotion portion of code) ! ------------------------------------------------------------------ currentCohort%prom_weight = currentCohort%hite**ED_val_comp_excln sumweights = sumweights + currentCohort%prom_weight else ! ------------------------------------------------------------------ ! Rank ordered deterministic method ! If there are cohorts that have the exact same height (which is possible, really) ! we don't want to unilaterally promote/demote one before the others. ! So we <>mote them as a unit ! now we need to go through and figure out how many equal-size cohorts there are. ! then we need to go through, add up the collective crown areas of all equal-sized ! and equal-canopy-layer cohorts, ! and then demote from each as if they were a single group ! ------------------------------------------------------------------ total_crownarea_of_tied_cohorts = currentCohort%c_area tied_size_with_neighbors = .false. nextc => currentCohort%shorter do while (associated(nextc)) if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if( nextc%canopy_layer .eq. currentCohort%canopy_layer ) then tied_size_with_neighbors = .true. total_crownarea_of_tied_cohorts = & total_crownarea_of_tied_cohorts + nextc%c_area end if else exit endif nextc => nextc%shorter end do if ( tied_size_with_neighbors ) then currentCohort%prom_weight = & max(0.0_r8,min(currentCohort%c_area, & (currentCohort%c_area/total_crownarea_of_tied_cohorts) * & (promote_area - sumweights) )) sumequal = currentCohort%prom_weight nextc => currentCohort%shorter do while (associated(nextc)) if ( abs(nextc%hite - currentCohort%hite) < similar_height_tol ) then if (nextc%canopy_layer .eq. currentCohort%canopy_layer ) then ! now we know the total crown area of all equal-sized, ! equal-canopy-layer cohorts nextc%prom_weight = & max(0.0_r8,min(nextc%c_area, & (nextc%c_area/total_crownarea_of_tied_cohorts) * & (promote_area - sumweights) )) sumequal = sumequal + nextc%prom_weight end if else exit endif nextc => nextc%shorter end do ! Update the current cohort pointer to the last similar cohort ! Its ok if this is not in the right layer if(associated(nextc))then currentCohort => nextc%taller else currentCohort => currentPatch%shortest end if sumweights = sumweights + sumequal else currentCohort%prom_weight = & max(min(currentCohort%c_area, promote_area - sumweights ), 0._r8) sumweights = sumweights + currentCohort%prom_weight end if endif endif currentCohort => currentCohort%shorter enddo !currentCohort ! If this is probabalistic promotion, we need to do a round of normalization. ! And then a few rounds where we pre-calculate the promotion areas ! and adjust things if the promoted area wants to be greater than ! what is available. if (ED_val_comp_excln .ge. 0.0_r8 ) then scale_factor_min = 1.e10_r8 scale_factor = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1) ) then currentCohort%prom_weight = currentCohort%prom_weight/sumweights if( 1._r8/currentCohort%prom_weight < scale_factor_min ) & scale_factor_min = 1._r8/currentCohort%prom_weight scale_factor = scale_factor + currentCohort%prom_weight * currentCohort%c_area endif currentCohort => currentCohort%shorter enddo ! This is the factor by which we need to multiply ! the demotion probabilities, so the sum result equals ! the total amount to demote scale_factor = promote_area/scale_factor if(scale_factor <= scale_factor_min) then ! Trivial case, all of the demotion fractions ! are less than 1. currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1) ) then currentCohort%prom_weight = currentCohort%c_area * & currentCohort%prom_weight * scale_factor if(debug)then if((currentCohort%prom_weight > & (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%prom_weight < 0._r8) ) then write(fates_log(),*) 'promotion area too big (1)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'currentCohort%prom_weight: ', & currentCohort%prom_weight write(fates_log(),*) 'excess: ', & currentCohort%prom_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if endif currentCohort => currentCohort%shorter enddo else ! Non-trivial case, at least 1 cohort's promotion ! rate would exceed its area, given the trivial scale factor area_res = 0._r8 scale_factor_res = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1) ) then area_res = area_res + & currentCohort%c_area*currentCohort%prom_weight*scale_factor_min scale_factor_res = scale_factor_res + & currentCohort%c_area * & (1._r8 - (currentCohort%prom_weight * scale_factor_min)) endif currentCohort => currentCohort%shorter enddo area_res = promote_area - area_res scale_factor_res = area_res / scale_factor_res currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1)) then currentCohort%prom_weight = currentCohort%c_area * & (currentCohort%prom_weight * scale_factor_min + & (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * & scale_factor_res) if(debug)then if((currentCohort%prom_weight > & (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%prom_weight < 0._r8) ) then write(fates_log(),*) 'promotion area error (2)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'currentCohort%prom_weight: ', & currentCohort%prom_weight write(fates_log(),*) 'excess: ', & currentCohort%prom_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if endif currentCohort => currentCohort%shorter enddo end if end if ! lets perform a check and see if the promotions meet the demand sumweights = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1)) then sumweights = sumweights + currentCohort%prom_weight end if currentCohort => currentCohort%shorter end do if(debug)then if (abs(sumweights - promote_area) > area_check_precision ) then write(fates_log(),*) 'promotions dont add up' write(fates_log(),*) 'sum promotions: ',sumweights write(fates_log(),*) 'area needed to be promoted: ',promote_area write(fates_log(),*) 'excess: ',sumweights - promote_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if currentCohort => currentPatch%tallest do while (associated(currentCohort)) !All the trees in this layer need to promote some area upwards... if( (currentCohort%canopy_layer == i_lyr+1) ) then cc_gain = currentCohort%prom_weight leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) if ( (cc_gain-currentCohort%c_area) > -nearzero .and. & (cc_gain-currentCohort%c_area) < area_target_precision ) then currentCohort%canopy_layer = i_lyr ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(currentCohort%size_class) = & currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + & (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * currentCohort%n elseif ( (cc_gain < currentCohort%c_area) .and. & (cc_gain > area_target_precision) ) then allocate(copyc) ! Initialize the PARTEH object and point to the ! correct boundary condition fields copyc%prt => null() call InitPRTObject(copyc%prt) call InitPRTBoundaryConditions(copyc) if( hlm_use_planthydro.eq.itrue ) then call InitHydrCohort(CurrentSite,copyc) endif call copy_cohort(currentCohort, copyc) !makes an identical copy... newarea = currentCohort%c_area - cc_gain !new area of existing cohort call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) ! number of individuals in promoted cohort. copyc%n = currentCohort%n*cc_gain/currentCohort%c_area ! number of individuals in cohort remaining in understorey currentCohort%n = currentCohort%n - copyc%n currentCohort%canopy_layer = i_lyr + 1 ! keep current cohort in the understory. copyc%canopy_layer = i_lyr ! promote copy to the higher canopy layer. ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(copyc%size_class) = & currentSite%promotion_rate(copyc%size_class) + copyc%n currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + & (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * copyc%n call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) call carea_allom(copyc%dbh,copyc%n,currentSite%spread,copyc%pft,copyc%c_area) !----------- Insert copy into linked list ------------------------! copyc%shorter => currentCohort if(associated(currentCohort%taller))then copyc%taller => currentCohort%taller currentCohort%taller%shorter => copyc else currentPatch%tallest => copyc copyc%taller => null() endif currentCohort%taller => copyc elseif(cc_gain > currentCohort%c_area)then write(fates_log(),*) 'more area than the cohort has is being promoted' write(fates_log(),*) 'loss:',cc_gain write(fates_log(),*) 'existing area:',currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) endif endif ! if(currentCohort%canopy_layer == i_lyr+1) then currentCohort => currentCohort%shorter enddo !currentCohort call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) if ((abs(arealayer_current - currentPatch%area)/arealayer_current > & area_check_rel_precision ) .or. & (abs(arealayer_current - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'promotion did not bring area within tolerance' write(fates_log(),*) 'arealayer:',arealayer_current write(fates_log(),*) 'patch%area:',currentPatch%area call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if end if return end subroutine PromoteIntoLayer ! ============================================================================ subroutine canopy_spread( currentSite ) ! ! !DESCRIPTION: ! Calculates the spatial spread of tree canopies based on canopy closure. ! ! !USES: use EDTypesMod , only : AREA use EDParamsMod, only : ED_val_canopy_closure_thresh ! ! !ARGUMENTS type (ed_site_type), intent(inout), target :: currentSite ! ! !LOCAL VARIABLES: type (ed_cohort_type), pointer :: currentCohort type (ed_patch_type) , pointer :: currentPatch real(r8) :: sitelevel_canopyarea ! Amount of canopy in top layer at the site level real(r8) :: inc ! Arbitrary daily incremental change in canopy area integer :: z !---------------------------------------------------------------------- inc = 0.05_r8 currentPatch => currentSite%oldest_patch sitelevel_canopyarea = 0.0_r8 do while (associated(currentPatch)) !calculate canopy area in each patch... currentCohort => currentPatch%tallest do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) if( (EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 ) .and. & (currentCohort%canopy_layer .eq. 1 ) ) then sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area endif currentCohort => currentCohort%shorter enddo currentPatch => currentPatch%younger enddo !currentPatch ! If the canopy area is approaching closure, ! squash the tree canopies and make them taller and thinner if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then currentSite%spread = currentSite%spread - inc else currentSite%spread = currentSite%spread + inc endif ! put within bounds to make sure it stays between 0 and 1 currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8) end subroutine canopy_spread ! ===================================================================================== subroutine canopy_summarization( nsites, sites, bc_in ) ! ---------------------------------------------------------------------------------- ! Much of this routine was once ed_clm_link minus all the IO and history stuff ! --------------------------------------------------------------------------------- use FatesInterfaceMod , only : bc_in_type use EDPatchDynamicsMod , only : set_patchno use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst ! !ARGUMENTS integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) type(bc_in_type) , intent(in) :: bc_in(nsites) ! ! !LOCAL VARIABLES: type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type) , pointer :: currentCohort integer :: s integer :: ft ! plant functional type integer :: ifp integer :: patchn ! identification number for each patch. real(r8) :: canopy_leaf_area ! total amount of leaf area in the vegetated area. m2. real(r8) :: leaf_c ! leaf carbon [kg] real(r8) :: fnrt_c ! fineroot carbon [kg] real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] !---------------------------------------------------------------------- if ( debug ) then write(fates_log(),*) 'in canopy_summarization' endif do s = 1,nsites ! -------------------------------------------------------------------------------- ! Set the patch indices (this is usefull mostly for communicating with a host or ! driving model. Loops through all patches and sets cpatch%patchno to the integer ! order of oldest to youngest where the oldest is 1. ! -------------------------------------------------------------------------------- call set_patchno( sites(s) ) currentPatch => sites(s)%oldest_patch do while(associated(currentPatch)) !zero cohort-summed variables. currentPatch%total_canopy_area = 0.0_r8 currentPatch%total_tree_area = 0.0_r8 canopy_leaf_area = 0.0_r8 !update cohort quantitie s currentCohort => currentPatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) ! Update the cohort's index within the size bin classes ! Update the cohort's index within the SCPF classification system call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,& currentCohort%pft,currentCohort%c_area) currentCohort%treelai = tree_lai(leaf_c, & currentCohort%pft, currentCohort%c_area, currentCohort%n, & currentCohort%canopy_layer, currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area if(currentCohort%canopy_layer==1)then currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area if(EDPftvarcon_inst%woody(ft)==1)then currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area endif endif ! Check for erroneous zero values. if(currentCohort%dbh <= 0._r8 .or. currentCohort%n == 0._r8)then write(fates_log(),*) 'FATES: dbh or n is zero in canopy_summarization', & currentCohort%dbh,currentCohort%n call endrun(msg=errMsg(sourcefile, __LINE__)) endif if(currentCohort%pft == 0.or.currentCohort%canopy_trim <= 0._r8)then write(fates_log(),*) 'FATES: PFT or trim is zero in canopy_summarization', & currentCohort%pft,currentCohort%canopy_trim call endrun(msg=errMsg(sourcefile, __LINE__)) endif if( (sapw_c + leaf_c + fnrt_c) <= 0._r8)then write(fates_log(),*) 'FATES: alive biomass is zero in canopy_summarization', & sapw_c + leaf_c + fnrt_c call endrun(msg=errMsg(sourcefile, __LINE__)) endif currentCohort => currentCohort%taller enddo ! ends 'do while(associated(currentCohort)) if ( currentPatch%total_canopy_area>currentPatch%area ) then if ( currentPatch%total_canopy_area-currentPatch%area > 0.001_r8 ) then write(fates_log(),*) 'FATES: canopy area bigger than area', & currentPatch%total_canopy_area ,currentPatch%area call endrun(msg=errMsg(sourcefile, __LINE__)) end if currentPatch%total_canopy_area = currentPatch%area endif currentPatch => currentPatch%younger end do !patch loop call leaf_area_profile(sites(s),bc_in(s)%snow_depth_si,bc_in(s)%frac_sno_eff_si) end do ! site loop return end subroutine canopy_summarization ! ===================================================================================== subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) ! ----------------------------------------------------------------------------------- ! This subroutine calculates how leaf and stem areas are distributed ! in vertical and horizontal space. ! ! The following cohort level diagnostics are updated here: ! ! currentCohort%treelai ! LAI per unit crown area (m2/m2) ! currentCohort%treesai ! SAI per unit crown area (m2/m2) ! currentCohort%lai ! LAI per unit canopy area (m2/m2) ! currentCohort%sai ! SAI per unit canopy area (m2/m2) ! currentCohort%NV ! The number of discrete vegetation ! ! layers needed to describe this crown ! ! The following patch level diagnostics are updated here: ! ! currentPatch%canopy_layer_tlai(cl) ! total leaf area index of canopy layer ! currentPatch%ncan(cl,ft) ! number of vegetation layers needed ! ! in this patch's pft/canopy-layer ! currentPatch%nrad(cl,ft) ! same as ncan, but does not include ! ! layers occluded by snow ! ! CURRENTLY SAME AS NCAN ! currentPatch%canopy_mask(cl,ft) ! are there canopy elements in this pft-layer? ! ! (This is redundant with nrad though...) ! currentPatch%tlai_profile(cl,ft,iv) ! m2 of leaves per m2 of the PFT's footprint ! currentPatch%elai_profile(cl,ft,iv) ! non-snow covered m2 of leaves per m2 of PFT footprint ! currentPatch%tsai_profile(cl,ft,iv) ! m2 of stems per m2 of PFT footprint ! currentPatch%esai_profile(cl,ft,iv) ! non-snow covered m2 of stems per m2 of PFT footprint ! currentPatch%canopy_area_profile(cl,ft,iv) ! Fractional area of leaf layer ! ! relative to vegetated area ! currentPatch%layer_height_profile(cl,ft,iv) ! Elevation of layer in m ! ! ----------------------------------------------------------------------------------- ! !USES: use EDtypesMod , only : area, dinc_ed, hitemax, n_hite_bins ! ! !ARGUMENTS type(ed_site_type) , intent(inout) :: currentSite real(r8) , intent(in) :: snow_depth_si real(r8) , intent(in) :: frac_sno_eff_si ! ! !LOCAL VARIABLES: type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type) , pointer :: currentCohort real(r8) :: remainder !Thickness of layer at bottom of canopy. real(r8) :: fleaf ! fraction of cohort incepting area that is leaves. integer :: ft ! Plant functional type index. integer :: iv ! Vertical leaf layer index integer :: cl ! Canopy layer index real(r8) :: fraction_exposed ! how much of this layer is not covered by snow? real(r8) :: layer_top_hite ! notional top height of this canopy layer (m) real(r8) :: layer_bottom_hite ! notional bottom height of this canopy layer (m) integer :: smooth_leaf_distribution ! is the leaf distribution this option (1) or not (0) real(r8) :: frac_canopy(N_HITE_BINS) ! amount of canopy in each height class real(r8) :: patch_lai ! LAI summed over the patch in m2/m2 of canopy area real(r8) :: minh(N_HITE_BINS) ! minimum height in height class (m) real(r8) :: maxh(N_HITE_BINS) ! maximum height in height class (m) real(r8) :: dh ! vertical detph of height class (m) real(r8) :: min_chite ! bottom of cohort canopy (m) real(r8) :: max_chite ! top of cohort canopy (m) real(r8) :: lai ! summed lai for checking m2 m-2 real(r8) :: snow_depth_avg ! avg snow over whole site real(r8) :: leaf_c ! leaf carbon [kg] !---------------------------------------------------------------------- smooth_leaf_distribution = 0 ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft ! We assume that each point in the canopy recieved the light attenuated by the average ! leaf area index above it, irrespective of PFT identity... ! Each leaf is defined by how deep in the canopy it is, in terms of LAI units. (FIX(RF,032414), GB) currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) ! -------------------------------------------------------------------------------- ! Calculate tree and canopy areas. ! calculate tree lai and sai. ! -------------------------------------------------------------------------------- currentPatch%canopy_layer_tlai(:) = 0._r8 currentPatch%ncan(:,:) = 0 currentPatch%nrad(:,:) = 0 patch_lai = 0._r8 currentPatch%tlai_profile(:,:,:) = 0._r8 currentPatch%tsai_profile(:,:,:) = 0._r8 currentPatch%elai_profile(:,:,:) = 0._r8 currentPatch%esai_profile(:,:,:) = 0._r8 currentPatch%layer_height_profile(:,:,:) = 0._r8 currentPatch%canopy_area_profile(:,:,:) = 0._r8 currentPatch%canopy_mask(:,:) = 0 ! ------------------------------------------------------------------------------ ! It is remotely possible that in deserts we will not have any canopy ! area, ie not plants at all... ! ------------------------------------------------------------------------------ if (currentPatch%total_canopy_area > nearzero ) then currentCohort => currentPatch%tallest do while(associated(currentCohort)) ft = currentCohort%pft cl = currentCohort%canopy_layer ! Calculate LAI of layers above ! Note that the canopy_layer_lai is also calculated in this loop ! but since we go top down in terms of plant size, we should be okay leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, & currentCohort%n, currentCohort%canopy_layer, & currentPatch%canopy_layer_tlai,currentCohort%vcmax25top ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & currentPatch%canopy_layer_tlai, currentCohort%treelai , & currentCohort%vcmax25top,4) currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area ! Number of actual vegetation layers in this cohort's crown currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV) patch_lai = patch_lai + currentCohort%lai currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + currentCohort%lai currentCohort => currentCohort%shorter enddo !currentCohort if(smooth_leaf_distribution == 1)then ! ----------------------------------------------------------------------------- ! we are going to ignore the concept of canopy layers, and put all of the leaf ! area into height banded bins. using the same domains as we had before, except ! that CL always = 1 ! ----------------------------------------------------------------------------- ! this is a crude way of dividing up the bins. Should it be a function of actual maximum height? dh = 1.0_r8*(HITEMAX/N_HITE_BINS) do iv = 1,N_HITE_BINS if (iv == 1) then minh(iv) = 0.0_r8 maxh(iv) = dh else minh(iv) = (iv-1)*dh maxh(iv) = (iv)*dh endif enddo currentCohort => currentPatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft min_chite = currentCohort%hite - currentCohort%hite * EDPftvarcon_inst%crown(ft) max_chite = currentCohort%hite do iv = 1,N_HITE_BINS frac_canopy(iv) = 0.0_r8 ! this layer is in the middle of the canopy if(max_chite > maxh(iv).and.min_chite < minh(iv))then frac_canopy(iv)= min(1.0_r8,dh / (currentCohort%hite*EDPftvarcon_inst%crown(ft))) ! this is the layer with the bottom of the canopy in it. elseif(min_chite < maxh(iv).and.min_chite > minh(iv).and.max_chite > maxh(iv))then frac_canopy(iv) = (maxh(iv) -min_chite ) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) ! this is the layer with the top of the canopy in it. elseif(max_chite > minh(iv).and.max_chite < maxh(iv).and.min_chite < minh(iv))then frac_canopy(iv) = (max_chite - minh(iv)) / (currentCohort%hite*EDPftvarcon_inst%crown(ft)) elseif(max_chite < maxh(iv).and.min_chite > minh(iv))then !the whole cohort is within this layer. frac_canopy(iv) = 1.0_r8 endif ! no m2 of leaf per m2 of ground in each height class currentPatch%tlai_profile(1,ft,iv) = currentPatch%tlai_profile(1,ft,iv) + frac_canopy(iv) * & currentCohort%lai currentPatch%tsai_profile(1,ft,iv) = currentPatch%tsai_profile(1,ft,iv) + frac_canopy(iv) * & currentCohort%sai !snow burial !write(fates_log(), *) 'calc snow' snow_depth_avg = snow_depth_si * frac_sno_eff_si if(snow_depth_avg > maxh(iv))then fraction_exposed = 0._r8 endif if(snow_depth_avg < minh(iv))then fraction_exposed = 1._r8 endif if(snow_depth_avg>= minh(iv).and.snow_depth_avg <= maxh(iv))then !only partly hidden... fraction_exposed = max(0._r8,(min(1.0_r8,(snow_depth_avg-minh(iv))/dh))) endif fraction_exposed = 1.0_r8 ! no m2 of leaf per m2 of ground in each height class ! FIX(SPM,032414) these should be uncommented this and double check if ( debug ) write(fates_log(), *) 'leaf_area_profile()', currentPatch%elai_profile(1,ft,iv) currentPatch%elai_profile(1,ft,iv) = currentPatch%tlai_profile(1,ft,iv) * fraction_exposed currentPatch%esai_profile(1,ft,iv) = currentPatch%tsai_profile(1,ft,iv) * fraction_exposed if ( debug ) write(fates_log(), *) 'leaf_area_profile()', currentPatch%elai_profile(1,ft,iv) enddo ! (iv) hite bins currentCohort => currentCohort%taller enddo !currentCohort ! ----------------------------------------------------------------------------- ! Perform a leaf area conservation check on the LAI profile lai = 0.0_r8 do ft = 1,numpft lai = lai+ sum(currentPatch%tlai_profile(1,ft,:)) enddo if(lai > patch_lai)then write(fates_log(), *) 'FATES: problem with lai assignments' call endrun(msg=errMsg(sourcefile, __LINE__)) endif else ! smooth leaf distribution ! ----------------------------------------------------------------------------- ! Standard canopy layering model. ! Go through all cohorts and add their leaf area ! and canopy area to the accumulators. ! ----------------------------------------------------------------------------- currentCohort => currentPatch%shortest do while(associated(currentCohort)) ft = currentCohort%pft cl = currentCohort%canopy_layer ! ---------------------------------------------------------------- ! How much of each tree is stem area index? Assuming that there is ! This may indeed be zero if there is a sensecent grass ! ---------------------------------------------------------------- if( (currentCohort%treelai+currentCohort%treesai) > 0._r8)then fleaf = currentCohort%lai / (currentCohort%lai + currentCohort%sai) else fleaf = 0._r8 endif ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX ! SNOW BURIAL IS CURRENTLY TURNED OFF ! WHEN IT IS TURNED ON, IT WILL HAVE TO BE COMPARED ! WITH SNOW HEIGHTS CALCULATED BELOW. ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX currentPatch%nrad(cl,ft) = currentPatch%ncan(cl,ft) if (currentPatch%nrad(cl,ft) > nlevleaf ) then write(fates_log(), *) 'Number of radiative leaf layers is larger' write(fates_log(), *) ' than the maximum allowed.' write(fates_log(), *) ' cl: ',cl write(fates_log(), *) ' ft: ',ft write(fates_log(), *) ' nlevleaf: ',nlevleaf write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft) call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! -------------------------------------------------------------------------- ! Whole layers. Make a weighted average of the leaf area in each layer ! before dividing it by the total area. Fill up layer for whole layers. ! -------------------------------------------------------------------------- do iv = 1,currentCohort%NV ! This loop builds the arrays that define the effective (not snow covered) ! and total (includes snow covered) area indices for leaves and stems ! We calculate the absolute elevation of each layer to help determine if the layer ! is obscured by snow. layer_top_hite = currentCohort%hite - & ( real(iv-1,r8)/currentCohort%NV * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) layer_bottom_hite = currentCohort%hite - & ( real(iv,r8)/currentCohort%NV * currentCohort%hite * & EDPftvarcon_inst%crown(currentCohort%pft) ) fraction_exposed = 1.0_r8 snow_depth_avg = snow_depth_si * frac_sno_eff_si if(snow_depth_avg > layer_top_hite)then fraction_exposed = 0._r8 endif if(snow_depth_avg < layer_bottom_hite)then fraction_exposed = 1._r8 endif if( snow_depth_avg>= layer_bottom_hite .and. & snow_depth_avg <= layer_top_hite) then !only partly hidden... fraction_exposed = max(0._r8,(min(1.0_r8,(snow_depth_avg-layer_bottom_hite)/ & (layer_top_hite-layer_bottom_hite )))) endif ! =========== OVER-WRITE ================= fraction_exposed= 1.0_r8 ! =========== OVER-WRITE ================= if(iv==currentCohort%NV) then remainder = (currentCohort%treelai + currentCohort%treesai) - & (dinc_ed*real(currentCohort%nv-1,r8)) if(remainder > dinc_ed )then write(fates_log(), *)'ED: issue with remainder', & currentCohort%treelai,currentCohort%treesai,dinc_ed, & currentCohort%NV,remainder call endrun(msg=errMsg(sourcefile, __LINE__)) endif else remainder = dinc_ed end if currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + & remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) + & remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & fraction_exposed currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) + & remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) + & remainder * (1._r8 - fleaf) * currentCohort%c_area/currentPatch%total_canopy_area * & fraction_exposed currentPatch%canopy_area_profile(cl,ft,iv) = currentPatch%canopy_area_profile(cl,ft,iv) + & currentCohort%c_area/currentPatch%total_canopy_area currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) + & (remainder * fleaf * currentCohort%c_area/currentPatch%total_canopy_area * & (layer_top_hite+layer_bottom_hite)/2.0_r8) !average height of layer. end do currentCohort => currentCohort%taller enddo !cohort ! -------------------------------------------------------------------------- ! If there is an upper-story, the top canopy layer ! should have a value of exactly 1.0 in its top leaf layer ! -------------------------------------------------------------------------- if ( (currentPatch%NCL_p > 1) .and. & (sum(currentPatch%canopy_area_profile(1,:,1)) < 0.9999 )) then write(fates_log(), *) 'FATES: canopy_area_profile was less than 1 at the canopy top' write(fates_log(), *) 'cl: ',1 write(fates_log(), *) 'iv: ',1 write(fates_log(), *) 'sum(cpatch%canopy_area_profile(1,:,1)): ', & sum(currentPatch%canopy_area_profile(1,:,1)) currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==1)then write(fates_log(), *) 'FATES: cohorts',currentCohort%dbh,currentCohort%c_area, & currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area endif currentCohort => currentCohort%taller enddo !currentCohort call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! -------------------------------------------------------------------------- ! In the following loop we are now normalizing the effective and ! total area profiles to convert from units of leaf/stem area per vegetated ! canopy area, into leaf/stem area per area of their own radiative column ! which is typically the footprint of all cohorts contained in the canopy ! layer x pft bins. ! Also perform some checks on area normalization. ! Check the area of each leaf layer, across pfts. ! It should never be larger than 1 or less than 0. ! -------------------------------------------------------------------------- do cl = 1,currentPatch%NCL_p do iv = 1,currentPatch%ncan(cl,ft) if( debug .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0' write(fates_log(), *) 'cl: ',cl write(fates_log(), *) 'iv: ',iv write(fates_log(), *) 'sum(cpatch%canopy_area_profile(cl,:,iv)): ', & sum(currentPatch%canopy_area_profile(cl,:,iv)) currentCohort => currentPatch%shortest do while(associated(currentCohort)) if(currentCohort%canopy_layer==cl)then write(fates_log(), *) 'FATES: cohorts in layer cl = ',cl, & currentCohort%dbh,currentCohort%c_area, & currentPatch%total_canopy_area,currentPatch%area write(fates_log(), *) 'ED: fracarea', currentCohort%pft, & currentCohort%c_area/currentPatch%total_canopy_area endif currentCohort => currentCohort%taller enddo !currentCohort call endrun(msg=errMsg(sourcefile, __LINE__)) end if end do do ft = 1,numpft do iv = 1,currentPatch%ncan(cl,ft) if( currentPatch%canopy_area_profile(cl,ft,iv) > nearzero ) then currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) / & currentPatch%canopy_area_profile(cl,ft,iv) currentPatch%tsai_profile(cl,ft,iv) = currentPatch%tsai_profile(cl,ft,iv) / & currentPatch%canopy_area_profile(cl,ft,iv) currentPatch%elai_profile(cl,ft,iv) = currentPatch%elai_profile(cl,ft,iv) / & currentPatch%canopy_area_profile(cl,ft,iv) currentPatch%esai_profile(cl,ft,iv) = currentPatch%esai_profile(cl,ft,iv) / & currentPatch%canopy_area_profile(cl,ft,iv) end if if(currentPatch%tlai_profile(cl,ft,iv)>nearzero )then currentPatch%layer_height_profile(cl,ft,iv) = currentPatch%layer_height_profile(cl,ft,iv) & /currentPatch%tlai_profile(cl,ft,iv) end if enddo enddo enddo ! -------------------------------------------------------------------------- ! Set the mask that identifies which PFT x can-layer combinations have ! scattering elements in them. ! -------------------------------------------------------------------------- do cl = 1,currentPatch%NCL_p do ft = 1,numpft do iv = 1, currentPatch%nrad(cl,ft) if(currentPatch%canopy_area_profile(cl,ft,iv) > 0._r8)then currentPatch%canopy_mask(cl,ft) = 1 endif end do !iv enddo !ft enddo ! loop over cl endif !leaf distribution end if currentPatch => currentPatch%younger enddo !patch return end subroutine leaf_area_profile ! ====================================================================================== subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out) ! ---------------------------------------------------------------------------------- ! The purpose of this routine is to package output boundary conditions related ! to vegetation coverage to the host land model. ! ---------------------------------------------------------------------------------- use EDTypesMod , only : ed_patch_type, ed_cohort_type, & ed_site_type, AREA use FatesInterfaceMod , only : bc_out_type use EDPftvarcon , only : EDPftvarcon_inst ! ! !ARGUMENTS integer, intent(in) :: nsites type(ed_site_type), intent(inout), target :: sites(nsites) integer, intent(in) :: fcolumn(nsites) type(bc_out_type), intent(inout) :: bc_out(nsites) ! Locals type (ed_cohort_type) , pointer :: currentCohort integer :: s, ifp, c, p type (ed_patch_type) , pointer :: currentPatch real(r8) :: bare_frac_area real(r8) :: total_patch_area real(r8) :: total_canopy_area real(r8) :: weight ! Weighting for cohort variables in patch do s = 1,nsites ifp = 0 total_patch_area = 0._r8 total_canopy_area = 0._r8 bc_out(s)%canopy_fraction_pa(:) = 0._r8 currentPatch => sites(s)%oldest_patch c = fcolumn(s) do while(associated(currentPatch)) ifp = ifp+1 if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then write(fates_log(),*) 'ED: canopy area bigger than area',currentPatch%total_canopy_area ,currentPatch%area currentPatch%total_canopy_area = currentPatch%area endif if (associated(currentPatch%tallest)) then bc_out(s)%htop_pa(ifp) = currentPatch%tallest%hite else ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height? bc_out(s)%htop_pa(ifp) = 0.1_r8 endif bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8)) ! Use leaf area weighting for all cohorts in the patch to define the characteristic ! leaf width used by the HLM ! ---------------------------------------------------------------------------- ! bc_out(s)%dleaf_pa(ifp) = 0.0_r8 ! if(currentPatch%lai>1.0e-9_r8) then ! currentCohort => currentPatch%shortest ! do while(associated(currentCohort)) ! weight = min(1.0_r8,currentCohort%lai/currentPatch%lai) ! bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + & ! EDPftvarcon_inst%dleaf(currentCohort%pft)*weight ! currentCohort => currentCohort%taller ! enddo ! end if ! Roughness length and displacement height are not PFT properties, they are ! properties of the canopy assemblage. Defining this needs an appropriate model. ! Right now z0 and d are pft level parameters. For the time being we will just ! use the 1st index until a suitable model is defined. (RGK 04-2017) ! ----------------------------------------------------------------------------- bc_out(s)%z0m_pa(ifp) = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp) bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp) bc_out(s)%dleaf_pa(ifp) = EDPftvarcon_inst%dleaf(1) ! We are assuming here that grass is all located underneath tree canopies. ! The alternative is to assume it is all spatial distinct from tree canopies. ! In which case, the bare area would have to be reduced by the grass area... ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants ! currentPatch%area/AREA is the fraction of the soil covered by this patch. bc_out(s)%canopy_fraction_pa(ifp) = & min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA) bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * & (currentPatch%area/AREA) total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp) ! Calculate area indices for output boundary to HLM ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called) bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai') bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai') bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai') bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai') ! Fraction of vegetation free of snow. This is used to flag those ! patches which shall under-go photosynthesis ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let ! FATES internal variables decide if photosynthesis is possible ! we are essentially calculating it inside FATES to tell the ! host to tell itself when to do things (circuitous). Just have ! to determine where else it is used if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8 else bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8 end if currentPatch => currentPatch%younger end do ! Apply patch and canopy area corrections ! If the difference is above reasonable math precision, apply a fix ! If the difference is way above reasonable math precision, gracefully exit if(abs(total_patch_area-1.0_r8) > rsnbl_math_prec ) then if(abs(total_patch_area-1.0_r8) > 1.0e-8_r8 )then write(fates_log(),*) 'total area is wrong in update_hlm_dynamics',total_patch_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if if(debug) then write(fates_log(),*) 'imprecise patch areas in update_hlm_dynamics',total_patch_area end if currentPatch => sites(s)%oldest_patch ifp = 0 do while(associated(currentPatch)) ifp = ifp+1 bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area currentPatch => currentPatch%younger end do endif end do ! If hydraulics is turned on, update the amount of water bound in vegetation if (hlm_use_planthydro.eq.itrue) then call RecruitWaterStorage(nsites,sites,bc_out) call UpdateH2OVeg(nsites,sites,bc_out) end if end subroutine update_hlm_dynamics ! ===================================================================================== function calc_areaindex(cpatch,ai_type) result(ai) ! ---------------------------------------------------------------------------------- ! This subroutine calculates the exposed leaf area index of a patch ! this is the square meters of leaf per square meter of ground area ! It does so by integrating over the depth and functional type profile of leaf area ! which are per area of crown. This value has to be scaled by crown area to convert ! to ground area. ! ---------------------------------------------------------------------------------- ! Arguments type(ed_patch_type),intent(in), target :: cpatch character(len=*),intent(in) :: ai_type integer :: cl,ft real(r8) :: ai ! TODO: THIS MIN LAI IS AN ARTIFACT FROM TESTING LONG-AGO AND SHOULD BE REMOVED ! THIS HAS BEEN KEPT THUS FAR TO MAINTAIN B4B IN TESTING OTHER COMMITS real(r8),parameter :: ai_min = 0.1_r8 real(r8),pointer :: ai_profile ai = 0._r8 if (trim(ai_type) == 'elai') then do cl = 1,cpatch%NCL_p do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%elai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tlai') then do cl = 1,cpatch%NCL_p do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tlai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'esai') then do cl = 1,cpatch%NCL_p do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%esai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo elseif (trim(ai_type) == 'tsai') then do cl = 1,cpatch%NCL_p do ft = 1,numpft ai = ai + sum(cpatch%canopy_area_profile(cl,ft,1:cpatch%nrad(cl,ft)) * & cpatch%tsai_profile(cl,ft,1:cpatch%nrad(cl,ft))) enddo enddo else write(fates_log(),*) 'Unsupported area index sent to calc_areaindex' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ai = max(ai_min,ai) return end function calc_areaindex ! =============================================================================================== subroutine CanopyLayerArea(currentPatch,site_spread,layer_index,layer_area) ! -------------------------------------------------------------------------------------------- ! This function calculates the total crown area footprint for a desired layer of the canopy ! within a patch. ! The return units are the same as patch%area, which is m2 ! --------------------------------------------------------------------------------------------- ! Arguments type(ed_patch_type),intent(inout), target :: currentPatch real(r8),intent(in) :: site_spread integer,intent(in) :: layer_index real(r8),intent(inout) :: layer_area type(ed_cohort_type), pointer :: currentCohort layer_area = 0.0_r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) call carea_allom(currentCohort%dbh,currentCohort%n,site_spread, & currentCohort%pft,currentCohort%c_area) if (currentCohort%canopy_layer .eq. layer_index) then layer_area = layer_area + currentCohort%c_area end if currentCohort => currentCohort%shorter enddo return end subroutine CanopyLayerArea ! =============================================================================================== function NumPotentialCanopyLayers(currentPatch,site_spread,include_substory) result(z) ! -------------------------------------------------------------------------------------------- ! Calculate the number of canopy layers in this patch. ! This simple call only determines total layering by querying the cohorts ! which layer they are in, it doesn't do any size evaluation. ! It may also, optionally, account for the temporary "substory", which is the imaginary ! layer below the understory which will be needed to temporarily accomodate demotions from ! the understory in the event the understory has reached maximum allowable area. ! -------------------------------------------------------------------------------------------- type(ed_patch_type),target :: currentPatch real(r8),intent(in) :: site_spread logical :: include_substory type(ed_cohort_type),pointer :: currentCohort integer :: z real(r8) :: c_area real(r8) :: arealayer z = 1 currentCohort => currentPatch%tallest do while (associated(currentCohort)) z = max(z,currentCohort%canopy_layer) currentCohort => currentCohort%shorter enddo if(include_substory)then arealayer = 0.0 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == z) then call carea_allom(currentCohort%dbh,currentCohort%n,site_spread,currentCohort%pft,c_area) arealayer = arealayer + c_area end if currentCohort => currentCohort%shorter enddo ! Does the bottom layer have more than a full canopy? ! If so we need to make another layer. if(arealayer > currentPatch%area)then z = z + 1 endif end if end function NumPotentialCanopyLayers end module EDCanopyStructureMod