From 10ee77bee0be3a20a8e45e29dc60180045ce6508 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 23 May 2022 17:27:13 -0600 Subject: [PATCH 1/8] first attempt to make disturbance rates happen in parallel --- biogeochem/EDPatchDynamicsMod.F90 | 129 ++++++++---------------------- 1 file changed, 34 insertions(+), 95 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 13552c899a..d2e612fd74 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -187,6 +187,7 @@ subroutine disturbance_rates( site_in, bc_in) integer :: i_dist real(r8) :: frac_site_primary real(r8) :: harvest_rate + real(r8) :: tempsum !---------------------------------------------------------------------------------------------- ! Calculate Mortality Rates (these were previously calculated during growth derivatives) @@ -343,90 +344,12 @@ subroutine disturbance_rates( site_in, bc_in) call FatesWarn(msg,index=2) endif - - - - ! ------------------------------------------------------------------------------------------ - ! Determine which disturbance is dominant, and force mortality diagnostics in the upper - ! canopy to be zero for the non-dominant mode. Note: upper-canopy tree-fall mortality is - ! not always disturbance generating, so when tree-fall mort is non-dominant, make sure - ! to still diagnose and track the non-disturbance rate - ! ------------------------------------------------------------------------------------------ - - ! DISTURBANCE IS LOGGING - if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. & - currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire) ) then - - currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ilog) - currentPatch%disturbance_mode = dtype_ilog - - ! Update diagnostics - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer == 1)then - currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction) - end if - currentCohort => currentCohort%taller - enddo !currentCohort - - ! DISTURBANCE IS FIRE - elseif (currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ifall) .and. & - currentPatch%disturbance_rates(dtype_ifire) > currentPatch%disturbance_rates(dtype_ilog) ) then - - currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ifire) - currentPatch%disturbance_mode = dtype_ifire - - ! Update diagnostics, zero non-fire mortality rates - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer == 1)then - currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction) - currentCohort%lmort_direct = 0.0_r8 - currentCohort%lmort_collateral = 0.0_r8 - currentCohort%lmort_infra = 0.0_r8 - currentCohort%l_degrad = 0.0_r8 - end if - - ! This may be counter-intuitive, but the diagnostic fire-mortality rate - ! will stay zero in the patch that undergoes fire, this is because - ! the actual cohorts who experience the fire are only those in the - ! newly created patch so currentCohort%fmort = 0.0_r8 - ! Don't worry, the cohorts in the newly created patch will reflect burn - - currentCohort => currentCohort%taller - enddo !currentCohort - - else ! If fire and logging are not greater than treefall, just set disturbance rate to tree-fall - ! which is most likely a 0.0 - - currentPatch%disturbance_rate = currentPatch%disturbance_rates(dtype_ifall) - currentPatch%disturbance_mode = dtype_ifall - - ! Update diagnostics, zero non-treefall mortality rates - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - if(currentCohort%canopy_layer == 1)then - currentCohort%lmort_direct = 0.0_r8 - currentCohort%lmort_collateral = 0.0_r8 - currentCohort%lmort_infra = 0.0_r8 - currentCohort%l_degrad = 0.0_r8 - end if - currentCohort => currentCohort%taller - enddo !currentCohort - - + ! if the sum of all disturbance rates is such that they will exceed total patch area on this day, then reduce them all proportionally. + if ( sum(currentPatch%disturbance_rates(:)) .gt. 1.0_r8 ) then + tempsum = sum(currentPatch%disturbance_rates(:)) + do i_dist = 1,N_DIST_TYPES + currentPatch%disturbance_rates(i_dist) = currentPatch%disturbance_rates(i_dist) / tempsum + end do endif currentPatch => currentPatch%younger @@ -491,6 +414,7 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_m ! leaf mass during partial burn calculations logical :: found_youngest_primary ! logical for finding the first primary forest patch integer :: min_nocomp_pft, max_nocomp_pft, i_nocomp_pft + integer :: i_dist, i_dist2 !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -515,6 +439,8 @@ subroutine spawn_patches( currentSite, bc_in) ! If nocomp is not enabled, then this is not much of a loop, it only passes through once. nocomp_pft_loop: do i_nocomp_pft = min_nocomp_pft,max_nocomp_pft + disturbance_type_loop: do i_dist = 1,N_DIST_TYPES + ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. currentPatch => currentSite%youngest_patch @@ -532,12 +458,8 @@ subroutine spawn_patches( currentSite, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - ! Check to make sure that the disturbance mode of the patch is set - if( .not.any(currentPatch%disturbance_mode == [dtype_ilog,dtype_ifall,dtype_ifire])) then - write(fates_log(),*) 'undefined disturbance mode? : ',currentPatch%disturbance_mode - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - + currentPatch%disturbance_mode = i_dist + currentPatch%disturbance_rate = currentPatch%disturbance_rates(i_dist) ! Only create new patches that have non-negligible amount of land if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then @@ -1118,8 +1040,18 @@ subroutine spawn_patches( currentSite, bc_in) call sort_cohorts(currentPatch) !update area of donor patch + oldarea = currentPatch%area currentPatch%area = currentPatch%area - patch_site_areadis + ! for all disturbance rates that haven't been resolved yet, increase their amount so that + ! they are the same amount of gridcell-scale disturbance relative to the original patch size + if (i_dist .ne. N_DIST_TYPES) then + do i_dist2 = i_dist+1,N_DIST_TYPES + currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & + * oldarea / currentPatch%area + end do + end if + ! sort out the cohorts, since some of them may be so small as to need removing. ! the first call to terminate cohorts removes sparse number densities, ! the second call removes for all other reasons (sparse culling must happen @@ -1131,11 +1063,6 @@ subroutine spawn_patches( currentSite, bc_in) end if ! if ( new_patch%area > nearzero ) then - !zero disturbance rate trackers - currentPatch%disturbance_rate = 0._r8 - currentPatch%disturbance_rates = 0._r8 - currentPatch%fract_ldist_not_harvested = 0._r8 - end if cp_nocomp_matches_2_if currentPatch => currentPatch%younger @@ -1218,7 +1145,19 @@ subroutine spawn_patches( currentSite, bc_in) call check_patch_area(currentSite) call set_patchno(currentSite) + end do disturbance_type_loop + end do nocomp_pft_loop + + !zero disturbance rate trackers on all patches + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) + currentPatch%disturbance_rate = 0._r8 + currentPatch%disturbance_rates = 0._r8 + currentPatch%fract_ldist_not_harvested = 0._r8 + currentPatch => currentPatch%younger + end do + return end subroutine spawn_patches From 6ac7604f1037fa863a3b2e10a624713fac9bd20c Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 23 May 2022 17:30:22 -0600 Subject: [PATCH 2/8] bugfix --- biogeochem/EDPatchDynamicsMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index d2e612fd74..66e01ca95c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -415,6 +415,7 @@ subroutine spawn_patches( currentSite, bc_in) logical :: found_youngest_primary ! logical for finding the first primary forest patch integer :: min_nocomp_pft, max_nocomp_pft, i_nocomp_pft integer :: i_dist, i_dist2 + real(r8) :: oldarea !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine From 4004259ee095bc2476fa708a94ee2173b7c7ed01 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 24 May 2022 23:52:29 -0600 Subject: [PATCH 3/8] removed patch%disturbance_rate and patch%disturbance_mode as no longer accurate --- biogeochem/EDPatchDynamicsMod.F90 | 89 +++++++++++++++---------------- main/EDMainMod.F90 | 1 - main/EDTypesMod.F90 | 4 -- 3 files changed, 43 insertions(+), 51 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 66e01ca95c..8f2b5e5fad 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -246,7 +246,6 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort => currentCohort%taller end do - currentPatch%disturbance_mode = fates_unset_int currentPatch => currentPatch%younger end do @@ -414,7 +413,8 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_m ! leaf mass during partial burn calculations logical :: found_youngest_primary ! logical for finding the first primary forest patch integer :: min_nocomp_pft, max_nocomp_pft, i_nocomp_pft - integer :: i_dist, i_dist2 + integer :: i_disturbance_type, i_dist2 + real(r8) :: disturbance_rate real(r8) :: oldarea !--------------------------------------------------------------------- @@ -440,7 +440,7 @@ subroutine spawn_patches( currentSite, bc_in) ! If nocomp is not enabled, then this is not much of a loop, it only passes through once. nocomp_pft_loop: do i_nocomp_pft = min_nocomp_pft,max_nocomp_pft - disturbance_type_loop: do i_dist = 1,N_DIST_TYPES + disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. currentPatch => currentSite%youngest_patch @@ -453,42 +453,41 @@ subroutine spawn_patches( currentSite, bc_in) cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - if(currentPatch%disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then - write(fates_log(),*) 'patch disturbance rate > 1 ?',currentPatch%disturbance_rate + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + + if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then + write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate call dump_patch(currentPatch) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - currentPatch%disturbance_mode = i_dist - currentPatch%disturbance_rate = currentPatch%disturbance_rates(i_dist) - ! Only create new patches that have non-negligible amount of land - if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then + if((currentPatch%area*disturbance_rate) > nearzero ) then ! figure out whether the receiver patch for disturbance from this patch will be ! primary or secondary land receiver patch is primary forest only if both the - ! donor patch is primary forest and the dominant disturbance type is not logging + ! donor patch is primary forest and the current disturbance type is not logging if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (currentPatch%disturbance_mode .ne. dtype_ilog) ) then + (i_disturbance_type .ne. dtype_ilog) ) then - site_areadis_primary = site_areadis_primary + currentPatch%area * currentPatch%disturbance_rate + site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate ! track disturbance rates to output to history - currentSite%disturbance_rates_primary_to_primary(currentPatch%disturbance_mode) = & - currentSite%disturbance_rates_primary_to_primary(currentPatch%disturbance_mode) + & - currentPatch%area * currentPatch%disturbance_rate * AREA_INV + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV else - site_areadis_secondary = site_areadis_secondary + currentPatch%area * currentPatch%disturbance_rate + site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate ! track disturbance rates to output to history if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then - currentSite%disturbance_rates_secondary_to_secondary(currentPatch%disturbance_mode) = & - currentSite%disturbance_rates_secondary_to_secondary(currentPatch%disturbance_mode) + & - currentPatch%area * currentPatch%disturbance_rate * AREA_INV + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV else - currentSite%disturbance_rates_primary_to_secondary(currentPatch%disturbance_mode) = & - currentSite%disturbance_rates_primary_to_secondary(currentPatch%disturbance_mode) + & - currentPatch%area * currentPatch%disturbance_rate * AREA_INV + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV endif endif @@ -562,17 +561,18 @@ subroutine spawn_patches( currentSite, bc_in) currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then ! This is the amount of patch area that is disturbed, and donated by the donor - patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + patch_site_areadis = currentPatch%area * disturbance_rate if ( patch_site_areadis > nearzero ) then ! figure out whether the receiver patch for disturbance from this patch ! will be primary or secondary land receiver patch is primary forest - ! only if both the donor patch is primary forest and the dominant + ! only if both the donor patch is primary forest and the current ! disturbance type is not logging if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (currentPatch%disturbance_mode .ne. dtype_ilog)) then + (i_disturbance_type .ne. dtype_ilog)) then new_patch => new_patch_primary else new_patch => new_patch_secondary @@ -585,11 +585,11 @@ subroutine spawn_patches( currentSite, bc_in) end if ! for the case where the donating patch is secondary forest, if - ! the dominant disturbance from this patch is non-anthropogenic, + ! the current disturbance from this patch is non-anthropogenic, ! we need to average in the time-since-anthropogenic-disturbance ! from the donor patch into that of the receiver patch if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & - (currentPatch%disturbance_mode .ne. dtype_ilog) ) then + (i_disturbance_type .ne. dtype_ilog) ) then new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) @@ -600,9 +600,9 @@ subroutine spawn_patches( currentSite, bc_in) ! Transfer the litter existing already in the donor patch to the new patch ! This call will only transfer non-burned litter to new patch ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when - ! fire is not the dominant disturbance regime. + ! fire is not the current disturbance regime. - if(currentPatch%disturbance_mode .ne. dtype_ifire) then + if(i_disturbance_type .ne. dtype_ifire) then currentPatch%burnt_frac_litter(:) = 0._r8 end if @@ -610,10 +610,10 @@ subroutine spawn_patches( currentSite, bc_in) ! Transfer in litter fluxes from plants in various contexts of death and destruction - if(currentPatch%disturbance_mode .eq. dtype_ilog) then + if(i_disturbance_type .eq. dtype_ilog) then call logging_litter_fluxes(currentSite, currentPatch, & new_patch, patch_site_areadis,bc_in) - elseif(currentPatch%disturbance_mode .eq. dtype_ifire) then + elseif(i_disturbance_type .eq. dtype_ifire) then call fire_litter_fluxes(currentSite, currentPatch, & new_patch, patch_site_areadis,bc_in) else @@ -672,8 +672,8 @@ subroutine spawn_patches( currentSite, bc_in) store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c - ! treefall mortality is the dominant disturbance - if(currentPatch%disturbance_mode .eq. dtype_ifall) then + ! treefall mortality is the current disturbance + if(i_disturbance_type .eq. dtype_ifall) then if(currentCohort%canopy_layer == 1)then @@ -782,8 +782,8 @@ subroutine spawn_patches( currentSite, bc_in) endif endif - ! Fire is the dominant disturbance - elseif (currentPatch%disturbance_mode .eq. dtype_ifire ) then + ! Fire is the current disturbance + elseif (i_disturbance_type .eq. dtype_ifire ) then ! Number of members in the new patch, before we impose fire survivorship nc%n = currentCohort%n * patch_site_areadis/currentPatch%area @@ -881,8 +881,8 @@ subroutine spawn_patches( currentSite, bc_in) - ! Logging is the dominant disturbance - elseif (currentPatch%disturbance_mode .eq. dtype_ilog ) then + ! Logging is the current disturbance + elseif (i_disturbance_type .eq. dtype_ilog ) then ! If this cohort is in the upper canopy. It generated if(currentCohort%canopy_layer == 1)then @@ -916,7 +916,7 @@ subroutine spawn_patches( currentSite, bc_in) else - ! WHat to do with cohorts in the understory of a logging generated + ! What to do with cohorts in the understory of a logging generated ! disturbance patch? if(int(prt_params%woody(currentCohort%pft)) == itrue)then @@ -1000,7 +1000,7 @@ subroutine spawn_patches( currentSite, bc_in) else write(fates_log(),*) 'unknown disturbance mode?' - write(fates_log(),*) 'disturbance_mode: ',currentPatch%disturbance_mode + write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Select disturbance mode @@ -1046,8 +1046,8 @@ subroutine spawn_patches( currentSite, bc_in) ! for all disturbance rates that haven't been resolved yet, increase their amount so that ! they are the same amount of gridcell-scale disturbance relative to the original patch size - if (i_dist .ne. N_DIST_TYPES) then - do i_dist2 = i_dist+1,N_DIST_TYPES + if (i_disturbance_type .ne. N_DIST_TYPES) then + do i_dist2 = i_disturbance_type+1,N_DIST_TYPES currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & * oldarea / currentPatch%area end do @@ -1153,8 +1153,7 @@ subroutine spawn_patches( currentSite, bc_in) !zero disturbance rate trackers on all patches currentPatch => currentSite%oldest_patch do while(associated(currentPatch)) - currentPatch%disturbance_rate = 0._r8 - currentPatch%disturbance_rates = 0._r8 + currentPatch%disturbance_rates(:) = 0._r8 currentPatch%fract_ldist_not_harvested = 0._r8 currentPatch => currentPatch%younger end do @@ -2047,7 +2046,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, label,nocomp_pft) ! This new value will be generated when the calculate disturbance ! rates routine is called. This does not need to be remembered or in the restart file. - new_patch%disturbance_mode = fates_unset_int new_patch%f_sun = 0._r8 new_patch%ed_laisun_z(:,:,:) = 0._r8 @@ -2150,8 +2148,7 @@ subroutine zero_patch(cp_p) currentPatch%pft_agb_profile(:,:) = nan ! DISTURBANCE - currentPatch%disturbance_rates = 0._r8 - currentPatch%disturbance_rate = 0._r8 + currentPatch%disturbance_rates(:) = 0._r8 currentPatch%fract_ldist_not_harvested = 0._r8 diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index d05966ef14..7b9d623413 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -847,7 +847,6 @@ subroutine TotalBalanceCheck (currentSite, call_index ) write(fates_log(),*) 'BG CWD (by layer): ', sum(litt%bg_cwd,dim=1) write(fates_log(),*) 'leaf litter:',sum(litt%leaf_fines) write(fates_log(),*) 'root litter (by layer): ',sum(litt%root_fines,dim=1) - write(fates_log(),*) 'dist mode: ',currentPatch%disturbance_mode write(fates_log(),*) 'anthro_disturbance_label: ',currentPatch%anthro_disturbance_label write(fates_log(),*) 'use_this_pft: ', currentSite%use_this_pft(:) if(print_cohorts)then diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 9ce4d5fe44..d108081045 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -537,9 +537,6 @@ module EDTypesMod real(r8) :: disturbance_rates(n_dist_types) ! disturbance rate from 1) mortality ! 2) fire: fraction/day ! 3) logging mortatliy - real(r8) :: disturbance_rate ! larger effective disturbance rate: fraction/day - integer :: disturbance_mode ! index identifying which disturbance was applied - ! can be one of: dtype_ifall, dtype_ilog or dtype_ifire real(r8) :: fract_ldist_not_harvested ! fraction of logged area that is canopy trees that weren't harvested @@ -1016,7 +1013,6 @@ subroutine dump_patch(cpatch) write(fates_log(),*) 'pa%gnd_alb_dir = ',cpatch%gnd_alb_dir(:) write(fates_log(),*) 'pa%c_stomata = ',cpatch%c_stomata write(fates_log(),*) 'pa%c_lblayer = ',cpatch%c_lblayer - write(fates_log(),*) 'pa%disturbance_rate = ',cpatch%disturbance_rate write(fates_log(),*) 'pa%disturbance_rates = ',cpatch%disturbance_rates(:) write(fates_log(),*) 'pa%anthro_disturbance_label = ',cpatch%anthro_disturbance_label write(fates_log(),*) '----------------------------------------' From 7dbe7538fead5d973cac2ab7ebaaf02d030b4063 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 24 May 2022 23:58:36 -0600 Subject: [PATCH 4/8] added variable descriptions --- biogeochem/EDPatchDynamicsMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 8f2b5e5fad..b3bab06bcb 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -413,9 +413,9 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_m ! leaf mass during partial burn calculations logical :: found_youngest_primary ! logical for finding the first primary forest patch integer :: min_nocomp_pft, max_nocomp_pft, i_nocomp_pft - integer :: i_disturbance_type, i_dist2 - real(r8) :: disturbance_rate - real(r8) :: oldarea + integer :: i_disturbance_type, i_dist2 ! iterators for looping over disturbance types + real(r8) :: disturbance_rate ! rate of disturbance being resolved [fraction of patch area / day] + real(r8) :: oldarea ! old patch area prior to disturbance !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine From 53ca295e84cdb6b0a0fe99ba937161dc19d76987 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 23 Jun 2022 14:55:45 -0600 Subject: [PATCH 5/8] auto-indented spawn_patches --- biogeochem/EDPatchDynamicsMod.F90 | 1326 ++++++++++++++--------------- 1 file changed, 663 insertions(+), 663 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index b3bab06bcb..dc3cfa373e 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -436,717 +436,717 @@ subroutine spawn_patches( currentSite, bc_in) ! in the nocomp cases, since every patch has a PFT identity, it can only receive patch area from patches ! that have the same identity. In order to allow this, we have this very high level loop over nocomp PFTs - ! and only do the disturbance for any patches that have that nocomp PFT identity. + ! and only do the disturbance for any patches that have that nocomp PFT identity. ! If nocomp is not enabled, then this is not much of a loop, it only passes through once. nocomp_pft_loop: do i_nocomp_pft = min_nocomp_pft,max_nocomp_pft - disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES + disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES - ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. - currentPatch => currentSite%youngest_patch + ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. + currentPatch => currentSite%youngest_patch - site_areadis_primary = 0.0_r8 - site_areadis_secondary = 0.0_r8 + site_areadis_primary = 0.0_r8 + site_areadis_secondary = 0.0_r8 - do while(associated(currentPatch)) + do while(associated(currentPatch)) - cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then - write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate - call dump_patch(currentPatch) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - ! Only create new patches that have non-negligible amount of land - if((currentPatch%area*disturbance_rate) > nearzero ) then - - ! figure out whether the receiver patch for disturbance from this patch will be - ! primary or secondary land receiver patch is primary forest only if both the - ! donor patch is primary forest and the current disturbance type is not logging - if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then - - site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate + if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then + write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate + call dump_patch(currentPatch) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! track disturbance rates to output to history - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate + ! Only create new patches that have non-negligible amount of land + if((currentPatch%area*disturbance_rate) > nearzero ) then - ! track disturbance rates to output to history - if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - endif + ! figure out whether the receiver patch for disturbance from this patch will be + ! primary or secondary land receiver patch is primary forest only if both the + ! donor patch is primary forest and the current disturbance type is not logging + if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then - endif - - end if - - end if cp_nocomp_matches_1_if - currentPatch => currentPatch%older - enddo ! end loop over patches. sum area disturbed for all patches. + site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate - ! It is possible that no disturbance area was generated - if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then - - age = 0.0_r8 + ! track disturbance rates to output to history + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate - ! create two empty patches, to absorb newly disturbed primary and secondary forest area - ! first create patch to receive primary forest area - if ( site_areadis_primary .gt. nearzero ) then - allocate(new_patch_primary) + ! track disturbance rates to output to history + if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + endif - call create_patch(currentSite, new_patch_primary, age, & - site_areadis_primary, primaryforest, i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_primary%tallest => null() - new_patch_primary%shortest => null() + endif - endif + end if - ! next create patch to receive secondary forest area - if ( site_areadis_secondary .gt. nearzero) then - allocate(new_patch_secondary) - call create_patch(currentSite, new_patch_secondary, age, & - site_areadis_secondary, secondaryforest,i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_secondary%tallest => null() - new_patch_secondary%shortest => null() + end if cp_nocomp_matches_1_if + currentPatch => currentPatch%older + enddo ! end loop over patches. sum area disturbed for all patches. + + ! It is possible that no disturbance area was generated + if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then + + age = 0.0_r8 + + ! create two empty patches, to absorb newly disturbed primary and secondary forest area + ! first create patch to receive primary forest area + if ( site_areadis_primary .gt. nearzero ) then + allocate(new_patch_primary) + + call create_patch(currentSite, new_patch_primary, age, & + site_areadis_primary, primaryforest, i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_primary%tallest => null() + new_patch_primary%shortest => null() - endif - - ! loop round all the patches that contribute surviving indivduals and litter - ! pools to the new patch. We only loop the pre-existing patches, so - ! quit the loop if the current patch is either null, or matches the - ! two new pointers. + endif - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + ! next create patch to receive secondary forest area + if ( site_areadis_secondary .gt. nearzero) then + allocate(new_patch_secondary) + call create_patch(currentSite, new_patch_secondary, age, & + site_areadis_secondary, secondaryforest,i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_secondary%tallest => null() + new_patch_secondary%shortest => null() - cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then + endif - ! This is the amount of patch area that is disturbed, and donated by the donor - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - patch_site_areadis = currentPatch%area * disturbance_rate + ! loop round all the patches that contribute surviving indivduals and litter + ! pools to the new patch. We only loop the pre-existing patches, so + ! quit the loop if the current patch is either null, or matches the + ! two new pointers. - - if ( patch_site_areadis > nearzero ) then - - ! figure out whether the receiver patch for disturbance from this patch - ! will be primary or secondary land receiver patch is primary forest - ! only if both the donor patch is primary forest and the current - ! disturbance type is not logging - if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog)) then - new_patch => new_patch_primary - else - new_patch => new_patch_secondary - endif - - if(.not.associated(new_patch))then - write(fates_log(),*) 'Patch spawning has attempted to point to' - write(fates_log(),*) 'an un-allocated patch' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - ! for the case where the donating patch is secondary forest, if - ! the current disturbance from this patch is non-anthropogenic, - ! we need to average in the time-since-anthropogenic-disturbance - ! from the donor patch into that of the receiver patch - if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then - - new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & - currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) - endif - - - ! Transfer the litter existing already in the donor patch to the new patch - ! This call will only transfer non-burned litter to new patch - ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when - ! fire is not the current disturbance regime. + cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - if(i_disturbance_type .ne. dtype_ifire) then - currentPatch%burnt_frac_litter(:) = 0._r8 - end if + ! This is the amount of patch area that is disturbed, and donated by the donor + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + patch_site_areadis = currentPatch%area * disturbance_rate - call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) - ! Transfer in litter fluxes from plants in various contexts of death and destruction + if ( patch_site_areadis > nearzero ) then - if(i_disturbance_type .eq. dtype_ilog) then - call logging_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - elseif(i_disturbance_type .eq. dtype_ifire) then - call fire_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - else - call mortality_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - endif + ! figure out whether the receiver patch for disturbance from this patch + ! will be primary or secondary land receiver patch is primary forest + ! only if both the donor patch is primary forest and the current + ! disturbance type is not logging + if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog)) then + new_patch => new_patch_primary + else + new_patch => new_patch_secondary + endif + if(.not.associated(new_patch))then + write(fates_log(),*) 'Patch spawning has attempted to point to' + write(fates_log(),*) 'an un-allocated patch' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! Copy any means or timers from the original patch to the new patch - ! These values will inherit all info from the original patch - ! -------------------------------------------------------------------------- - call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) - call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) - - - ! -------------------------------------------------------------------------- - ! The newly formed patch from disturbance (new_patch), has now been given - ! some litter from dead plants and pre-existing litter from the donor patches. - ! - ! Next, we loop through the cohorts in the donor patch, copy them with - ! area modified number density into the new-patch, and apply survivorship. - ! ------------------------------------------------------------------------- - - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - - allocate(nc) - if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) - - ! Initialize the PARTEH object and point to the - ! correct boundary condition fields - nc%prt => null() - call InitPRTObject(nc%prt) - call InitPRTBoundaryConditions(nc) - - ! (Keeping as an example) - ! Allocate running mean functions - !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) - - call zero_cohort(nc) - - ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort - ! is the curent cohort that stays in the donor patch (currentPatch) - call copy_cohort(currentCohort, nc) - - !this is the case as the new patch probably doesn't have a closed canopy, and - ! even if it does, that will be sorted out in canopy_structure. - nc%canopy_layer = 1 - nc%canopy_layer_yesterday = 1._r8 - - sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) - store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) - total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c - - ! treefall mortality is the current disturbance - if(i_disturbance_type .eq. dtype_ifall) then - - if(currentCohort%canopy_layer == 1)then - - ! In the donor patch we are left with fewer trees because the area has decreased - ! the plant density for large trees does not actually decrease in the donor patch - ! because this is the part of the original patch where no trees have actually fallen - ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved - - currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & - min(1.0_r8,currentCohort%dmort * hlm_freq_day)) - - nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. - - nc%cmort = nan ! The mortality diagnostics are set to nan - ! because the cohort should dissappear - nc%hmort = nan - nc%bmort = nan - nc%frmort = nan - nc%smort = nan - nc%asmort = nan - nc%lmort_direct = nan - nc%lmort_collateral = nan - nc%lmort_infra = nan - nc%l_degrad = nan - - else - ! small trees - if( int(prt_params%woody(currentCohort%pft)) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller and cohort counts - ! are absolute, reduce this number. - - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which - ! had been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to history - ! results in large numerical errors, on account of the sharply - ! reduced number densities. so instead pass this info via a - ! site-level diagnostic variable before reducing the number density. - - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * ED_val_understorey_death / hlm_freq_day - - - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & - (nc%n * ED_val_understorey_death / hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - ! Step 2: Apply survivor ship function based on the understory death fraction - ! remaining of understory plants of those that are knocked over - ! by the overstorey trees dying... - nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) - - ! since the donor patch split and sent a fraction of its members - ! to the new patch and a fraction to be preserved in itself, - ! when reporting diagnostic rates, we must carry over the mortality rates from - ! the donor that were applied before the patch split. Remember this is only - ! for diagnostics. But think of it this way, the rates are weighted by - ! number density in EDCLMLink, and the number density of this new patch is donated - ! so with the number density must come the effective mortality rates. - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - ! understory trees that might potentially be knocked over in the disturbance. - ! The existing (donor) patch should not have any impact mortality, it should - ! only lose cohorts due to the decrease in area. This is not mortality. - ! Besides, the current and newly created patch sum to unity - - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - else - ! grass is not killed by mortality disturbance events. Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif - endif - - ! Fire is the current disturbance - elseif (i_disturbance_type .eq. dtype_ifire ) then - - ! Number of members in the new patch, before we impose fire survivorship - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! loss of individuals from source patch due to area shrinking - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - levcan = currentCohort%canopy_layer - - if(levcan==ican_upper) then - - ! before changing number densities, track total rate of trees that died - ! due to fire, as well as from each fire mortality term - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - - else - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - end if - - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%cambial_mort / hlm_freq_day - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%crownfire_mort / hlm_freq_day - - ! loss of individual from fire in new patch. - nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - - ! Some of of the leaf mass from living plants has been - ! burned off. Here, we remove that mass, and - ! tally it in the flux we sent to the atmosphere - - if(int(prt_params%woody(currentCohort%pft)) == itrue)then - leaf_burn_frac = currentCohort%fraction_crown_burned - else + ! for the case where the donating patch is secondary forest, if + ! the current disturbance from this patch is non-anthropogenic, + ! we need to average in the time-since-anthropogenic-disturbance + ! from the donor patch into that of the receiver patch + if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then - ! Grasses determine their fraction of leaves burned here + new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & + currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) - leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) - endif - - ! Perform a check to make sure that spitfire gave - ! us reasonable mortality and burn fraction rates - - if( (leaf_burn_frac < 0._r8) .or. & - (leaf_burn_frac > 1._r8) .or. & - (currentCohort%fire_mort < 0._r8) .or. & - (currentCohort%fire_mort > 1._r8)) then - write(fates_log(),*) 'unexpected fire fractions' - write(fates_log(),*) prt_params%woody(currentCohort%pft) - write(fates_log(),*) leaf_burn_frac - write(fates_log(),*) currentCohort%fire_mort - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - do el = 1,num_elements - - leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) - - currentSite%mass_balance(el)%burn_flux_to_atm = & - currentSite%mass_balance(el)%burn_flux_to_atm + & - leaf_burn_frac * leaf_m * nc%n - end do + endif - ! Here the mass is removed from the plant - call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) - currentCohort%fraction_crown_burned = 0.0_r8 - nc%fraction_crown_burned = 0.0_r8 + ! Transfer the litter existing already in the donor patch to the new patch + ! This call will only transfer non-burned litter to new patch + ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when + ! fire is not the current disturbance regime. + if(i_disturbance_type .ne. dtype_ifire) then + currentPatch%burnt_frac_litter(:) = 0._r8 + end if + call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) - ! Logging is the current disturbance - elseif (i_disturbance_type .eq. dtype_ilog ) then - - ! If this cohort is in the upper canopy. It generated - if(currentCohort%canopy_layer == 1)then - - ! calculate the survivorship of disturbed trees because non-harvested - nc%n = currentCohort%n * currentCohort%l_degrad - ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & - ! currentCohort%lmort_direct + currentCohort%lmort_collateral + - ! currentCohort%lmort_infra) ) * & - ! currentCohort%n * patch_site_areadis/currentPatch%area - - ! Reduce counts in the existing/donor patch according to the logging rate - currentCohort%n = currentCohort%n * & - (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & - currentCohort%lmort_collateral + & - currentCohort%lmort_infra + currentCohort%l_degrad))) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - - ! since these are the ones that weren't logged, - ! set the logging mortality rates as zero - nc%lmort_direct = 0._r8 - nc%lmort_collateral = 0._r8 - nc%lmort_infra = 0._r8 - - else - - ! What to do with cohorts in the understory of a logging generated - ! disturbance patch? - - if(int(prt_params%woody(currentCohort%pft)) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller - ! and cohort counts are absolute, reduce this number. - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which had - ! been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to - ! history results in large numerical errors, on account - ! of the sharply reduced number densities. so instead pass this info - ! via a site-level diagnostic variable before reducing - ! the number density. - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac / hlm_freq_day - - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & - (nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac/ hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - - ! Step 2: Apply survivor ship function based on the understory death fraction - - ! remaining of understory plants of those that are knocked - ! over by the overstorey trees dying... - ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER - ! in the fatesparameter files - nc%n = nc%n * (1.0_r8 - & - (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) - - ! Step 3: Reduce the number count of cohorts in the - ! original/donor/non-disturbed patch to reflect the area change - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - + ! Transfer in litter fluxes from plants in various contexts of death and destruction + + if(i_disturbance_type .eq. dtype_ilog) then + call logging_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + elseif(i_disturbance_type .eq. dtype_ifire) then + call fire_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) else - - ! grass is not killed by mortality disturbance events. - ! Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - ! No grass impact mortality imposed on the newly created patch - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif ! is/is-not woody - - endif ! Select canopy layer - - else - write(fates_log(),*) 'unknown disturbance mode?' - write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if ! Select disturbance mode - - if (nc%n > 0.0_r8) then - storebigcohort => new_patch%tallest - storesmallcohort => new_patch%shortest - if(associated(new_patch%tallest))then - tnull = 0 - else - tnull = 1 - new_patch%tallest => nc - nc%taller => null() - endif - - if(associated(new_patch%shortest))then - snull = 0 + call mortality_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + endif + + + ! Copy any means or timers from the original patch to the new patch + ! These values will inherit all info from the original patch + ! -------------------------------------------------------------------------- + call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) + call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) + + + ! -------------------------------------------------------------------------- + ! The newly formed patch from disturbance (new_patch), has now been given + ! some litter from dead plants and pre-existing litter from the donor patches. + ! + ! Next, we loop through the cohorts in the donor patch, copy them with + ! area modified number density into the new-patch, and apply survivorship. + ! ------------------------------------------------------------------------- + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + + allocate(nc) + if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) + + ! Initialize the PARTEH object and point to the + ! correct boundary condition fields + nc%prt => null() + call InitPRTObject(nc%prt) + call InitPRTBoundaryConditions(nc) + + ! (Keeping as an example) + ! Allocate running mean functions + !allocate(nc%tveg_lpa) + !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + + call zero_cohort(nc) + + ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort + ! is the curent cohort that stays in the donor patch (currentPatch) + call copy_cohort(currentCohort, nc) + + !this is the case as the new patch probably doesn't have a closed canopy, and + ! even if it does, that will be sorted out in canopy_structure. + nc%canopy_layer = 1 + nc%canopy_layer_yesterday = 1._r8 + + sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) + total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c + + ! treefall mortality is the current disturbance + if(i_disturbance_type .eq. dtype_ifall) then + + if(currentCohort%canopy_layer == 1)then + + ! In the donor patch we are left with fewer trees because the area has decreased + ! the plant density for large trees does not actually decrease in the donor patch + ! because this is the part of the original patch where no trees have actually fallen + ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved + + currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & + min(1.0_r8,currentCohort%dmort * hlm_freq_day)) + + nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + + nc%cmort = nan ! The mortality diagnostics are set to nan + ! because the cohort should dissappear + nc%hmort = nan + nc%bmort = nan + nc%frmort = nan + nc%smort = nan + nc%asmort = nan + nc%lmort_direct = nan + nc%lmort_collateral = nan + nc%lmort_infra = nan + nc%l_degrad = nan + + else + ! small trees + if( int(prt_params%woody(currentCohort%pft)) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller and cohort counts + ! are absolute, reduce this number. + + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which + ! had been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to history + ! results in large numerical errors, on account of the sharply + ! reduced number densities. so instead pass this info via a + ! site-level diagnostic variable before reducing the number density. + + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * ED_val_understorey_death / hlm_freq_day + + + currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + (nc%n * ED_val_understorey_death / hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + ! Step 2: Apply survivor ship function based on the understory death fraction + ! remaining of understory plants of those that are knocked over + ! by the overstorey trees dying... + nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) + + ! since the donor patch split and sent a fraction of its members + ! to the new patch and a fraction to be preserved in itself, + ! when reporting diagnostic rates, we must carry over the mortality rates from + ! the donor that were applied before the patch split. Remember this is only + ! for diagnostics. But think of it this way, the rates are weighted by + ! number density in EDCLMLink, and the number density of this new patch is donated + ! so with the number density must come the effective mortality rates. + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + ! understory trees that might potentially be knocked over in the disturbance. + ! The existing (donor) patch should not have any impact mortality, it should + ! only lose cohorts due to the decrease in area. This is not mortality. + ! Besides, the current and newly created patch sum to unity + + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + else + ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif + endif + + ! Fire is the current disturbance + elseif (i_disturbance_type .eq. dtype_ifire ) then + + ! Number of members in the new patch, before we impose fire survivorship + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! loss of individuals from source patch due to area shrinking + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + levcan = currentCohort%canopy_layer + + if(levcan==ican_upper) then + + ! before changing number densities, track total rate of trees that died + ! due to fire, as well as from each fire mortality term + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + + else + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + end if + + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%cambial_mort / hlm_freq_day + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%crownfire_mort / hlm_freq_day + + ! loss of individual from fire in new patch. + nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + + ! Some of of the leaf mass from living plants has been + ! burned off. Here, we remove that mass, and + ! tally it in the flux we sent to the atmosphere + + if(int(prt_params%woody(currentCohort%pft)) == itrue)then + leaf_burn_frac = currentCohort%fraction_crown_burned + else + + ! Grasses determine their fraction of leaves burned here + + leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) + endif + + ! Perform a check to make sure that spitfire gave + ! us reasonable mortality and burn fraction rates + + if( (leaf_burn_frac < 0._r8) .or. & + (leaf_burn_frac > 1._r8) .or. & + (currentCohort%fire_mort < 0._r8) .or. & + (currentCohort%fire_mort > 1._r8)) then + write(fates_log(),*) 'unexpected fire fractions' + write(fates_log(),*) prt_params%woody(currentCohort%pft) + write(fates_log(),*) leaf_burn_frac + write(fates_log(),*) currentCohort%fire_mort + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + do el = 1,num_elements + + leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) + + currentSite%mass_balance(el)%burn_flux_to_atm = & + currentSite%mass_balance(el)%burn_flux_to_atm + & + leaf_burn_frac * leaf_m * nc%n + end do + + ! Here the mass is removed from the plant + + call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) + currentCohort%fraction_crown_burned = 0.0_r8 + nc%fraction_crown_burned = 0.0_r8 + + + + ! Logging is the current disturbance + elseif (i_disturbance_type .eq. dtype_ilog ) then + + ! If this cohort is in the upper canopy. It generated + if(currentCohort%canopy_layer == 1)then + + ! calculate the survivorship of disturbed trees because non-harvested + nc%n = currentCohort%n * currentCohort%l_degrad + ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & + ! currentCohort%lmort_direct + currentCohort%lmort_collateral + + ! currentCohort%lmort_infra) ) * & + ! currentCohort%n * patch_site_areadis/currentPatch%area + + ! Reduce counts in the existing/donor patch according to the logging rate + currentCohort%n = currentCohort%n * & + (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra + currentCohort%l_degrad))) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + + ! since these are the ones that weren't logged, + ! set the logging mortality rates as zero + nc%lmort_direct = 0._r8 + nc%lmort_collateral = 0._r8 + nc%lmort_infra = 0._r8 + + else + + ! What to do with cohorts in the understory of a logging generated + ! disturbance patch? + + if(int(prt_params%woody(currentCohort%pft)) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller + ! and cohort counts are absolute, reduce this number. + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which had + ! been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to + ! history results in large numerical errors, on account + ! of the sharply reduced number densities. so instead pass this info + ! via a site-level diagnostic variable before reducing + ! the number density. + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac / hlm_freq_day + + currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + (nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac/ hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + + ! Step 2: Apply survivor ship function based on the understory death fraction + + ! remaining of understory plants of those that are knocked + ! over by the overstorey trees dying... + ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER + ! in the fatesparameter files + nc%n = nc%n * (1.0_r8 - & + (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) + + ! Step 3: Reduce the number count of cohorts in the + ! original/donor/non-disturbed patch to reflect the area change + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + else + + ! grass is not killed by mortality disturbance events. + ! Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + ! No grass impact mortality imposed on the newly created patch + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif ! is/is-not woody + + endif ! Select canopy layer + + else + write(fates_log(),*) 'unknown disturbance mode?' + write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Select disturbance mode + + if (nc%n > 0.0_r8) then + storebigcohort => new_patch%tallest + storesmallcohort => new_patch%shortest + if(associated(new_patch%tallest))then + tnull = 0 + else + tnull = 1 + new_patch%tallest => nc + nc%taller => null() + endif + + if(associated(new_patch%shortest))then + snull = 0 + else + snull = 1 + new_patch%shortest => nc + nc%shorter => null() + endif + nc%patchptr => new_patch + call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & + tnull, snull, storebigcohort, storesmallcohort) + + new_patch%tallest => storebigcohort + new_patch%shortest => storesmallcohort + else + + ! Get rid of the new temporary cohort + call DeallocateCohort(nc) + deallocate(nc) + + endif + + currentCohort => currentCohort%taller + enddo ! currentCohort + call sort_cohorts(currentPatch) + + !update area of donor patch + oldarea = currentPatch%area + currentPatch%area = currentPatch%area - patch_site_areadis + + ! for all disturbance rates that haven't been resolved yet, increase their amount so that + ! they are the same amount of gridcell-scale disturbance relative to the original patch size + if (i_disturbance_type .ne. N_DIST_TYPES) then + do i_dist2 = i_disturbance_type+1,N_DIST_TYPES + currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & + * oldarea / currentPatch%area + end do + end if + + ! sort out the cohorts, since some of them may be so small as to need removing. + ! the first call to terminate cohorts removes sparse number densities, + ! the second call removes for all other reasons (sparse culling must happen + ! before fusion) + call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) + call fuse_cohorts(currentSite,currentPatch, bc_in) + call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) + call sort_cohorts(currentPatch) + + end if ! if ( new_patch%area > nearzero ) then + + end if cp_nocomp_matches_2_if + currentPatch => currentPatch%younger + + enddo ! currentPatch patch loop. + + !*************************/ + !** INSERT NEW PATCH(ES) INTO LINKED LIST + !*************************/ + + if ( site_areadis_primary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + ! insert new youngest primary patch after all the secondary patches, if there are any. + ! this requires first finding the current youngest primary to insert the new one ahead of + if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then + found_youngest_primary = .false. + do while(associated(currentPatch) .and. .not. found_youngest_primary) + currentPatch => currentPatch%older + if (associated(currentPatch)) then + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + found_youngest_primary = .true. + endif + endif + end do + if (associated(currentPatch)) then + ! the case where we've found a youngest primary patch + new_patch_primary%older => currentPatch + new_patch_primary%younger => currentPatch%younger + currentPatch%younger%older => new_patch_primary + currentPatch%younger => new_patch_primary else - snull = 1 - new_patch%shortest => nc - nc%shorter => null() + ! the case where we haven't, because the patches are all secondaary, + ! and are putting a primary patch at the oldest end of the + ! linked list (not sure how this could happen, but who knows...) + new_patch_primary%older => null() + new_patch_primary%younger => currentSite%oldest_patch + currentSite%oldest_patch%older => new_patch_primary + currentSite%oldest_patch => new_patch_primary endif - nc%patchptr => new_patch - call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & - tnull, snull, storebigcohort, storesmallcohort) - - new_patch%tallest => storebigcohort - new_patch%shortest => storesmallcohort else - - ! Get rid of the new temporary cohort - call DeallocateCohort(nc) - deallocate(nc) - + ! the case where there are no secondary patches at the start of the linked list (prior logic) + new_patch_primary%older => currentPatch + new_patch_primary%younger => null() + currentPatch%younger => new_patch_primary + currentSite%youngest_patch => new_patch_primary endif - - currentCohort => currentCohort%taller - enddo ! currentCohort - call sort_cohorts(currentPatch) - - !update area of donor patch - oldarea = currentPatch%area - currentPatch%area = currentPatch%area - patch_site_areadis - - ! for all disturbance rates that haven't been resolved yet, increase their amount so that - ! they are the same amount of gridcell-scale disturbance relative to the original patch size - if (i_disturbance_type .ne. N_DIST_TYPES) then - do i_dist2 = i_disturbance_type+1,N_DIST_TYPES - currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & - * oldarea / currentPatch%area - end do - end if + endif - ! sort out the cohorts, since some of them may be so small as to need removing. + ! insert first secondary at the start of the list + if ( site_areadis_secondary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + new_patch_secondary%older => currentPatch + new_patch_secondary%younger=> null() + currentPatch%younger => new_patch_secondary + currentSite%youngest_patch => new_patch_secondary + endif + + + ! sort out the cohorts, since some of them may be so small as to need removing. ! the first call to terminate cohorts removes sparse number densities, ! the second call removes for all other reasons (sparse culling must happen ! before fusion) - call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) - call fuse_cohorts(currentSite,currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) - call sort_cohorts(currentPatch) - end if ! if ( new_patch%area > nearzero ) then - - end if cp_nocomp_matches_2_if - currentPatch => currentPatch%younger - - enddo ! currentPatch patch loop. + if ( site_areadis_primary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) + call fuse_cohorts(currentSite,new_patch_primary, bc_in) + call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) + call sort_cohorts(new_patch_primary) + endif - !*************************/ - !** INSERT NEW PATCH(ES) INTO LINKED LIST - !*************************/ - - if ( site_areadis_primary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - ! insert new youngest primary patch after all the secondary patches, if there are any. - ! this requires first finding the current youngest primary to insert the new one ahead of - if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then - found_youngest_primary = .false. - do while(associated(currentPatch) .and. .not. found_youngest_primary) - currentPatch => currentPatch%older - if (associated(currentPatch)) then - if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - found_youngest_primary = .true. - endif - endif - end do - if (associated(currentPatch)) then - ! the case where we've found a youngest primary patch - new_patch_primary%older => currentPatch - new_patch_primary%younger => currentPatch%younger - currentPatch%younger%older => new_patch_primary - currentPatch%younger => new_patch_primary - else - ! the case where we haven't, because the patches are all secondaary, - ! and are putting a primary patch at the oldest end of the - ! linked list (not sure how this could happen, but who knows...) - new_patch_primary%older => null() - new_patch_primary%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => new_patch_primary - currentSite%oldest_patch => new_patch_primary + if ( site_areadis_secondary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) + call fuse_cohorts(currentSite,new_patch_secondary, bc_in) + call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) + call sort_cohorts(new_patch_secondary) endif - else - ! the case where there are no secondary patches at the start of the linked list (prior logic) - new_patch_primary%older => currentPatch - new_patch_primary%younger => null() - currentPatch%younger => new_patch_primary - currentSite%youngest_patch => new_patch_primary - endif - endif - - ! insert first secondary at the start of the list - if ( site_areadis_secondary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - new_patch_secondary%older => currentPatch - new_patch_secondary%younger=> null() - currentPatch%younger => new_patch_secondary - currentSite%youngest_patch => new_patch_secondary - endif - - - ! sort out the cohorts, since some of them may be so small as to need removing. - ! the first call to terminate cohorts removes sparse number densities, - ! the second call removes for all other reasons (sparse culling must happen - ! before fusion) - - if ( site_areadis_primary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) - call fuse_cohorts(currentSite,new_patch_primary, bc_in) - call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) - call sort_cohorts(new_patch_primary) - endif - - if ( site_areadis_secondary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) - call fuse_cohorts(currentSite,new_patch_secondary, bc_in) - call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) - call sort_cohorts(new_patch_secondary) - endif - - endif !end new_patch area - - - call check_patch_area(currentSite) - call set_patchno(currentSite) - - end do disturbance_type_loop + + endif !end new_patch area + + + call check_patch_area(currentSite) + call set_patchno(currentSite) + + end do disturbance_type_loop end do nocomp_pft_loop From 485fe8373256d356db22e62f9b74f8f460af9128 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 23 Jun 2022 15:06:01 -0600 Subject: [PATCH 6/8] Revert "auto-indented spawn_patches" This reverts commit 53ca295e84cdb6b0a0fe99ba937161dc19d76987. --- biogeochem/EDPatchDynamicsMod.F90 | 1326 ++++++++++++++--------------- 1 file changed, 663 insertions(+), 663 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index dc3cfa373e..b3bab06bcb 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -436,717 +436,717 @@ subroutine spawn_patches( currentSite, bc_in) ! in the nocomp cases, since every patch has a PFT identity, it can only receive patch area from patches ! that have the same identity. In order to allow this, we have this very high level loop over nocomp PFTs - ! and only do the disturbance for any patches that have that nocomp PFT identity. + ! and only do the disturbance for any patches that have that nocomp PFT identity. ! If nocomp is not enabled, then this is not much of a loop, it only passes through once. nocomp_pft_loop: do i_nocomp_pft = min_nocomp_pft,max_nocomp_pft - disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES + disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES - ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. - currentPatch => currentSite%youngest_patch - - site_areadis_primary = 0.0_r8 - site_areadis_secondary = 0.0_r8 - - do while(associated(currentPatch)) + ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. + currentPatch => currentSite%youngest_patch - cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then + site_areadis_primary = 0.0_r8 + site_areadis_secondary = 0.0_r8 - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + do while(associated(currentPatch)) - if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then - write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate - call dump_patch(currentPatch) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then + + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - ! Only create new patches that have non-negligible amount of land - if((currentPatch%area*disturbance_rate) > nearzero ) then + if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then + write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate + call dump_patch(currentPatch) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! figure out whether the receiver patch for disturbance from this patch will be - ! primary or secondary land receiver patch is primary forest only if both the - ! donor patch is primary forest and the current disturbance type is not logging - if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then + ! Only create new patches that have non-negligible amount of land + if((currentPatch%area*disturbance_rate) > nearzero ) then + + ! figure out whether the receiver patch for disturbance from this patch will be + ! primary or secondary land receiver patch is primary forest only if both the + ! donor patch is primary forest and the current disturbance type is not logging + if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then + + site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate - site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate + ! track disturbance rates to output to history + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate - ! track disturbance rates to output to history - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate + ! track disturbance rates to output to history + if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + endif - ! track disturbance rates to output to history - if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - endif + endif + + end if + + end if cp_nocomp_matches_1_if + currentPatch => currentPatch%older + enddo ! end loop over patches. sum area disturbed for all patches. - endif + ! It is possible that no disturbance area was generated + if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then + + age = 0.0_r8 - end if + ! create two empty patches, to absorb newly disturbed primary and secondary forest area + ! first create patch to receive primary forest area + if ( site_areadis_primary .gt. nearzero ) then + allocate(new_patch_primary) - end if cp_nocomp_matches_1_if - currentPatch => currentPatch%older - enddo ! end loop over patches. sum area disturbed for all patches. - - ! It is possible that no disturbance area was generated - if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then - - age = 0.0_r8 - - ! create two empty patches, to absorb newly disturbed primary and secondary forest area - ! first create patch to receive primary forest area - if ( site_areadis_primary .gt. nearzero ) then - allocate(new_patch_primary) - - call create_patch(currentSite, new_patch_primary, age, & - site_areadis_primary, primaryforest, i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_primary%tallest => null() - new_patch_primary%shortest => null() + call create_patch(currentSite, new_patch_primary, age, & + site_areadis_primary, primaryforest, i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_primary%tallest => null() + new_patch_primary%shortest => null() - endif + endif - ! next create patch to receive secondary forest area - if ( site_areadis_secondary .gt. nearzero) then - allocate(new_patch_secondary) - call create_patch(currentSite, new_patch_secondary, age, & - site_areadis_secondary, secondaryforest,i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_secondary%tallest => null() - new_patch_secondary%shortest => null() + ! next create patch to receive secondary forest area + if ( site_areadis_secondary .gt. nearzero) then + allocate(new_patch_secondary) + call create_patch(currentSite, new_patch_secondary, age, & + site_areadis_secondary, secondaryforest,i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_secondary%tallest => null() + new_patch_secondary%shortest => null() - endif + endif + + ! loop round all the patches that contribute surviving indivduals and litter + ! pools to the new patch. We only loop the pre-existing patches, so + ! quit the loop if the current patch is either null, or matches the + ! two new pointers. - ! loop round all the patches that contribute surviving indivduals and litter - ! pools to the new patch. We only loop the pre-existing patches, so - ! quit the loop if the current patch is either null, or matches the - ! two new pointers. + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then + ! This is the amount of patch area that is disturbed, and donated by the donor + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + patch_site_areadis = currentPatch%area * disturbance_rate - ! This is the amount of patch area that is disturbed, and donated by the donor - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - patch_site_areadis = currentPatch%area * disturbance_rate + + if ( patch_site_areadis > nearzero ) then + + ! figure out whether the receiver patch for disturbance from this patch + ! will be primary or secondary land receiver patch is primary forest + ! only if both the donor patch is primary forest and the current + ! disturbance type is not logging + if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog)) then + new_patch => new_patch_primary + else + new_patch => new_patch_secondary + endif + + if(.not.associated(new_patch))then + write(fates_log(),*) 'Patch spawning has attempted to point to' + write(fates_log(),*) 'an un-allocated patch' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! for the case where the donating patch is secondary forest, if + ! the current disturbance from this patch is non-anthropogenic, + ! we need to average in the time-since-anthropogenic-disturbance + ! from the donor patch into that of the receiver patch + if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then + + new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & + currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) + endif + + + ! Transfer the litter existing already in the donor patch to the new patch + ! This call will only transfer non-burned litter to new patch + ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when + ! fire is not the current disturbance regime. - if ( patch_site_areadis > nearzero ) then + if(i_disturbance_type .ne. dtype_ifire) then + currentPatch%burnt_frac_litter(:) = 0._r8 + end if - ! figure out whether the receiver patch for disturbance from this patch - ! will be primary or secondary land receiver patch is primary forest - ! only if both the donor patch is primary forest and the current - ! disturbance type is not logging - if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog)) then - new_patch => new_patch_primary - else - new_patch => new_patch_secondary - endif + call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) - if(.not.associated(new_patch))then - write(fates_log(),*) 'Patch spawning has attempted to point to' - write(fates_log(),*) 'an un-allocated patch' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + ! Transfer in litter fluxes from plants in various contexts of death and destruction - ! for the case where the donating patch is secondary forest, if - ! the current disturbance from this patch is non-anthropogenic, - ! we need to average in the time-since-anthropogenic-disturbance - ! from the donor patch into that of the receiver patch - if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then + if(i_disturbance_type .eq. dtype_ilog) then + call logging_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + elseif(i_disturbance_type .eq. dtype_ifire) then + call fire_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + else + call mortality_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + endif - new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & - currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) + ! Copy any means or timers from the original patch to the new patch + ! These values will inherit all info from the original patch + ! -------------------------------------------------------------------------- + call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) + call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) + + + ! -------------------------------------------------------------------------- + ! The newly formed patch from disturbance (new_patch), has now been given + ! some litter from dead plants and pre-existing litter from the donor patches. + ! + ! Next, we loop through the cohorts in the donor patch, copy them with + ! area modified number density into the new-patch, and apply survivorship. + ! ------------------------------------------------------------------------- + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + + allocate(nc) + if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) + + ! Initialize the PARTEH object and point to the + ! correct boundary condition fields + nc%prt => null() + call InitPRTObject(nc%prt) + call InitPRTBoundaryConditions(nc) + + ! (Keeping as an example) + ! Allocate running mean functions + !allocate(nc%tveg_lpa) + !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + + call zero_cohort(nc) + + ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort + ! is the curent cohort that stays in the donor patch (currentPatch) + call copy_cohort(currentCohort, nc) + + !this is the case as the new patch probably doesn't have a closed canopy, and + ! even if it does, that will be sorted out in canopy_structure. + nc%canopy_layer = 1 + nc%canopy_layer_yesterday = 1._r8 + + sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) + total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c + + ! treefall mortality is the current disturbance + if(i_disturbance_type .eq. dtype_ifall) then + + if(currentCohort%canopy_layer == 1)then + + ! In the donor patch we are left with fewer trees because the area has decreased + ! the plant density for large trees does not actually decrease in the donor patch + ! because this is the part of the original patch where no trees have actually fallen + ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved + + currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & + min(1.0_r8,currentCohort%dmort * hlm_freq_day)) + + nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + + nc%cmort = nan ! The mortality diagnostics are set to nan + ! because the cohort should dissappear + nc%hmort = nan + nc%bmort = nan + nc%frmort = nan + nc%smort = nan + nc%asmort = nan + nc%lmort_direct = nan + nc%lmort_collateral = nan + nc%lmort_infra = nan + nc%l_degrad = nan + + else + ! small trees + if( int(prt_params%woody(currentCohort%pft)) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller and cohort counts + ! are absolute, reduce this number. + + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which + ! had been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to history + ! results in large numerical errors, on account of the sharply + ! reduced number densities. so instead pass this info via a + ! site-level diagnostic variable before reducing the number density. + + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * ED_val_understorey_death / hlm_freq_day + + + currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + (nc%n * ED_val_understorey_death / hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + ! Step 2: Apply survivor ship function based on the understory death fraction + ! remaining of understory plants of those that are knocked over + ! by the overstorey trees dying... + nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) + + ! since the donor patch split and sent a fraction of its members + ! to the new patch and a fraction to be preserved in itself, + ! when reporting diagnostic rates, we must carry over the mortality rates from + ! the donor that were applied before the patch split. Remember this is only + ! for diagnostics. But think of it this way, the rates are weighted by + ! number density in EDCLMLink, and the number density of this new patch is donated + ! so with the number density must come the effective mortality rates. + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + ! understory trees that might potentially be knocked over in the disturbance. + ! The existing (donor) patch should not have any impact mortality, it should + ! only lose cohorts due to the decrease in area. This is not mortality. + ! Besides, the current and newly created patch sum to unity + + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + else + ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + endif + endif + + ! Fire is the current disturbance + elseif (i_disturbance_type .eq. dtype_ifire ) then + + ! Number of members in the new patch, before we impose fire survivorship + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! loss of individuals from source patch due to area shrinking + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + levcan = currentCohort%canopy_layer + + if(levcan==ican_upper) then + + ! before changing number densities, track total rate of trees that died + ! due to fire, as well as from each fire mortality term + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + + else + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + end if + + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%cambial_mort / hlm_freq_day + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%crownfire_mort / hlm_freq_day + + ! loss of individual from fire in new patch. + nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + + ! Some of of the leaf mass from living plants has been + ! burned off. Here, we remove that mass, and + ! tally it in the flux we sent to the atmosphere + + if(int(prt_params%woody(currentCohort%pft)) == itrue)then + leaf_burn_frac = currentCohort%fraction_crown_burned + else + ! Grasses determine their fraction of leaves burned here - ! Transfer the litter existing already in the donor patch to the new patch - ! This call will only transfer non-burned litter to new patch - ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when - ! fire is not the current disturbance regime. - - if(i_disturbance_type .ne. dtype_ifire) then - currentPatch%burnt_frac_litter(:) = 0._r8 - end if - - call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) + leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) + endif + + ! Perform a check to make sure that spitfire gave + ! us reasonable mortality and burn fraction rates + + if( (leaf_burn_frac < 0._r8) .or. & + (leaf_burn_frac > 1._r8) .or. & + (currentCohort%fire_mort < 0._r8) .or. & + (currentCohort%fire_mort > 1._r8)) then + write(fates_log(),*) 'unexpected fire fractions' + write(fates_log(),*) prt_params%woody(currentCohort%pft) + write(fates_log(),*) leaf_burn_frac + write(fates_log(),*) currentCohort%fire_mort + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + do el = 1,num_elements + + leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) + + currentSite%mass_balance(el)%burn_flux_to_atm = & + currentSite%mass_balance(el)%burn_flux_to_atm + & + leaf_burn_frac * leaf_m * nc%n + end do - ! Transfer in litter fluxes from plants in various contexts of death and destruction + ! Here the mass is removed from the plant - if(i_disturbance_type .eq. dtype_ilog) then - call logging_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - elseif(i_disturbance_type .eq. dtype_ifire) then - call fire_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - else - call mortality_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - endif + call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) + currentCohort%fraction_crown_burned = 0.0_r8 + nc%fraction_crown_burned = 0.0_r8 - ! Copy any means or timers from the original patch to the new patch - ! These values will inherit all info from the original patch - ! -------------------------------------------------------------------------- - call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) - call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) - - ! -------------------------------------------------------------------------- - ! The newly formed patch from disturbance (new_patch), has now been given - ! some litter from dead plants and pre-existing litter from the donor patches. - ! - ! Next, we loop through the cohorts in the donor patch, copy them with - ! area modified number density into the new-patch, and apply survivorship. - ! ------------------------------------------------------------------------- - - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - - allocate(nc) - if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) - - ! Initialize the PARTEH object and point to the - ! correct boundary condition fields - nc%prt => null() - call InitPRTObject(nc%prt) - call InitPRTBoundaryConditions(nc) - - ! (Keeping as an example) - ! Allocate running mean functions - !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) - - call zero_cohort(nc) - - ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort - ! is the curent cohort that stays in the donor patch (currentPatch) - call copy_cohort(currentCohort, nc) - - !this is the case as the new patch probably doesn't have a closed canopy, and - ! even if it does, that will be sorted out in canopy_structure. - nc%canopy_layer = 1 - nc%canopy_layer_yesterday = 1._r8 - - sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) - store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) - total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c - - ! treefall mortality is the current disturbance - if(i_disturbance_type .eq. dtype_ifall) then - - if(currentCohort%canopy_layer == 1)then - - ! In the donor patch we are left with fewer trees because the area has decreased - ! the plant density for large trees does not actually decrease in the donor patch - ! because this is the part of the original patch where no trees have actually fallen - ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved - - currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & - min(1.0_r8,currentCohort%dmort * hlm_freq_day)) - - nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. - - nc%cmort = nan ! The mortality diagnostics are set to nan - ! because the cohort should dissappear - nc%hmort = nan - nc%bmort = nan - nc%frmort = nan - nc%smort = nan - nc%asmort = nan - nc%lmort_direct = nan - nc%lmort_collateral = nan - nc%lmort_infra = nan - nc%l_degrad = nan - - else - ! small trees - if( int(prt_params%woody(currentCohort%pft)) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller and cohort counts - ! are absolute, reduce this number. - - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which - ! had been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to history - ! results in large numerical errors, on account of the sharply - ! reduced number densities. so instead pass this info via a - ! site-level diagnostic variable before reducing the number density. - - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * ED_val_understorey_death / hlm_freq_day - - - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & - (nc%n * ED_val_understorey_death / hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - ! Step 2: Apply survivor ship function based on the understory death fraction - ! remaining of understory plants of those that are knocked over - ! by the overstorey trees dying... - nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) - - ! since the donor patch split and sent a fraction of its members - ! to the new patch and a fraction to be preserved in itself, - ! when reporting diagnostic rates, we must carry over the mortality rates from - ! the donor that were applied before the patch split. Remember this is only - ! for diagnostics. But think of it this way, the rates are weighted by - ! number density in EDCLMLink, and the number density of this new patch is donated - ! so with the number density must come the effective mortality rates. - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - ! understory trees that might potentially be knocked over in the disturbance. - ! The existing (donor) patch should not have any impact mortality, it should - ! only lose cohorts due to the decrease in area. This is not mortality. - ! Besides, the current and newly created patch sum to unity - - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - else - ! grass is not killed by mortality disturbance events. Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif - endif - - ! Fire is the current disturbance - elseif (i_disturbance_type .eq. dtype_ifire ) then - - ! Number of members in the new patch, before we impose fire survivorship - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! loss of individuals from source patch due to area shrinking - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - levcan = currentCohort%canopy_layer - - if(levcan==ican_upper) then - - ! before changing number densities, track total rate of trees that died - ! due to fire, as well as from each fire mortality term - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - - else - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - end if - - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%cambial_mort / hlm_freq_day - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%crownfire_mort / hlm_freq_day - - ! loss of individual from fire in new patch. - nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - - ! Some of of the leaf mass from living plants has been - ! burned off. Here, we remove that mass, and - ! tally it in the flux we sent to the atmosphere - - if(int(prt_params%woody(currentCohort%pft)) == itrue)then - leaf_burn_frac = currentCohort%fraction_crown_burned - else - - ! Grasses determine their fraction of leaves burned here - - leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) - endif - - ! Perform a check to make sure that spitfire gave - ! us reasonable mortality and burn fraction rates - - if( (leaf_burn_frac < 0._r8) .or. & - (leaf_burn_frac > 1._r8) .or. & - (currentCohort%fire_mort < 0._r8) .or. & - (currentCohort%fire_mort > 1._r8)) then - write(fates_log(),*) 'unexpected fire fractions' - write(fates_log(),*) prt_params%woody(currentCohort%pft) - write(fates_log(),*) leaf_burn_frac - write(fates_log(),*) currentCohort%fire_mort - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - do el = 1,num_elements - - leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) - - currentSite%mass_balance(el)%burn_flux_to_atm = & - currentSite%mass_balance(el)%burn_flux_to_atm + & - leaf_burn_frac * leaf_m * nc%n - end do - - ! Here the mass is removed from the plant - - call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) - currentCohort%fraction_crown_burned = 0.0_r8 - nc%fraction_crown_burned = 0.0_r8 - - - - ! Logging is the current disturbance - elseif (i_disturbance_type .eq. dtype_ilog ) then - - ! If this cohort is in the upper canopy. It generated - if(currentCohort%canopy_layer == 1)then - - ! calculate the survivorship of disturbed trees because non-harvested - nc%n = currentCohort%n * currentCohort%l_degrad - ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & - ! currentCohort%lmort_direct + currentCohort%lmort_collateral + - ! currentCohort%lmort_infra) ) * & - ! currentCohort%n * patch_site_areadis/currentPatch%area - - ! Reduce counts in the existing/donor patch according to the logging rate - currentCohort%n = currentCohort%n * & - (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & - currentCohort%lmort_collateral + & - currentCohort%lmort_infra + currentCohort%l_degrad))) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - - ! since these are the ones that weren't logged, - ! set the logging mortality rates as zero - nc%lmort_direct = 0._r8 - nc%lmort_collateral = 0._r8 - nc%lmort_infra = 0._r8 - - else - - ! What to do with cohorts in the understory of a logging generated - ! disturbance patch? - - if(int(prt_params%woody(currentCohort%pft)) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller - ! and cohort counts are absolute, reduce this number. - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which had - ! been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to - ! history results in large numerical errors, on account - ! of the sharply reduced number densities. so instead pass this info - ! via a site-level diagnostic variable before reducing - ! the number density. - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac / hlm_freq_day - - currentSite%imort_carbonflux = currentSite%imort_carbonflux + & - (nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac/ hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - - ! Step 2: Apply survivor ship function based on the understory death fraction - - ! remaining of understory plants of those that are knocked - ! over by the overstorey trees dying... - ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER - ! in the fatesparameter files - nc%n = nc%n * (1.0_r8 - & - (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) - - ! Step 3: Reduce the number count of cohorts in the - ! original/donor/non-disturbed patch to reflect the area change - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - else - - ! grass is not killed by mortality disturbance events. - ! Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - ! No grass impact mortality imposed on the newly created patch - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif ! is/is-not woody - - endif ! Select canopy layer - - else - write(fates_log(),*) 'unknown disturbance mode?' - write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if ! Select disturbance mode - - if (nc%n > 0.0_r8) then - storebigcohort => new_patch%tallest - storesmallcohort => new_patch%shortest - if(associated(new_patch%tallest))then - tnull = 0 - else - tnull = 1 - new_patch%tallest => nc - nc%taller => null() - endif - - if(associated(new_patch%shortest))then - snull = 0 - else - snull = 1 - new_patch%shortest => nc - nc%shorter => null() - endif - nc%patchptr => new_patch - call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & - tnull, snull, storebigcohort, storesmallcohort) - - new_patch%tallest => storebigcohort - new_patch%shortest => storesmallcohort - else - - ! Get rid of the new temporary cohort - call DeallocateCohort(nc) - deallocate(nc) - - endif - - currentCohort => currentCohort%taller - enddo ! currentCohort - call sort_cohorts(currentPatch) - - !update area of donor patch - oldarea = currentPatch%area - currentPatch%area = currentPatch%area - patch_site_areadis - - ! for all disturbance rates that haven't been resolved yet, increase their amount so that - ! they are the same amount of gridcell-scale disturbance relative to the original patch size - if (i_disturbance_type .ne. N_DIST_TYPES) then - do i_dist2 = i_disturbance_type+1,N_DIST_TYPES - currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & - * oldarea / currentPatch%area - end do - end if - - ! sort out the cohorts, since some of them may be so small as to need removing. - ! the first call to terminate cohorts removes sparse number densities, - ! the second call removes for all other reasons (sparse culling must happen - ! before fusion) - call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) - call fuse_cohorts(currentSite,currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) - call sort_cohorts(currentPatch) - - end if ! if ( new_patch%area > nearzero ) then - - end if cp_nocomp_matches_2_if - currentPatch => currentPatch%younger - - enddo ! currentPatch patch loop. - - !*************************/ - !** INSERT NEW PATCH(ES) INTO LINKED LIST - !*************************/ - - if ( site_areadis_primary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - ! insert new youngest primary patch after all the secondary patches, if there are any. - ! this requires first finding the current youngest primary to insert the new one ahead of - if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then - found_youngest_primary = .false. - do while(associated(currentPatch) .and. .not. found_youngest_primary) - currentPatch => currentPatch%older - if (associated(currentPatch)) then - if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - found_youngest_primary = .true. - endif - endif - end do - if (associated(currentPatch)) then - ! the case where we've found a youngest primary patch - new_patch_primary%older => currentPatch - new_patch_primary%younger => currentPatch%younger - currentPatch%younger%older => new_patch_primary - currentPatch%younger => new_patch_primary + ! Logging is the current disturbance + elseif (i_disturbance_type .eq. dtype_ilog ) then + + ! If this cohort is in the upper canopy. It generated + if(currentCohort%canopy_layer == 1)then + + ! calculate the survivorship of disturbed trees because non-harvested + nc%n = currentCohort%n * currentCohort%l_degrad + ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & + ! currentCohort%lmort_direct + currentCohort%lmort_collateral + + ! currentCohort%lmort_infra) ) * & + ! currentCohort%n * patch_site_areadis/currentPatch%area + + ! Reduce counts in the existing/donor patch according to the logging rate + currentCohort%n = currentCohort%n * & + (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra + currentCohort%l_degrad))) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + + ! since these are the ones that weren't logged, + ! set the logging mortality rates as zero + nc%lmort_direct = 0._r8 + nc%lmort_collateral = 0._r8 + nc%lmort_infra = 0._r8 + else - ! the case where we haven't, because the patches are all secondaary, - ! and are putting a primary patch at the oldest end of the - ! linked list (not sure how this could happen, but who knows...) - new_patch_primary%older => null() - new_patch_primary%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => new_patch_primary - currentSite%oldest_patch => new_patch_primary + + ! What to do with cohorts in the understory of a logging generated + ! disturbance patch? + + if(int(prt_params%woody(currentCohort%pft)) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller + ! and cohort counts are absolute, reduce this number. + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which had + ! been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to + ! history results in large numerical errors, on account + ! of the sharply reduced number densities. so instead pass this info + ! via a site-level diagnostic variable before reducing + ! the number density. + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac / hlm_freq_day + + currentSite%imort_carbonflux = currentSite%imort_carbonflux + & + (nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac/ hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + + ! Step 2: Apply survivor ship function based on the understory death fraction + + ! remaining of understory plants of those that are knocked + ! over by the overstorey trees dying... + ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER + ! in the fatesparameter files + nc%n = nc%n * (1.0_r8 - & + (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) + + ! Step 3: Reduce the number count of cohorts in the + ! original/donor/non-disturbed patch to reflect the area change + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + else + + ! grass is not killed by mortality disturbance events. + ! Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + ! No grass impact mortality imposed on the newly created patch + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif ! is/is-not woody + + endif ! Select canopy layer + + else + write(fates_log(),*) 'unknown disturbance mode?' + write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Select disturbance mode + + if (nc%n > 0.0_r8) then + storebigcohort => new_patch%tallest + storesmallcohort => new_patch%shortest + if(associated(new_patch%tallest))then + tnull = 0 + else + tnull = 1 + new_patch%tallest => nc + nc%taller => null() endif + + if(associated(new_patch%shortest))then + snull = 0 + else + snull = 1 + new_patch%shortest => nc + nc%shorter => null() + endif + nc%patchptr => new_patch + call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & + tnull, snull, storebigcohort, storesmallcohort) + + new_patch%tallest => storebigcohort + new_patch%shortest => storesmallcohort else - ! the case where there are no secondary patches at the start of the linked list (prior logic) - new_patch_primary%older => currentPatch - new_patch_primary%younger => null() - currentPatch%younger => new_patch_primary - currentSite%youngest_patch => new_patch_primary + + ! Get rid of the new temporary cohort + call DeallocateCohort(nc) + deallocate(nc) + endif - endif - - ! insert first secondary at the start of the list - if ( site_areadis_secondary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - new_patch_secondary%older => currentPatch - new_patch_secondary%younger=> null() - currentPatch%younger => new_patch_secondary - currentSite%youngest_patch => new_patch_secondary - endif - + + currentCohort => currentCohort%taller + enddo ! currentCohort + call sort_cohorts(currentPatch) + + !update area of donor patch + oldarea = currentPatch%area + currentPatch%area = currentPatch%area - patch_site_areadis + + ! for all disturbance rates that haven't been resolved yet, increase their amount so that + ! they are the same amount of gridcell-scale disturbance relative to the original patch size + if (i_disturbance_type .ne. N_DIST_TYPES) then + do i_dist2 = i_disturbance_type+1,N_DIST_TYPES + currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & + * oldarea / currentPatch%area + end do + end if - ! sort out the cohorts, since some of them may be so small as to need removing. + ! sort out the cohorts, since some of them may be so small as to need removing. ! the first call to terminate cohorts removes sparse number densities, ! the second call removes for all other reasons (sparse culling must happen ! before fusion) + call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) + call fuse_cohorts(currentSite,currentPatch, bc_in) + call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) + call sort_cohorts(currentPatch) - if ( site_areadis_primary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) - call fuse_cohorts(currentSite,new_patch_primary, bc_in) - call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) - call sort_cohorts(new_patch_primary) - endif + end if ! if ( new_patch%area > nearzero ) then + + end if cp_nocomp_matches_2_if + currentPatch => currentPatch%younger + + enddo ! currentPatch patch loop. - if ( site_areadis_secondary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) - call fuse_cohorts(currentSite,new_patch_secondary, bc_in) - call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) - call sort_cohorts(new_patch_secondary) + !*************************/ + !** INSERT NEW PATCH(ES) INTO LINKED LIST + !*************************/ + + if ( site_areadis_primary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + ! insert new youngest primary patch after all the secondary patches, if there are any. + ! this requires first finding the current youngest primary to insert the new one ahead of + if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then + found_youngest_primary = .false. + do while(associated(currentPatch) .and. .not. found_youngest_primary) + currentPatch => currentPatch%older + if (associated(currentPatch)) then + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + found_youngest_primary = .true. + endif + endif + end do + if (associated(currentPatch)) then + ! the case where we've found a youngest primary patch + new_patch_primary%older => currentPatch + new_patch_primary%younger => currentPatch%younger + currentPatch%younger%older => new_patch_primary + currentPatch%younger => new_patch_primary + else + ! the case where we haven't, because the patches are all secondaary, + ! and are putting a primary patch at the oldest end of the + ! linked list (not sure how this could happen, but who knows...) + new_patch_primary%older => null() + new_patch_primary%younger => currentSite%oldest_patch + currentSite%oldest_patch%older => new_patch_primary + currentSite%oldest_patch => new_patch_primary endif - - endif !end new_patch area - - - call check_patch_area(currentSite) - call set_patchno(currentSite) - - end do disturbance_type_loop + else + ! the case where there are no secondary patches at the start of the linked list (prior logic) + new_patch_primary%older => currentPatch + new_patch_primary%younger => null() + currentPatch%younger => new_patch_primary + currentSite%youngest_patch => new_patch_primary + endif + endif + + ! insert first secondary at the start of the list + if ( site_areadis_secondary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + new_patch_secondary%older => currentPatch + new_patch_secondary%younger=> null() + currentPatch%younger => new_patch_secondary + currentSite%youngest_patch => new_patch_secondary + endif + + + ! sort out the cohorts, since some of them may be so small as to need removing. + ! the first call to terminate cohorts removes sparse number densities, + ! the second call removes for all other reasons (sparse culling must happen + ! before fusion) + + if ( site_areadis_primary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) + call fuse_cohorts(currentSite,new_patch_primary, bc_in) + call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) + call sort_cohorts(new_patch_primary) + endif + + if ( site_areadis_secondary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) + call fuse_cohorts(currentSite,new_patch_secondary, bc_in) + call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) + call sort_cohorts(new_patch_secondary) + endif + + endif !end new_patch area + + + call check_patch_area(currentSite) + call set_patchno(currentSite) + + end do disturbance_type_loop end do nocomp_pft_loop From bd994a1c548d9c13e1bcc9a40657de2f670515c4 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 23 Jun 2022 15:09:28 -0600 Subject: [PATCH 7/8] re-auto-indented spawn_patches (after merging main branch) --- biogeochem/EDPatchDynamicsMod.F90 | 1334 ++++++++++++++--------------- 1 file changed, 667 insertions(+), 667 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 584494de8b..d837ed737a 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -437,721 +437,721 @@ subroutine spawn_patches( currentSite, bc_in) ! in the nocomp cases, since every patch has a PFT identity, it can only receive patch area from patches ! that have the same identity. In order to allow this, we have this very high level loop over nocomp PFTs - ! and only do the disturbance for any patches that have that nocomp PFT identity. + ! and only do the disturbance for any patches that have that nocomp PFT identity. ! If nocomp is not enabled, then this is not much of a loop, it only passes through once. nocomp_pft_loop: do i_nocomp_pft = min_nocomp_pft,max_nocomp_pft - disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES + disturbance_type_loop: do i_disturbance_type = 1,N_DIST_TYPES - ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. - currentPatch => currentSite%youngest_patch + ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. + currentPatch => currentSite%youngest_patch - site_areadis_primary = 0.0_r8 - site_areadis_secondary = 0.0_r8 + site_areadis_primary = 0.0_r8 + site_areadis_secondary = 0.0_r8 - do while(associated(currentPatch)) + do while(associated(currentPatch)) - cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + cp_nocomp_matches_1_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then - write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate - call dump_patch(currentPatch) - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - ! Only create new patches that have non-negligible amount of land - if((currentPatch%area*disturbance_rate) > nearzero ) then - - ! figure out whether the receiver patch for disturbance from this patch will be - ! primary or secondary land receiver patch is primary forest only if both the - ! donor patch is primary forest and the current disturbance type is not logging - if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then - - site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate + if(disturbance_rate > (1.0_r8 + rsnbl_math_prec)) then + write(fates_log(),*) 'patch disturbance rate > 1 ?',disturbance_rate + call dump_patch(currentPatch) + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! track disturbance rates to output to history - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate + ! Only create new patches that have non-negligible amount of land + if((currentPatch%area*disturbance_rate) > nearzero ) then - ! track disturbance rates to output to history - if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - else - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & - currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & - currentPatch%area * disturbance_rate * AREA_INV - endif + ! figure out whether the receiver patch for disturbance from this patch will be + ! primary or secondary land receiver patch is primary forest only if both the + ! donor patch is primary forest and the current disturbance type is not logging + if ( currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then - endif - - end if - - end if cp_nocomp_matches_1_if - currentPatch => currentPatch%older - enddo ! end loop over patches. sum area disturbed for all patches. + site_areadis_primary = site_areadis_primary + currentPatch%area * disturbance_rate - ! It is possible that no disturbance area was generated - if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then - - age = 0.0_r8 + ! track disturbance rates to output to history + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_primary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + site_areadis_secondary = site_areadis_secondary + currentPatch%area * disturbance_rate - ! create two empty patches, to absorb newly disturbed primary and secondary forest area - ! first create patch to receive primary forest area - if ( site_areadis_primary .gt. nearzero ) then - allocate(new_patch_primary) + ! track disturbance rates to output to history + if (currentPatch%anthro_disturbance_label .eq. secondaryforest) then + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_secondary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + else + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) = & + currentSite%disturbance_rates_primary_to_secondary(i_disturbance_type) + & + currentPatch%area * disturbance_rate * AREA_INV + endif - call create_patch(currentSite, new_patch_primary, age, & - site_areadis_primary, primaryforest, i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_primary%tallest => null() - new_patch_primary%shortest => null() + endif - endif + end if - ! next create patch to receive secondary forest area - if ( site_areadis_secondary .gt. nearzero) then - allocate(new_patch_secondary) - call create_patch(currentSite, new_patch_secondary, age, & - site_areadis_secondary, secondaryforest,i_nocomp_pft) - - ! Initialize the litter pools to zero, these - ! pools will be populated by looping over the existing patches - ! and transfering in mass - do el=1,num_elements - call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & - init_root_fines=0._r8, & - init_ag_cwd=0._r8, & - init_bg_cwd=0._r8, & - init_seed=0._r8, & - init_seed_germ=0._r8) - end do - new_patch_secondary%tallest => null() - new_patch_secondary%shortest => null() + end if cp_nocomp_matches_1_if + currentPatch => currentPatch%older + enddo ! end loop over patches. sum area disturbed for all patches. + + ! It is possible that no disturbance area was generated + if ( (site_areadis_primary + site_areadis_secondary) > nearzero) then + + age = 0.0_r8 + + ! create two empty patches, to absorb newly disturbed primary and secondary forest area + ! first create patch to receive primary forest area + if ( site_areadis_primary .gt. nearzero ) then + allocate(new_patch_primary) + + call create_patch(currentSite, new_patch_primary, age, & + site_areadis_primary, primaryforest, i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_primary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_primary%tallest => null() + new_patch_primary%shortest => null() - endif - - ! loop round all the patches that contribute surviving indivduals and litter - ! pools to the new patch. We only loop the pre-existing patches, so - ! quit the loop if the current patch is either null, or matches the - ! two new pointers. + endif - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) + ! next create patch to receive secondary forest area + if ( site_areadis_secondary .gt. nearzero) then + allocate(new_patch_secondary) + call create_patch(currentSite, new_patch_secondary, age, & + site_areadis_secondary, secondaryforest,i_nocomp_pft) + + ! Initialize the litter pools to zero, these + ! pools will be populated by looping over the existing patches + ! and transfering in mass + do el=1,num_elements + call new_patch_secondary%litter(el)%InitConditions(init_leaf_fines=0._r8, & + init_root_fines=0._r8, & + init_ag_cwd=0._r8, & + init_bg_cwd=0._r8, & + init_seed=0._r8, & + init_seed_germ=0._r8) + end do + new_patch_secondary%tallest => null() + new_patch_secondary%shortest => null() - cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & - currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then + endif - ! This is the amount of patch area that is disturbed, and donated by the donor - disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) - patch_site_areadis = currentPatch%area * disturbance_rate + ! loop round all the patches that contribute surviving indivduals and litter + ! pools to the new patch. We only loop the pre-existing patches, so + ! quit the loop if the current patch is either null, or matches the + ! two new pointers. - - if ( patch_site_areadis > nearzero ) then - - ! figure out whether the receiver patch for disturbance from this patch - ! will be primary or secondary land receiver patch is primary forest - ! only if both the donor patch is primary forest and the current - ! disturbance type is not logging - if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & - (i_disturbance_type .ne. dtype_ilog)) then - new_patch => new_patch_primary - else - new_patch => new_patch_secondary - endif - - if(.not.associated(new_patch))then - write(fates_log(),*) 'Patch spawning has attempted to point to' - write(fates_log(),*) 'an un-allocated patch' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - ! for the case where the donating patch is secondary forest, if - ! the current disturbance from this patch is non-anthropogenic, - ! we need to average in the time-since-anthropogenic-disturbance - ! from the donor patch into that of the receiver patch - if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & - (i_disturbance_type .ne. dtype_ilog) ) then - - new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & - currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) + currentPatch => currentSite%oldest_patch + do while(associated(currentPatch)) - endif - - - ! Transfer the litter existing already in the donor patch to the new patch - ! This call will only transfer non-burned litter to new patch - ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when - ! fire is not the current disturbance regime. + cp_nocomp_matches_2_if: if ( hlm_use_nocomp .eq. ifalse .or. & + currentPatch%nocomp_pft_label .eq. i_nocomp_pft ) then - if(i_disturbance_type .ne. dtype_ifire) then - currentPatch%burnt_frac_litter(:) = 0._r8 - end if + ! This is the amount of patch area that is disturbed, and donated by the donor + disturbance_rate = currentPatch%disturbance_rates(i_disturbance_type) + patch_site_areadis = currentPatch%area * disturbance_rate - call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) - ! Transfer in litter fluxes from plants in various contexts of death and destruction + if ( patch_site_areadis > nearzero ) then - if(i_disturbance_type .eq. dtype_ilog) then - call logging_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - elseif(i_disturbance_type .eq. dtype_ifire) then - call fire_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - else - call mortality_litter_fluxes(currentSite, currentPatch, & - new_patch, patch_site_areadis,bc_in) - endif + ! figure out whether the receiver patch for disturbance from this patch + ! will be primary or secondary land receiver patch is primary forest + ! only if both the donor patch is primary forest and the current + ! disturbance type is not logging + if (currentPatch%anthro_disturbance_label .eq. primaryforest .and. & + (i_disturbance_type .ne. dtype_ilog)) then + new_patch => new_patch_primary + else + new_patch => new_patch_secondary + endif + if(.not.associated(new_patch))then + write(fates_log(),*) 'Patch spawning has attempted to point to' + write(fates_log(),*) 'an un-allocated patch' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! Copy any means or timers from the original patch to the new patch - ! These values will inherit all info from the original patch - ! -------------------------------------------------------------------------- - call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) - call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) - - - ! -------------------------------------------------------------------------- - ! The newly formed patch from disturbance (new_patch), has now been given - ! some litter from dead plants and pre-existing litter from the donor patches. - ! - ! Next, we loop through the cohorts in the donor patch, copy them with - ! area modified number density into the new-patch, and apply survivorship. - ! ------------------------------------------------------------------------- - - currentCohort => currentPatch%shortest - do while(associated(currentCohort)) - - allocate(nc) - if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) - - ! Initialize the PARTEH object and point to the - ! correct boundary condition fields - nc%prt => null() - call InitPRTObject(nc%prt) - call InitPRTBoundaryConditions(nc) - - ! (Keeping as an example) - ! Allocate running mean functions - !allocate(nc%tveg_lpa) - !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) - - call zero_cohort(nc) - - ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort - ! is the curent cohort that stays in the donor patch (currentPatch) - call copy_cohort(currentCohort, nc) - - !this is the case as the new patch probably doesn't have a closed canopy, and - ! even if it does, that will be sorted out in canopy_structure. - nc%canopy_layer = 1 - nc%canopy_layer_yesterday = 1._r8 - - sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) - fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) - store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) - total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c - - ! treefall mortality is the current disturbance - if(i_disturbance_type .eq. dtype_ifall) then - - if(currentCohort%canopy_layer == 1)then - - ! In the donor patch we are left with fewer trees because the area has decreased - ! the plant density for large trees does not actually decrease in the donor patch - ! because this is the part of the original patch where no trees have actually fallen - ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved - - currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & - min(1.0_r8,currentCohort%dmort * hlm_freq_day)) - - nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. - - nc%cmort = nan ! The mortality diagnostics are set to nan - ! because the cohort should dissappear - nc%hmort = nan - nc%bmort = nan - nc%frmort = nan - nc%smort = nan - nc%asmort = nan - nc%lmort_direct = nan - nc%lmort_collateral = nan - nc%lmort_infra = nan - nc%l_degrad = nan - - else - ! small trees - if( prt_params%woody(currentCohort%pft) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller and cohort counts - ! are absolute, reduce this number. - - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which - ! had been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to history - ! results in large numerical errors, on account of the sharply - ! reduced number densities. so instead pass this info via a - ! site-level diagnostic variable before reducing the number density. - - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * ED_val_understorey_death / hlm_freq_day - - - currentSite%imort_carbonflux(currentCohort%pft) = & - currentSite%imort_carbonflux(currentCohort%pft) + & - (nc%n * ED_val_understorey_death / hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - ! Step 2: Apply survivor ship function based on the understory death fraction - ! remaining of understory plants of those that are knocked over - ! by the overstorey trees dying... - nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) - - ! since the donor patch split and sent a fraction of its members - ! to the new patch and a fraction to be preserved in itself, - ! when reporting diagnostic rates, we must carry over the mortality rates from - ! the donor that were applied before the patch split. Remember this is only - ! for diagnostics. But think of it this way, the rates are weighted by - ! number density in EDCLMLink, and the number density of this new patch is donated - ! so with the number density must come the effective mortality rates. - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - ! understory trees that might potentially be knocked over in the disturbance. - ! The existing (donor) patch should not have any impact mortality, it should - ! only lose cohorts due to the decrease in area. This is not mortality. - ! Besides, the current and newly created patch sum to unity - - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - else - ! grass is not killed by mortality disturbance events. Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif - endif - - ! Fire is the current disturbance - elseif (i_disturbance_type .eq. dtype_ifire ) then - - ! Number of members in the new patch, before we impose fire survivorship - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! loss of individuals from source patch due to area shrinking - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - levcan = currentCohort%canopy_layer - - if(levcan==ican_upper) then - - ! before changing number densities, track total rate of trees that died - ! due to fire, as well as from each fire mortality term - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_canopy(currentCohort%pft) = & - currentSite%fmort_carbonflux_canopy(currentCohort%pft) + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - - else - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%fire_mort / hlm_freq_day - - currentSite%fmort_carbonflux_ustory(currentCohort%pft) = & - currentSite%fmort_carbonflux_ustory(currentCohort%pft) + & - (nc%n * currentCohort%fire_mort) * & - total_c * g_per_kg * days_per_sec * ha_per_m2 - end if - - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%cambial_mort / hlm_freq_day - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & - currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentCohort%crownfire_mort / hlm_freq_day - - ! loss of individual from fire in new patch. - nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - - ! Some of of the leaf mass from living plants has been - ! burned off. Here, we remove that mass, and - ! tally it in the flux we sent to the atmosphere - - if(prt_params%woody(currentCohort%pft) == itrue)then - leaf_burn_frac = currentCohort%fraction_crown_burned - else + ! for the case where the donating patch is secondary forest, if + ! the current disturbance from this patch is non-anthropogenic, + ! we need to average in the time-since-anthropogenic-disturbance + ! from the donor patch into that of the receiver patch + if ( currentPatch%anthro_disturbance_label .eq. secondaryforest .and. & + (i_disturbance_type .ne. dtype_ilog) ) then - ! Grasses determine their fraction of leaves burned here + new_patch%age_since_anthro_disturbance = new_patch%age_since_anthro_disturbance + & + currentPatch%age_since_anthro_disturbance * (patch_site_areadis / site_areadis_secondary) - leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) - endif - - ! Perform a check to make sure that spitfire gave - ! us reasonable mortality and burn fraction rates - - if( (leaf_burn_frac < 0._r8) .or. & - (leaf_burn_frac > 1._r8) .or. & - (currentCohort%fire_mort < 0._r8) .or. & - (currentCohort%fire_mort > 1._r8)) then - write(fates_log(),*) 'unexpected fire fractions' - write(fates_log(),*) prt_params%woody(currentCohort%pft) - write(fates_log(),*) leaf_burn_frac - write(fates_log(),*) currentCohort%fire_mort - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - do el = 1,num_elements - - leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) - - currentSite%mass_balance(el)%burn_flux_to_atm = & - currentSite%mass_balance(el)%burn_flux_to_atm + & - leaf_burn_frac * leaf_m * nc%n - end do + endif - ! Here the mass is removed from the plant - call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) - currentCohort%fraction_crown_burned = 0.0_r8 - nc%fraction_crown_burned = 0.0_r8 + ! Transfer the litter existing already in the donor patch to the new patch + ! This call will only transfer non-burned litter to new patch + ! and burned litter to atmosphere. Thus it is important to zero burnt_frac_litter when + ! fire is not the current disturbance regime. + if(i_disturbance_type .ne. dtype_ifire) then + currentPatch%burnt_frac_litter(:) = 0._r8 + end if + call TransLitterNewPatch( currentSite, currentPatch, new_patch, patch_site_areadis) - ! Logging is the current disturbance - elseif (i_disturbance_type .eq. dtype_ilog ) then - - ! If this cohort is in the upper canopy. It generated - if(currentCohort%canopy_layer == 1)then - - ! calculate the survivorship of disturbed trees because non-harvested - nc%n = currentCohort%n * currentCohort%l_degrad - ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & - ! currentCohort%lmort_direct + currentCohort%lmort_collateral + - ! currentCohort%lmort_infra) ) * & - ! currentCohort%n * patch_site_areadis/currentPatch%area - - ! Reduce counts in the existing/donor patch according to the logging rate - currentCohort%n = currentCohort%n * & - (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & - currentCohort%lmort_collateral + & - currentCohort%lmort_infra + currentCohort%l_degrad))) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - - ! since these are the ones that weren't logged, - ! set the logging mortality rates as zero - nc%lmort_direct = 0._r8 - nc%lmort_collateral = 0._r8 - nc%lmort_infra = 0._r8 - - else - - ! What to do with cohorts in the understory of a logging generated - ! disturbance patch? - - if(prt_params%woody(currentCohort%pft) == itrue)then - - - ! Survivorship of undestory woody plants. Two step process. - ! Step 1: Reduce current number of plants to reflect the - ! change in area. - ! The number density per square are doesn't change, - ! but since the patch is smaller - ! and cohort counts are absolute, reduce this number. - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! because the mortality rate due to impact for the cohorts which had - ! been in the understory and are now in the newly- - ! disturbed patch is very high, passing the imort directly to - ! history results in large numerical errors, on account - ! of the sharply reduced number densities. so instead pass this info - ! via a site-level diagnostic variable before reducing - ! the number density. - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & - currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & - nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac / hlm_freq_day - - currentSite%imort_carbonflux(currentCohort%pft) = & - currentSite%imort_carbonflux(currentCohort%pft) + & - (nc%n * currentPatch%fract_ldist_not_harvested * & - logging_coll_under_frac/ hlm_freq_day ) * & - total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 - - - ! Step 2: Apply survivor ship function based on the understory death fraction - - ! remaining of understory plants of those that are knocked - ! over by the overstorey trees dying... - ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER - ! in the fatesparameter files - nc%n = nc%n * (1.0_r8 - & - (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) - - ! Step 3: Reduce the number count of cohorts in the - ! original/donor/non-disturbed patch to reflect the area change - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - + ! Transfer in litter fluxes from plants in various contexts of death and destruction + + if(i_disturbance_type .eq. dtype_ilog) then + call logging_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + elseif(i_disturbance_type .eq. dtype_ifire) then + call fire_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) else - - ! grass is not killed by mortality disturbance events. - ! Just move it into the new patch area. - ! Just split the grass into the existing and new patch structures - nc%n = currentCohort%n * patch_site_areadis/currentPatch%area - - ! Those remaining in the existing - currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) - - ! No grass impact mortality imposed on the newly created patch - nc%cmort = currentCohort%cmort - nc%hmort = currentCohort%hmort - nc%bmort = currentCohort%bmort - nc%frmort = currentCohort%frmort - nc%smort = currentCohort%smort - nc%asmort = currentCohort%asmort - nc%dmort = currentCohort%dmort - nc%lmort_direct = currentCohort%lmort_direct - nc%lmort_collateral = currentCohort%lmort_collateral - nc%lmort_infra = currentCohort%lmort_infra - - endif ! is/is-not woody - - endif ! Select canopy layer - - else - write(fates_log(),*) 'unknown disturbance mode?' - write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if ! Select disturbance mode - - if (nc%n > 0.0_r8) then - storebigcohort => new_patch%tallest - storesmallcohort => new_patch%shortest - if(associated(new_patch%tallest))then - tnull = 0 - else - tnull = 1 - new_patch%tallest => nc - nc%taller => null() - endif - - if(associated(new_patch%shortest))then - snull = 0 + call mortality_litter_fluxes(currentSite, currentPatch, & + new_patch, patch_site_areadis,bc_in) + endif + + + ! Copy any means or timers from the original patch to the new patch + ! These values will inherit all info from the original patch + ! -------------------------------------------------------------------------- + call new_patch%tveg24%CopyFromDonor(currentPatch%tveg24) + call new_patch%tveg_lpa%CopyFromDonor(currentPatch%tveg_lpa) + + + ! -------------------------------------------------------------------------- + ! The newly formed patch from disturbance (new_patch), has now been given + ! some litter from dead plants and pre-existing litter from the donor patches. + ! + ! Next, we loop through the cohorts in the donor patch, copy them with + ! area modified number density into the new-patch, and apply survivorship. + ! ------------------------------------------------------------------------- + + currentCohort => currentPatch%shortest + do while(associated(currentCohort)) + + allocate(nc) + if(hlm_use_planthydro.eq.itrue) call InitHydrCohort(CurrentSite,nc) + + ! Initialize the PARTEH object and point to the + ! correct boundary condition fields + nc%prt => null() + call InitPRTObject(nc%prt) + call InitPRTBoundaryConditions(nc) + + ! (Keeping as an example) + ! Allocate running mean functions + !allocate(nc%tveg_lpa) + !call nc%tveg_lpa%InitRMean(ema_lpa,init_value=new_patch%tveg_lpa%GetMean()) + + call zero_cohort(nc) + + ! nc is the new cohort that goes in the disturbed patch (new_patch)... currentCohort + ! is the curent cohort that stays in the donor patch (currentPatch) + call copy_cohort(currentCohort, nc) + + !this is the case as the new patch probably doesn't have a closed canopy, and + ! even if it does, that will be sorted out in canopy_structure. + nc%canopy_layer = 1 + nc%canopy_layer_yesterday = 1._r8 + + sapw_c = currentCohort%prt%GetState(sapw_organ, all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ, all_carbon_elements) + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ, all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ, all_carbon_elements) + total_c = sapw_c + struct_c + leaf_c + fnrt_c + store_c + + ! treefall mortality is the current disturbance + if(i_disturbance_type .eq. dtype_ifall) then + + if(currentCohort%canopy_layer == 1)then + + ! In the donor patch we are left with fewer trees because the area has decreased + ! the plant density for large trees does not actually decrease in the donor patch + ! because this is the part of the original patch where no trees have actually fallen + ! The diagnostic cmort,bmort,hmort, and frmort rates have already been saved + + currentCohort%n = currentCohort%n * (1.0_r8 - fates_mortality_disturbance_fraction * & + min(1.0_r8,currentCohort%dmort * hlm_freq_day)) + + nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. + + nc%cmort = nan ! The mortality diagnostics are set to nan + ! because the cohort should dissappear + nc%hmort = nan + nc%bmort = nan + nc%frmort = nan + nc%smort = nan + nc%asmort = nan + nc%lmort_direct = nan + nc%lmort_collateral = nan + nc%lmort_infra = nan + nc%l_degrad = nan + + else + ! small trees + if( prt_params%woody(currentCohort%pft) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller and cohort counts + ! are absolute, reduce this number. + + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which + ! had been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to history + ! results in large numerical errors, on account of the sharply + ! reduced number densities. so instead pass this info via a + ! site-level diagnostic variable before reducing the number density. + + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * ED_val_understorey_death / hlm_freq_day + + + currentSite%imort_carbonflux(currentCohort%pft) = & + currentSite%imort_carbonflux(currentCohort%pft) + & + (nc%n * ED_val_understorey_death / hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + ! Step 2: Apply survivor ship function based on the understory death fraction + ! remaining of understory plants of those that are knocked over + ! by the overstorey trees dying... + nc%n = nc%n * (1.0_r8 - ED_val_understorey_death) + + ! since the donor patch split and sent a fraction of its members + ! to the new patch and a fraction to be preserved in itself, + ! when reporting diagnostic rates, we must carry over the mortality rates from + ! the donor that were applied before the patch split. Remember this is only + ! for diagnostics. But think of it this way, the rates are weighted by + ! number density in EDCLMLink, and the number density of this new patch is donated + ! so with the number density must come the effective mortality rates. + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + ! understory trees that might potentially be knocked over in the disturbance. + ! The existing (donor) patch should not have any impact mortality, it should + ! only lose cohorts due to the decrease in area. This is not mortality. + ! Besides, the current and newly created patch sum to unity + + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + else + ! grass is not killed by mortality disturbance events. Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif + endif + + ! Fire is the current disturbance + elseif (i_disturbance_type .eq. dtype_ifire ) then + + ! Number of members in the new patch, before we impose fire survivorship + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! loss of individuals from source patch due to area shrinking + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + levcan = currentCohort%canopy_layer + + if(levcan==ican_upper) then + + ! before changing number densities, track total rate of trees that died + ! due to fire, as well as from each fire mortality term + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_canopy(currentCohort%pft) = & + currentSite%fmort_carbonflux_canopy(currentCohort%pft) + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + + else + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%fire_mort / hlm_freq_day + + currentSite%fmort_carbonflux_ustory(currentCohort%pft) = & + currentSite%fmort_carbonflux_ustory(currentCohort%pft) + & + (nc%n * currentCohort%fire_mort) * & + total_c * g_per_kg * days_per_sec * ha_per_m2 + end if + + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%cambial_mort / hlm_freq_day + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = & + currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentCohort%crownfire_mort / hlm_freq_day + + ! loss of individual from fire in new patch. + nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + + ! Some of of the leaf mass from living plants has been + ! burned off. Here, we remove that mass, and + ! tally it in the flux we sent to the atmosphere + + if(prt_params%woody(currentCohort%pft) == itrue)then + leaf_burn_frac = currentCohort%fraction_crown_burned + else + + ! Grasses determine their fraction of leaves burned here + + leaf_burn_frac = currentPatch%burnt_frac_litter(lg_sf) + endif + + ! Perform a check to make sure that spitfire gave + ! us reasonable mortality and burn fraction rates + + if( (leaf_burn_frac < 0._r8) .or. & + (leaf_burn_frac > 1._r8) .or. & + (currentCohort%fire_mort < 0._r8) .or. & + (currentCohort%fire_mort > 1._r8)) then + write(fates_log(),*) 'unexpected fire fractions' + write(fates_log(),*) prt_params%woody(currentCohort%pft) + write(fates_log(),*) leaf_burn_frac + write(fates_log(),*) currentCohort%fire_mort + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + do el = 1,num_elements + + leaf_m = nc%prt%GetState(leaf_organ, element_list(el)) + + currentSite%mass_balance(el)%burn_flux_to_atm = & + currentSite%mass_balance(el)%burn_flux_to_atm + & + leaf_burn_frac * leaf_m * nc%n + end do + + ! Here the mass is removed from the plant + + call PRTBurnLosses(nc%prt, leaf_organ, leaf_burn_frac) + currentCohort%fraction_crown_burned = 0.0_r8 + nc%fraction_crown_burned = 0.0_r8 + + + + ! Logging is the current disturbance + elseif (i_disturbance_type .eq. dtype_ilog ) then + + ! If this cohort is in the upper canopy. It generated + if(currentCohort%canopy_layer == 1)then + + ! calculate the survivorship of disturbed trees because non-harvested + nc%n = currentCohort%n * currentCohort%l_degrad + ! nc%n = (currentCohort%l_degrad / (currentCohort%l_degrad + & + ! currentCohort%lmort_direct + currentCohort%lmort_collateral + + ! currentCohort%lmort_infra) ) * & + ! currentCohort%n * patch_site_areadis/currentPatch%area + + ! Reduce counts in the existing/donor patch according to the logging rate + currentCohort%n = currentCohort%n * & + (1.0_r8 - min(1.0_r8,(currentCohort%lmort_direct + & + currentCohort%lmort_collateral + & + currentCohort%lmort_infra + currentCohort%l_degrad))) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + + ! since these are the ones that weren't logged, + ! set the logging mortality rates as zero + nc%lmort_direct = 0._r8 + nc%lmort_collateral = 0._r8 + nc%lmort_infra = 0._r8 + + else + + ! What to do with cohorts in the understory of a logging generated + ! disturbance patch? + + if(prt_params%woody(currentCohort%pft) == itrue)then + + + ! Survivorship of undestory woody plants. Two step process. + ! Step 1: Reduce current number of plants to reflect the + ! change in area. + ! The number density per square are doesn't change, + ! but since the patch is smaller + ! and cohort counts are absolute, reduce this number. + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! because the mortality rate due to impact for the cohorts which had + ! been in the understory and are now in the newly- + ! disturbed patch is very high, passing the imort directly to + ! history results in large numerical errors, on account + ! of the sharply reduced number densities. so instead pass this info + ! via a site-level diagnostic variable before reducing + ! the number density. + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) = & + currentSite%imort_rate(currentCohort%size_class, currentCohort%pft) + & + nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac / hlm_freq_day + + currentSite%imort_carbonflux(currentCohort%pft) = & + currentSite%imort_carbonflux(currentCohort%pft) + & + (nc%n * currentPatch%fract_ldist_not_harvested * & + logging_coll_under_frac/ hlm_freq_day ) * & + total_c * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + + + ! Step 2: Apply survivor ship function based on the understory death fraction + + ! remaining of understory plants of those that are knocked + ! over by the overstorey trees dying... + ! LOGGING SURVIVORSHIP OF UNDERSTORY PLANTS IS SET AS A NEW PARAMETER + ! in the fatesparameter files + nc%n = nc%n * (1.0_r8 - & + (1.0_r8-currentPatch%fract_ldist_not_harvested) * logging_coll_under_frac) + + ! Step 3: Reduce the number count of cohorts in the + ! original/donor/non-disturbed patch to reflect the area change + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + else + + ! grass is not killed by mortality disturbance events. + ! Just move it into the new patch area. + ! Just split the grass into the existing and new patch structures + nc%n = currentCohort%n * patch_site_areadis/currentPatch%area + + ! Those remaining in the existing + currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area) + + ! No grass impact mortality imposed on the newly created patch + nc%cmort = currentCohort%cmort + nc%hmort = currentCohort%hmort + nc%bmort = currentCohort%bmort + nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort + nc%dmort = currentCohort%dmort + nc%lmort_direct = currentCohort%lmort_direct + nc%lmort_collateral = currentCohort%lmort_collateral + nc%lmort_infra = currentCohort%lmort_infra + + endif ! is/is-not woody + + endif ! Select canopy layer + + else + write(fates_log(),*) 'unknown disturbance mode?' + write(fates_log(),*) 'i_disturbance_type: ',i_disturbance_type + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if ! Select disturbance mode + + if (nc%n > 0.0_r8) then + storebigcohort => new_patch%tallest + storesmallcohort => new_patch%shortest + if(associated(new_patch%tallest))then + tnull = 0 + else + tnull = 1 + new_patch%tallest => nc + nc%taller => null() + endif + + if(associated(new_patch%shortest))then + snull = 0 + else + snull = 1 + new_patch%shortest => nc + nc%shorter => null() + endif + nc%patchptr => new_patch + call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & + tnull, snull, storebigcohort, storesmallcohort) + + new_patch%tallest => storebigcohort + new_patch%shortest => storesmallcohort + else + + ! Get rid of the new temporary cohort + call DeallocateCohort(nc) + deallocate(nc) + + endif + + currentCohort => currentCohort%taller + enddo ! currentCohort + call sort_cohorts(currentPatch) + + !update area of donor patch + oldarea = currentPatch%area + currentPatch%area = currentPatch%area - patch_site_areadis + + ! for all disturbance rates that haven't been resolved yet, increase their amount so that + ! they are the same amount of gridcell-scale disturbance relative to the original patch size + if (i_disturbance_type .ne. N_DIST_TYPES) then + do i_dist2 = i_disturbance_type+1,N_DIST_TYPES + currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & + * oldarea / currentPatch%area + end do + end if + + ! sort out the cohorts, since some of them may be so small as to need removing. + ! the first call to terminate cohorts removes sparse number densities, + ! the second call removes for all other reasons (sparse culling must happen + ! before fusion) + call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) + call fuse_cohorts(currentSite,currentPatch, bc_in) + call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) + call sort_cohorts(currentPatch) + + end if ! if ( new_patch%area > nearzero ) then + + end if cp_nocomp_matches_2_if + currentPatch => currentPatch%younger + + enddo ! currentPatch patch loop. + + !*************************/ + !** INSERT NEW PATCH(ES) INTO LINKED LIST + !*************************/ + + if ( site_areadis_primary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + ! insert new youngest primary patch after all the secondary patches, if there are any. + ! this requires first finding the current youngest primary to insert the new one ahead of + if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then + found_youngest_primary = .false. + do while(associated(currentPatch) .and. .not. found_youngest_primary) + currentPatch => currentPatch%older + if (associated(currentPatch)) then + if (currentPatch%anthro_disturbance_label .eq. primaryforest) then + found_youngest_primary = .true. + endif + endif + end do + if (associated(currentPatch)) then + ! the case where we've found a youngest primary patch + new_patch_primary%older => currentPatch + new_patch_primary%younger => currentPatch%younger + currentPatch%younger%older => new_patch_primary + currentPatch%younger => new_patch_primary else - snull = 1 - new_patch%shortest => nc - nc%shorter => null() + ! the case where we haven't, because the patches are all secondaary, + ! and are putting a primary patch at the oldest end of the + ! linked list (not sure how this could happen, but who knows...) + new_patch_primary%older => null() + new_patch_primary%younger => currentSite%oldest_patch + currentSite%oldest_patch%older => new_patch_primary + currentSite%oldest_patch => new_patch_primary endif - nc%patchptr => new_patch - call insert_cohort(nc, new_patch%tallest, new_patch%shortest, & - tnull, snull, storebigcohort, storesmallcohort) - - new_patch%tallest => storebigcohort - new_patch%shortest => storesmallcohort else - - ! Get rid of the new temporary cohort - call DeallocateCohort(nc) - deallocate(nc) - + ! the case where there are no secondary patches at the start of the linked list (prior logic) + new_patch_primary%older => currentPatch + new_patch_primary%younger => null() + currentPatch%younger => new_patch_primary + currentSite%youngest_patch => new_patch_primary endif - - currentCohort => currentCohort%taller - enddo ! currentCohort - call sort_cohorts(currentPatch) - - !update area of donor patch - oldarea = currentPatch%area - currentPatch%area = currentPatch%area - patch_site_areadis - - ! for all disturbance rates that haven't been resolved yet, increase their amount so that - ! they are the same amount of gridcell-scale disturbance relative to the original patch size - if (i_disturbance_type .ne. N_DIST_TYPES) then - do i_dist2 = i_disturbance_type+1,N_DIST_TYPES - currentPatch%disturbance_rates(i_dist2) = currentPatch%disturbance_rates(i_dist2) & - * oldarea / currentPatch%area - end do - end if + endif - ! sort out the cohorts, since some of them may be so small as to need removing. + ! insert first secondary at the start of the list + if ( site_areadis_secondary .gt. nearzero) then + currentPatch => currentSite%youngest_patch + new_patch_secondary%older => currentPatch + new_patch_secondary%younger=> null() + currentPatch%younger => new_patch_secondary + currentSite%youngest_patch => new_patch_secondary + endif + + + ! sort out the cohorts, since some of them may be so small as to need removing. ! the first call to terminate cohorts removes sparse number densities, ! the second call removes for all other reasons (sparse culling must happen ! before fusion) - call terminate_cohorts(currentSite, currentPatch, 1,16,bc_in) - call fuse_cohorts(currentSite,currentPatch, bc_in) - call terminate_cohorts(currentSite, currentPatch, 2,16,bc_in) - call sort_cohorts(currentPatch) - end if ! if ( new_patch%area > nearzero ) then - - end if cp_nocomp_matches_2_if - currentPatch => currentPatch%younger - - enddo ! currentPatch patch loop. + if ( site_areadis_primary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) + call fuse_cohorts(currentSite,new_patch_primary, bc_in) + call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) + call sort_cohorts(new_patch_primary) + endif - !*************************/ - !** INSERT NEW PATCH(ES) INTO LINKED LIST - !*************************/ - - if ( site_areadis_primary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - ! insert new youngest primary patch after all the secondary patches, if there are any. - ! this requires first finding the current youngest primary to insert the new one ahead of - if (currentPatch%anthro_disturbance_label .eq. secondaryforest ) then - found_youngest_primary = .false. - do while(associated(currentPatch) .and. .not. found_youngest_primary) - currentPatch => currentPatch%older - if (associated(currentPatch)) then - if (currentPatch%anthro_disturbance_label .eq. primaryforest) then - found_youngest_primary = .true. - endif - endif - end do - if (associated(currentPatch)) then - ! the case where we've found a youngest primary patch - new_patch_primary%older => currentPatch - new_patch_primary%younger => currentPatch%younger - currentPatch%younger%older => new_patch_primary - currentPatch%younger => new_patch_primary - else - ! the case where we haven't, because the patches are all secondaary, - ! and are putting a primary patch at the oldest end of the - ! linked list (not sure how this could happen, but who knows...) - new_patch_primary%older => null() - new_patch_primary%younger => currentSite%oldest_patch - currentSite%oldest_patch%older => new_patch_primary - currentSite%oldest_patch => new_patch_primary + if ( site_areadis_secondary .gt. nearzero) then + call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) + call fuse_cohorts(currentSite,new_patch_secondary, bc_in) + call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) + call sort_cohorts(new_patch_secondary) endif - else - ! the case where there are no secondary patches at the start of the linked list (prior logic) - new_patch_primary%older => currentPatch - new_patch_primary%younger => null() - currentPatch%younger => new_patch_primary - currentSite%youngest_patch => new_patch_primary - endif - endif - - ! insert first secondary at the start of the list - if ( site_areadis_secondary .gt. nearzero) then - currentPatch => currentSite%youngest_patch - new_patch_secondary%older => currentPatch - new_patch_secondary%younger=> null() - currentPatch%younger => new_patch_secondary - currentSite%youngest_patch => new_patch_secondary - endif - - - ! sort out the cohorts, since some of them may be so small as to need removing. - ! the first call to terminate cohorts removes sparse number densities, - ! the second call removes for all other reasons (sparse culling must happen - ! before fusion) - - if ( site_areadis_primary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_primary, 1,17, bc_in) - call fuse_cohorts(currentSite,new_patch_primary, bc_in) - call terminate_cohorts(currentSite, new_patch_primary, 2,17, bc_in) - call sort_cohorts(new_patch_primary) - endif - - if ( site_areadis_secondary .gt. nearzero) then - call terminate_cohorts(currentSite, new_patch_secondary, 1,18,bc_in) - call fuse_cohorts(currentSite,new_patch_secondary, bc_in) - call terminate_cohorts(currentSite, new_patch_secondary, 2,18,bc_in) - call sort_cohorts(new_patch_secondary) - endif - - endif !end new_patch area - - - call check_patch_area(currentSite) - call set_patchno(currentSite) - - end do disturbance_type_loop + + endif !end new_patch area + + + call check_patch_area(currentSite) + call set_patchno(currentSite) + + end do disturbance_type_loop end do nocomp_pft_loop From 6e00ae0178f81b563aa6c287b069bab6f8278951 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 12 Jul 2022 14:30:00 -0400 Subject: [PATCH 8/8] Remove endrun --- biogeochem/EDPatchDynamicsMod.F90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index d837ed737a..a25045a58c 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -2330,11 +2330,6 @@ subroutine fuse_patches( csite, bc_in ) tpp => currentSite%youngest_patch tpp_loop: do while(associated(tpp)) - if(.not.associated(currentPatch))then - write(fates_log(),*) 'FATES fuse_patches(): currentPatch is not associated?' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif - both_associated_if: if(associated(tpp).and.associated(currentPatch))then !--------------------------------------------------------------------! ! only fuse patches whose anthropogenic disturbance category matches !