From 7cdca1363307e3657c893d67623ba8b37535cb6e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 13 Dec 2018 12:30:43 -0800 Subject: [PATCH] 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