From 80099339ce25878b0bc625b2ee1db5b80d5a5519 Mon Sep 17 00:00:00 2001 From: Joshua Rady Date: Sun, 8 Aug 2021 17:22:09 -0400 Subject: [PATCH] Clean up thinning routines from recent edits. --- biogeochem/FatesVegetationManagementMod.F90 | 484 ++------------------ 1 file changed, 38 insertions(+), 446 deletions(-) diff --git a/biogeochem/FatesVegetationManagementMod.F90 b/biogeochem/FatesVegetationManagementMod.F90 index 25b643696a..d98114628b 100644 --- a/biogeochem/FatesVegetationManagementMod.F90 +++ b/biogeochem/FatesVegetationManagementMod.F90 @@ -3274,7 +3274,6 @@ subroutine hardwood_control(site, efficiency) ! method, where, area/patch_fracti type(ed_patch_type), pointer :: current_patch type(ed_cohort_type), pointer :: current_cohort - ! Locals: real(r8) :: the_efficiency ! ---------------------------------------------------------------------------------------------- @@ -3347,7 +3346,7 @@ subroutine thin_patch_proportional(patch, pfts, patch_fraction, thin_fraction, f ! !----------------------------------------------------------------------------------------------- ! Operation Type: Thinnning - ! Operation Level: Site + ! Operation Level: Patch ! Regime: Multiple ! ---------------------------------------------------------------------------------------------- @@ -3405,8 +3404,6 @@ subroutine thin_patch_proportional(patch, pfts, patch_fraction, thin_fraction, f thin_pfts => tree_pfts endif ! present(pfts)... - if (vm_debug_verbose) write(fates_log(), *) 'Done with thin_pfts:', thin_pfts ! Very temporary. - ! Validate the patch fraction: if (present(patch_fraction)) then if (patch_fraction > 1.0_r8 .or. patch_fraction <= 0.0_r8) then @@ -3419,8 +3416,6 @@ subroutine thin_patch_proportional(patch, pfts, patch_fraction, thin_fraction, f the_patch_fraction = 1.0_r8 end if ! present(patch_fraction) - if (vm_debug_verbose) write(fates_log(), *) 'Done with the_patch_fraction:', the_patch_fraction ! Very temporary. - ! Only thin_fraction, final_basal_area, or final_stem_density should be provided: ! Since we thin proportionally we can convert stem densities and basal areas to a fraction. ! Since the thinning specifications are per area the patch fraction does not effect them. @@ -3461,25 +3456,17 @@ subroutine thin_patch_proportional(patch, pfts, patch_fraction, thin_fraction, f call endrun(msg = errMsg(__FILE__, __LINE__)) endif ! count(...) - if (vm_debug_verbose) write(fates_log(), *) 'Done with the_thin_fraction:', the_thin_fraction ! Very temporary. - if (the_thin_fraction > 0) then ! Remove a fraction of all trees from all matching cohorts, adjusting for the patch fraction: current_cohort => patch%shortest do while(associated(current_cohort)) - if (vm_debug_verbose) write(fates_log(), *) 'Cohort...' ! Very temporary. - if (any(thin_pfts == current_cohort%pft)) then - if (vm_debug_verbose) write(fates_log(), *) 'kill()...' - call kill(cohort = current_cohort, flux_profile = bole_harvest, & kill_fraction = (the_thin_fraction * the_patch_fraction), & area_fraction = the_patch_fraction) - if (vm_debug_verbose) write(fates_log(), *) 'harvest calc...' ! Very temporary. - ! Accumulate harvest estimate: harvest = harvest + (cohort_harvestable_biomass(current_cohort) * & the_thin_fraction * the_patch_fraction) @@ -3612,26 +3599,20 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi ! Locals: integer(i4), pointer, dimension(:) :: thin_pfts ! Holds the PFTs to thin, computed from arguments. real(r8) :: the_patch_fraction ! Patch fraction after application of defaults. - !real(r8) :: the_final_stem_density ! Stem density goal possibly converted from thin_fraction. + real(r8) :: harvest ! Accumulator logical :: use_bai ! If true use basal area index as the criteria, otherwise use stem density. - !real(r8) :: pfts_bai ! Basal area of the patch, including only PFTs in pfts. real(r8) :: pfts_sd ! Stem density of the patch, including only PFTs in pfts. - !real(r8) :: cohort_ba ! Basal area of a cohort - !real(r8) :: cohort_stems - !real(r8) :: thin_ba_remaining - !real(r8) :: thin_sd_remaining + ! Variables used for iterative thinning: real(r8) :: current_state real(r8) :: goal_state real(r8) :: cohort_state real(r8) :: thin_amount_remaining - real(r8) :: cohort_fraction - real(r8) :: harvest ! Accumulator + real(r8) :: cohort_fraction ! Fraction of the last cohort needing to be thinned. type(ed_cohort_type), pointer :: current_cohort - integer :: i ! Iterator ! ---------------------------------------------------------------------------------------------- @@ -3639,7 +3620,6 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi ! Initialize variables: harvest = 0.0_r8 - !thin_ba_remaining = 0.0_r8 ! Determine which arguments were passed in, validity check them, and supply default values: @@ -3674,173 +3654,7 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi else the_patch_fraction = 1.0_r8 end if - ! the_patch_fraction is not actually used below!!!!! - - ! The following code block works, at least if the patch area is assumed to be 1.---------------- - ! I'm preserving it for reference and to revert if needed while I rewrite it below. -! -! ! Determine the basal area and stem densities for the relavant PFTs: -! pfts_bai = disturbed_basal_area(patch, thin_pfts) -! pfts_sd = patch_disturbed_n(patch, thin_pfts) !disturbed_stem_density(patch, thin_pfts) -! -! ! Only thin_fraction, final_basal_area, or final_stem_density should be provided: -! ! Update logic using provided()!!!!! -! if (count([(present(final_basal_area) .and. final_basal_area /= vm_empty_real), & -! (present(final_stem_density) .and. final_stem_density /= vm_empty_real), & -! (present(thin_fraction) .and. thin_fraction /= vm_empty_real)]) > 1) then -! write(fates_log(),*) 'thin_patch_low_perfect(): Provide thinning fraction, basal area index, or stem density, not more than one.' -! call endrun(msg = errMsg(__FILE__, __LINE__)) -! else if (present(final_basal_area) .and. final_basal_area /= vm_empty_real) then -! use_bai = .true. -! else if (present(final_stem_density) .and. final_stem_density /= vm_empty_real) then -! use_bai = .false. -! the_final_stem_density = final_stem_density -! else if (present(thin_fraction) .and. thin_fraction /= vm_empty_real) then -! ! Validity checking: -! if (thin_fraction <= 0.0_r8 .or. thin_fraction > 1.0_r8) then -! write(fates_log(),*) 'Invalid value for thin_fraction argument.' -! call endrun(msg = errMsg(__FILE__, __LINE__)) -! endif -! -! ! Convert the thinning fraction to a stem density goal: -! use_bai = .false. -! the_final_stem_density = pfts_sd * (1.0_r8 - thin_fraction) -! else -! write(fates_log(),*) 'thin_patch_low_perfect(): Must provide thinning fraction, basal area index, or stem density.' -! call endrun(msg = errMsg(__FILE__, __LINE__)) -! endif -! -! ! Temporary reporting: -! ! if (debug) then -! ! write(fates_log(), *) 'thin_pfts: ', thin_pfts -! ! write(fates_log(), *) 'the_patch_fraction: ', the_patch_fraction -! ! if (present(thin_fraction)) then -! ! write(fates_log(), *) 'thin_fraction: ', thin_fraction -! ! end if -! ! write(fates_log(), *) 'pfts_bai: ', pfts_bai -! ! if (present(final_basal_area)) then -! ! write(fates_log(), *) 'final_basal_area: ', final_basal_area -! ! endif -! ! write(fates_log(), *) 'pfts_sd: ', pfts_sd -! ! write(fates_log(), *) 'the_final_stem_density: ', the_final_stem_density -! ! write(fates_log(), *) 'use_bai: ', use_bai -! ! endif -! -! ! If the stand is above the goal value (BAI or stem density) start thinning: -! if ((use_bai .and. (pfts_bai > final_basal_area)) .or. & -! ((.not. use_bai) .and. (pfts_sd > the_final_stem_density))) then -! -! ! There may be previously staged mortality. Therefore we need to keep track of the changes to -! ! BAI, density, and plant number's: -! pfts_bai = disturbed_basal_area(patch, thin_pfts) -! pfts_sd = patch_disturbed_n(patch, thin_pfts) !disturbed_stem_density(patch, thin_pfts) -! -! ! Cull the smallest trees (low thinning) recursively until the goal has been reached: -! -! if (use_bai) then ! Thin to a goal basal area: -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() starting low thinning by BAI.' -! -! ! Given the fixed allometry this will also give us cohorts from the lowest DBH: -! current_cohort => patch%shortest -! do while(associated(current_cohort) .and. pfts_bai > final_basal_area) -! -! if (any(thin_pfts == current_cohort%pft)) then -! ! Get the effective basal area (after any staged removals) of the cohort and determine -! ! if it can be removed in part or in whole: -! cohort_ba = disturbed_basal_area(current_cohort) -! -! ! Remaining basal area: -! thin_ba_remaining = pfts_bai - final_basal_area -! -! ! If the cohort basal area is less that what still needs to be removed kill all of it: -! if (cohort_ba <= thin_ba_remaining) then -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() cut whole cohort.' -! -! call kill_disturbed(cohort = current_cohort, flux_profile = bole_harvest, & -! kill_fraction = 1.0_r8) -! -! else ! Otherwise only take part of the cohort: -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() cut part of cohort.' -! -! cohort_fraction = thin_ba_remaining / cohort_ba -! -! if (debug) then -! write(fates_log(), *) 'thin_ba_remaining = ', thin_ba_remaining -! write(fates_log(), *) 'cohort_ba = ', cohort_ba -! write(fates_log(), *) 'cohort_fraction = ', cohort_fraction -! end if -! -! call kill_disturbed(cohort = current_cohort, flux_profile = bole_harvest, & -! kill_fraction = cohort_fraction) -! -! end if ! (cohort_ba <= thin_ba_remaining) -! -! ! The harvest estimate and BAI only needs to be updated if we harvested something: -! ! Accumulate harvest estimate: -! ! This is wrong isn't it! This will count all cohort regardless of harvest extent!!!!! -! harvest = harvest + cohort_harvestable_biomass(current_cohort) ! staged = true!!!! -! pfts_bai = disturbed_basal_area(patch, thin_pfts) ! Update the BAI calculation. -! end if ! (any(thin_pfts == current_cohort%pft)) -! -! current_cohort => current_cohort%taller -! end do ! Cohort loop. -! -! else ! Thin to a goal stem density: -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() starting low thinning by stem density.' -! -! ! This loop is identically structured to the above so the two could be combined by -! ! generalizing the comparator variables. -! -! current_cohort => patch%shortest -! do while(pfts_sd > the_final_stem_density) -! -! !if (debug) call dump_cohort(current_cohort) ! Temporary. -! -! if (any(thin_pfts == current_cohort%pft)) then -! ! Get the effective number of stems in cohort and determine if they can be removed in -! ! part or in whole. Because n is per nominal hectare n = stem density (n/ha): -! cohort_stems = cohort_disturbed_n(current_cohort) -! !if (debug) write(fates_log(), *) "cohort_stems: ", cohort_stems ! Temporary. -! -! ! Remaining stems to remove: -! thin_sd_remaining = pfts_sd - the_final_stem_density -! !if (debug) write(fates_log(), *) "thin_sd_remaining: ", thin_sd_remaining ! Temporary. -! -! ! If the cohort stem count is less that what still needs to be removed kill all of it: -! if (cohort_stems <= thin_sd_remaining) then -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() cut whole cohort.' -! -! call kill_disturbed(cohort = current_cohort, flux_profile = bole_harvest, & -! kill_fraction = 1.0_r8) -! -! else ! Otherwise only take part of the cohort: -! if (debug) write(fates_log(), *) 'thin_patch_low_perfect() cut part of cohort.' -! -! cohort_fraction = thin_sd_remaining / cohort_stems -! call kill_disturbed(cohort = current_cohort, flux_profile = bole_harvest, & -! kill_fraction = cohort_fraction) -! -! end if -! -! ! The harvest estimate and stem density only need to updated if we harvested something: -! ! Accumulate harvest estimate: -! harvest = harvest + cohort_harvestable_biomass(current_cohort) ! staged = true!!!! -! pfts_sd = patch_disturbed_n(patch, thin_pfts) ! disturbed_stem_density(patch, thin_pfts) ! Update the stem density. -! !if (debug) write(fates_log(), *) "pfts_sd: ", pfts_sd ! Temporary. -! end if ! (any(thin_pfts == current_cohort%pft)) -! -! current_cohort => current_cohort%taller -! end do ! Cohort loop. -! endif ! (use_bai) -! endif ! Stand is above goal loop. -! ! Consider adding warning here if patch is below goal? -! -! if (present(harvest_estimate)) then -! harvest_estimate = harvest -! endif -! ! Should anything be done or reported if no thinning was needed? - ! Rewrite using provided(), addressing the patch area, and making the loop logic generic:------- ! Only thin_fraction, final_basal_area, or final_stem_density should be provided: if (count([provided(final_basal_area), provided(final_stem_density), & @@ -3852,7 +3666,6 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi goal_state = final_basal_area else if (provided(final_stem_density)) then use_bai = .false. - !the_final_stem_density = final_stem_density goal_state = final_stem_density else if (provided(thin_fraction)) then ! Validity checking: @@ -3864,7 +3677,6 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi ! Convert the thinning fraction to a stem density goal: use_bai = .false. pfts_sd = disturbed_stem_density(patch, thin_pfts) - !the_final_stem_density = pfts_sd * (1.0_r8 - thin_fraction) goal_state = pfts_sd * (1.0_r8 - thin_fraction) else write(fates_log(),*) 'thin_patch_low_perfect(): Must provide thinning fraction, basal area index, or stem density.' @@ -3946,8 +3758,6 @@ subroutine thin_patch_low_perfect(patch, pfts, patch_fraction, thin_fraction, fi if (debug) write(fates_log(), *) 'thin_patch_low_perfect() harvest estimate: ', harvest_estimate endif - ! End Rewrite. - if (debug) write(fates_log(), *) 'thin_patch_low_perfect() exiting.' end subroutine thin_patch_low_perfect @@ -3963,8 +3773,6 @@ subroutine thin_low_perfect(site, pfts, thin_fraction, final_basal_area, final_s ! In the future we will (may) allow targeting of thinning behavior to specific patches via the ! 'where' parameter. ! - ! This routine does not yet return a harvest estimate. - ! !----------------------------------------------------------------------------------------------- ! Operation Type: Thinnning ! Operation Level: Site @@ -4084,17 +3892,14 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti ! In reality some very small trees may be left on site. ! ! Notes: - ! - This routine does use effective values in its calculations! - ! - Instead of looping through all the cohorts multiple times we could make a list of the valid cohorts once and loop through those... - ! - ! Use fraction or specific number of trees????? + ! - This routine does use effective values in all its calculations. + ! - Instead of looping through all the cohorts multiple times we could make a list of the valid + ! cohorts once and loop through those. ! !----------------------------------------------------------------------------------------------- ! Operation Type: Thinnning ! Operation Level: Patch ! Regime: Multiple - ! - ! Being debugged!!!!! ! ---------------------------------------------------------------------------------------------- ! Uses: NA @@ -4122,25 +3927,17 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti integer(i4), pointer, dimension(:) :: thin_pfts ! Holds the PFTs to thin, computed from arguments. integer :: i ! PFT iterator. real(r8) :: the_patch_fraction ! Patch fraction after application of defaults. - real(r8) :: the_thin_fraction ! The calculated actual or estimated fraction of trees to thin. - real(r8) :: sum_dbh ! Accumulator - real(r8) :: num_trees ! Accumulator + real(r8) :: sum_dbh ! Accumulator. + real(r8) :: num_trees ! Accumulator. real(r8) :: mean_dbh ! Mean (weighted) DBH of the patch. - !real(r8) :: sum_sq_dev - !real(r8) :: std_dev ! Standard deviation of DBH for the patch. - real(r8) :: model_midpoint ! Thinning probability function midpoint term in centered DBH space. real(r8) :: midpoint_step ! Step for linear search phase of solving the thinning weights. - !real(r8) :: thin_goal ! The number of trees that should be thinned. - !real(r8) :: thin_total ! Accumulator for thinning calculation - !real(r8) :: thin_last ! The number of trees to thin calculated by the last iteration. real(r8) :: thin_prob ! The probability of thinning for a cohort. real(r8) :: dbh_tranformed ! Cohort DBH adjusted so the mean DBH is at 0. - logical :: use_bai ! If true use basal area index as the criteria, otherwise use stem density. ! Could reuse current_state rather than pfts_sd and pfts_ba... real(r8) :: pfts_sd ! Stem density of the patch, including only PFTs in pfts. @@ -4149,30 +3946,35 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti real(r8) :: last_sd ! The stem density for the last midpoint value assessed. real(r8) :: goal_state ! The SD or BA we are trying to thin to. real(r8) :: current_state ! The 'state' (SD or BA) resulting from the current midpoint. - !real(r8) :: last_state ! last_state The state with the previous midpoint. real(r8) :: best_state ! The best state (closest to goal) achieved so far. - real(r8) :: best_midpoint ! The best midpoint tested so far. - !real(r8) :: next_step ! Midpoint adjustment for the next round. [May remove] + real(r8) :: best_midpoint ! The best midpoint value tested so far. real(r8) :: midpoint_min ! Midpoint search range minimum. real(r8) :: midpoint_max ! Midpoint search range maximum. integer :: cycles ! Loop counter for debugging purposes. - real(r8) :: harvest ! Accumulator + real(r8) :: harvest ! Accumulator. type(ed_cohort_type), pointer :: current_cohort ! ---------------------------------------------------------------------------------------------- if (debug) write(fates_log(), *) 'thin_patch_low_probabilistic() entering.' + ! ---------------------------------------------------------------------------------------------- ! Initialize variables: + ! ---------------------------------------------------------------------------------------------- sum_dbh = 0.0_r8 num_trees = 0.0_r8 - !thin_total = 0.0_r8 - !thin_last = 0.0_r8 - !midpoint_step = 0.05_r8 ! The initial step size is arbitrary. We assume we are pretty close when we start. - !midpoint_step = 0.5_r8 ! Start by adjusting the midpoint by 0.5 cm if we are not within the tolerance. midpoint_step = 1.0_r8 ! Start by adjusting the midpoint by 1 cm if we are not within the tolerance. + current_sd = 0.0_r8 + last_sd = 100.0_r8 ! Set to a large initial value to get through the first loop check... Or add cycle check to loop? + current_state = 0.0_r8 + best_state = 0.0_r8 + ! The maximum search range for the midpoint is ill-defined because at high or low thinning + ! fractions the midpoint may extend beyond the range of actual sizes. We set them to initial + ! unreasonably extreem values and let the linear search phase find bounding values: + midpoint_min = -1000_r8 + midpoint_max = 1000_r8 cycles = 0 harvest = 0 @@ -4181,6 +3983,10 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti harvest_estimate = 0.0_r8 endif + ! ---------------------------------------------------------------------------------------------- + ! Validate arguments and calculate thinning goals: + ! ---------------------------------------------------------------------------------------------- + ! Validate the PFTs: [This block is repeated several times elsewhere.] ! Only allow trees and default to all trees. if (present(pfts) .and. any(pfts /= vm_empty_integer)) then @@ -4267,6 +4073,10 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti call endrun(msg = errMsg(__FILE__, __LINE__)) endif + ! ---------------------------------------------------------------------------------------------- + ! Solve for a set of thinning probabilities / fractions that yield the requested thinning goal. + ! ---------------------------------------------------------------------------------------------- + ! Calculate the starting midpoint parameter: model_midpoint = midpoint_slope * the_thin_fraction + midpoint_intercept !if (vm_debug_verbose) write(fates_log(), *) 'Initial model_midpoint:', model_midpoint ! Temporary reporting. @@ -4275,6 +4085,8 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti write(fates_log(), *) 'Initial model_midpoint:', model_midpoint endif + best_midpoint = model_midpoint + ! Calculate the mean DBH of the PFTs (cohort DBH weighted by number) to be thinned: current_cohort => patch%shortest do while(associated(current_cohort)) @@ -4287,219 +4099,11 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti end do ! Cohort loop. mean_dbh = sum_dbh / num_trees - ! Calculate the number of trees to be removed based on the thinning fraction: - !thin_goal = num_trees * thin_fraction - - ! Calculate the starting midpoint parameter based on the fraction of trees to be thinned: - !model_midpoint = midpoint_slope * thin_fraction + midpoint_intercept - - - ! V1: This works for a thinning fraction but needs to be modified for other specifications. - - ! Starting with the initial midpoint value calculate thinning weights and repeat the process - ! with adjusted values until we are within the tolerance: -! do while(abs(thin_total - thin_goal) > thin_tolerance) -! cycles = cycles + 1 ! Counter ... -! thin_total = 0.0_r8 -! -! ! For each valid cohort determine the probability of thinning and number of trees to thin: -! current_cohort => patch%shortest -! do while(associated(current_cohort)) -! -! if (any(thin_pfts == current_cohort%pft)) then -! dbh_tranformed = current_cohort%dbh - mean_dbh -! !thin_prob = 1 / (1 + exp(-1.0r8 * model_steepness * (dbh_tranformed - model_midpoint))) ! -model_steepness? -! thin_prob = 1 / (1 + exp(-model_steepness * (dbh_tranformed - model_midpoint))) -! thin_total = thin_total + (thin_prob * current_cohort%n) -! endif -! -! current_cohort => current_cohort%taller -! end do ! Cohort loop. -! -! ! Adjust the midpoint for the next cycle (may not be used if we are close enough already): -! if (cycles > 1) then -! ! If we passed over the goal since the last step decrease the step size: -! if ((thin_last < thin_goal .and. thin_total > thin_goal) .or. & -! (thin_last > thin_goal .and. thin_total < thin_goal)) then -! midpoint_step = midpoint_step / 2.0_r8 -! endif -! endif -! -! ! Step the midpoint up or down: -! if (thin_total < thin_goal) then -! model_midpoint = model_midpoint + midpoint_step -! elseif (thin_total > thin_goal) then -! model_midpoint = model_midpoint - midpoint_step -! endif -! -! ! Are there invalid midpoint values that will let us know the search has gone wrong? -! -! thin_last = thin_total ! Record the trees thinned for comparison during the next cycle. -! ! cycles = cycles + 1 ! Counter for debugging purposes only. -! -! if (vm_debug_verbose) then ! Temporary reporting: -! write(fates_log(), *) 'cycles:', cycles -! write(fates_log(), *) 'model_midpoint:', model_midpoint -! write(fates_log(), *) 'midpoint_step:', midpoint_step -! write(fates_log(), *) 'thin_total:', thin_total -! endif -! -! end do ! Thinning calculation -! -! if (vm_debug_verbose) write(fates_log(), *) 'Thinning calculation completed.' ! Temporary reporting. -! -! ! Once the proper thinning weights have been solved apply the thinning mortality: -! ! This is a bit repetetive but there is little avoiding that this must be done in a loop. -! current_cohort => patch%shortest -! do while(associated(current_cohort)) -! if (vm_debug_verbose) write(fates_log(), *) 'Starting cohort.' ! Temporary reporting. -! -! if (any(thin_pfts == current_cohort%pft)) then -! dbh_tranformed = current_cohort%dbh - mean_dbh -! thin_prob = 1 / (1 + exp(-model_steepness * (dbh_tranformed - model_midpoint))) -! -! if (vm_debug_verbose) write(fates_log(), *) 'thin_prob:', thin_prob! Temporary reporting. -! -! ! Some or all trees could be left in place but we assume a commercial harvest: -! call kill(cohort = current_cohort, flux_profile = bole_harvest, kill_fraction = thin_prob) -! endif -! -! current_cohort => current_cohort%taller -! end do ! Cohort loop. - - - ! V2: - ! Add current_state, goal_state, state_last / last_state -! -! current_sd = 0 -! current_state = 0.0_r8 -! last_state = 0.0_r8 ! last_state -! best_state = = 0.0_r8 -! best_midpoint = model_midpoint -! next_step = 0.0_r8 -! -! ! Starting with the initial midpoint value calculate thinning weights and repeat the process -! ! with adjusted values until we are within the tolerance: -! ! Here the tolerance is applied to the amount by which the stem density is changing, with the -! ! idea that the value decreases as we get closer. -! do while (abs(current_sd - last_sd) > thin_tolerance) -! cycles = cycles + 1 ! Counter for debugging purposes only. (Actually used below), -! current_state = 0.0_r8 -! -! model_midpoint = model_midpoint + next_step ! Adjust the midpoint if appropriate. -! -! ! Calculate the result of thinning with the current midpoint value: -! current_cohort => patch%shortest -! do while(associated(current_cohort)) -! -! if (any(thin_pfts == current_cohort%pft)) then -! dbh_tranformed = current_cohort%dbh - mean_dbh -! thin_prob = 1 / (1 + exp(-model_steepness * (dbh_tranformed - model_midpoint))) -! -! ! Calculate the stem density in all cases for the master loop conditional: -! current_sd = current_sd + ((1 - thin_prob) * current_cohort%n) -! -! ! Calculate the impact on the state: -! if (use_bai) then -! current_state = current_state + (cohort_basal_area(NOT IMPLEMENTED) * (1 - thin_prob)) !!!!!!!!!!!!!!!!!! -! else -! !current_state = current_state + ((1 - thin_prob) * current_cohort%n) -! current_state = current_sd -! endif -! endif -! -! current_cohort => current_cohort%taller -! end do ! Cohort loop. -! -! ! If this is the closest we have approached the goal store it: -! if (abs(goal_state - current_state) < abs(goal_state - best_state)) then -! best_state = current_state ! or state_closest = abs(goal_state - best_state) -! best_midpoint = model_midpoint -! endif -! -! ! Adjust the midpoint for the next cycle (may not be used if we are close enough already): -! if (cycles > 1) then -! ! If we passed over the goal since the last step decrease the step size: -! if ((last_state < goal_state .and. current_state > goal_state) .or. & -! (last_state > goal_state .and. current_state < goal_state)) then -! midpoint_step = midpoint_step / 2.0_r8 -! endif -! endif -! -! ! Step the midpoint up or down (but wait until we start the next loop): -! ! Delaying the adjustment is probably not necessary with the addition of best_midpoint? -! if (current_state < goal_state) then -! !model_midpoint = model_midpoint + midpoint_step -! next_step = midpoint_step -! elseif (current_state > goal_state) then -! !model_midpoint = model_midpoint - midpoint_step -! next_step = -midpoint_step -! endif -! -! ! Record information about this thinning attempt for comparison during the next cycle: -! last_sd = current_sd -! last_state = current_state -! -! if (vm_debug_verbose) then ! Temporary reporting: -! write(fates_log(), *) 'cycles:', cycles -! write(fates_log(), *) 'model_midpoint:', model_midpoint -! write(fates_log(), *) 'midpoint_step:', midpoint_step -! write(fates_log(), *) 'current_sd:', current_sd -! write(fates_log(), *) 'current_state:', current_state -! endif -! -! end do ! Thinning calculation -! -! ! Note: I think it is possible we could get to this point without reversing direction in the -! ! search. If so the final value of the midpoint could actually yield result further from the -! ! goal than the last one value. I think this is unlikely so we will just monitor for it. -! ! if (debug .and. (abs(goal_state - last_state) < abs(goal_state - state_current))) then -! ! write(fates_log(), *) 'Did not end on best midpoint value.' -! ! end if -! ! I think I handled this above with best_state... -! -! model_midpoint = best_midpoint ! The last value might not have been the closest one. -! -! if (vm_debug_verbose) write(fates_log(), *) 'Thinning calculation completed.' ! Temporary reporting. - - ! V3: - ! - !sum_sq_dev = 0.0_r8 - current_sd = 0.0_r8 - last_sd = 100.0_r8 ! Set to a large initial to get through the first loop check... Or add cycle check to loop? - current_state = 0.0_r8 - !last_state = 100.0_r8 ! Set to a large initial to get through the first loop check... Or add cycle check to loop? - best_state = 0.0_r8 - best_midpoint = model_midpoint - - ! The maximum search range for the midpoint is ill-defined because at high or low thinning - ! fractions the midpoint may extend beyond the range of actual sizes. - - ! Approach 1: Using the standard deviation of the DBH seems sensible. - ! Calculate the standard deviation of the DBHs: -! current_cohort => patch%shortest -! do while(associated(current_cohort)) -! if (any(thin_pfts == current_cohort%pft)) then -! sum_sq_dev = sum_sq_dev + (current_cohort%dbh - mean_dbh)**2 * current_cohort%n -! endif -! current_cohort => current_cohort%taller -! end do ! Cohort loop. -! std_dev = sum_sq_dev / num_trees -! -! ! Make the range 3 SD wide (transformed DBH is centereded at 0): -! midpoint_min = -1.5_r8 * std_dev -! midpoint_max = 1.5_r8 * std_dev - - ! Approach 2: - ! Alternatively, we can just set them to initial unreasonable values and let the linear search - ! phase find bounding values: (this argues for a larger steps size maybe?) - midpoint_min = -1000_r8 - midpoint_max = 1000_r8 - ! Starting with the initial midpoint value calculate thinning weights and repeat the process ! with adjusted values until we are within the tolerance: ! Here the tolerance is applied to the amount by which the stem density is changing, with the - ! idea that the value decreases as we get closer. + ! idea that the value decreases as we get closer. SD is used for all specification modes to + ! give consistency. do while (abs(current_sd - last_sd) > thin_tolerance) cycles = cycles + 1 ! Counter for debugging purposes only. current_state = 0.0_r8 @@ -4520,7 +4124,6 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti ! Calculate the impact on the state: if (use_bai) then - ! The basal area doesn't need to be effective here (there is no previous mortality) but it will work fine: current_state = current_state + (effective_basal_area(current_cohort) * (1 - thin_prob)) else current_state = current_sd @@ -4572,9 +4175,7 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti model_midpoint = best_midpoint ! The last value might not have been the closest one. - ! Once the proper thinning weights have been solved apply the thinning mortality: - ! This is a bit repetetive but there is little avoiding that this must be done in a loop. current_cohort => patch%shortest do while(associated(current_cohort)) if (vm_debug_verbose) write(fates_log(), *) 'Starting cohort.' ! Temporary reporting. @@ -4598,10 +4199,6 @@ subroutine thin_patch_low_probabilistic(patch, pfts, patch_fraction, thin_fracti if (debug) then ! Report algorithm stats: - !write(fates_log(), *) 'Thinning fraction specified: ', thin_fraction - !write(fates_log(), *) 'As trees: ', thin_goal - !write(fates_log(), *) 'Number of trees thinned: ', thin_total - if (use_bai) then write(fates_log(), *) 'BA specified: ', goal_state else @@ -4642,7 +4239,7 @@ subroutine thin_low_probabilistic(site, pfts, thin_fraction, final_basal_area, & ! Regime: Multiple ! ! VM Event Interface: - ! thin_low_probabilistic([pfts], thin_fraction OR final_basal_area OR final_stem_density) Being debugged!!!!! + ! thin_low_probabilistic([pfts], thin_fraction OR final_basal_area OR final_stem_density) ! ---------------------------------------------------------------------------------------------- ! Uses: NA @@ -4650,29 +4247,25 @@ subroutine thin_low_probabilistic(site, pfts, thin_fraction, final_basal_area, & ! Arguments: type(ed_site_type), intent(in), target :: site ! The current site object. integer(i4), dimension(:), intent(in), optional :: pfts ! An array of PFT IDs to thin. - !real(r8), intent(in) :: thin_fraction ! Fraction of trees to thin. ! Three ways to specify how much to thin: real(r8), intent(in), optional :: thin_fraction ! Fraction of trees (in pfts) to thin. - real(r8), intent(in), optional :: final_basal_area ! Goal final basal area index (m^2 / ha) - real(r8), intent(in), optional :: final_stem_density ! Goal final stem density (trees / ha) + real(r8), intent(in), optional :: final_basal_area ! Goal final basal area index (m^2 / ha). + real(r8), intent(in), optional :: final_stem_density ! Goal final stem density (trees / ha). real(r8), intent(out), optional :: harvest_estimate ! The wood harvested by this operation. ! Locals: type(ed_patch_type), pointer :: current_patch real(r8) :: patch_harvest - real(r8) :: total_harvest ! Accumulator + real(r8) :: total_harvest ! Accumulator. ! ---------------------------------------------------------------------------------------------- if (debug) write(fates_log(), *) 'thin_low_probabilistic() beginning.' - ! Add handling for 'empty' pfts value. - current_patch => site%oldest_patch do while (associated(current_patch)) - !call thin_patch_low_probabilistic(current_patch, pfts, thin_fraction) call thin_patch_low_probabilistic(patch = current_patch, pfts = pfts, & thin_fraction = thin_fraction, & final_basal_area = final_basal_area, & @@ -4730,7 +4323,6 @@ subroutine thin_patch_row_low(patch, pfts, row_fraction, patch_fraction, & ! Uses: use FatesInterfaceTypesMod, only : numpft - !use EDPftvarcon, only : EDPftvarcon_inst ! Will change to prt_params soon when merged to master. ! Arguments: type(ed_patch_type), intent(inout), target :: patch ! Or pointer?????