From e48065bcd435f7c91b446b10cbdcad750f83df34 Mon Sep 17 00:00:00 2001 From: Charles Koven Date: Fri, 16 Nov 2018 16:28:33 -0800 Subject: [PATCH 01/28] changed cohort fusion to conserve crown area --- biogeochem/EDCohortDynamicsMod.F90 | 37 +++++++++++-------- biogeochem/FatesAllometryMod.F90 | 59 ++++++++++++++++++++++-------- main/FatesConstantsMod.F90 | 2 +- 3 files changed, 66 insertions(+), 32 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 71a1416cfc..6bef2c35ad 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -11,6 +11,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : itrue,ifalse + use FatesConstantsMod , only : fates_unset_r8 use FatesInterfaceMod , only : hlm_days_per_year use EDPftvarcon , only : EDPftvarcon_inst use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type @@ -729,6 +730,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) real(r8) :: newn real(r8) :: diff real(r8) :: dynamic_fusion_tolerance + real(r8) :: dbh integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion real(r8) :: larger_n, smaller_n @@ -794,6 +796,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) if( currentCohort%isnew.eqv.nextc%isnew ) then newn = currentCohort%n + nextc%n + fusion_took_place = 1 if ( fuse_debug .and. currentCohort%isnew ) then @@ -822,26 +825,28 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & + nextc%n*nextc%laimemory)/newn - currentCohort%dbh = (currentCohort%n*currentCohort%dbh & - + nextc%n*nextc%dbh)/newn - - call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite) - currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + nextc%n*nextc%canopy_trim)/newn - ! ----------------------------------------------------------------- - ! If fusion pushed structural biomass to be larger than - ! the allometric target value derived by diameter, we - ! then increase diameter and height until the allometric - ! target matches actual bdead. (if it is the other way around - ! we then just let the carbon pools grow to fill-out allometry) - ! ----------------------------------------------------------------- + ! conserve total crown area during the fusion step, and then calculate dbh of the + ! fused cohort as that which conserves both crown area and the dbh to crown area allometry. + ! dbh will be updated in the next growth step in the (likely) event that dbh to structural + ! biomass allometry is exceeded. if using a capped crown area allometry and above the cap, + ! then calculate as the weighted average of fusing cohorts' dbh + + currentCohort%c_area = currentCohort%c_area + nextc%c_area + + call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& + currentCohort%c_area,inverse=.true.) - if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then - call StructureResetOfDH( currentCohort%prt%GetState(struct_organ,all_carbon_elements), currentCohort%pft, & - currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite ) - end if + if (dbh .eq. fates_unset_r8) then + currentCohort%dbh = (currentCohort%n*currentCohort%dbh & + + nextc%n*nextc%dbh)/newn + else + currentCohort%dbh = dbh + endif + + call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 81881a190a..74f08b6072 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -90,6 +90,7 @@ module FatesAllometryMod use FatesConstantsMod, only : cm2_per_m2 use FatesConstantsMod, only : kg_per_Megag use FatesConstantsMod, only : calloc_abs_error + use FatesConstantsMod, only : fates_unset_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -448,15 +449,17 @@ end subroutine blmax_allom ! Generic crown area allometry wrapper ! ============================================================================ - subroutine carea_allom(d,nplant,site_spread,ipft,c_area) + subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse) - real(r8),intent(in) :: d ! plant diameter [cm] + real(r8),intent(inout) :: d ! plant diameter [cm] real(r8),intent(in) :: site_spread ! site level spread factor (crowdedness) real(r8),intent(in) :: nplant ! number of plants [1/ha] integer(i4),intent(in) :: ipft ! PFT index - real(r8),intent(out) :: c_area ! crown area per cohort (m2) + real(r8),intent(inout) :: c_area ! crown area per cohort (m2) + logical,optional,intent(in) :: inverse! if true, calculate dbh from crown area instead of crown area from dbh real(r8) :: d_eff ! Effective diameter (cm) + logical :: do_inverse, capped_allom associate( dbh_maxh => EDPftvarcon_inst%allom_dbh_maxheight(ipft), & allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), & @@ -465,24 +468,45 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area) d2ca_min => EDPftvarcon_inst%allom_d2ca_coefficient_min(ipft), & d2ca_max => EDPftvarcon_inst%allom_d2ca_coefficient_max(ipft)) + if( .not. present(inverse) ) then + do_inverse = .false. + else + do_inverse = inverse + endif + + if (do_inverse) then + c_area = c_area / nplant + endif + select case(int(allom_lmode)) case(1) d_eff = min(d,dbh_maxh) - call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + capped_allom = .true. case(2) ! "2par_pwr") - call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + capped_allom = .false. case(3) d_eff = min(d,dbh_maxh) - call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + capped_allom = .true. case DEFAULT - write(fates_log(),*) 'An undefined leaf allometry was specified: ', & + write(fates_log(),*) 'An undefined leaf allometry was specified: ', & allom_lmode - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + - c_area = c_area * nplant + if (capped_allom .and. do_inverse) then + if (d_eff .lt. dbh_maxh) then + d = d_eff + else + d = fates_unset_r8 + endif + endif + c_area = c_area * nplant end associate return @@ -1848,20 +1872,21 @@ end subroutine h2d_martcano ! ============================================================================= - subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) + subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse) ! ============================================================================ ! Calculate area of ground covered by entire cohort. (m2) ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. ! ============================================================================ - real(r8),intent(in) :: d ! diameter [cm] + real(r8),intent(inout) :: d ! diameter [cm] real(r8),intent(in) :: spread ! site level relative spread score [0-1] real(r8),intent(in) :: d2bl_p2 ! parameter 2 in the diameter->bleaf allometry (exponent) real(r8),intent(in) :: d2bl_ediff ! area difference factor in the diameter-bleaf allometry (exponent) real(r8),intent(in) :: d2ca_min ! minimum diameter to crown area scaling factor real(r8),intent(in) :: d2ca_max ! maximum diameter to crown area scaling factor - real(r8),intent(out) :: c_area ! crown area for one plant [m2] + real(r8),intent(inout) :: c_area ! crown area for one plant [m2] + logical,intent(in) :: inverse ! if true, calculate dbh from crown area rather than its reverse real(r8) :: crown_area_to_dbh_exponent real(r8) :: spreadterm ! Effective 2bh to crown area scaling factor @@ -1885,7 +1910,11 @@ subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area) spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min - c_area = spreadterm * d ** crown_area_to_dbh_exponent + if ( .not. inverse) then + c_area = spreadterm * d ** crown_area_to_dbh_exponent + else + d = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent) + endif end subroutine carea_2pwr diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index f585c53915..eadc1df658 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -17,7 +17,7 @@ module FatesConstantsMod ! Unset and various other 'special' values integer, parameter :: fates_unset_int = -9999 - + real(fates_r8), parameter :: fates_unset_r8 = -9999._fates_r8 ! Integer equivalent of true (in case some compilers dont auto convert) integer, parameter :: itrue = 1 From 7cdca1363307e3657c893d67623ba8b37535cb6e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 13 Dec 2018 12:30:43 -0800 Subject: [PATCH 02/28] Some minor syntax updates to inverse crown area allometry. Minor optimization to inverse crown area allometry. Implemented safer logic against a real. --- biogeochem/EDCohortDynamicsMod.F90 | 18 +++++---- biogeochem/FatesAllometryMod.F90 | 63 +++++++++++++++++------------- 2 files changed, 47 insertions(+), 34 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 6bef2c35ad..30c6abd8e7 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -12,6 +12,7 @@ module EDCohortDynamicsMod use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : fates_unset_r8 + use FatesConstantsMod , only : nearzero use FatesInterfaceMod , only : hlm_days_per_year use EDPftvarcon , only : EDPftvarcon_inst use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type @@ -820,7 +821,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! Fuse all mass pools - call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, currentCohort%n/newn ) + call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, & + currentCohort%n/newn ) currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & + nextc%n*nextc%laimemory)/newn @@ -828,18 +830,20 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + nextc%n*nextc%canopy_trim)/newn - ! conserve total crown area during the fusion step, and then calculate dbh of the - ! fused cohort as that which conserves both crown area and the dbh to crown area allometry. - ! dbh will be updated in the next growth step in the (likely) event that dbh to structural - ! biomass allometry is exceeded. if using a capped crown area allometry and above the cap, - ! then calculate as the weighted average of fusing cohorts' dbh + ! conserve total crown area during the fusion step, and then calculate + ! dbh of the fused cohort as that which conserves both crown area and + ! the dbh to crown area allometry. dbh will be updated in the next + ! growth step in the (likely) event that dbh to structural iomass + ! allometry is exceeded. if using a capped crown area allometry and + ! above the cap, then calculate as the weighted average of fusing + ! cohorts' dbh currentCohort%c_area = currentCohort%c_area + nextc%c_area call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& currentCohort%c_area,inverse=.true.) - if (dbh .eq. fates_unset_r8) then + if (abs(dbh-fates_unset_r8) EDPftvarcon_inst%allom_dbh_maxheight(ipft), & allom_lmode => EDPftvarcon_inst%allom_lmode(ipft), & @@ -472,23 +477,22 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse) do_inverse = .false. else do_inverse = inverse - endif - - if (do_inverse) then - c_area = c_area / nplant + if (do_inverse) then + c_area = c_area / nplant + endif endif select case(int(allom_lmode)) case(1) - d_eff = min(d,dbh_maxh) - call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + dbh_eff = min(dbh,dbh_maxh) + call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) capped_allom = .true. case(2) ! "2par_pwr") - call carea_2pwr(d,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + call carea_2pwr(dbh,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) capped_allom = .false. case(3) - d_eff = min(d,dbh_maxh) - call carea_2pwr(d_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) + dbh_eff = min(dbh,dbh_maxh) + call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,do_inverse) capped_allom = .true. case DEFAULT write(fates_log(),*) 'An undefined leaf allometry was specified: ', & @@ -499,10 +503,15 @@ subroutine carea_allom(d,nplant,site_spread,ipft,c_area,inverse) if (capped_allom .and. do_inverse) then - if (d_eff .lt. dbh_maxh) then - d = d_eff + if (dbh_eff .lt. dbh_maxh) then + dbh = dbh_eff else - d = fates_unset_r8 + ! In this situation, we are signaling to the + ! calling routine that we we cannot calculate + ! dbh from crown area, because we have already + ! hit the area cap, and the two are not proportional + ! anymore. hopefully, the calling routine has an alternative + dbh = fates_unset_r8 endif endif @@ -1872,20 +1881,20 @@ end subroutine h2d_martcano ! ============================================================================= - subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse) + subroutine carea_2pwr(dbh,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inverse) ! ============================================================================ ! Calculate area of ground covered by entire cohort. (m2) ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. ! ============================================================================ - real(r8),intent(inout) :: d ! diameter [cm] + real(r8),intent(inout) :: dbh ! diameter at breast (refernce) height [cm] real(r8),intent(in) :: spread ! site level relative spread score [0-1] real(r8),intent(in) :: d2bl_p2 ! parameter 2 in the diameter->bleaf allometry (exponent) real(r8),intent(in) :: d2bl_ediff ! area difference factor in the diameter-bleaf allometry (exponent) real(r8),intent(in) :: d2ca_min ! minimum diameter to crown area scaling factor real(r8),intent(in) :: d2ca_max ! maximum diameter to crown area scaling factor - real(r8),intent(inout) :: c_area ! crown area for one plant [m2] + real(r8),intent(inout) :: c_area ! crown area for one plant [m2] logical,intent(in) :: inverse ! if true, calculate dbh from crown area rather than its reverse real(r8) :: crown_area_to_dbh_exponent @@ -1911,9 +1920,9 @@ subroutine carea_2pwr(d,spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max,c_area,inver spreadterm = spread * d2ca_max + (1._r8 - spread) * d2ca_min if ( .not. inverse) then - c_area = spreadterm * d ** crown_area_to_dbh_exponent + c_area = spreadterm * dbh ** crown_area_to_dbh_exponent else - d = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent) + dbh = (c_area / spreadterm) ** (1./crown_area_to_dbh_exponent) endif end subroutine carea_2pwr From 7d49bef8337204fdeca11ca08a938283b8c50bdb Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Sat, 16 Feb 2019 14:32:24 -0800 Subject: [PATCH 03/28] reincorporated old conservation-of-dbh logic as an option in a case select --- biogeochem/EDCohortDynamicsMod.F90 | 96 +++++++++++++++++++++++------- 1 file changed, 76 insertions(+), 20 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 006090dc6d..5dcb3b02a0 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -91,6 +91,12 @@ module EDCohortDynamicsMod character(len=*), parameter, private :: sourcefile = & __FILE__ + + integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1 + integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2 + + integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh + ! 10/30/09: Created by Rosie Fisher !-------------------------------------------------------------------------------------! @@ -870,27 +876,77 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + nextc%n*nextc%canopy_trim)/newn - ! conserve total crown area during the fusion step, and then calculate - ! dbh of the fused cohort as that which conserves both crown area and - ! the dbh to crown area allometry. dbh will be updated in the next - ! growth step in the (likely) event that dbh to structural iomass - ! allometry is exceeded. if using a capped crown area allometry and - ! above the cap, then calculate as the weighted average of fusing - ! cohorts' dbh - - currentCohort%c_area = currentCohort%c_area + nextc%c_area - - call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& - currentCohort%c_area,inverse=.true.) - - if (abs(dbh-fates_unset_r8) Date: Mon, 25 Feb 2019 09:49:41 -0800 Subject: [PATCH 04/28] Increased the number of maximum leaf layers back up, to 30. Was creating crashes. --- main/EDTypesMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index ee2e574310..682df90aae 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -24,7 +24,7 @@ module EDTypesMod ! to understory layers (all layers that ! are not the top canopy layer) - integer, parameter :: nlevleaf = 20 ! number of leaf layers in canopy layer + integer, parameter :: nlevleaf = 30 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed ! the parameter file may determine that fewer ! are used, but this helps allocate scratch From e3a3b75db37c6dde149260a6e88005e1674048ff Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 14:43:18 -0700 Subject: [PATCH 05/28] Updating some error diagnostics and the precision thresholds for demotion/promotion --- biogeochem/EDCanopyStructureMod.F90 | 43 ++++++++++++++++++++--------- biogeochem/EDCohortDynamicsMod.F90 | 4 +-- 2 files changed, 32 insertions(+), 15 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a1326d482c..099d63e3bf 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -56,8 +56,14 @@ module EDCanopyStructureMod character(len=*), parameter, private :: sourcefile = & __FILE__ - real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation must be within this tolerance - real(r8), parameter :: area_check_precision = 1.0E-9_r8 ! Area conservation checks must be within this tolerance + real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation + ! will attempt to reduce errors + ! below this level + + real(r8), parameter :: area_check_precision = 1.0E-6_r8 ! Area conservation checks must + ! be within this absolute tolerance + real(r8), parameter :: area_check_rel_precision = 1.0E-3_r8 ! Area conservation checks must + ! be within this relative tolerance real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm ! can be roughly considered the same right? @@ -242,7 +248,8 @@ subroutine canopy_structure( currentSite , bc_in ) area_not_balanced = .false. do i_lyr = 1,z call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer(i_lyr)) - if( (arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_precision )then + if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. & + ((arealayer(i_lyr)-currentPatch%area) > area_check_precision ) ) then area_not_balanced = .true. endif enddo @@ -258,6 +265,7 @@ subroutine canopy_structure( currentSite , bc_in ) do i_lyr = 1,z write(fates_log(),*) 'layer: ',i_lyr,' area: ',arealayer(i_lyr) write(fates_log(),*) 'rel error: ',(arealayer(i_lyr)-currentPatch%area)/currentPatch%area + write(fates_log(),*) 'abs error: ',arealayer(i_lyr)-currentPatch%area enddo write(fates_log(),*) 'lat:',currentSite%lat write(fates_log(),*) 'lon:',currentSite%lon @@ -363,7 +371,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) demote_area = arealayer - currentPatch%area - if ( demote_area/currentPatch%area > area_target_precision ) then + if ( demote_area > nearzero ) then ! Is this layer currently over-occupied? ! In that case, we need to work out which cohorts to demote. @@ -493,7 +501,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) enddo exceedance_counter = 0 - do while(remainder_area/demote_area > area_target_precision ) + do while(remainder_area > area_target_precision ) ! Keep attempting to add exceedance to members that have not ! lost more area than they started with.. @@ -527,6 +535,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if( currentCohort%excl_weight < (currentCohort%c_area-nearzero) ) then sumweights = sumweights + currentCohort%excl_weight end if + + elseif(currentCohort%excl_weight > (currentCohort%c_area+area_target_precision) ) then + write(fates_log(),*) 'more area than the cohort wants to be domoted' + write(fates_log(),*) 'loss:',currentCohort%excl_weight + write(fates_log(),*) 'existing area:',currentCohort%c_area + write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if currentCohort => currentCohort%shorter @@ -542,7 +557,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'remainder_area:',remainder_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + end do end if @@ -629,6 +644,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'more area than the cohort has is being demoted' write(fates_log(),*) 'loss:',cc_loss write(fates_log(),*) 'existing area:',currentCohort%c_area + write(fates_log(),*) 'excess: ',cc_loss - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -702,7 +718,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer) - if ( abs(arealayer - currentPatch%area)/arealayer > area_target_precision ) then + if ( (abs(arealayer - currentPatch%area)/arealayer > area_check_rel_precision ) .or. & + (abs(arealayer - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'demotion did not trim area within tolerance' write(fates_log(),*) 'arealayer:',arealayer write(fates_log(),*) 'patch%area:',currentPatch%area @@ -712,7 +729,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) end if - + return end subroutine DemoteFromLayer @@ -770,7 +787,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! how much do we need to gain? promote_area = currentPatch%area - arealayer_current - if( promote_area/currentPatch%area > area_target_precision ) then + if( promote_area > nearzero ) then if(arealayer_below <= promote_area ) then @@ -930,7 +947,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) enddo exceedance_counter = 0 - do while(remainder_area/promote_area > area_target_precision ) + do while(remainder_area > area_target_precision ) sumweights = 0.0_r8 remainder_area_hold = remainder_area @@ -1076,14 +1093,14 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) - if ( abs(arealayer_current - currentPatch%area)/arealayer_current & - > area_target_precision ) then + if ((abs(arealayer_current - currentPatch%area)/arealayer_current > area_check_rel_precision ) .or. & + (abs(arealayer_current - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'promotion did not bring area within tolerance' write(fates_log(),*) 'arealayer:',arealayer_current write(fates_log(),*) 'patch%area:',currentPatch%area call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + end if end if diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0f9faaf5d6..d7b82a0530 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -103,7 +103,7 @@ module EDCohortDynamicsMod integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1 integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2 - integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh + integer, parameter, private :: cohort_fusion_conservation_method = conserve_dbh_and_number_not_crownarea ! 10/30/09: Created by Rosie Fisher !-------------------------------------------------------------------------------------! @@ -1003,7 +1003,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft,& currentCohort%c_area,inverse=.false.) ! - case(default) + case default write(fates_log(),*) 'FATES: Invalid choice for cohort_fusion_conservation_method' call endrun(msg=errMsg(sourcefile, __LINE__)) end select From 4b6a363e9040ecde23e660a42495aee9aada7bf0 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 14:43:18 -0700 Subject: [PATCH 06/28] Added some updates to area threshold checks in promotion/demotion, updated error messages --- biogeochem/EDCanopyStructureMod.F90 | 43 ++++++++++++++++++++--------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a1326d482c..099d63e3bf 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -56,8 +56,14 @@ module EDCanopyStructureMod character(len=*), parameter, private :: sourcefile = & __FILE__ - real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation must be within this tolerance - real(r8), parameter :: area_check_precision = 1.0E-9_r8 ! Area conservation checks must be within this tolerance + real(r8), parameter :: area_target_precision = 1.0E-11_r8 ! Area conservation + ! will attempt to reduce errors + ! below this level + + real(r8), parameter :: area_check_precision = 1.0E-6_r8 ! Area conservation checks must + ! be within this absolute tolerance + real(r8), parameter :: area_check_rel_precision = 1.0E-3_r8 ! Area conservation checks must + ! be within this relative tolerance real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm ! can be roughly considered the same right? @@ -242,7 +248,8 @@ subroutine canopy_structure( currentSite , bc_in ) area_not_balanced = .false. do i_lyr = 1,z call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer(i_lyr)) - if( (arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_precision )then + if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. & + ((arealayer(i_lyr)-currentPatch%area) > area_check_precision ) ) then area_not_balanced = .true. endif enddo @@ -258,6 +265,7 @@ subroutine canopy_structure( currentSite , bc_in ) do i_lyr = 1,z write(fates_log(),*) 'layer: ',i_lyr,' area: ',arealayer(i_lyr) write(fates_log(),*) 'rel error: ',(arealayer(i_lyr)-currentPatch%area)/currentPatch%area + write(fates_log(),*) 'abs error: ',arealayer(i_lyr)-currentPatch%area enddo write(fates_log(),*) 'lat:',currentSite%lat write(fates_log(),*) 'lon:',currentSite%lon @@ -363,7 +371,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) demote_area = arealayer - currentPatch%area - if ( demote_area/currentPatch%area > area_target_precision ) then + if ( demote_area > nearzero ) then ! Is this layer currently over-occupied? ! In that case, we need to work out which cohorts to demote. @@ -493,7 +501,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) enddo exceedance_counter = 0 - do while(remainder_area/demote_area > area_target_precision ) + do while(remainder_area > area_target_precision ) ! Keep attempting to add exceedance to members that have not ! lost more area than they started with.. @@ -527,6 +535,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if( currentCohort%excl_weight < (currentCohort%c_area-nearzero) ) then sumweights = sumweights + currentCohort%excl_weight end if + + elseif(currentCohort%excl_weight > (currentCohort%c_area+area_target_precision) ) then + write(fates_log(),*) 'more area than the cohort wants to be domoted' + write(fates_log(),*) 'loss:',currentCohort%excl_weight + write(fates_log(),*) 'existing area:',currentCohort%c_area + write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if currentCohort => currentCohort%shorter @@ -542,7 +557,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'remainder_area:',remainder_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + end do end if @@ -629,6 +644,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'more area than the cohort has is being demoted' write(fates_log(),*) 'loss:',cc_loss write(fates_log(),*) 'existing area:',currentCohort%c_area + write(fates_log(),*) 'excess: ',cc_loss - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -702,7 +718,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer) - if ( abs(arealayer - currentPatch%area)/arealayer > area_target_precision ) then + if ( (abs(arealayer - currentPatch%area)/arealayer > area_check_rel_precision ) .or. & + (abs(arealayer - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'demotion did not trim area within tolerance' write(fates_log(),*) 'arealayer:',arealayer write(fates_log(),*) 'patch%area:',currentPatch%area @@ -712,7 +729,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) end if - + return end subroutine DemoteFromLayer @@ -770,7 +787,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! how much do we need to gain? promote_area = currentPatch%area - arealayer_current - if( promote_area/currentPatch%area > area_target_precision ) then + if( promote_area > nearzero ) then if(arealayer_below <= promote_area ) then @@ -930,7 +947,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) enddo exceedance_counter = 0 - do while(remainder_area/promote_area > area_target_precision ) + do while(remainder_area > area_target_precision ) sumweights = 0.0_r8 remainder_area_hold = remainder_area @@ -1076,14 +1093,14 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) - if ( abs(arealayer_current - currentPatch%area)/arealayer_current & - > area_target_precision ) then + if ((abs(arealayer_current - currentPatch%area)/arealayer_current > area_check_rel_precision ) .or. & + (abs(arealayer_current - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'promotion did not bring area within tolerance' write(fates_log(),*) 'arealayer:',arealayer_current write(fates_log(),*) 'patch%area:',currentPatch%area call endrun(msg=errMsg(sourcefile, __LINE__)) end if - + end if end if From fb9aaf3cc54eaa22dcbe5b968406493ff9c59b99 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 15:35:42 -0800 Subject: [PATCH 07/28] Added new algorithm for scaling promotion/demotion weights to within cohort crown areas --- biogeochem/EDCanopyStructureMod.F90 | 303 +++++++++++++++------------- 1 file changed, 167 insertions(+), 136 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 099d63e3bf..261b890771 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -349,20 +349,17 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] - real(r8) :: lossarea + real(r8) :: scale_factor ! for prob. exclusion - scales weight to a fraction + real(r8) :: scale_factor_min ! "" minimum before exeedance of 1 + real(r8) :: scale_factor_res ! "" applied to residual areas + real(r8) :: area_res ! residual area to demote after weakest cohort hits max real(r8) :: newarea real(r8) :: demote_area - real(r8) :: remainder_area - real(r8) :: remainder_area_hold real(r8) :: sumweights - real(r8) :: sumweights_old real(r8) :: sumequal ! for rank-ordered same-size cohorts ! this tallies their excluded area real(r8) :: arealayer ! the area of the current canopy layer - integer :: exceedance_counter ! when seeking to rebalance demotion exceedance - ! keep a loop counter to check for hangs logical :: tied_size_with_neighbors - type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts ! First, determine how much total canopy area we have in this layer @@ -393,8 +390,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! inverse size to a constant power ! ---------------------------------------------------------- - currentCohort%excl_weight = & - currentCohort%n/(currentCohort%hite**ED_val_comp_excln) + currentCohort%excl_weight = 1._r8 / (currentCohort%hite**ED_val_comp_excln) sumweights = sumweights + currentCohort%excl_weight else @@ -481,87 +477,104 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! area of those cohorts that have not reached their total canopy area ! yet. + + + if (ED_val_comp_excln .ge. 0.0_r8 ) then - remainder_area = 0.0_r8 - sumweights_old = 0.0_r8 + scale_factor_min = 1.e10_r8 + scale_factor = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then - currentCohort%excl_weight = demote_area*currentCohort%excl_weight/sumweights - if( currentCohort%excl_weight > currentCohort%c_area ) then - remainder_area = remainder_area + currentCohort%excl_weight-currentCohort%c_area - currentCohort%excl_weight = currentCohort%c_area - else - sumweights_old = sumweights_old + currentCohort%excl_weight - end if + + currentCohort%excl_weight = currentCohort%excl_weight/sumweights + if( 1._r8/currentCohort%excl_weight < scale_factor_min ) & + scale_factor_min = 1._r8/currentCohort%excl_weight + + scale_factor = scale_factor + currentCohort%excl_weight * currentCohort%c_area + endif currentCohort => currentCohort%shorter enddo - exceedance_counter = 0 - do while(remainder_area > area_target_precision ) - - ! Keep attempting to add exceedance to members that have not - ! lost more area than they started with.. - - sumweights = 0.0_r8 - remainder_area_hold = remainder_area - currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - if(currentCohort%canopy_layer == i_lyr)then - if ( currentCohort%excl_weight < (currentCohort%c_area-nearzero) ) then + ! This is the factor by which we need to multiply + ! the demotion probabilities, so the sum result equals + ! the total amount to demote + scale_factor = demote_area/scale_factor + - ! Calculate how much exceeded demotion from filled - ! cohorts must be transferred to this cohort - ! Two requirements: 1) the tacked-on demotion can not also - ! exceed this cohort's area. And, 2) the tacked-on - ! demotion can't exceed the amount left + if(scale_factor <= scale_factor_min) then - cc_loss = min(remainder_area, & - remainder_area_hold * currentCohort%excl_weight/sumweights_old, & - currentCohort%c_area-currentCohort%excl_weight) + ! Trivial case, all of the demotion fractions + ! are less than 1. + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == i_lyr) then + currentCohort%excl_weight = currentCohort%c_area * currentCohort%excl_weight * scale_factor + + if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%excl_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area too big (1)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight + write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if - ! Reduce the remainder_area - remainder_area = remainder_area - cc_loss + endif + currentCohort => currentCohort%shorter + enddo + + else - ! Update this cohorts exclusion - currentCohort%excl_weight = currentCohort%excl_weight + cc_loss - - ! Update the sum of weights for cohorts where exceedance can - ! still be spread to - if( currentCohort%excl_weight < (currentCohort%c_area-nearzero) ) then - sumweights = sumweights + currentCohort%excl_weight - end if + ! Non-trivial case, at least 1 cohort's demotion + ! rate would exceed its area, given the trivial scale factor + + area_res = 0._r8 + scale_factor_res = 0._r8 + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == i_lyr) then + area_res = area_res + & + currentCohort%c_area*currentCohort%excl_weight*scale_factor_min + scale_factor_res = scale_factor_res + & + currentCohort%c_area * (1._r8 - (currentCohort%excl_weight * scale_factor_min)) + endif + currentCohort => currentCohort%shorter + enddo + + area_res = demote_area - area_res + + scale_factor_res = area_residual / scale_factor_res + + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == i_lyr) then + + currentCohort%excl_weight = currentCohort%c_area * & + (currentCohort%excl_weight * scale_factor_min + & + (1._r8 - (currentCohort%excl_weight*scale_factor_min) ) * scale_factor_res) - elseif(currentCohort%excl_weight > (currentCohort%c_area+area_target_precision) ) then - write(fates_log(),*) 'more area than the cohort wants to be domoted' - write(fates_log(),*) 'loss:',currentCohort%excl_weight - write(fates_log(),*) 'existing area:',currentCohort%c_area + if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%excl_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if - currentCohort => currentCohort%shorter - end do - ! Update the sum of weights for the next loop - sumweights_old = sumweights - - exceedance_counter = exceedance_counter + 1 - if( exceedance_counter > 100 ) then - write(fates_log(),*) 'It is taking a long time to distribute demotion exceedance' - write(fates_log(),*) '(ie greater than c_area) to neighbors... exiting' - write(fates_log(),*) 'sumweights:',sumweights - write(fates_log(),*) 'remainder_area:',remainder_area - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - end do - end if + endif + currentCohort => currentCohort%shorter + enddo + + end if + end if + ! Weights have been calculated. Now move them to the lower layer currentCohort => currentPatch%tallest @@ -757,16 +770,15 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) type(ed_cohort_type), pointer :: nextc ! the next cohort, or used for looping ! cohorts against the current - + real(r8) :: scale_factor ! for prob. exclusion - scales weight to a fraction + real(r8) :: scale_factor_min ! "" minimum before exeedance of 1 + real(r8) :: scale_factor_res ! "" applied to residual areas + real(r8) :: area_res ! residual area to demote after weakest cohort hits max real(r8) :: promote_area real(r8) :: newarea real(r8) :: sumweights - real(r8) :: sumweights_old real(r8) :: sumequal ! for tied cohorts, the sum of weights in ! their group - integer :: exceedance_counter - real(r8) :: remainder_area - real(r8) :: remainder_area_hold real(r8) :: cc_gain ! cohort crown area gain in promotion (m2) real(r8) :: arealayer_current ! area (m2) of the current canopy layer real(r8) :: arealayer_below ! area (m2) of the layer below the current layer @@ -777,7 +789,6 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) real(r8) :: struct_c ! structure carbon [kg] logical :: tied_size_with_neighbors - type(ed_cohort_type), pointer :: cohort_tosearch_relative_to, cohort_tocompare_to real(r8) :: total_crownarea_of_tied_cohorts call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) @@ -840,7 +851,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if (ED_val_comp_excln .ge. 0.0_r8 ) then ! normal (stochastic) case, as above. - currentCohort%prom_weight = currentCohort%n*currentCohort%hite**ED_val_comp_excln + currentCohort%prom_weight = currentCohort%hite**ED_val_comp_excln sumweights = sumweights + currentCohort%prom_weight else @@ -923,84 +934,104 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! And then a few rounds where we pre-calculate the promotion areas ! and adjust things if the promoted area wants to be greater than ! what is available. - if (ED_val_comp_excln .ge. 0.0_r8 ) then - remainder_area = 0.0_r8 - sumweights_old = 0.0_r8 - - currentCohort => currentPatch%tallest !start from the tallest cohort + scale_factor_min = 1.e10_r8 + scale_factor = 0._r8 + currentCohort => currentPatch%tallest do while (associated(currentCohort)) - if(currentCohort%canopy_layer == i_lyr+1) then !still looking at the layer beneath. - - currentCohort%prom_weight = promote_area*currentCohort%prom_weight/sumweights + + if(currentCohort%canopy_layer == (i_lyr+1) ) then - if( currentCohort%prom_weight > currentCohort%c_area ) then - remainder_area = remainder_area + currentCohort%prom_weight-currentCohort%c_area - currentCohort%prom_weight = currentCohort%c_area - else - sumweights_old = sumweights_old + currentCohort%prom_weight - end if + currentCohort%prom_weight = currentCohort%prom_weight/sumweights + if( 1._r8/currentCohort%prom_weight < scale_factor_min ) & + scale_factor_min = 1._r8/currentCohort%prom_weight + + scale_factor = scale_factor + currentCohort%prom_weight * currentCohort%c_area endif currentCohort => currentCohort%shorter enddo - - exceedance_counter = 0 - do while(remainder_area > area_target_precision ) + + ! This is the factor by which we need to multiply + ! the demotion probabilities, so the sum result equals + ! the total amount to demote + scale_factor = promote_area/scale_factor + + + if(scale_factor <= scale_factor_min) then + + ! Trivial case, all of the demotion fractions + ! are less than 1. - sumweights = 0.0_r8 - remainder_area_hold = remainder_area currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - if(currentCohort%canopy_layer == i_lyr+1)then - - if ( currentCohort%prom_weight < (currentCohort%c_area-nearzero) ) then - - ! Calculate how much exceeded promotion from filled - ! cohorts must be transferred to this cohort - ! Two requirements: 1) the tacked-on promotion can not also - ! exceed this cohort's area. And, 2) the tacked-on - ! promotion can't exceed the amount left - ! This is how promotion is transferred from this cohort - - cc_gain = min(remainder_area, & - remainder_area_hold * currentCohort%prom_weight/sumweights_old, & - currentCohort%c_area-currentCohort%prom_weight) - - - ! Reduce the remainder_area - remainder_area = remainder_area - cc_gain - - ! Update this cohort's promotion - currentCohort%prom_weight = currentCohort%prom_weight + cc_gain - - ! Update the sum of weights for cohorts where exceedance can - ! still be spread to - if( currentCohort%prom_weight < (currentCohort%c_area-nearzero) ) then - sumweights = sumweights + currentCohort%prom_weight - end if - + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == (i_lyr+1) ) then + currentCohort%prom_weight = currentCohort%c_area * currentCohort%prom_weight * scale_factor + + if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area too big (1)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight + write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if - currentCohort => currentCohort%shorter - end do - sumweights_old = sumweights + + endif + currentCohort => currentCohort%shorter + enddo - exceedance_counter = exceedance_counter + 1 - if( exceedance_counter > 100 ) then - write(fates_log(),*) 'It is taking a long time to distribute promotion exceedance' - write(fates_log(),*) '(ie greater than c_area) to neighbors... exiting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if + else - end do + ! Non-trivial case, at least 1 cohort's demotion + ! rate would exceed its area, given the trivial scale factor + + area_res = 0._r8 + scale_factor_res = 0._r8 + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == (i_lyr+1) ) then + area_res = area_res + & + currentCohort%c_area*currentCohort%prom_weight*scale_factor_min + scale_factor_res = scale_factor_res + & + currentCohort%c_area * (1._r8 - (currentCohort%prom_weight * scale_factor_min)) + endif + currentCohort => currentCohort%shorter + enddo + + area_res = promote_area - area_res + + scale_factor_res = area_residual / scale_factor_res + + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + if(currentCohort%canopy_layer == (i_lyr+1)) then + + currentCohort%prom_weight = currentCohort%c_area * & + (currentCohort%prom_weight * scale_factor_min + & + (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * scale_factor_res) + + if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight + write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + endif + currentCohort => currentCohort%shorter + enddo + + end if end if + currentCohort => currentPatch%tallest do while (associated(currentCohort)) - !All the trees in this layer need to promote some area upwards... From 7ff6e9d281ad4a14ab4e520b4e377eb463b48ff6 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 15:49:11 -0800 Subject: [PATCH 08/28] Small fixes and checks to the updated promotion/demotion logic --- biogeochem/EDCanopyStructureMod.F90 | 40 ++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 261b890771..a03902fbd3 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -548,7 +548,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) area_res = demote_area - area_res - scale_factor_res = area_residual / scale_factor_res + scale_factor_res = area_res / scale_factor_res currentCohort => currentPatch%tallest do while (associated(currentCohort)) @@ -575,6 +575,24 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) end if + + ! lets perform a check and see if the demotions meet the demand + sumweights = 0._r8 + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + sumweights = sumweights + currentCohort%excl_weight + currentCohort => currentCohort%shorter + end do + + if (abs(sumweights - demote_area) > area_check_precision ) then + write(fates_log(),*) 'demotions dont add up' + write(fates_log(),*) 'sum demotions: ',sumweights + write(fates_log(),*) 'area needed to be demoted: ',demote_area + write(fates_log(),*) 'excess: ',sumweights - demote_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Weights have been calculated. Now move them to the lower layer currentCohort => currentPatch%tallest @@ -1002,7 +1020,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) area_res = promote_area - area_res - scale_factor_res = area_residual / scale_factor_res + scale_factor_res = area_res / scale_factor_res currentCohort => currentPatch%tallest do while (associated(currentCohort)) @@ -1029,7 +1047,23 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) end if - + + ! lets perform a check and see if the promotions meet the demand + sumweights = 0._r8 + currentCohort => currentPatch%tallest + do while (associated(currentCohort)) + sumweights = sumweights + currentCohort%prom_weight + currentCohort => currentCohort%shorter + end do + + if (abs(sumweights - promote_area) > area_check_precision ) then + write(fates_log(),*) 'promotions dont add up' + write(fates_log(),*) 'sum promotions: ',sumweights + write(fates_log(),*) 'area needed to be promoted: ',promote_area + write(fates_log(),*) 'excess: ',sumweights - promote_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + currentCohort => currentPatch%tallest do while (associated(currentCohort)) From f2896038819f0845e08eb443f0b52bf3eefd1329 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 20:15:26 -0700 Subject: [PATCH 09/28] Various fixes to promotion/demotion. Notably, fixed filter on zeroing out cohorts that surpass maxlayers during demotion. --- biogeochem/EDCanopyStructureMod.F90 | 126 ++++++++++++++++++---------- 1 file changed, 84 insertions(+), 42 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a03902fbd3..b53cfb0074 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -195,6 +195,10 @@ subroutine canopy_structure( currentSite , bc_in ) ! the layers below. ! --------------------------------------------------------------------------- + ! Its possible that before we even enter this scheme + ! some cohort numbers are very low. Terminate them. + call terminate_cohorts(currentSite, currentPatch, 1) + ! Calculate how many layers we have in this canopy ! This also checks the understory to see if its crown ! area is large enough to warrant a temporary sub-understory layer @@ -204,7 +208,8 @@ subroutine canopy_structure( currentSite , bc_in ) call DemoteFromLayer(currentSite, currentPatch, i_lyr) end do - ! Remove cohorts that are incredibly sparse + ! After demotions, we may then again have cohorts that are very very + ! very sparse, remove them call terminate_cohorts(currentSite, currentPatch, 1) call fuse_cohorts(currentSite, currentPatch, bc_in) @@ -368,7 +373,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) demote_area = arealayer - currentPatch%area - if ( demote_area > nearzero ) then + if ( demote_area > area_target_precision ) then ! Is this layer currently over-occupied? ! In that case, we need to work out which cohorts to demote. @@ -380,6 +385,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) + + if(currentCohort%c_area<0._r8)then + write(fates_log(),*) 'negative c_area stage 1d: ',currentCohort%dbh,i_lyr,currentCohort%n, & + currentSite%spread,currentCohort%pft,currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + if( currentCohort%canopy_layer == i_lyr)then @@ -473,12 +485,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! If this is probabalistic demotion, we need to do a round of normalization. ! And then a few rounds where we pre-calculate the demotion areas ! and adjust things if the demoted area wants to be greater than - ! what is available. sumweights_old is used to sum up the exclusion - ! area of those cohorts that have not reached their total canopy area - ! yet. - - - + ! what is available. The math is too hard to explain here, see + ! the tech note section on promotion/demotion. if (ED_val_comp_excln .ge. 0.0_r8 ) then @@ -519,6 +527,10 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) (currentCohort%excl_weight < 0._r8) ) then write(fates_log(),*) 'exclusion area too big (1)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'dbh: ',currentCohort%dbh + write(fates_log(),*) 'n: ',currentCohort%n + write(fates_log(),*) 'spread: ',currentSite%spread + write(fates_log(),*) 'pft: ',currentCohort%pft write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area call endrun(msg=errMsg(sourcefile, __LINE__)) @@ -580,7 +592,9 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) sumweights = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) - sumweights = sumweights + currentCohort%excl_weight + if(currentCohort%canopy_layer == i_lyr) then + sumweights = sumweights + currentCohort%excl_weight + end if currentCohort => currentCohort%shorter end do @@ -596,17 +610,16 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! Weights have been calculated. Now move them to the lower layer currentCohort => currentPatch%tallest - do while (associated(currentCohort)) - - cc_loss = currentCohort%excl_weight - - leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) - store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) - fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) - sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) + do while (associated(currentCohort)) - if(currentCohort%canopy_layer == i_lyr .and. cc_loss>nearzero )then + if(currentCohort%canopy_layer == i_lyr )then + + cc_loss = currentCohort%excl_weight + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) if ( (cc_loss-currentCohort%c_area) > -nearzero .and. & (cc_loss-currentCohort%c_area) < area_target_precision ) then @@ -622,7 +635,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentSite%demotion_carbonflux = currentSite%demotion_carbonflux + & (leaf_c + store_c + fnrt_c + sapw_c + struct_c) * currentCohort%n - elseif(cc_loss > nearzero .and. cc_loss < currentCohort%c_area )then + elseif( (cc_loss < currentCohort%c_area) .and. & + (cc_loss > area_target_precision) ) then ! If only part of the cohort is demoted ! then it must be split (little more complicated) @@ -631,6 +645,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! conserve total number density of the original. The copy ! remains in the upper-story. The original is the one ! demoted to the understory + allocate(copyc) @@ -639,10 +654,21 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call InitHydrCohort(currentSite,copyc) endif call copy_cohort(currentCohort, copyc) - + + if(currentCohort%n < 0._r8) then + write(fates_log(),*) 'negatives (0_?): ',currentCohort%n,currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + newarea = currentCohort%c_area - cc_loss copyc%n = currentCohort%n*newarea/currentCohort%c_area currentCohort%n = currentCohort%n - copyc%n + + if(copyc%n < 0._r8 .or. currentCohort%n < 0._r8) then + write(fates_log(),*) 'negatives?: ',newarea,cc_loss,copyc%n,currentCohort%n,currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + copyc%canopy_layer = i_lyr !the taller cohort is the copy @@ -659,6 +685,16 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) + ! Calculate how much area loss from recalculation + if( abs(copyc%c_area-newarea)>area_target_precision .or. & + abs(currentCohort%c_area-cc_loss)>area_target_precision) then + write(fates_log(),*) 'recalculation losses?',newarea-copyc%c_area,newarea, & + copyc%c_area,currentCohort%c_area-cc_loss,currentCohort%c_area,& + cc_loss + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + !----------- Insert copy into linked list ------------------------! copyc%shorter => currentCohort if(associated(currentCohort%taller))then @@ -681,7 +717,8 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) end if ! kill the ones which go into canopy layers that are not allowed - if(i_lyr+1 > nclmax)then + + if(currentCohort%canopy_layer>nclmax )then ! put the litter from the terminated cohorts into the fragmenting pools do i_cwd=1,ncwd @@ -754,7 +791,9 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'demotion did not trim area within tolerance' write(fates_log(),*) 'arealayer:',arealayer write(fates_log(),*) 'patch%area:',currentPatch%area - write(fates_log(),*) 'error:',abs(arealayer - currentPatch%area) + write(fates_log(),*) 'ilayer: ',i_lyr + write(fates_log(),*) 'bias:',arealayer - currentPatch%area + write(fates_log(),*) 'demote_area:',demote_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -816,7 +855,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! how much do we need to gain? promote_area = currentPatch%area - arealayer_current - if( promote_area > nearzero ) then + if( promote_area > area_target_precision ) then if(arealayer_below <= promote_area ) then @@ -1052,7 +1091,9 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) sumweights = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) - sumweights = sumweights + currentCohort%prom_weight + if(currentCohort%canopy_layer == (i_lyr+1)) then + sumweights = sumweights + currentCohort%prom_weight + end if currentCohort => currentCohort%shorter end do @@ -1069,29 +1110,29 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) !All the trees in this layer need to promote some area upwards... - if(currentCohort%canopy_layer == i_lyr+1)then - - cc_gain = currentCohort%prom_weight + if( (currentCohort%canopy_layer == i_lyr+1) ) then + + cc_gain = currentCohort%prom_weight + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) + fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) + sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) + struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) if ( (cc_gain-currentCohort%c_area) > -nearzero .and. & (cc_gain-currentCohort%c_area) < area_target_precision ) then - + currentCohort%canopy_layer = i_lyr ! keep track of number and biomass of promoted cohort currentSite%promotion_rate(currentCohort%size_class) = & currentSite%promotion_rate(currentCohort%size_class) + currentCohort%n - leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) - store_c = currentCohort%prt%GetState(store_organ,all_carbon_elements) - fnrt_c = currentCohort%prt%GetState(fnrt_organ,all_carbon_elements) - sapw_c = currentCohort%prt%GetState(sapw_organ,all_carbon_elements) - struct_c = currentCohort%prt%GetState(struct_organ,all_carbon_elements) - currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + & (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * currentCohort%n - elseif ( cc_gain > nearzero .and. cc_gain < currentCohort%c_area) then + elseif ( (cc_gain < currentCohort%c_area) .and. & + (cc_gain > area_target_precision) ) then allocate(copyc) @@ -1112,6 +1153,13 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! number of individuals in cohort remaining in understorey currentCohort%n = currentCohort%n - copyc%n + + if(copyc%n < 0._r8 .or. currentCohort%n < 0._r8) then + write(fates_log(),*) 'negatives (2)?: ',newarea,cc_gain,copyc%n,currentCohort%n + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + currentCohort%canopy_layer = i_lyr + 1 ! keep current cohort in the understory. copyc%canopy_layer = i_lyr ! promote copy to the higher canopy layer. @@ -1119,12 +1167,6 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) currentSite%promotion_rate(copyc%size_class) = & currentSite%promotion_rate(copyc%size_class) + copyc%n - leaf_c = copyc%prt%GetState(leaf_organ,all_carbon_elements) - store_c = copyc%prt%GetState(store_organ,all_carbon_elements) - fnrt_c = copyc%prt%GetState(fnrt_organ,all_carbon_elements) - sapw_c = copyc%prt%GetState(sapw_organ,all_carbon_elements) - struct_c = copyc%prt%GetState(struct_organ,all_carbon_elements) - currentSite%promotion_carbonflux = currentSite%promotion_carbonflux + & (leaf_c + fnrt_c + store_c + sapw_c + struct_c) * copyc%n From bcc719d38b50fdbcb4569b07aa83118c50e052fb Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 23:18:36 -0700 Subject: [PATCH 10/28] Expanded leaf layers to 30 --- main/EDTypesMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index c6a22ad6b1..4fc7ee7d71 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -24,7 +24,7 @@ module EDTypesMod ! to understory layers (all layers that ! are not the top canopy layer) - integer, parameter :: nlevleaf = 20 ! number of leaf layers in canopy layer + integer, parameter :: nlevleaf = 30 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed ! the parameter file may determine that fewer ! are used, but this helps allocate scratch From 999122741844beb905c3c4695dfee24349bef240 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 28 Feb 2019 23:22:07 -0700 Subject: [PATCH 11/28] Changed flag back to conservation of crown area --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 7066224beb..ccde161014 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -103,7 +103,7 @@ module EDCohortDynamicsMod integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1 integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2 - integer, parameter, private :: cohort_fusion_conservation_method = conserve_dbh_and_number_not_crownarea + integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh ! 10/30/09: Created by Rosie Fisher !-------------------------------------------------------------------------------------! From fdebaa755c77cf68e1ac63019b8b0b596203e54e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 1 Mar 2019 13:49:00 -0700 Subject: [PATCH 12/28] increased number of leaf levels to 40 after overflows --- main/EDTypesMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 4fc7ee7d71..b04cea860e 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -24,7 +24,7 @@ module EDTypesMod ! to understory layers (all layers that ! are not the top canopy layer) - integer, parameter :: nlevleaf = 30 ! number of leaf layers in canopy layer + integer, parameter :: nlevleaf = 40 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed ! the parameter file may determine that fewer ! are used, but this helps allocate scratch From 53775ba321ebd7379133a4bb2a951c0e4c38359f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 1 Mar 2019 15:30:29 -0800 Subject: [PATCH 13/28] Minor comment updates --- biogeochem/EDCanopyStructureMod.F90 | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index b53cfb0074..da5ea89ec2 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -398,8 +398,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if (ED_val_comp_excln .ge. 0.0_r8 ) then ! ---------------------------------------------------------- - ! normal (stochastic) case. weight cohort demotion by - ! inverse size to a constant power + ! Stochastic method. + ! Weight cohort demotion by inverse size to a constant power. + ! In this hypothesis, it is assumed that even the tallest + ! cohorts have a chance (although smaller) of being forced + ! to the understory. ! ---------------------------------------------------------- currentCohort%excl_weight = 1._r8 / (currentCohort%hite**ED_val_comp_excln) @@ -510,13 +513,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) ! This is the factor by which we need to multiply ! the demotion probabilities, so the sum result equals ! the total amount to demote + scale_factor = demote_area/scale_factor if(scale_factor <= scale_factor_min) then - ! Trivial case, all of the demotion fractions - ! are less than 1. + ! Trivial case, all of the demotion fractions are less than 1. currentCohort => currentPatch%tallest do while (associated(currentCohort)) @@ -588,7 +591,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) end if - ! lets perform a check and see if the demotions meet the demand + ! perform a check and see if the demotions meet the demand sumweights = 0._r8 currentCohort => currentPatch%tallest do while (associated(currentCohort)) @@ -907,7 +910,11 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if(currentCohort%canopy_layer == i_lyr+1)then !look at the cohorts in the canopy layer below... if (ED_val_comp_excln .ge. 0.0_r8 ) then - ! normal (stochastic) case, as above. + + ! ------------------------------------------------------------------ + ! Stochastic case, as above (in demotion portion of code) + ! ------------------------------------------------------------------ + currentCohort%prom_weight = currentCohort%hite**ED_val_comp_excln sumweights = sumweights + currentCohort%prom_weight else @@ -991,6 +998,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! And then a few rounds where we pre-calculate the promotion areas ! and adjust things if the promoted area wants to be greater than ! what is available. + if (ED_val_comp_excln .ge. 0.0_r8 ) then scale_factor_min = 1.e10_r8 @@ -1028,7 +1036,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'exclusion area too big (1)' + write(fates_log(),*) 'promotion area too big (1)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area @@ -1071,7 +1079,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'exclusion area error (2)' + write(fates_log(),*) 'promotion area error (2)' write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area From 9c282067775bb44ba67748b78c016a2f78aac3b3 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 4 Mar 2019 12:15:32 -0700 Subject: [PATCH 14/28] Updated diagnostics on LAI exceedance check --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- biogeochem/FatesAllometryMod.F90 | 10 ++++++++++ main/EDTypesMod.F90 | 2 +- 3 files changed, 12 insertions(+), 2 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index ccde161014..7066224beb 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -103,7 +103,7 @@ module EDCohortDynamicsMod integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1 integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2 - integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh + integer, parameter, private :: cohort_fusion_conservation_method = conserve_dbh_and_number_not_crownarea ! 10/30/09: Created by Rosie Fisher !-------------------------------------------------------------------------------------! diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index ea11b292ad..94bdc6c715 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -744,8 +744,18 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_la write(fates_log(),*) 'The leaf and stem are predicted for a cohort, maxed out the array size' write(fates_log(),*) 'lai: ',treelai write(fates_log(),*) 'sai: ',tree_sai + write(fates_log(),*) 'target_lai: ',target_lai write(fates_log(),*) 'lai+sai: ',treelai+tree_sai write(fates_log(),*) 'nlevleaf,dinc_ed,nlevleaf*dinc_ed :',nlevleaf,dinc_ed,nlevleaf*dinc_ed + write(fates_log(),*) 'pft: ',pft + write(fates_log(),*) 'n: ',nplant + write(fates_log(),*) 'c_area: ',c_area + write(fates_log(),*) 'dbh: ',dbh + write(fates_log(),*) 'canopy_trim: ',canopy_trim + write(fates_log(),*) 'target_bleaf: ',target_bleaf + write(fates_log(),*) 'canopy layer: ',cl + write(fates_log(),*) 'canopy_tlai: ',canopy_lai(:) + write(fates_log(),*) 'vcmax25top: ',vcmax25top call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index b04cea860e..4fc7ee7d71 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -24,7 +24,7 @@ module EDTypesMod ! to understory layers (all layers that ! are not the top canopy layer) - integer, parameter :: nlevleaf = 40 ! number of leaf layers in canopy layer + integer, parameter :: nlevleaf = 30 ! number of leaf layers in canopy layer integer, parameter :: maxpft = 15 ! maximum number of PFTs allowed ! the parameter file may determine that fewer ! are used, but this helps allocate scratch From 7e1b7d41b348cbfb76295c64926522c024175320 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 4 Mar 2019 17:22:50 -0700 Subject: [PATCH 15/28] adding a lai check function --- biogeochem/EDCohortDynamicsMod.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 7066224beb..2b736ada24 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1022,6 +1022,13 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end select + + ! Perform check on LAI (specifically, see if LAI is greater than target) + call check_lai(leaf_c, currentCohort%dbh,currentCohort%pft, currentCohort%c_area, & + currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai,currentCohort%vcmax25top, 0 ) + + call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) From 5ad22a1b636b741a2d87084816f9eb1cf232223d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 4 Mar 2019 16:41:35 -0800 Subject: [PATCH 16/28] Updating diagnostics in tree_sai --- biogeochem/EDCohortDynamicsMod.F90 | 15 ++++++++++----- biogeochem/EDPhysiologyMod.F90 | 3 ++- biogeochem/FatesAllometryMod.F90 | 7 +++++-- main/EDMainMod.F90 | 12 ++++++++++-- 4 files changed, 27 insertions(+), 10 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 2b736ada24..577e25d1b1 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -277,7 +277,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine new_cohort%treesai = tree_sai(new_cohort%pft, new_cohort%dbh, new_cohort%canopy_trim, & new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, & - patchptr%canopy_layer_tlai, new_cohort%treelai,new_cohort%vcmax25top ) + patchptr%canopy_layer_tlai, new_cohort%treelai,new_cohort%vcmax25top,2 ) new_cohort%lai = new_cohort%treelai * new_cohort%c_area/patchptr%area @@ -1022,11 +1022,16 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end select + + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + + currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & + currentCohort%vcmax25top) + currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,1 ) - ! Perform check on LAI (specifically, see if LAI is greater than target) - call check_lai(leaf_c, currentCohort%dbh,currentCohort%pft, currentCohort%c_area, & - currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai,currentCohort%vcmax25top, 0 ) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 80738a1bd8..1e4a69d0cd 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -259,7 +259,8 @@ subroutine trim_canopy( currentSite ) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top ) + currentPatch%canopy_layer_tlai, currentCohort%treelai, & + currentCohort%vcmax25top,0 ) currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 94bdc6c715..0a6edd9909 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -712,7 +712,8 @@ end function tree_lai ! ============================================================================ - real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_lai, treelai, vcmax25top ) + real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, & + canopy_lai, treelai, vcmax25top, call_id ) ! ============================================================================ ! SAI of individual trees is a function of the LAI of individual trees @@ -728,7 +729,8 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_la ! each canopy layer real(r8), intent(in) :: treelai ! tree LAI for checking purposes only real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown - + integer,intent(in) :: call_id ! flag specifying where this is called + ! from real(r8) :: target_bleaf real(r8) :: target_lai @@ -748,6 +750,7 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, canopy_la write(fates_log(),*) 'lai+sai: ',treelai+tree_sai write(fates_log(),*) 'nlevleaf,dinc_ed,nlevleaf*dinc_ed :',nlevleaf,dinc_ed,nlevleaf*dinc_ed write(fates_log(),*) 'pft: ',pft + write(fates_log(),*) 'call id: ',call_id write(fates_log(),*) 'n: ',nplant write(fates_log(),*) 'c_area: ',c_area write(fates_log(),*) 'dbh: ',dbh diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index cc366cc0fe..dae9699c30 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -48,7 +48,7 @@ module EDMainMod use FatesPlantHydraulicsMod , only : updateSizeDepTreeHydStates use FatesPlantHydraulicsMod , only : initTreeHydStates use FatesPlantHydraulicsMod , only : updateSizeDepRhizHydProps - use FatesAllometryMod , only : h_allom + use FatesAllometryMod , only : h_allom,tree_sai,tree_lai use FatesPlantHydraulicsMod , only : updateSizeDepRhizHydStates use EDLoggingMortalityMod , only : IsItLoggingTime use FatesGlobals , only : endrun => fates_endrun @@ -265,6 +265,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) real(r8) :: cohort_biomass_store ! remembers the biomass in the cohort for balance checking real(r8) :: dbh_old ! dbh of plant before daily PRT [cm] real(r8) :: hite_old ! height of plant before daily PRT [m] + real(r8) :: leaf_c !----------------------------------------------------------------------- @@ -353,7 +354,14 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! routine is also called following fusion call UpdateCohortBioPhysRates(currentCohort) - + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & + currentCohort%vcmax25top) + currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,3 ) + ! Transfer all reproductive tissues into seed production call PRTReproRelease(currentCohort%prt,repro_organ,carbon12_element, & 1.0_r8, currentCohort%seed_prod) From 3204b39ee2e98c6866674b234395809198ffdf0e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 6 Mar 2019 13:42:37 -0700 Subject: [PATCH 17/28] Now tying allometry in grasses to root biomass. ie if root biomass exceeds its target, dbh will increase to match. This will allow grasses to not force growth, and should reduce divergence from allometries. --- biogeochem/EDCanopyStructureMod.F90 | 3 +- biogeochem/EDCohortDynamicsMod.F90 | 10 +- biogeochem/EDPhysiologyMod.F90 | 1 - biogeochem/FatesAllometryMod.F90 | 151 ++++++++++++++++++--------- main/EDMainMod.F90 | 19 ++++ parteh/PRTAllometricCarbonMod.F90 | 156 +++++++++++++++++----------- 6 files changed, 223 insertions(+), 117 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index c7ffec0b06..16108301f9 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1553,7 +1553,8 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & - currentPatch%canopy_layer_tlai, currentCohort%treelai ,currentCohort%vcmax25top) + currentPatch%canopy_layer_tlai, currentCohort%treelai , & + currentCohort%vcmax25top,4) currentCohort%lai = currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area currentCohort%sai = currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 577e25d1b1..5bdd971c97 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -49,7 +49,7 @@ module EDCohortDynamicsMod use FatesAllometryMod , only : bfineroot use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : carea_allom - use FatesAllometryMod , only : StructureResetOfDH + use FatesAllometryMod , only : ForceDBH use FatesAllometryMod , only : tree_lai, tree_sai use PRTGenericMod, only : prt_carbon_allom_hyp @@ -103,7 +103,7 @@ module EDCohortDynamicsMod integer, parameter, private :: conserve_crownarea_and_number_not_dbh = 1 integer, parameter, private :: conserve_dbh_and_number_not_crownarea = 2 - integer, parameter, private :: cohort_fusion_conservation_method = conserve_dbh_and_number_not_crownarea + integer, parameter, private :: cohort_fusion_conservation_method = conserve_crownarea_and_number_not_dbh ! 10/30/09: Created by Rosie Fisher !-------------------------------------------------------------------------------------! @@ -1010,8 +1010,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! ----------------------------------------------------------------- ! if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then - call StructureResetOfDH( currentCohort%prt%GetState(struct_organ,all_carbon_elements), currentCohort%pft, & - currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite ) + call ForceDBH( currentCohort%pft, currentCohort%canopy_trim, & + currentCohort%dbh, currentCohort%hite, & + bdead = currentCohort%prt%GetState(struct_organ,all_carbon_elements)) + end if ! call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft,& diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 1e4a69d0cd..50973ec892 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -60,7 +60,6 @@ module EDPhysiologyMod use FatesAllometryMod , only : bbgw_allom use FatesAllometryMod , only : carea_allom use FatesAllometryMod , only : CheckIntegratedAllometries - use FatesAllometryMod , only : StructureResetOfDH use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : leaf_organ diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 0a6edd9909..2fefc3fd18 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -91,6 +91,7 @@ module FatesAllometryMod use FatesConstantsMod, only : kg_per_Megag use FatesConstantsMod, only : calloc_abs_error use FatesConstantsMod, only : fates_unset_r8 + use FatesConstantsMod, only : itrue use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -116,7 +117,8 @@ module FatesAllometryMod public :: carea_allom ! Generic crown area wrapper public :: bstore_allom ! Generic maximum storage carbon wrapper public :: decay_coeff_kn - public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass + public :: ForceDBH ! Method to set DBH to sync with structure + ! or fineroot biomass public :: CheckIntegratedAllometries public :: CrownDepth public :: set_root_fraction ! Generic wrapper to calculate normalized @@ -731,7 +733,7 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, & real(r8), intent(in) :: vcmax25top ! maximum carboxylation rate at top of crown integer,intent(in) :: call_id ! flag specifying where this is called ! from - + real(r8) :: h real(r8) :: target_bleaf real(r8) :: target_lai @@ -741,8 +743,10 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, & tree_sai = EDPftvarcon_inst%allom_sai_scaler(pft) * target_lai - if( (treelai + tree_sai) > (nlevleaf*dinc_ed) )then + + call h_allom(dbh,pft,h) + write(fates_log(),*) 'The leaf and stem are predicted for a cohort, maxed out the array size' write(fates_log(),*) 'lai: ',treelai write(fates_log(),*) 'sai: ',tree_sai @@ -753,7 +757,8 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, & write(fates_log(),*) 'call id: ',call_id write(fates_log(),*) 'n: ',nplant write(fates_log(),*) 'c_area: ',c_area - write(fates_log(),*) 'dbh: ',dbh + write(fates_log(),*) 'dbh: ',dbh,' dbh_max: ',EDPftvarcon_inst%allom_dbh_maxheight(pft) + write(fates_log(),*) 'h: ',h write(fates_log(),*) 'canopy_trim: ',canopy_trim write(fates_log(),*) 'target_bleaf: ',target_bleaf write(fates_log(),*) 'canopy layer: ',cl @@ -2195,90 +2200,140 @@ end function decay_coeff_kn ! ===================================================================================== - subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h ) + subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bfnrt ) ! ========================================================================= - ! This subroutine estimates the diameter based on the structural biomass - ! using the allometric functions. Since allometry is specified with diameter + ! This subroutine estimates the diameter based on either the structural biomass + ! (if woody) or the fine root biomass (if a grass) using the allometric + ! functions. Since allometry is specified with diameter ! as the independant variable, we must do this through a search algorithm. ! Here, we keep searching until the difference between actual structure and ! the predicted structure based on the searched diameter is within a tolerance. - ! T ! ============================================================================ use FatesConstantsMod , only : calloc_abs_error ! Arguments - real(r8),intent(in) :: bdead ! actual bdead [kgC] + integer(i4),intent(in) :: ipft ! PFT index real(r8),intent(in) :: canopy_trim real(r8),intent(inout) :: d ! plant diameter [cm] real(r8),intent(out) :: h ! plant height + real(r8),intent(in),optional :: bdead ! Structural biomass + real(r8),intent(in),optional :: bfnrt ! Fine root biomass + ! Locals real(r8) :: bt_sap,dbt_sap_dd ! target sap wood at current d real(r8) :: bt_agw,dbt_agw_dd ! target AG wood at current d real(r8) :: bt_bgw,dbt_bgw_dd ! target BG wood at current d real(r8) :: bt_dead,dbt_dead_dd ! target struct wood at current d + real(r8) :: bt_fnrt,dbt_fnrt_dd ! target fineroot at current d real(r8) :: at_sap ! sapwood area (dummy) m2 real(r8) :: dd ! diameter increment for each step real(r8) :: d_try ! trial diameter real(r8) :: bt_dead_try ! trial structure biomasss real(r8) :: dbt_dead_dd_try ! trial structural derivative + real(r8) :: bt_fnrt_try ! trial fineroot biomass + real(r8) :: dbt_fnrt_dd_try ! trial fineroot derivative real(r8) :: step_frac ! step fraction integer :: counter real(r8), parameter :: step_frac0 = 0.9_r8 integer, parameter :: max_counter = 200 + + ! Do reduce "if" calls, we break this call into two parts + if ( EDPftvarcon_inst%woody(ipft) == itrue ) then - call bsap_allom(d,ipft,canopy_trim,at_sap,bt_sap,dbt_sap_dd) - call bagw_allom(d,ipft,bt_agw,dbt_agw_dd) - call bbgw_allom(d,ipft,bt_bgw,dbt_bgw_dd) - call bdead_allom(bt_agw,bt_bgw, bt_sap, ipft, bt_dead, dbt_agw_dd, & - dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd) - - ! This calculates a diameter increment based on the difference - ! in structural mass and the target mass, and sets it to a fraction - ! of the diameter increment - counter = 0 - step_frac = step_frac0 - do while( (bdead-bt_dead) > calloc_abs_error .and. dbt_dead_dd>0.0_r8) - - ! vulnerable to div0 - dd = step_frac*(bdead-bt_dead)/dbt_dead_dd - d_try = d + dd - - call bsap_allom(d_try,ipft,canopy_trim,at_sap,bt_sap,dbt_sap_dd) - call bagw_allom(d_try,ipft,bt_agw,dbt_agw_dd) - call bbgw_allom(d_try,ipft,bt_bgw,dbt_bgw_dd) - call bdead_allom(bt_agw,bt_bgw, bt_sap, ipft, bt_dead_try, dbt_agw_dd, & - dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd_try) - - ! Prevent overshooting - if(bt_dead_try > (bdead+calloc_abs_error)) then - step_frac = step_frac*0.5_r8 - else - step_frac = step_frac0 - d = d_try - bt_dead = bt_dead_try - dbt_dead_dd = dbt_dead_dd_try + if(.not.present(bdead)) then + write(fates_log(),*) 'woody plants must use structure for dbh reset' + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - counter = counter + 1 - if (counter>max_counter) then - write(fates_log(),*) 'Having trouble converging on dbh reset' + + call bsap_allom(d,ipft,canopy_trim,at_sap,bt_sap,dbt_sap_dd) + call bagw_allom(d,ipft,bt_agw,dbt_agw_dd) + call bbgw_allom(d,ipft,bt_bgw,dbt_bgw_dd) + call bdead_allom(bt_agw,bt_bgw, bt_sap, ipft, bt_dead, dbt_agw_dd, & + dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd) + + ! This calculates a diameter increment based on the difference + ! in structural mass and the target mass, and sets it to a fraction + ! of the diameter increment + counter = 0 + step_frac = step_frac0 + do while( (bdead-bt_dead) > calloc_abs_error .and. dbt_dead_dd>0.0_r8) + + ! vulnerable to div0 + dd = step_frac*(bdead-bt_dead)/dbt_dead_dd + d_try = d + dd + + call bsap_allom(d_try,ipft,canopy_trim,at_sap,bt_sap,dbt_sap_dd) + call bagw_allom(d_try,ipft,bt_agw,dbt_agw_dd) + call bbgw_allom(d_try,ipft,bt_bgw,dbt_bgw_dd) + call bdead_allom(bt_agw,bt_bgw, bt_sap, ipft, bt_dead_try, dbt_agw_dd, & + dbt_bgw_dd, dbt_sap_dd, dbt_dead_dd_try) + + ! Prevent overshooting + if(bt_dead_try > (bdead+calloc_abs_error)) then + step_frac = step_frac*0.5_r8 + else + step_frac = step_frac0 + d = d_try + bt_dead = bt_dead_try + dbt_dead_dd = dbt_dead_dd_try + end if + counter = counter + 1 + if (counter>max_counter) then + write(fates_log(),*) 'Having trouble converging on dbh reset' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + + + else + + if(.not.present(bfnrt)) then + write(fates_log(),*) 'grasses must use fine-root for dbh reset' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end do + + call bfineroot(d,ipft,canopy_trim,bt_fnrt,dbt_fnrt_dd) + + counter = 0 + step_frac = step_frac0 + do while( (bfnrt-bt_fnrt) > calloc_abs_error .and. dbt_fnrt_dd>0.0_r8) + + dd = step_frac*(bfnrt-bt_fnrt)/dbt_fnrt_dd + d_try = d + dd + + call bfineroot(d_try,ipft,canopy_trim,bt_fnrt_try,dbt_fnrt_dd_try) + + ! Prevent overshooting + if(bt_fnrt_try > (bfnrt+calloc_abs_error)) then + step_frac = step_frac*0.5_r8 + else + step_frac = step_frac0 + d = d_try + bt_fnrt = bt_fnrt_try + dbt_fnrt_dd = dbt_fnrt_dd_try + end if + counter = counter + 1 + if (counter>max_counter) then + write(fates_log(),*) 'Having trouble converging on dbh reset' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do + + end if call h_allom(d,ipft,h) if(counter>10)then - write(fates_log(),*) 'dbh counter: ',counter + write(fates_log(),*) 'dbh counter: ',counter,EDPftvarcon_inst%woody(ipft) end if - ! At this point, the diameter, height and their target structural biomass - ! should be pretty close to and greater than actual + return - end subroutine StructureResetOfDH + end subroutine ForceDBH ! ========================================================================= diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index dae9699c30..f9ee74af22 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -339,11 +339,30 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n + + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & + currentCohort%vcmax25top) + currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,6 ) + + ! Conducte Maintenance Turnover (parteh) call currentCohort%prt%CheckMassConservation(ft,3) call PRTMaintTurnover(currentCohort%prt,ft,currentSite%is_drought) call currentCohort%prt%CheckMassConservation(ft,4) + leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & + currentCohort%vcmax25top) + currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & + currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,7 ) + + ! Conduct Growth (parteh) call currentCohort%prt%DailyPRT() call currentCohort%prt%CheckMassConservation(ft,5) diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index 43f250183f..cffd717d40 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -34,7 +34,7 @@ module PRTAllometricCarbonMod use FatesAllometryMod , only : bagw_allom use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : CheckIntegratedAllometries - use FatesAllometryMod , only : StructureResetOfDH + use FatesAllometryMod , only : ForceDBH use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log @@ -336,6 +336,7 @@ subroutine DailyPRTAllometricCarbon(this) real(r8) :: repro_c0 ! "" real(r8) :: struct_c0 ! "" + logical :: grow_struct logical :: grow_leaf ! Are leaves at allometric target and should be grown? logical :: grow_fnrt ! Are fine-roots at allometric target and should be grown? logical :: grow_sapw ! Is sapwood at allometric target and should be grown? @@ -441,55 +442,84 @@ subroutine DailyPRTAllometricCarbon(this) ! ----------------------------------------------------------------------------------- ! II. Calculate target size of the biomass compartment for a given dbh. ! ----------------------------------------------------------------------------------- - - ! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm] - call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) - - ! Target total above ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] - call bagw_allom(dbh,ipft,target_agw_c) - ! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] - call bbgw_allom(dbh,ipft,target_bgw_c) + if( EDPftvarcon_inst%woody(ipft) == itrue) then - ! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm] - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + ! Target total above ground biomass in woody/fibrous tissues [kgC] + call bagw_allom(dbh,ipft,target_agw_c) + + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) + + ! Target total dead (structrual) biomass [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! ------------------------------------------------------------------------------------ + ! If structure is larger than target, then we need to correct some integration errors + ! by slightly increasing dbh to match it. + ! For grasses, if leaf biomass is larger than target, then we reset dbh to match + ! ----------------------------------------------------------------------------------- + + if( (struct_c - target_struct_c ) > calloc_abs_error ) then + + call ForceDBH( ipft, canopy_trim, dbh, hite_out, bdead=struct_c ) + + ! Set the structural target biomass to the current structural boimass [kgC] + target_struct_c = struct_c + ! Target sapwood biomass according to allometry and trimming [kgC] + call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) + + end if + + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + + ! Target storage carbon [kgC,kgC/cm] + call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + + else - ! ------------------------------------------------------------------------------------ - ! If structure is larger than target, then we need to correct some integration errors - ! by slightly increasing dbh to match it. - ! For grasses, if leaf biomass is larger than target, then we reset dbh to match - ! ----------------------------------------------------------------------------------- - if( (( struct_c - target_struct_c ) > calloc_abs_error) .and. & - (EDPftvarcon_inst%woody(ipft) == itrue) ) then + ! Target fine-root biomass and deriv. according to + ! allometry and trimming [kgC] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) - call StructureResetOfDH( struct_c, ipft, & - canopy_trim, dbh, hite_out ) + if( (fnrt_c - target_fnrt_c ) > calloc_abs_error ) then + + call ForceDBH( ipft, canopy_trim, dbh, hite_out, bfnrt=fnrt_c ) + + target_fnrt_c = fnrt_c + + end if - ! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + ! Target sapwood biomass according to allometry and trimming [kgC] call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) - ! Target total above ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] + ! Target total above ground biomass in woody/fibrous tissues [kgC] call bagw_allom(dbh,ipft,target_agw_c) - - ! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] - call bbgw_allom(dbh,ipft,target_bgw_c) - - ! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm] - call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) - end if - - ! Target leaf biomass according to allometry and trimming - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) + ! Target total below ground biomass in woody/fibrous tissues [kgC] + call bbgw_allom(dbh,ipft,target_bgw_c) - ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] - call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + ! Target total dead (structrual) biomass and [kgC] + call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) + + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - ! Target storage carbon [kgC,kgC/cm] - call bstore_allom(dbh,ipft,canopy_trim,target_store_c) + ! Target storage carbon [kgC] + call bstore_allom(dbh,ipft,canopy_trim,target_store_c) - + + end if ! ----------------------------------------------------------------------------------- @@ -667,32 +697,24 @@ subroutine DailyPRTAllometricCarbon(this) ! allow actual pools to be above the target, and in these cases, it sends ! a false on the "grow_<>" flag, allowing the plant to grow into these pools. ! It also checks to make sure that structural biomass is not above the target. - if ( EDPftvarcon_inst%woody(ipft) == itrue ) then - - - if( (target_store_c - store_c)>calloc_abs_error) then - write(fates_log(),*) 'storage is not on-allometry at the growth step' - write(fates_log(),*) 'exiting' - write(fates_log(),*) 'cbal: ',carbon_balance - write(fates_log(),*) 'near-zero',nearzero - write(fates_log(),*) 'store_c: ',store_c - write(fates_log(),*) 'target c: ',target_store_c - write(fates_log(),*) 'store_c0:', store_c0 - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, & - store_c, struct_c, & - target_leaf_c, target_fnrt_c, & - target_sapw_c, target_store_c, target_struct_c, & - grow_leaf, grow_fnrt, grow_sapw, grow_store) - else ! for grasses - grow_leaf = .true. - grow_fnrt = .true. - grow_sapw = .true. - grow_store = .true. + if( (target_store_c - store_c)>calloc_abs_error) then + write(fates_log(),*) 'storage is not on-allometry at the growth step' + write(fates_log(),*) 'exiting' + write(fates_log(),*) 'cbal: ',carbon_balance + write(fates_log(),*) 'near-zero',nearzero + write(fates_log(),*) 'store_c: ',store_c + write(fates_log(),*) 'target c: ',target_store_c + write(fates_log(),*) 'store_c0:', store_c0 + call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + + call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, & + store_c, struct_c, & + target_leaf_c, target_fnrt_c, & + target_sapw_c, target_store_c, target_struct_c, & + grow_struct, grow_leaf, grow_fnrt, grow_sapw, grow_store) ! -------------------------------------------------------------------------------- ! The numerical integration of growth requires that the instantaneous state @@ -729,7 +751,7 @@ subroutine DailyPRTAllometricCarbon(this) c_mask(fnrt_c_id) = grow_fnrt c_mask(sapw_c_id) = grow_sapw c_mask(store_c_id) = grow_store - c_mask(struct_c_id) = .true. ! Always increment dead on growth step + c_mask(struct_c_id) = grow_struct c_mask(repro_c_id) = .true. ! Always calculate reproduction on growth c_mask(dbh_id) = .true. ! Always increment dbh on growth step @@ -1045,7 +1067,7 @@ end function AllomCGrowthDeriv subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, & - grow_leaf,grow_froot,grow_sapw,grow_store) + grow_dead,grow_leaf,grow_froot,grow_sapw,grow_store) ! Arguments real(r8),intent(in) :: bleaf !actual @@ -1062,6 +1084,7 @@ subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & logical,intent(out) :: grow_froot logical,intent(out) :: grow_sapw logical,intent(out) :: grow_store + logical,intent(out) :: grow_dead if( (bt_leaf - bleaf)>calloc_abs_error) then write(fates_log(),*) 'leaves are not on-allometry at the growth step' @@ -1108,7 +1131,14 @@ subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, & write(fates_log(),*) 'structure not on-allometry at the growth step' write(fates_log(),*) 'exiting',bdead,bt_dead call endrun(msg=errMsg(sourcefile, __LINE__)) + elseif( (bdead-bt_dead)> calloc_abs_error) then + grow_dead = .false. + else + grow_dead = .true. end if + + + return end subroutine TargetAllometryCheck ! ===================================================================================== From 6b22ae78a9db2f88a6cc3da28a51b77a4de92244 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 7 Mar 2019 17:19:21 -0700 Subject: [PATCH 18/28] Reduced number of canopy layers back to 2 in attempts to avoid canopy radiation errors. Also, made provisions to forcefully terminate patches that have super small areas, even if they want to fuse with the youngest patch. --- biogeochem/EDCanopyStructureMod.F90 | 5 +++-- biogeochem/EDPatchDynamicsMod.F90 | 24 +++++++++++++----------- main/EDTypesMod.F90 | 7 ++++++- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 16108301f9..59b0173ce8 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -60,9 +60,9 @@ module EDCanopyStructureMod ! will attempt to reduce errors ! below this level - real(r8), parameter :: area_check_precision = 1.0E-6_r8 ! Area conservation checks must + real(r8), parameter :: area_check_precision = 1.0E-7_r8 ! Area conservation checks must ! be within this absolute tolerance - real(r8), parameter :: area_check_rel_precision = 1.0E-3_r8 ! Area conservation checks must + real(r8), parameter :: area_check_rel_precision = 1.0E-4_r8 ! Area conservation checks must ! be within this relative tolerance real(r8), parameter :: similar_height_tol = 1.0E-3_r8 ! I think trees that differ by 1mm @@ -796,6 +796,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) write(fates_log(),*) 'patch%area:',currentPatch%area write(fates_log(),*) 'ilayer: ',i_lyr write(fates_log(),*) 'bias:',arealayer - currentPatch%area + write(fates_log(),*) 'rel bias:',(arealayer - currentPatch%area)/arealayer write(fates_log(),*) 'demote_area:',demote_area call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 7803d8c7b6..174d97d159 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -12,6 +12,7 @@ module EDPatchDynamicsMod use EDTypesMod , only : maxPatchesPerSite use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : min_patch_area + use EDTypesMod , only : min_patch_area_forced use EDTypesMod , only : nclmax use EDTypesMod , only : maxpft use EDTypesMod , only : dtype_ifall @@ -1887,7 +1888,15 @@ subroutine terminate_patches(currentSite) if(currentPatch%area <= min_patch_area)then - if ( .not.associated(currentPatch,currentSite%youngest_patch) ) then + ! Even if the patch area is small, avoid fusing it into its neighbor + ! if it is the youngest of all patches. We do this in attempts to maintain + ! a discrete patch for very young patches + ! However, if the patch to be fused is excessivlely small, then fuse + ! at all costs. If it is not fused, it will make + + + if ( .not.associated(currentPatch,currentSite%youngest_patch) .or. & + currentPatch%area <= min_patch_area_forced ) then if(associated(currentPatch%older) )then @@ -1908,7 +1917,7 @@ subroutine terminate_patches(currentSite) ! This logic checks to make sure that the younger patch is not the youngest ! patch. As mentioned earlier, we try not to fuse it. - elseif( .not. associated(currentPatch%younger,currentSite%youngest_patch) ) then + elseif( associated(currentPatch%younger) ) then if(debug) & write(fates_log(),*) 'fusing to younger patch because oldest one is too small', & @@ -1928,16 +1937,9 @@ subroutine terminate_patches(currentSite) enddo !check area is not exceeded - areatot = 0._r8 - currentPatch => currentSite%oldest_patch - do while(associated(currentPatch)) - areatot = areatot + currentPatch%area - currentPatch => currentPatch%younger - if((areatot-area) > 0.0000001_r8)then - write(fates_log(),*) 'ED: areatot too large. end terminate', areatot - endif - enddo + call check_patch_area( currentSite ) + return end subroutine terminate_patches ! ===================================================================================== diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 4fc7ee7d71..834873813e 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -18,7 +18,7 @@ module EDTypesMod integer, parameter :: maxPatchesPerSite = 10 ! maximum number of patches to live on a site integer, parameter :: maxCohortsPerPatch = 100 ! maximum number of cohorts per patch - integer, parameter :: nclmax = 3 ! Maximum number of canopy layers + integer, parameter :: nclmax = 2 ! Maximum number of canopy layers integer, parameter :: ican_upper = 1 ! Nominal index for the upper canopy integer, parameter :: ican_ustory = 2 ! Nominal index for diagnostics that refer ! to understory layers (all layers that @@ -151,6 +151,11 @@ module EDTypesMod real(r8), parameter :: min_npm2 = 1.0E-7_r8 ! minimum cohort number density per m2 before termination real(r8), parameter :: min_patch_area = 0.01_r8 ! smallest allowable patch area before termination + real(r8), parameter :: min_patch_area_forced = 0.0001_r8 ! patch termination will not fuse the youngest patch + ! if the area is less than min_patch_area. + ! however, it is allowed to fuse the youngest patch + ! if the fusion area is less than min_patch_area_forced + real(r8), parameter :: min_nppatch = min_npm2*min_patch_area ! minimum number of cohorts per patch (min_npm2*min_patch_area) real(r8), parameter :: min_n_safemath = 1.0E-12_r8 ! in some cases, we want to immediately remove super small ! number densities of cohorts to prevent FPEs From bc933893864c80d8bdc3e0133cabc0469723237c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 8 Mar 2019 12:50:32 -0700 Subject: [PATCH 19/28] Changed logic to patch termination to ensure no super small patches escape detection and fusion --- biogeochem/EDPatchDynamicsMod.F90 | 36 ++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 174d97d159..a7403a17d8 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1538,7 +1538,7 @@ subroutine fuse_patches( csite, bc_in ) ! Keep doing this until nopatches >= maxPatchesPerSite ! !---------------------------------------------------------------------! - do while(iterate == 1) + do while(iterate == 1 .and. nopatches>1) !---------------------------------------------------------------------! ! Calculate the biomass profile of each patch ! !---------------------------------------------------------------------! @@ -1556,10 +1556,6 @@ subroutine fuse_patches( csite, bc_in ) tpp => currentSite%youngest_patch do while(associated(tpp)) - if(.not.associated(currentPatch))then - write(fates_log(),*) 'ED: issue with currentPatch' - endif - if(associated(tpp).and.associated(currentPatch))then !-------------------------------------------------------------------------------------------- @@ -1893,11 +1889,15 @@ subroutine terminate_patches(currentSite) ! a discrete patch for very young patches ! However, if the patch to be fused is excessivlely small, then fuse ! at all costs. If it is not fused, it will make - + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area: ',currentPatch%area,currentSite%lat,currentSite%lon if ( .not.associated(currentPatch,currentSite%youngest_patch) .or. & currentPatch%area <= min_patch_area_forced ) then + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area(2): ',currentPatch%area,currentSite%lat,currentSite%lon + if(associated(currentPatch%older) )then if(debug) & @@ -1905,12 +1905,18 @@ subroutine terminate_patches(currentSite) currentPatch%area, & currentPatch%older%area + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area(3): ',currentPatch%area,currentSite%lat,currentSite%lon + ! We set a pointer to this patch, because ! it will be returned by the subroutine as de-referenced olderPatch => currentPatch%older call fuse_2_patches(currentSite, olderPatch, currentPatch) + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area(4): ',currentPatch%area,currentSite%lat,currentSite%lon + ! The fusion process has updated the "older" pointer on currentPatch ! for us. @@ -1923,17 +1929,31 @@ subroutine terminate_patches(currentSite) write(fates_log(),*) 'fusing to younger patch because oldest one is too small', & currentPatch%area + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area(5): ',currentPatch%area,currentSite%lat,currentSite%lon + youngerPatch => currentPatch%younger call fuse_2_patches(currentSite, youngerPatch, currentPatch) + if ( currentPatch%area <= min_patch_area_forced ) & + write(fates_log(),*) 'small area(6): ',currentPatch%area,currentSite%lat,currentSite%lon + ! The fusion process has updated the "younger" pointer on currentPatch endif endif endif - - currentPatch => currentPatch%older + ! It is possible that an incredibly small patch just fused into another incredibly + ! small patch, resulting in an incredibly small patch. It is also possible that this + ! resulting incredibly small patch is the oldest patch. If this was true than + ! we would had been at the end of the loop, and left with an incredibly small patch. + ! Think this is impossible? No, this really happens, especially when we have fires. + ! So, we don't move forward until we have merged enough area into this thing. + + if(currentPatch%area > min_patch_area_forced)then + currentPatch => currentPatch%older + end if enddo !check area is not exceeded From 5b41ed61baaabf637c2340ed13a815a2e7baaf0d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 11 Mar 2019 17:07:56 -0700 Subject: [PATCH 20/28] Syntax fix on select case --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0f9faaf5d6..cd88956da7 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1003,7 +1003,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft,& currentCohort%c_area,inverse=.false.) ! - case(default) + case default write(fates_log(),*) 'FATES: Invalid choice for cohort_fusion_conservation_method' call endrun(msg=errMsg(sourcefile, __LINE__)) end select From 668e3de6f4d1c43d9b0a677f4a12d9bc9b0e8f14 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 13 Mar 2019 10:59:19 -0700 Subject: [PATCH 21/28] Added a re-calculation of crown areas prior to crown fusion --- biogeochem/EDCohortDynamicsMod.F90 | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index cd88956da7..58f9689c97 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -844,7 +844,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) logical, parameter :: fuse_debug = .false. ! This debug is over-verbose ! and gets its own flag - + real(r8) :: next_c_area + real(r8) :: curr_c_area !---------------------------------------------------------------------- !set initial fusion tolerance @@ -961,7 +962,16 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! cohorts' dbh ! ----------------------------------------------------------------- ! - currentCohort%c_area = currentCohort%c_area + nextc%c_area + call carea_allom(currentCohort%dbh,currentCohort%n, & + currentSite%spread,currentCohort%pft,& + curr_c_area,inverse=.false.) + + call carea_allom(nextc%dbh,nextc%n, & + currentSite%spread,nextc%pft,& + next_c_area,inverse=.false.) + + !currentCohort%c_area = currentCohort%c_area + nextc%c_area + currentCohort%c_area = curr_c_area + next_c_area ! call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& currentCohort%c_area,inverse=.true.) @@ -969,6 +979,17 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) if (abs(dbh-fates_unset_r8) Date: Wed, 13 Mar 2019 11:07:47 -0700 Subject: [PATCH 22/28] Cleaned up some temporary diagnostics related to crown fusion --- biogeochem/EDCohortDynamicsMod.F90 | 31 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index d589888484..69770aa7ba 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -849,8 +849,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) logical, parameter :: fuse_debug = .false. ! This debug is over-verbose ! and gets its own flag - real(r8) :: next_c_area - real(r8) :: curr_c_area + !---------------------------------------------------------------------- !set initial fusion tolerance @@ -942,15 +941,15 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & + nextc%n*nextc%canopy_trim)/newn - ! c13disc_acc calculation; weighted mean by GPP - if ((currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) .eq. 0.0_r8) then - currentCohort%c13disc_acc = 0.0_r8 - else - currentCohort%c13disc_acc = (currentCohort%n * currentCohort%gpp_acc * currentCohort%c13disc_acc + & - nextc%n * nextc%gpp_acc * nextc%c13disc_acc)/ & - (currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) - endif - + ! c13disc_acc calculation; weighted mean by GPP + if ((currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) .eq. 0.0_r8) then + currentCohort%c13disc_acc = 0.0_r8 + else + currentCohort%c13disc_acc = (currentCohort%n * currentCohort%gpp_acc * currentCohort%c13disc_acc + & + nextc%n * nextc%gpp_acc * nextc%c13disc_acc)/ & + (currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) + endif + select case(cohort_fusion_conservation_method) ! ! ----------------------------------------------------------------- @@ -978,14 +977,14 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,& - curr_c_area,inverse=.false.) - + currentCohort%c_area,inverse=.false.) + call carea_allom(nextc%dbh,nextc%n, & currentSite%spread,nextc%pft,& - next_c_area,inverse=.false.) + nextc%c_area,inverse=.false.) - !currentCohort%c_area = currentCohort%c_area + nextc%c_area - currentCohort%c_area = curr_c_area + next_c_area + currentCohort%c_area = currentCohort%c_area + nextc%c_area + ! call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& currentCohort%c_area,inverse=.true.) From 86a3b28cb47e0be1e93ada3b41fb56549c7fbad9 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 13 Mar 2019 11:19:20 -0700 Subject: [PATCH 23/28] Updated structure reset to use the new caling convention of forcedbh --- biogeochem/EDCohortDynamicsMod.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 69770aa7ba..50c046ba33 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -994,10 +994,11 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) + nextc%n*nextc%dbh)/newn if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then - call StructureResetOfDH( & - currentCohort%prt%GetState(struct_organ,all_carbon_elements), & - currentCohort%pft, & - currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite ) + + call ForceDBH( currentCohort%pft, currentCohort%canopy_trim, & + currentCohort%dbh, currentCohort%hite, & + bdead = currentCohort%prt%GetState(struct_organ,all_carbon_elements)) + end if ! call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft,& From d6c81e044203102d08872ec77f04db6b19465858 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 15 Mar 2019 13:58:37 -0600 Subject: [PATCH 24/28] Fixed a diagnostic in LAI during fusion and changed forced dbh resets to use leaves instead of fineroots for grasses --- biogeochem/EDCohortDynamicsMod.F90 | 38 ++++++++++++++++++++++-------- biogeochem/FatesAllometryMod.F90 | 30 +++++++++++------------ parteh/PRTAllometricCarbonMod.F90 | 18 +++++++------- 3 files changed, 52 insertions(+), 34 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 50c046ba33..a92deeb5b6 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -838,8 +838,9 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) integer :: nocohorts real(r8) :: newn real(r8) :: diff - real(r8) :: leaf_c_next ! Leaf carbon * plant density of current (for weighting) - real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting) + real(r8) :: leaf_c_next ! Leaf carbon * plant density of current (for weighting) + real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting) + real(r8) :: leaf_c_target real(r8) :: dynamic_fusion_tolerance real(r8) :: dbh real(r8) :: leaf_c ! leaf carbon [kg] @@ -935,6 +936,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, & currentCohort%n/newn ) + ! Leaf biophysical rates (use leaf mass weighting) + ! ----------------------------------------------------------------- + call UpdateCohortBioPhysRates(currentCohort) + currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & + nextc%n*nextc%laimemory)/newn @@ -996,8 +1001,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) if( EDPftvarcon_inst%woody(currentCohort%pft) == itrue ) then call ForceDBH( currentCohort%pft, currentCohort%canopy_trim, & - currentCohort%dbh, currentCohort%hite, & - bdead = currentCohort%prt%GetState(struct_organ,all_carbon_elements)) + currentCohort%dbh, currentCohort%hite, & + bdead = currentCohort%prt%GetState(struct_organ,all_carbon_elements)) end if ! @@ -1007,6 +1012,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) else currentCohort%dbh = dbh endif + ! call h_allom(currentCohort%dbh,currentCohort%pft,currentCohort%hite) ! @@ -1045,14 +1051,29 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call endrun(msg=errMsg(sourcefile, __LINE__)) end select + + ! If fusion forces the actual leaf biomass to be unreasonably + ! greater than the target (ie 25%), reset the DBH + +! call bleaf(currentCohort%dbh,currentCohort%pft, & +! currentCohort%canopy_trim,leaf_c_target) - leaf_c = currentCohort%prt%GetState(leaf_organ, all_carbon_elements) + leaf_c = currentCohort%prt%GetState(leaf_organ,all_carbon_elements) + +! if (leaf_c > leaf_c_target*1.25_r8) then +! call ForceDBH( currentCohort%pft, currentCohort%canopy_trim, & +! currentCohort%dbh, currentCohort%hite, & +! bl = leaf_c) +! call carea_allom(currentCohort%dbh,newn,currentSite%spread,currentCohort%pft, & +! currentCohort%c_area,inverse=.false.) +! end if + - currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%treelai = tree_lai(leaf_c, currentCohort%pft, currentCohort%c_area, newn, & currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & currentCohort%vcmax25top) currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, & - currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, & + currentCohort%c_area, newn, currentCohort%canopy_layer, & currentPatch%canopy_layer_tlai, currentCohort%treelai,currentCohort%vcmax25top,1 ) @@ -1069,9 +1090,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%canopy_layer_yesterday = (currentCohort%n*currentCohort%canopy_layer_yesterday + & nextc%n*nextc%canopy_layer_yesterday)/newn - ! Leaf biophysical rates (use leaf mass weighting) - ! ----------------------------------------------------------------- - call UpdateCohortBioPhysRates(currentCohort) ! keep track of the size class bins so that we can monitor growth fluxes ! compare the values. if they are the same, then nothing needs to be done. if not, track the diagnostic flux diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 2fefc3fd18..e7470b076c 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -2200,11 +2200,11 @@ end function decay_coeff_kn ! ===================================================================================== - subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bfnrt ) + subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bl ) ! ========================================================================= ! This subroutine estimates the diameter based on either the structural biomass - ! (if woody) or the fine root biomass (if a grass) using the allometric + ! (if woody) or the leaf biomass using the allometric ! functions. Since allometry is specified with diameter ! as the independant variable, we must do this through a search algorithm. ! Here, we keep searching until the difference between actual structure and @@ -2220,7 +2220,7 @@ subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bfnrt ) real(r8),intent(inout) :: d ! plant diameter [cm] real(r8),intent(out) :: h ! plant height real(r8),intent(in),optional :: bdead ! Structural biomass - real(r8),intent(in),optional :: bfnrt ! Fine root biomass + real(r8),intent(in),optional :: bl ! Leaf biomass ! Locals @@ -2228,14 +2228,14 @@ subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bfnrt ) real(r8) :: bt_agw,dbt_agw_dd ! target AG wood at current d real(r8) :: bt_bgw,dbt_bgw_dd ! target BG wood at current d real(r8) :: bt_dead,dbt_dead_dd ! target struct wood at current d - real(r8) :: bt_fnrt,dbt_fnrt_dd ! target fineroot at current d + real(r8) :: bt_leaf,dbt_leaf_dd ! target leaf at current d real(r8) :: at_sap ! sapwood area (dummy) m2 real(r8) :: dd ! diameter increment for each step real(r8) :: d_try ! trial diameter real(r8) :: bt_dead_try ! trial structure biomasss real(r8) :: dbt_dead_dd_try ! trial structural derivative - real(r8) :: bt_fnrt_try ! trial fineroot biomass - real(r8) :: dbt_fnrt_dd_try ! trial fineroot derivative + real(r8) :: bt_leaf_try ! trial leaf biomass + real(r8) :: dbt_leaf_dd_try ! trial leaf derivative real(r8) :: step_frac ! step fraction integer :: counter real(r8), parameter :: step_frac0 = 0.9_r8 @@ -2291,30 +2291,30 @@ subroutine ForceDBH( ipft, canopy_trim, d, h, bdead, bfnrt ) else - if(.not.present(bfnrt)) then - write(fates_log(),*) 'grasses must use fine-root for dbh reset' + if(.not.present(bl)) then + write(fates_log(),*) 'grasses must use leaf for dbh reset' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - call bfineroot(d,ipft,canopy_trim,bt_fnrt,dbt_fnrt_dd) + call bleaf(d,ipft,canopy_trim,bt_leaf,dbt_leaf_dd) counter = 0 step_frac = step_frac0 - do while( (bfnrt-bt_fnrt) > calloc_abs_error .and. dbt_fnrt_dd>0.0_r8) + do while( (bl-bt_leaf) > calloc_abs_error .and. dbt_leaf_dd>0.0_r8) - dd = step_frac*(bfnrt-bt_fnrt)/dbt_fnrt_dd + dd = step_frac*(bl-bt_leaf)/dbt_leaf_dd d_try = d + dd - call bfineroot(d_try,ipft,canopy_trim,bt_fnrt_try,dbt_fnrt_dd_try) + call bleaf(d_try,ipft,canopy_trim,bt_leaf_try,dbt_leaf_dd_try) ! Prevent overshooting - if(bt_fnrt_try > (bfnrt+calloc_abs_error)) then + if(bt_leaf_try > (bl+calloc_abs_error)) then step_frac = step_frac*0.5_r8 else step_frac = step_frac0 d = d_try - bt_fnrt = bt_fnrt_try - dbt_fnrt_dd = dbt_fnrt_dd_try + bt_leaf = bt_leaf_try + dbt_leaf_dd = dbt_leaf_dd_try end if counter = counter + 1 if (counter>max_counter) then diff --git a/parteh/PRTAllometricCarbonMod.F90 b/parteh/PRTAllometricCarbonMod.F90 index cffd717d40..3a443c9050 100644 --- a/parteh/PRTAllometricCarbonMod.F90 +++ b/parteh/PRTAllometricCarbonMod.F90 @@ -488,18 +488,21 @@ subroutine DailyPRTAllometricCarbon(this) else - ! Target fine-root biomass and deriv. according to - ! allometry and trimming [kgC] - call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + ! Target leaf biomass according to allometry and trimming + call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - if( (fnrt_c - target_fnrt_c ) > calloc_abs_error ) then + + if( ( sum(leaf_c) - target_leaf_c ) > calloc_abs_error ) then - call ForceDBH( ipft, canopy_trim, dbh, hite_out, bfnrt=fnrt_c ) + call ForceDBH( ipft, canopy_trim, dbh, hite_out, bl=sum(leaf_c) ) - target_fnrt_c = fnrt_c + target_leaf_c = sum(leaf_c) end if + ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm] + call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c) + ! Target sapwood biomass according to allometry and trimming [kgC] call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c) @@ -512,9 +515,6 @@ subroutine DailyPRTAllometricCarbon(this) ! Target total dead (structrual) biomass and [kgC] call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c) - ! Target leaf biomass according to allometry and trimming - call bleaf(dbh,ipft,canopy_trim,target_leaf_c) - ! Target storage carbon [kgC] call bstore_allom(dbh,ipft,canopy_trim,target_store_c) From 2bfdb1d99e864ec52d1aec3c680638b2deefea8e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 26 Mar 2019 17:58:15 -0600 Subject: [PATCH 25/28] Added a check on the number of loops allowed to catch super-small patches --- biogeochem/EDPatchDynamicsMod.F90 | 43 +++++++++++++++++++------------ 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index a7403a17d8..9834e9da72 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1875,10 +1875,15 @@ subroutine terminate_patches(currentSite) type(ed_patch_type), pointer :: currentPatch type(ed_patch_type), pointer :: olderPatch type(ed_patch_type), pointer :: youngerPatch + integer, parameter :: max_cycles = 10 ! After 10 loops through + ! You should had fused + integer :: count_cycles real(r8) areatot ! variable for checking whether the total patch area is wrong. !--------------------------------------------------------------------- + count_cycles = 0 + currentPatch => currentSite%youngest_patch do while(associated(currentPatch)) @@ -1889,15 +1894,10 @@ subroutine terminate_patches(currentSite) ! a discrete patch for very young patches ! However, if the patch to be fused is excessivlely small, then fuse ! at all costs. If it is not fused, it will make - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area: ',currentPatch%area,currentSite%lat,currentSite%lon if ( .not.associated(currentPatch,currentSite%youngest_patch) .or. & currentPatch%area <= min_patch_area_forced ) then - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area(2): ',currentPatch%area,currentSite%lat,currentSite%lon - if(associated(currentPatch%older) )then if(debug) & @@ -1905,18 +1905,12 @@ subroutine terminate_patches(currentSite) currentPatch%area, & currentPatch%older%area - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area(3): ',currentPatch%area,currentSite%lat,currentSite%lon - ! We set a pointer to this patch, because ! it will be returned by the subroutine as de-referenced olderPatch => currentPatch%older call fuse_2_patches(currentSite, olderPatch, currentPatch) - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area(4): ',currentPatch%area,currentSite%lat,currentSite%lon - ! The fusion process has updated the "older" pointer on currentPatch ! for us. @@ -1928,16 +1922,10 @@ subroutine terminate_patches(currentSite) if(debug) & write(fates_log(),*) 'fusing to younger patch because oldest one is too small', & currentPatch%area - - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area(5): ',currentPatch%area,currentSite%lat,currentSite%lon youngerPatch => currentPatch%younger call fuse_2_patches(currentSite, youngerPatch, currentPatch) - if ( currentPatch%area <= min_patch_area_forced ) & - write(fates_log(),*) 'small area(6): ',currentPatch%area,currentSite%lat,currentSite%lon - ! The fusion process has updated the "younger" pointer on currentPatch endif @@ -1953,7 +1941,28 @@ subroutine terminate_patches(currentSite) if(currentPatch%area > min_patch_area_forced)then currentPatch => currentPatch%older + count_cycles = 0 + else + count_cycles = count_cycles + 1 end if + + if(count_cycles > max_cycles) then + write(fates_log(),*) 'FATES is having difficulties fusing very small patches.' + write(fates_log(),*) 'It is possible that a either a secondary or primary' + write(fates_log(),*) 'patch has become the only patch of its kind, and it is' + write(fates_log(),*) 'is very very small. You can test your luck by' + write(fates_log(),*) 'disabling the endrun statement following this message.' + write(fates_log(),*) 'FATES may or may not continue to operate within error' + write(fates_log(),*) 'tolerances, but will generate another fail if it does not.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + + ! Note to user. If you DO decide to remove the end-run above this line + ! Make sure that you keep the pointer below this line, or you will get + ! an infinite loop. + currentPatch => currentPatch%older + count_cycles = 0 + end if + enddo !check area is not exceeded From 19bcfb3f0906880b556681461e7a275410e27dca Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 27 Mar 2019 16:44:16 -0600 Subject: [PATCH 26/28] Removed unnecessary error checks, and bound others behind logical constants --- biogeochem/EDCanopyStructureMod.F90 | 131 ++++++++++++---------------- 1 file changed, 57 insertions(+), 74 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 59b0173ce8..02eb8ab563 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -386,12 +386,13 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call carea_allom(currentCohort%dbh,currentCohort%n, & currentSite%spread,currentCohort%pft,currentCohort%c_area) - if(currentCohort%c_area<0._r8)then - write(fates_log(),*) 'negative c_area stage 1d: ',currentCohort%dbh,i_lyr,currentCohort%n, & - currentSite%spread,currentCohort%pft,currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(debug) then + if(currentCohort%c_area<0._r8)then + write(fates_log(),*) 'negative c_area stage 1d: ',currentCohort%dbh,i_lyr,currentCohort%n, & + currentSite%spread,currentCohort%pft,currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if - if( currentCohort%canopy_layer == i_lyr)then @@ -525,17 +526,19 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) if(currentCohort%canopy_layer == i_lyr) then currentCohort%excl_weight = currentCohort%c_area * currentCohort%excl_weight * scale_factor - if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%excl_weight < 0._r8) ) then - write(fates_log(),*) 'exclusion area too big (1)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'dbh: ',currentCohort%dbh - write(fates_log(),*) 'n: ',currentCohort%n - write(fates_log(),*) 'spread: ',currentSite%spread - write(fates_log(),*) 'pft: ',currentCohort%pft - write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight - write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(debug) then + if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%excl_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area too big (1)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'dbh: ',currentCohort%dbh + write(fates_log(),*) 'n: ',currentCohort%n + write(fates_log(),*) 'spread: ',currentSite%spread + write(fates_log(),*) 'pft: ',currentCohort%pft + write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight + write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if endif @@ -572,14 +575,16 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) currentCohort%excl_weight = currentCohort%c_area * & (currentCohort%excl_weight * scale_factor_min + & (1._r8 - (currentCohort%excl_weight*scale_factor_min) ) * scale_factor_res) - - if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%excl_weight < 0._r8) ) then - write(fates_log(),*) 'exclusion area error (2)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight - write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + + if(debug)then + if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%excl_weight < 0._r8) ) then + write(fates_log(),*) 'exclusion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight + write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if endif @@ -658,21 +663,10 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) endif call copy_cohort(currentCohort, copyc) - if(currentCohort%n < 0._r8) then - write(fates_log(),*) 'negatives (0_?): ',currentCohort%n,currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - newarea = currentCohort%c_area - cc_loss copyc%n = currentCohort%n*newarea/currentCohort%c_area currentCohort%n = currentCohort%n - copyc%n - if(copyc%n < 0._r8 .or. currentCohort%n < 0._r8) then - write(fates_log(),*) 'negatives?: ',newarea,cc_loss,copyc%n,currentCohort%n,currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - copyc%canopy_layer = i_lyr !the taller cohort is the copy ! Demote the current cohort to the understory. @@ -688,16 +682,6 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & currentCohort%pft,currentCohort%c_area) - ! Calculate how much area loss from recalculation - if( abs(copyc%c_area-newarea)>area_target_precision .or. & - abs(currentCohort%c_area-cc_loss)>area_target_precision) then - write(fates_log(),*) 'recalculation losses?',newarea-copyc%c_area,newarea, & - copyc%c_area,currentCohort%c_area-cc_loss,currentCohort%c_area,& - cc_loss - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - !----------- Insert copy into linked list ------------------------! copyc%shorter => currentCohort if(associated(currentCohort%taller))then @@ -1034,14 +1018,16 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1) ) then currentCohort%prom_weight = currentCohort%c_area * currentCohort%prom_weight * scale_factor - - if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'promotion area too big (1)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight - write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + + if(debug)then + if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'promotion area too big (1)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight + write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if endif @@ -1078,15 +1064,17 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) (currentCohort%prom_weight * scale_factor_min + & (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * scale_factor_res) - if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'promotion area error (2)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight - write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(debug)then + if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'promotion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight + write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if - + endif currentCohort => currentCohort%shorter enddo @@ -1106,12 +1094,14 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) currentCohort => currentCohort%shorter end do - if (abs(sumweights - promote_area) > area_check_precision ) then - write(fates_log(),*) 'promotions dont add up' - write(fates_log(),*) 'sum promotions: ',sumweights - write(fates_log(),*) 'area needed to be promoted: ',promote_area - write(fates_log(),*) 'excess: ',sumweights - promote_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if(debug)then + if (abs(sumweights - promote_area) > area_check_precision ) then + write(fates_log(),*) 'promotions dont add up' + write(fates_log(),*) 'sum promotions: ',sumweights + write(fates_log(),*) 'area needed to be promoted: ',promote_area + write(fates_log(),*) 'excess: ',sumweights - promote_area + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if end if currentCohort => currentPatch%tallest @@ -1162,13 +1152,6 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) ! number of individuals in cohort remaining in understorey currentCohort%n = currentCohort%n - copyc%n - - if(copyc%n < 0._r8 .or. currentCohort%n < 0._r8) then - write(fates_log(),*) 'negatives (2)?: ',newarea,cc_gain,copyc%n,currentCohort%n - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - currentCohort%canopy_layer = i_lyr + 1 ! keep current cohort in the understory. copyc%canopy_layer = i_lyr ! promote copy to the higher canopy layer. From a1eff334855d017ce0290b51a7000ee6adcb20ec Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 27 Mar 2019 15:56:07 -0700 Subject: [PATCH 27/28] Reduced some line lengths --- biogeochem/EDCanopyStructureMod.F90 | 96 ++++++++++++++++------------- 1 file changed, 54 insertions(+), 42 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 02eb8ab563..c027cfdc16 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -167,10 +167,6 @@ subroutine canopy_structure( currentSite , bc_in ) ! Perform numerical checks on some cohort and patch structures ! ------------------------------------------------------------------------------ -! call val_check_ed_vars(currentPatch,'co_n:co_dbh:pa_area',return_code) -! ! No need to make error message, already generated in math_check_ed_vars -! if(return_code>0) call endrun(msg=errMsg(sourcefile, __LINE__)) - ! canopy layer has a special bounds check currentCohort => currentPatch%tallest do while (associated(currentCohort)) @@ -253,7 +249,7 @@ subroutine canopy_structure( currentSite , bc_in ) area_not_balanced = .false. do i_lyr = 1,z call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer(i_lyr)) - if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. & + if( ((arealayer(i_lyr)-currentPatch%area)/currentPatch%area > area_check_rel_precision) .or. & ((arealayer(i_lyr)-currentPatch%area) > area_check_precision ) ) then area_not_balanced = .true. endif @@ -557,9 +553,11 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) do while (associated(currentCohort)) if(currentCohort%canopy_layer == i_lyr) then area_res = area_res + & - currentCohort%c_area*currentCohort%excl_weight*scale_factor_min + currentCohort%c_area * currentCohort%excl_weight * & + scale_factor_min scale_factor_res = scale_factor_res + & - currentCohort%c_area * (1._r8 - (currentCohort%excl_weight * scale_factor_min)) + currentCohort%c_area * & + (1._r8 - (currentCohort%excl_weight * scale_factor_min)) endif currentCohort => currentCohort%shorter enddo @@ -577,15 +575,18 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr) (1._r8 - (currentCohort%excl_weight*scale_factor_min) ) * scale_factor_res) if(debug)then - if((currentCohort%excl_weight > (currentCohort%c_area+area_target_precision)) .or. & + if((currentCohort%excl_weight > & + (currentCohort%c_area+area_target_precision)) .or. & (currentCohort%excl_weight < 0._r8) ) then - write(fates_log(),*) 'exclusion area error (2)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%excl_weight: ',currentCohort%excl_weight - write(fates_log(),*) 'excess: ',currentCohort%excl_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) 'exclusion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%excl_weight: ', & + currentCohort%excl_weight + write(fates_log(),*) 'excess: ', & + currentCohort%excl_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end if endif currentCohort => currentCohort%shorter @@ -1017,18 +1018,22 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) currentCohort => currentPatch%tallest do while (associated(currentCohort)) if(currentCohort%canopy_layer == (i_lyr+1) ) then - currentCohort%prom_weight = currentCohort%c_area * currentCohort%prom_weight * scale_factor + currentCohort%prom_weight = currentCohort%c_area * & + currentCohort%prom_weight * scale_factor if(debug)then - if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'promotion area too big (1)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight - write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if((currentCohort%prom_weight > & + (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'promotion area too big (1)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ', & + currentCohort%prom_weight + write(fates_log(),*) 'excess: ', & + currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if - end if + end if endif currentCohort => currentCohort%shorter @@ -1047,7 +1052,8 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) area_res = area_res + & currentCohort%c_area*currentCohort%prom_weight*scale_factor_min scale_factor_res = scale_factor_res + & - currentCohort%c_area * (1._r8 - (currentCohort%prom_weight * scale_factor_min)) + currentCohort%c_area * & + (1._r8 - (currentCohort%prom_weight * scale_factor_min)) endif currentCohort => currentCohort%shorter enddo @@ -1062,16 +1068,20 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) currentCohort%prom_weight = currentCohort%c_area * & (currentCohort%prom_weight * scale_factor_min + & - (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * scale_factor_res) + (1._r8 - (currentCohort%prom_weight*scale_factor_min) ) * & + scale_factor_res) if(debug)then - if((currentCohort%prom_weight > (currentCohort%c_area+area_target_precision)) .or. & - (currentCohort%prom_weight < 0._r8) ) then - write(fates_log(),*) 'promotion area error (2)' - write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area - write(fates_log(),*) 'currentCohort%prom_weight: ',currentCohort%prom_weight - write(fates_log(),*) 'excess: ',currentCohort%prom_weight - currentCohort%c_area - call endrun(msg=errMsg(sourcefile, __LINE__)) + if((currentCohort%prom_weight > & + (currentCohort%c_area+area_target_precision)) .or. & + (currentCohort%prom_weight < 0._r8) ) then + write(fates_log(),*) 'promotion area error (2)' + write(fates_log(),*) 'currentCohort%c_area: ',currentCohort%c_area + write(fates_log(),*) 'currentCohort%prom_weight: ', & + currentCohort%prom_weight + write(fates_log(),*) 'excess: ', & + currentCohort%prom_weight - currentCohort%c_area + call endrun(msg=errMsg(sourcefile, __LINE__)) end if end if @@ -1192,7 +1202,8 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) call CanopyLayerArea(currentPatch,currentSite%spread,i_lyr,arealayer_current) - if ((abs(arealayer_current - currentPatch%area)/arealayer_current > area_check_rel_precision ) .or. & + if ((abs(arealayer_current - currentPatch%area)/arealayer_current > & + area_check_rel_precision ) .or. & (abs(arealayer_current - currentPatch%area) > area_check_precision) ) then write(fates_log(),*) 'promotion did not bring area within tolerance' write(fates_log(),*) 'arealayer:',arealayer_current @@ -1252,7 +1263,8 @@ subroutine canopy_spread( currentSite ) enddo !currentPatch - !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner + ! If the canopy area is approaching closure, + ! squash the tree canopies and make them taller and thinner if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then currentSite%spread = currentSite%spread - inc else @@ -1290,15 +1302,15 @@ subroutine canopy_summarization( nsites, sites, bc_in ) type (ed_patch_type) , pointer :: currentPatch type (ed_cohort_type) , pointer :: currentCohort integer :: s - integer :: ft ! plant functional type + integer :: ft ! plant functional type integer :: ifp - integer :: patchn ! identification number for each patch. - real(r8) :: canopy_leaf_area ! total amount of leaf area in the vegetated area. m2. - real(r8) :: leaf_c ! leaf carbon [kg] - real(r8) :: fnrt_c ! fineroot carbon [kg] - real(r8) :: sapw_c ! sapwood carbon [kg] - real(r8) :: store_c ! storage carbon [kg] - real(r8) :: struct_c ! structure carbon [kg] + integer :: patchn ! identification number for each patch. + real(r8) :: canopy_leaf_area ! total amount of leaf area in the vegetated area. m2. + real(r8) :: leaf_c ! leaf carbon [kg] + real(r8) :: fnrt_c ! fineroot carbon [kg] + real(r8) :: sapw_c ! sapwood carbon [kg] + real(r8) :: store_c ! storage carbon [kg] + real(r8) :: struct_c ! structure carbon [kg] !---------------------------------------------------------------------- if ( debug ) then From f6734abc87077334a0c912e844905e11b24b4d2f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 5 Apr 2019 13:29:21 -0700 Subject: [PATCH 28/28] Fixed comment typo --- biogeochem/EDCanopyStructureMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index c027cfdc16..8dd5e139c1 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1041,7 +1041,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr) else - ! Non-trivial case, at least 1 cohort's demotion + ! Non-trivial case, at least 1 cohort's promotion ! rate would exceed its area, given the trivial scale factor area_res = 0._r8