diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index a52be9a30e..7845bf72c8 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -19,6 +19,7 @@ module EDCanopyStructureMod use EDCohortDynamicsMod , only : InitPRTObject use FatesAllometryMod , only : tree_lai use FatesAllometryMod , only : tree_sai + use FatesAllometryMod , only : CrownDepth use EDtypesMod , only : ed_site_type use FatesPatchMod, only : fates_patch_type use FatesCohortMod, only : fates_cohort_type @@ -1527,6 +1528,7 @@ subroutine leaf_area_profile( currentSite ) real(r8) :: max_cheight ! top of cohort canopy (m) real(r8) :: lai ! leaf area per canopy area real(r8) :: sai ! stem area per canopy area + real(r8) :: crown_depth ! Current cohort's crown depth !---------------------------------------------------------------------- @@ -1603,6 +1605,13 @@ subroutine leaf_area_profile( currentSite ) end if + !---~--- + ! Find current crown depth using the allometric function. + !---~--- + call CrownDepth(currentCohort%height,currentCohort%pft,crown_depth) + !---~--- + + ! -------------------------------------------------------------------------- ! Whole layers. Make a weighted average of the leaf area in each layer ! before dividing it by the total area. Fill up layer for whole layers. @@ -1616,12 +1625,10 @@ subroutine leaf_area_profile( currentSite ) ! is obscured by snow. layer_top_height = currentCohort%height - & - ( real(iv-1,r8)/currentCohort%NV * currentCohort%height * & - prt_params%crown_depth_frac(currentCohort%pft) ) + ( real(iv-1,r8)/currentCohort%NV * crown_depth ) layer_bottom_height = currentCohort%height - & - ( real(iv,r8)/currentCohort%NV * currentCohort%height * & - prt_params%crown_depth_frac(currentCohort%pft) ) + ( real(iv,r8)/currentCohort%NV * crown_depth ) fraction_exposed = 1.0_r8 if(currentSite%snow_depth > layer_top_height)then diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90 index 431219cfda..45b8217518 100644 --- a/biogeochem/FatesAllometryMod.F90 +++ b/biogeochem/FatesAllometryMod.F90 @@ -92,6 +92,7 @@ module FatesAllometryMod use FatesConstantsMod, only : calloc_abs_error use FatesConstantsMod, only : fates_unset_r8 use FatesConstantsMod, only : itrue + use FatesConstantsMod, only : nearzero use shr_log_mod , only : errMsg => shr_log_errMsg use FatesGlobals , only : fates_log use FatesGlobals , only : endrun => fates_endrun @@ -305,7 +306,7 @@ subroutine h2d_allom(h,ipft,d,dddh) p3 => prt_params%allom_d2h3(ipft), & allom_hmode => prt_params%allom_hmode(ipft)) - select case(int(allom_hmode)) + select case(allom_hmode) case (1) ! O'Brien et al 1995, BCI call h2d_obrien(h,p1,p2,d,dddh) case (2) ! poorter 2006 @@ -344,7 +345,7 @@ subroutine h_allom(d,ipft,h,dhdd) p3 => prt_params%allom_d2h3(ipft), & allom_hmode => prt_params%allom_hmode(ipft)) - select case(int(allom_hmode)) + select case(allom_hmode) case (1) ! "obrien" call d2h_obrien(d,p1,p2,dbh_maxh,h,dhdd) case (2) ! "poorter06" @@ -397,7 +398,7 @@ subroutine bagw_allom(d,ipft,crowndamage, elongf_stem, bagw,dbagwdd) branch_frac = param_derived%branch_frac(ipft) - select case(int(allom_amode)) + select case(allom_amode) case (1) !"salda") call h_allom(d,ipft,h,dhdd) call dh2bagw_salda(d,h,dhdd,p1,p2,p3,p4,wood_density,c2b,agb_frac,bagw,dbagwdd) @@ -407,6 +408,9 @@ subroutine bagw_allom(d,ipft,crowndamage, elongf_stem, bagw,dbagwdd) case (3) !"chave14") call h_allom(d,ipft,h,dhdd) call dh2bagw_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bagw,dbagwdd) + case (4) ! 3par_pwr + call h_allom(d,ipft,h,dhdd) + call dh2bagw_3pwr(d,h,dhdd,p1,p2,p3,wood_density,c2b,bagw,dbagwdd) case DEFAULT write(fates_log(),*) 'An undefined AGB allometry was specified: ',allom_amode write(fates_log(),*) 'Aborting' @@ -444,21 +448,28 @@ subroutine blmax_allom(d,ipft,blmax,dblmaxdd) real(r8),intent(out) :: blmax ! plant leaf biomass [kg] real(r8),intent(out),optional :: dblmaxdd ! change leaf bio per diameter [kgC/cm] + real(r8) :: h ! height + real(r8) :: dhdd ! change in height wrt d + associate( dbh_maxh => prt_params%allom_dbh_maxheight(ipft), & rho => prt_params%wood_density(ipft), & + slatop => prt_params%slatop(ipft), & c2b => prt_params%c2b(ipft), & allom_lmode => prt_params%allom_lmode(ipft), & p1 => prt_params%allom_d2bl1(ipft), & p2 => prt_params%allom_d2bl2(ipft), & p3 => prt_params%allom_d2bl3(ipft)) - select case(int(allom_lmode)) + select case(allom_lmode) case(1) !"salda") call d2blmax_salda(d,p1,p2,p3,rho,dbh_maxh,c2b,blmax,dblmaxdd) case(2) !"2par_pwr") call d2blmax_2pwr(d,p1,p2,c2b,blmax,dblmaxdd) case(3) ! dh2blmax_2pwr call dh2blmax_2pwr(d,p1,p2,dbh_maxh,c2b,blmax,dblmaxdd) + case(4) ! dh2blmax_3pwr + call h_allom(d,ipft,h,dhdd) + call dh2blmax_3pwr(d,h,dhdd,p1,p2,p3,slatop,dbh_maxh,c2b,blmax,dblmaxdd) case DEFAULT write(fates_log(),*) 'An undefined leaf allometry was specified: ', & allom_lmode @@ -485,6 +496,7 @@ subroutine carea_allom(dbh,nplant,site_spread,ipft,crowndamage,c_area,inverse) ! instead of crown area from dbh real(r8) :: dbh_eff ! Effective diameter (cm) + real(r8) :: height ! height logical :: do_inverse ! local copy of the inverse argument ! defaults to false logical :: capped_allom ! if we are using an allometry that caps @@ -507,7 +519,7 @@ subroutine carea_allom(dbh,nplant,site_spread,ipft,crowndamage,c_area,inverse) endif endif - select case(int(allom_lmode)) + select case(allom_lmode) case(1) dbh_eff = min(dbh,dbh_maxh) call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max, & @@ -522,6 +534,12 @@ subroutine carea_allom(dbh,nplant,site_spread,ipft,crowndamage,c_area,inverse) call carea_2pwr(dbh_eff,site_spread,d2bl_p2,d2bl_ediff,d2ca_min,d2ca_max, & crowndamage, c_area, do_inverse) capped_allom = .true. + case (4) + dbh_eff = min(dbh,dbh_maxh) + call h_allom(dbh,ipft,height) + call carea_3pwr(dbh_eff,height,ipft,dbh_maxh, site_spread,d2bl_p2, & + d2bl_ediff, d2ca_min,d2ca_max,crowndamage, c_area, do_inverse) + capped_allom = .true. case DEFAULT write(fates_log(),*) 'An undefined leaf allometry was specified: ', & allom_lmode @@ -959,7 +977,7 @@ subroutine bsap_allom(d,ipft,crowndamage,canopy_trim,elongf_stem, sapw_area,bsap branch_frac = param_derived%branch_frac(ipft) - select case(int(prt_params%allom_smode(ipft))) + select case(prt_params%allom_smode(ipft)) ! --------------------------------------------------------------------- ! Currently only one sapwood allometry model. the slope ! of the la:sa to diameter line is zero. @@ -1033,7 +1051,7 @@ subroutine bbgw_allom(d,ipft,elongf_stem,bbgw,dbbgwdd) real(r8) :: bagw ! above ground biomass [kgC] real(r8) :: dbagwdd ! change in agb per diameter [kgC/cm] - select case(int(prt_params%allom_cmode(ipft))) + select case(prt_params%allom_cmode(ipft)) case(1) !"constant") ! bbgw not affected by damage so use target allometry no damage. But note that bbgw ! is affected by stem phenology (typically applied only to grasses). We do not need @@ -1078,7 +1096,7 @@ subroutine bfineroot(d,ipft,canopy_trim,l2fr,elongf_fnrt,bfr,dbfrdd) real(r8) :: dbfrmaxdd real(r8) :: slascaler - select case(int(prt_params%allom_fmode(ipft))) + select case(prt_params%allom_fmode(ipft)) case(1) ! "constant proportionality with TRIMMED target bleaf" call blmax_allom(d,ipft,blmax,dblmaxdd) @@ -1139,7 +1157,7 @@ subroutine bstore_allom(d,ipft,crowndamage, canopy_trim,bstore,dbstoredd) associate( allom_stmode => prt_params%allom_stmode(ipft), & cushion => prt_params%cushion(ipft) ) - select case(int(allom_stmode)) + select case(allom_stmode) case(1) ! Storage is constant proportionality of trimmed maximum leaf ! biomass (ie cushion * bleaf), and thus leaf phenology is ignored. call bleaf(d,ipft, crowndamage, canopy_trim, 1.0_r8, bl, dbldd) @@ -1190,7 +1208,7 @@ subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeadd associate( agb_fraction => prt_params%allom_agb_frac(ipft)) - select case(int(prt_params%allom_amode(ipft))) + select case(prt_params%allom_amode(ipft)) case(1) ! Saldariagga mass allometry originally calculated bdead directly. ! we assume proportionality between bdead and bagw @@ -1199,7 +1217,7 @@ subroutine bdead_allom(bagw,bbgw,bsap,ipft,bdead,dbagwdd,dbbgwdd,dbsapdd,dbdeadd dbdeaddd = dbagwdd/agb_fraction end if - case(2,3) + case(2,3,4) bdead = bagw + bbgw - bsap if(present(dbagwdd) .and. present(dbbgwdd) .and. & @@ -1318,7 +1336,7 @@ subroutine bsap_ltarg_slatop(d,h,dhdd,bleaf,dbleafdd,ipft, & hbl2bsap = slatop * g_per_kg * wood_density * kg_per_Megag / (c2b*cm2_per_m2 ) ! Calculate area. Note that no c2b conversion here, because it is - ! wood density that is in biomass units, SLA is in units [m2/gC. + ! wood density that is in biomass units, SLA is in units [m2/gC]. ! [m2] = [m2/gC] * [kgC] * [gC/kgC] / ( [m2/cm2] * [cm2/m2]) la_per_sa = la_per_sa_int + h*la_per_sa_slp sapw_area = slatop * bleaf * g_per_kg / (la_per_sa*cm2_per_m2 ) @@ -1489,7 +1507,109 @@ subroutine dh2blmax_2pwr(d,p1,p2,dbh_maxh,c2b,blmax,dblmaxdd) return end subroutine dh2blmax_2pwr - + + ! =========================================================================== + + subroutine dh2blmax_3pwr(d,h,dhdd,p1,p2,p3,slatop,c2b,dbh_maxh,blmax,dblmaxdd) + !-------------------------------------------------------------------------- + ! + ! This function calculates the maximum leaf biomass from reference + ! diameter, plant height and top-of-the-canopy (fully sunlit) SLA. This + ! functional form is similar to Lescure et al. (1983) and Longo et al. + ! (2020), except that it uses SLA as an additional scaler for the + ! allometric equation that can have a different exponent from + ! (DBH^2 * Height). + ! + ! ----------------- + ! References + ! ----------------- + ! Lescure JP, Puig H, Riera B, Leclerc D, Beekman A , Beneteau A. 1983. + ! La phytomasse epigee d'une foret dense en Guiane Francaise + ! Acta Oecol.-Oec. Gen. 4: 237-251. + ! URL http://www.documentation.ird.fr/hor/fdi:010005089 + ! + ! Longo M, Saatchi SS, Keller M, Bowman KW, Ferraz A, Moorcroft PR, + ! Morton D, Bonal D, Brando P, Burban B et al. 2020. Impacts of + ! degradation on water, energy, and carbon cycling of the Amazon + ! tropical forests. J. Geophys. Res.-Biogeosci. 125: + ! e2020JG005677. doi:10.1029/2020JG005677. + ! + ! ----------------- + ! Input arguments + ! ----------------- + ! d -- Diameter at breast height [ cm] + ! h -- Total tree height [ m] + ! dhdd -- Height derivative with dbh [ m/cm] + ! p1 -- Parameter 1 (log-intercept) [ --] + ! p2 -- Parameter 2 (power, or log-slope) [ --] + ! p3 -- Parameter 3 (power, or log-slope) [ --] + ! slatop -- Top-of-canopy specific leaf area [ m2/gC] + ! c2b -- Carbon to biomass multiplier ~ 2 [ kg/kgC] + ! dbh_maxh -- DBH at maximum height [ cm] + ! + ! ------------------ + ! Output arguments + ! ------------------ + ! blmax -- Plant leaf biomass [ kgC] + ! dblmaxdd -- Plant leaf biomass derivative [ kgC/cm] + ! + !-------------------------------------------------------------------------- + + + !--- Arguments + real(r8), intent(in) :: d ! plant diameter [ cm] + real(r8), intent(in) :: h ! plant height [ m] + real(r8), intent(in) :: dhdd ! Height derivative wrt diameter [ m/cm] + real(r8), intent(in) :: p1 ! Log-intercept parameter [ -] + real(r8), intent(in) :: p2 ! Log-slope parameter for size [ -] + real(r8), intent(in) :: p3 ! Log-slope parameter for SLA [ -] + real(r8), intent(in) :: slatop ! Top canopy specific leaf area [ m2/gC] + real(r8), intent(in) :: c2b ! Carbon to biomass multiplier [kg/kgC] + real(r8), intent(in) :: dbh_maxh ! dbh at maximum height [ cm] + real(r8), intent(out) :: blmax ! Leaf biomass [ kgC] + real(r8), intent(out), optional :: dblmaxdd ! Leaf biomass derivative [kgC/cm] + !--- Local variables + real(r8) :: duse + !---~--- + + + + !--- Cap DBH + duse = min(d,dbh_maxh) + !---~--- + + + !--- Find leaf biomass + blmax = p1 * (duse*duse*h)**p2 * slatop**p3 / c2b + !---~--- + + + !---~--- + ! Compute the leaf biomass derivative with DBH if needed. + !---~--- + if (present(dblmaxdd)) then + if (d >= dbh_maxh) then + !---~--- + ! Leaf area is capped at the maximum DBH. This may be removed in the + ! future. + !---~--- + dblmaxdd = 0._r8 + !---~--- + else + !---~--- + ! Find the leaf biomass derivative, noting that height is actually + ! a function of DBH. + !---~--- + dblmaxdd = p2 * blmax * ( 2._r8 / duse + dhdd / h ) + !---~--- + end if + end if + !---~--- + + return + end subroutine dh2blmax_3pwr + + ! ========================================================================= ! Diameter to height (D2H) functions ! ========================================================================= @@ -1500,7 +1620,7 @@ subroutine d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd) ! "d to height via Chave et al. 2014" ! This function calculates tree height based on tree diameter and the - ! environmental stress factor "E", as per Chave et al. 2015 GCB + ! environmental stress factor "E", as per Chave et al. 2014 GCB ! As opposed to previous allometric models in ED, in this formulation ! we do not impose a hard cap on tree height. But, maximum_height ! is an important parameter, but instead of imposing a hard limit, in @@ -1509,7 +1629,7 @@ subroutine d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd) ! begin to route available NPP into seed and defense respiration. ! ! The stress function is based on the geographic location of the site. If - ! a user decides to use Chave2015 allometry, the E factor will be read in + ! a user decides to use Chave2014 allometry, the E factor will be read in ! from a global gridded dataset and assigned for each ED patch (note it ! will be the same for each ED patch, but this distinction will help in ! porting ED into different models (patches are pure ED). It @@ -1747,7 +1867,7 @@ subroutine dh2bagw_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bagw,dbagwdd) ! height and wood density. ! ! Chave et al. Improved allometric models to estimate the abovegroud - ! biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2015. + ! biomass of tropical trees. Global Change Biology. V20, p3177-3190. 2014. ! ! Input arguments: ! d: Diameter at breast height [cm] @@ -1791,6 +1911,90 @@ subroutine dh2bagw_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bagw,dbagwdd) return end subroutine dh2bagw_chave2014 + ! ============================================================================ + + + + subroutine dh2bagw_3pwr(d,h,dhdd,p1,p2,p3,wood_density,c2b,bagw,dbagwdd) + !-------------------------------------------------------------------------- + ! + ! This function calculates the maximum above-ground biomass from + ! reference diameter, plant height and wood density. This functional form + ! is an intermediate between Saldarriaga et al. (1988) and Chave et al. + ! (2014), because the wood-density exponent is independent on the + ! plant size (DBH^2 * Height) but the diameter and height are scaled + ! together through the size function. + ! + ! ----------------- + ! References + ! ----------------- + ! Chave J, Rejou-Mechain M, Burquez A, Chidumayo E, Colgan MS, Delitti WB, + ! Duque A, Eid T, Fearnside PM, Goodman RC et al. 2014. Improved + ! allometric models to estimate the aboveground biomass of tropical + ! trees. Glob. Change Biol. 20: 3177-3190. doi:10.1111/gcb.12629. + ! + ! Saldarriaga JG, West DC, Tharp ML , Uhl C. 1988. Long-term + ! chronosequence of forest succession in the upper Rio Negro of + ! Colombia and Venezuela. J. Ecol. 76: 938-958. + ! doi:10.2307/2260625. + ! + ! ----------------- + ! Input arguments + ! ----------------- + ! d -- Diameter at breast height [ cm] + ! h -- Total tree height [ m] + ! dhdd -- Height derivative with dbh [ m/cm] + ! p1 -- Parameter 1 (log-intercept) [ --] + ! p2 -- Parameter 2 (power, or log-slope) [ --] + ! p3 -- Parameter 3 (power, or log-slope) [ --] + ! wood_density -- Wood density [ g/cm3] + ! c2b -- Carbon to biomass multiplier ~ 2 [ kg/kgC] + ! dbh_maxh -- DBH at maximum height [ cm] + ! + ! ------------------ + ! Output arguments + ! ------------------ + ! bagw -- Above-ground biomass per individual [ kgC] + ! dbagwdd -- Above-ground biomass derivative [ kgC/cm] + ! + !-------------------------------------------------------------------------- + + + !--- Arguments + real(r8), intent(in) :: d ! plant diameter [ cm] + real(r8), intent(in) :: h ! plant height [ m] + real(r8), intent(in) :: dhdd ! Height deriv. wrt diameter [ m/cm] + real(r8), intent(in) :: p1 ! Log-intercept parameter [ -] + real(r8), intent(in) :: p2 ! Log-slope parameter (size) [ -] + real(r8), intent(in) :: p3 ! Log-slope parameter (WD) [ -] + real(r8), intent(in) :: wood_density ! Wood density [ g/cm3] + real(r8), intent(in) :: c2b ! Carbon to biomass factor [kg/kgC] + real(r8), intent(out) :: bagw ! Above-ground biomass [ kgC] + real(r8), intent(out), optional :: dbagwdd ! AG biomass derivative [kgC/cm] + !---~--- + + + !--- Find above-ground biomass + bagw = p1 * (d*d*h)**p2 * wood_density**p3 / c2b + !---~--- + + + !---~--- + ! Compute the above-ground biomass derivative with DBH if needed, noting that + ! height is actually a function of DBH. + !---~--- + if (present(dbagwdd)) then + dbagwdd = p2 * bagw * ( 2._r8 / d + dhdd / h ) + end if + !---~--- + + return + end subroutine dh2bagw_3pwr + + + ! ============================================================================ + + subroutine d2bagw_2pwr(d,p1,p2,c2b,bagw,dbagwdd) ! ========================================================================= @@ -2077,36 +2281,62 @@ subroutine h2d_martcano(h,p1,p2,p3,d,dddh) return end subroutine h2d_martcano - ! ===================================================================================== + ! ===================================================================================== - subroutine CrownDepth(height,ft,crown_depth) - ! ----------------------------------------------------------------------------------- - ! This routine returns the depth of a plant's crown. Which is the length - ! from the bottom of the crown to the top in the vertical dimension. - ! - ! This code may be used as a wrapper if different hypotheses are wished to be - ! optioned. - ! ----------------------------------------------------------------------------------- - - real(r8),intent(in) :: height ! The height of the plant [m] - integer,intent(in) :: ft ! functional type index - real(r8),intent(out) :: crown_depth ! The depth of the crown [m] + subroutine CrownDepth(height,ipft,crown_depth) - ! Alternative Hypothesis: - ! crown depth from Poorter, Bongers & Bongers - ! crown_depth = exp(-1.169_r8)*cCohort%height**1.098_r8 + !-------------------------------------------------------------------------- + ! This routine returns the depth of a plant's crown. Which is the + ! length from the bottom of the crown to the top in the vertical dimension. + ! + ! The original mode (now allom_dmode = 1) used only a fraction of the + ! plant's height. Alternatively, allom_dmode = 2 uses the same functional + ! form as Poorter et al. (2006), with an additional constraint to prevent + ! crown depth to exceed the plant's height. + ! + ! ----------------- + ! References + ! ----------------- + ! Poorter L, Bongers L , Bongers F. 2006. Architecture of 54 moist-forest + ! tree species: traits, trade-offs, and functional groups. Ecology 87: + ! 1289-1301. doi:10.1890/0012-9658(2006)87[1289:AOMTST]2.0.CO;2. + ! + !-------------------------------------------------------------------------- + + real(r8),intent(in) :: height ! The height of the plant [m] + integer ,intent(in) :: ipft ! functional type index + real(r8),intent(out) :: crown_depth ! The depth of the crown [m] + + associate( p1 => prt_params%allom_h2cd1(ipft), & + p2 => prt_params%allom_h2cd2(ipft), & + allom_dmode => prt_params%allom_dmode(ipft)) + + select case (allom_dmode) + case (1) ! Default, linear relationship with height + crown_depth = p1 * height + case (2) ! Power law, akin to Poorter et al. (2006). + !---~--- + ! Apply the two coefficients, but make sure crown depth does not exceed + ! the plant's height. + !---~--- + crown_depth = min(height, p1 * height ** p2) + !---~--- + case default + write(fates_log(),*) 'Invalid settings for crown depth mode for PFT ',ipft,'.' + write(fates_log(),*) 'Current allom_dmode: ',allom_dmode,'.' + write(fates_log(),*) 'Valid allom_dmode values: 1 or 2.' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + end associate - ! Alternative Hypothesis: - ! Original FATES crown depth heigh used for hydraulics - ! crown_depth = min(height,0.1_r8) + return + end subroutine CrownDepth + ! ===================================================================================== - crown_depth = prt_params%crown_depth_frac(ft) * height - - return - end subroutine CrownDepth @@ -2173,7 +2403,99 @@ subroutine carea_2pwr(dbh,spread,d2bl_p2,d2bl_ediff,d2ca_min, & endif end subroutine carea_2pwr - + + + ! ============================================================================= + + + subroutine carea_3pwr(dbh,height,ipft,dbh_maxh,spread,dh2bl_p2,dh2bl_ediff, & + dh2ca_min,dh2ca_max,crowndamage,c_area,inverse) + !---~--- + ! Calculate area of ground covered by entire cohort. (m2) + ! Function of DBH (cm), height (m), canopy spread (m/cm) and number of + ! individuals. + !---~--- + + !--- List of arguments + real(r8) , intent(inout) :: dbh ! Diameter at breast/ref/ height [ cm] + real(r8) , intent(inout) :: height ! Height [ m] + integer(i4), intent(in) :: ipft ! PFT index + real(r8) , intent(in) :: dbh_maxh ! Minimum DBH at maximum height [ m] + real(r8) , intent(in) :: spread ! site level relative spread score [ 0-1] + real(r8) , intent(in) :: dh2bl_p2 ! Exponent for size (bleaf) [ -] + real(r8) , intent(in) :: dh2bl_ediff ! Difference in size exponent [ -] + ! between crown area and bleaf + real(r8) , intent(in) :: dh2ca_min ! Minimum (closed forest) scaling [ -] + ! coefficient for crown area + real(r8) , intent(in) :: dh2ca_max ! Maximum (savannah) scaling [ -] + ! coefficient for crown area + integer , intent(in) :: crowndamage ! Crown damage class [ -] + ! [1: undamaged, >1: damaged] + 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 + !--- Local variables + real(r8) :: size ! Size (Diameter^2 * Height) [cm2 m] + real(r8) :: dh2ca_p1 ! Effective scaling factor (crown area) [ -] + real(r8) :: dh2ca_p2 ! Effective exponent (crown area) [ -] + real(r8) :: crown_reduction ! Crown area reduction due to damage. [ -] + !---~--- + + + !---~--- + ! Define the scaling (log-intercept) and exponent (log-slope) parameters for + ! crown area. The scaling parameter accounts for the site-level spread elasticity. + ! The exponent is defined in terms of the leaf biomass exponent plus an offset + ! parameter (allom_blca_expnt_diff). This is done because the default in FATES is + ! for both exponents to be same (i.e., allom_blca_expnt_diff = 0.) so the per-plant + ! canopy area remains invariant during growth. However, allometric models in general + ! predict that leaf area grows faster than crown area. + !---~--- + dh2ca_p1 = spread * dh2ca_max + (1._r8 - spread) * dh2ca_min + dh2ca_p2 = dh2bl_p2 + dh2bl_ediff + !---~--- + + + + !---~--- + ! Decide whether to use DBH and height to find crown area (default) or the + ! other way round. + !---~--- + select case (inverse) + case (.false.) + !--- Find the maximum area + size = dbh * dbh * height + c_area = dh2ca_p1 * size ** dh2ca_p2 + !---~--- + + !--- Reduce area if the crown is damaged. + if (crowndamage > 1) then + call GetCrownReduction(crowndamage, crown_reduction) + c_area = c_area * (1.0_r8 - crown_reduction) + end if + !---~--- + case (.true.) + !--- Reduce area if the crown is damaged. + if (crowndamage > 1) then + call GetCrownReduction(crowndamage, crown_reduction) + c_area = c_area * (1.0_r8 - crown_reduction) + end if + !---~--- + + + !---~--- + ! Find the size, then use a root-finding algorithm to find DBH. + !---~--- + size = ( c_area / dh2ca_p1 ) ** ( 1.0_r8 / dh2ca_p2 ) + call size2dbh(size,ipft,dbh,dbh_maxh) + !---~--- + end select + !---~--- + + return + end subroutine carea_3pwr + + ! ========================================================================= subroutine set_root_fraction(root_fraction, ft, zi, max_nlevroot) @@ -2625,7 +2947,265 @@ subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx) dydx = (y2-y1)/(x2-x1) + (1.0_r8-2.0_r8*t)*(a*(1.0_r8-t)+b*t)/(x2-x1) + t*(1.0_r8-t)*(b-a)/(x2-x1) return end subroutine cspline - + ! ============================================================================ + + + + ! ============================================================================ + ! This function finds the DBH when size (DBH^2 * Height) is known but we + ! cannot find DBH analytically due to the non-linear relationship between DBH + ! and height. This is borrowed from the same approach applied in ED2 for + ! root finding. It starts with the Newton's method, which should quickly + ! converge to the solution. In the unlikely case of failure, we use the + ! Regula Falsi (Illinois) method as a back-up. + ! ============================================================================ + subroutine size2dbh(size,ipft,dbh,dbh_maxh) + !--- Arguments. + real(r8) , intent(in) :: size ! Size (DBH^2 * Height) [cm2 m] + integer(i4), intent(in) :: ipft ! PFT index [ -] + real(r8) , intent(inout) :: dbh ! Diameter at breast height [ cm] + real(r8) , intent(in) :: dbh_maxh ! Minimum DBH at maximum height [ cm] + !--- Local variables + real(r8) :: hgt ! Height [ m] + real(r8) :: dhgtddbh ! Height derivative [ m/cm] + real(r8) :: size_maxh ! Minimum size at maximum height [cm2 m] + real(r8) :: deriv ! Function derivative [ cm m] + real(r8) :: afun ! Function value (lower guess) [cm2 m] + real(r8) :: rfun ! Function value (RF new guess) [cm2 m] + real(r8) :: zfun ! Function value (upper guess) [cm2 m] + real(r8) :: adbh ! DBH: lower guess [ cm] + real(r8) :: rdbh ! DBH: updated guess (Reg. Falsi) [ cm] + real(r8) :: zdbh ! DBH: upper guess [ cm] + real(r8) :: delta ! Second guess for the RF method [ cm] + integer :: itn ! Iteration counter -- Newton [ -] + integer :: iti ! Iteration counter -- Reg. Falsi [ -] + logical :: converged ! Has the solution converged? [ T|F] + logical :: zside ! Converging on the upper size? [ T|F] + !--- Local constants. + real(r8) , parameter :: toler =1.0e-12_r8 ! Relative tolerance [ --] + integer , parameter :: maxit_newt = 10 ! Cap in iterations -- Newton [ --] + integer , parameter :: maxit_rf = 100 ! Cap in iterations -- Reg. Falsi [ --] + !---~--- + + + !---~--- + ! Find the maximum size beyond which the height is assumed constant. In this + ! case, DBH can be determined without the iterative approach. + !---~--- + call h_allom(dbh_maxh,ipft,hgt) + size_maxh = dbh_maxh * dbh_maxh * hgt + if (size >= size_maxh) then + dbh = sqrt(size/hgt) + return + end if + !---~--- + + + !--- First guess: use current DBH. + adbh = dbh + call h_allom(adbh,ipft,hgt,dhgtddbh) + afun = adbh * adbh * hgt - size + deriv = 2.0_r8 * adbh * hgt + adbh * adbh * dhgtddbh + !---~--- + + + !--- Copy just in case it fails at the first iteration. + zdbh = adbh + zfun = afun + !---~--- + + + !---~--- + ! Enter the Newton's method loop + !---~--- + converged = .false. + newton_loop: do itn = 1, maxit_newt + !--- If derivative is too flat, go to Regula Falsi + if ( abs(deriv) < toler) exit newton_loop + !---~--- + + + !--- Copy the previous guess. + adbh = zdbh + afun = zfun + !---~--- + + + !--- Find the new guess, and evaluate the function and derivative. + zdbh = adbh - afun / deriv + call h_allom(zdbh,ipft,hgt,dhgtddbh) + zfun = zdbh * zdbh * hgt - size + deriv = 2.0_r8 * zdbh * hgt + zdbh * zdbh * dhgtddbh + !---~--- + + !--- Check convergence. + converged = abs(adbh - zdbh) < toler * zdbh + if (converged) then + !--- Convergence by iterations. + dbh = 0.5_r8 * (adbh + zdbh) + return + !---~--- + else if (abs(zfun) < nearzero) then + !--- Convergence by luck. + dbh = zdbh + return + !---~--- + end if + !---~--- + end do newton_loop + !---~--- + + + + !---~--- + ! If we have reached this point, then Newton's method has failed. Use Regula + ! Falsi instead. For this, we must have two guesses whose function evaluation has + ! opposite signs. + !---~--- + if (afun * zfun <= -nearzero) then + !--- We already have two guesses with opposite signs. + zside = .true. + !---~--- + else + !--- Look for another guess with opposite sign. + if (abs(zfun-afun) < 100._r8 * toler * adbh) then + delta = 100._r8 * toler * adbh + else + delta = max( abs( afun * (zdbh-adbh) / (zfun-afun) ),100._r8 * toler * adbh ) + end if + !---~--- + + + !---~--- + ! Try guesses on both sides of the first guess, sending guesses increasingly + ! further away until we find a good guess. + !---~--- + zdbh = adbh + delta + zside = .false. + zguess_loop: do iti=1,maxit_rf + zdbh = adbh + real((-1)**iti * (iti+3)/2,r8) * delta + call h_allom(zdbh,ipft,hgt) + zfun = zdbh * zdbh * hgt - size + zside = afun * zfun < -nearzero + if (zside) exit zguess_loop + end do zguess_loop + + !---~--- + ! Issue an error in case the function failed finding a second guess. + !---~--- + if (.not. zside) then + write (unit=*,fmt='(a)') '---~---' + write (unit=*,fmt='(a)') ' Failed finding the second guess:' + write (unit=*,fmt='(a)') '---~---' + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') ' Input: ' + write (unit=*,fmt='(a,1x,es14.7)') ' + size =',size + write (unit=*,fmt='(a,1x,es14.7)') ' + dbh =',dbh + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') ' Current guesses and evaluations:' + write (unit=*,fmt='(a,1x,es14.7)') ' + adbh =',adbh + write (unit=*,fmt='(a,1x,es14.7)') ' + afun =',afun + write (unit=*,fmt='(a,1x,es14.7)') ' + zdbh =',zdbh + write (unit=*,fmt='(a,1x,es14.7)') ' + zfun =',zfun + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') '---~---' + write(fates_log(),*) 'Second guess for Regula Falsi method not found.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + !---~--- + end if + !---~--- + + !---~--- + ! Proceed to the regula falsi loop. + !---~--- + regfalsi_loop: do iti=1,maxit_rf + + !--- Update solution. + rdbh = ( zfun * adbh - afun * zdbh ) / ( zfun - afun) + !---~--- + + + !---~--- + ! Check for convergence. In case it converged, we can exit the sub-routine. + !---~--- + converged = abs(rdbh - adbh) < toler * max(rdbh,adbh) + if (converged) exit regfalsi_loop + !---~--- + + + !--- Find the new function evaluation. + call h_allom(rdbh,ipft,hgt) + rfun = rdbh * rdbh * hgt - size + !---~--- + + + !---~--- + ! Define the new searching interval based on the intermediate value theorem. + !---~--- + if (abs(rfun) < nearzero) then + !--- Converged by luck. + converged = .true. + exit regfalsi_loop + !---~--- + else if (rfun * afun <= -nearzero ) then + !--- Guess is between lower and current guess. + zdbh = rdbh + zfun = rfun + !--- If we are updating the upper side again, halve afun (Regula Falsi method). + if (zside) afun = afun * 0.5_r8 + !--- Flag that we have just updated the upper side. + zside = .true. + !---~--- + else + !--- Guess is between current and upper guess. + adbh = rdbh + afun = rfun + !--- If we are updating the lower side again, halve zfun (Regula Falsi method). + if (.not. zside) zfun = zfun * 0.5_r8 + !--- Flag that we have just updated the lower side. + zside = .false. + end if + end do regfalsi_loop + !---~--- + + + !---~--- + ! Check that the Regula Falsi method indeed converged. + !---~--- + if (converged) then + !--- Yes, return the last guess + dbh = rdbh + !---~--- + else + !--- No, report the bad news + write (unit=*,fmt='(a)') '---~---' + write (unit=*,fmt='(a)') ' Failed finding the solution:' + write (unit=*,fmt='(a)') '---~---' + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') ' Input: ' + write (unit=*,fmt='(a,1x,es14.7)') ' + size =',size + write (unit=*,fmt='(a,1x,es14.7)') ' + dbh =',dbh + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') ' Current guesses and evaluations:' + write (unit=*,fmt='(a,1x,es14.7)') ' + adbh =',adbh + write (unit=*,fmt='(a,1x,es14.7)') ' + afun =',afun + write (unit=*,fmt='(a,1x,es14.7)') ' + rdbh =',rdbh + write (unit=*,fmt='(a,1x,es14.7)') ' + rfun =',rfun + write (unit=*,fmt='(a,1x,es14.7)') ' + zdbh =',zdbh + write (unit=*,fmt='(a,1x,es14.7)') ' + zfun =',zfun + write (unit=*,fmt='(a)') ' ' + write (unit=*,fmt='(a)') '---~---' + write(fates_log(),*) 'Size to DBH routine failed to converge!' + call endrun(msg=errMsg(sourcefile, __LINE__)) + !---~--- + end if + !---~--- + + return + end subroutine size2dbh + ! ============================================================================ + end module FatesAllometryMod diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index f170fe2275..1a3fde13e9 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -80,9 +80,6 @@ variables: double fates_allom_cmode(fates_pft) ; fates_allom_cmode:units = "index" ; fates_allom_cmode:long_name = "coarse root biomass allometry function index." ; - double fates_allom_crown_depth_frac(fates_pft) ; - fates_allom_crown_depth_frac:units = "fraction" ; - fates_allom_crown_depth_frac:long_name = "the depth of a cohort crown as a fraction of its height" ; double fates_allom_d2bl1(fates_pft) ; fates_allom_d2bl1:units = "variable" ; fates_allom_d2bl1:long_name = "Parameter 1 for d2bl allometry" ; @@ -110,6 +107,9 @@ variables: double fates_allom_dbh_maxheight(fates_pft) ; fates_allom_dbh_maxheight:units = "cm" ; fates_allom_dbh_maxheight:long_name = "the diameter (if any) corresponding to maximum height, diameters may increase beyond this" ; + double fates_allom_dmode(fates_pft) ; + fates_allom_dmode:units = "index" ; + fates_allom_dmode:long_name = "crown depth allometry function index." ; double fates_allom_fmode(fates_pft) ; fates_allom_fmode:units = "index" ; fates_allom_fmode:long_name = "fine root biomass allometry function index." ; @@ -125,6 +125,12 @@ variables: double fates_allom_frbstor_repro(fates_pft) ; fates_allom_frbstor_repro:units = "fraction" ; fates_allom_frbstor_repro:long_name = "fraction of bstore goes to reproduction after plant dies" ; + double fates_allom_h2cd1(fates_pft) ; + fates_allom_h2cd1:units = "variable" ; + fates_allom_h2cd1:long_name = "Parameter 1 for h2cd allometry (exp(log-intercept) or scaling). If allom_dmode=1; this is the same as former crown_depth_frac parameter" ; + double fates_allom_h2cd2(fates_pft) ; + fates_allom_h2cd2:units = "variable" ; + fates_allom_h2cd2:long_name = "Parameter 2 for h2cd allometry (log-slope or exponent). If allom_dmode=1; this is not needed (as exponent is assumed 1)" ; double fates_allom_hmode(fates_pft) ; fates_allom_hmode:units = "index" ; fates_allom_hmode:long_name = "height allometry function index." ; @@ -971,9 +977,6 @@ data: fates_allom_cmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; - fates_allom_crown_depth_frac = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.95, 0.95, - 0.95, 1, 1, 1 ; - fates_allom_d2bl1 = 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07 ; @@ -1002,6 +1005,8 @@ data: fates_allom_dbh_maxheight = 90, 80, 80, 80, 90, 80, 3, 3, 2, 0.35, 0.35, 0.35 ; + fates_allom_dmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; + fates_allom_fmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; fates_allom_fnrt_prof_a = 7, 7, 7, 7, 6, 6, 7, 7, 7, 11, 11, 11 ; @@ -1012,6 +1017,11 @@ data: fates_allom_frbstor_repro = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_allom_h2cd1 = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.95, 0.95, + 0.95, 1, 1, 1 ; + + fates_allom_h2cd2 = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; + fates_allom_hmode = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; fates_allom_l2fr = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90 index f2419c3334..bb3f98f6be 100644 --- a/parteh/PRTParametersMod.F90 +++ b/parteh/PRTParametersMod.F90 @@ -117,24 +117,24 @@ module PRTParametersMod real(r8), allocatable :: c2b(:) ! Carbon to biomass multiplier [kg/kgC] real(r8), allocatable :: wood_density(:) ! wood density g cm^-3 ... - real(r8), allocatable :: crown_depth_frac(:) ! fraction of the height of the plant integer , allocatable :: woody(:) ! Does the plant have wood? (1=yes, 0=no) ! that is occupied by crown real(r8), allocatable :: slamax(:) ! Maximum specific leaf area of plant (at bottom) [m2/gC] real(r8), allocatable :: slatop(:) ! Specific leaf area at canopy top [m2/gC] real(r8), allocatable :: allom_sai_scaler(:) ! real(r8), allocatable :: allom_dbh_maxheight(:) ! dbh at which height growth ceases - real(r8), allocatable :: allom_hmode(:) ! height allometry function type - real(r8), allocatable :: allom_lmode(:) ! maximum leaf allometry function type - real(r8), allocatable :: allom_fmode(:) ! maximum root allometry function type - real(r8), allocatable :: allom_amode(:) ! AGB allometry function type - real(r8), allocatable :: allom_cmode(:) ! Coarse root allometry function type - real(r8), allocatable :: allom_smode(:) ! sapwood allometry function type - real(r8), allocatable :: allom_stmode(:) ! storage allometry functional type + integer , allocatable :: allom_hmode(:) ! height allometry function type + integer , allocatable :: allom_lmode(:) ! maximum leaf allometry function type + integer , allocatable :: allom_fmode(:) ! maximum root allometry function type + integer , allocatable :: allom_amode(:) ! AGB allometry function type + integer , allocatable :: allom_cmode(:) ! Coarse root allometry function type + integer , allocatable :: allom_smode(:) ! sapwood allometry function type + integer , allocatable :: allom_stmode(:) ! storage allometry functional type ! 0 - storage is proportional to maximum leaf biomass ! (considering trimmed) ! 1 - storage is proportional to maximum leaf biomass ! (untrimmed) + integer , allocatable :: allom_dmode(:) ! crown depth allometry function type ! (HARD-CODED FOR TIME BEING, RGK 11-2017) real(r8), allocatable :: allom_la_per_sa_int(:) ! Leaf area to sap area conversion, intercept ! (sapwood area / leaf area) [cm2/m2] @@ -161,6 +161,12 @@ module PRTParametersMod real(r8), allocatable :: allom_agb3(:) ! Parameter 3 for agb allometry real(r8), allocatable :: allom_agb4(:) ! Parameter 3 for agb allometry + real(r8), allocatable :: allom_h2cd1(:) ! Parameter 1 for crown depth allometry. When allom_dmode=1 + ! this is fraction of the height of the plant that is + ! considered crown (former parameter crown_depth_frac). + real(r8), allocatable :: allom_h2cd2(:) ! Exponent for crown depth allometry. Used only when + ! allom_dmode /= 1. + real(r8), allocatable :: allom_zroot_max_dbh(:) ! dbh at which maximum rooting depth saturates (largest possible) [cm] real(r8), allocatable :: allom_zroot_max_z(:) ! the maximum rooting depth defined at dbh = fates_allom_zroot_max_dbh [m] real(r8), allocatable :: allom_zroot_min_dbh(:) ! dbh at which the maximum rooting depth for a recruit is defined [cm] diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90 index 098b1e3738..381546f1d9 100644 --- a/parteh/PRTParamsFATESMod.F90 +++ b/parteh/PRTParamsFATESMod.F90 @@ -201,10 +201,6 @@ subroutine PRTRegisterPFT(fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_allom_crown_depth_frac' - call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names, lower_bounds=dim_lower_bound) - name = 'fates_wood_density' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -321,6 +317,10 @@ subroutine PRTRegisterPFT(fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_dmode' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_la_per_sa_int' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -385,6 +385,14 @@ subroutine PRTRegisterPFT(fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_h2cd1' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_allom_h2cd2' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_zroot_max_dbh' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -501,10 +509,6 @@ subroutine PRTReceivePFT(fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=prt_params%fnrt_prof_mode) - name = 'fates_allom_crown_depth_frac' - call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%crown_depth_frac) - name = 'fates_woody' call fates_params%RetrieveParameterAllocate(name=name, & data=tmpreal) @@ -566,31 +570,59 @@ subroutine PRTReceivePFT(fates_params) name = 'fates_allom_hmode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_hmode) + data=tmpreal) + allocate(prt_params%allom_hmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_hmode) + deallocate(tmpreal) name = 'fates_allom_lmode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_lmode) + data=tmpreal) + allocate(prt_params%allom_lmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_lmode) + deallocate(tmpreal) name = 'fates_allom_fmode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_fmode) + data=tmpreal) + allocate(prt_params%allom_fmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_fmode) + deallocate(tmpreal) name = 'fates_allom_amode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_amode) + data=tmpreal) + allocate(prt_params%allom_amode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_amode) + deallocate(tmpreal) name = 'fates_allom_stmode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_stmode) + data=tmpreal) + allocate(prt_params%allom_stmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_stmode) + deallocate(tmpreal) name = 'fates_allom_cmode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_cmode) + data=tmpreal) + allocate(prt_params%allom_cmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_cmode) + deallocate(tmpreal) name = 'fates_allom_smode' call fates_params%RetrieveParameterAllocate(name=name, & - data=prt_params%allom_smode) + data=tmpreal) + allocate(prt_params%allom_smode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_smode) + deallocate(tmpreal) + + name = 'fates_allom_dmode' + call fates_params%RetrieveParameterAllocate(name=name, & + data=tmpreal) + allocate(prt_params%allom_dmode(size(tmpreal,dim=1))) + call ArrayNint(tmpreal,prt_params%allom_dmode) + deallocate(tmpreal) name = 'fates_allom_la_per_sa_int' call fates_params%RetrieveParameterAllocate(name=name, & @@ -680,6 +712,14 @@ subroutine PRTReceivePFT(fates_params) call fates_params%RetrieveParameterAllocate(name=name, & data=prt_params%allom_agb4) + name = 'fates_allom_h2cd1' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%allom_h2cd1) + + name = 'fates_allom_h2cd2' + call fates_params%RetrieveParameterAllocate(name=name, & + data=prt_params%allom_h2cd2) + name = 'fates_allom_zroot_max_dbh' call fates_params%RetrieveParameterAllocate(name=name, & data=prt_params%allom_zroot_max_dbh) @@ -979,6 +1019,8 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'allom_agb2 = ',prt_params%allom_agb2 write(fates_log(),fmt0) 'allom_agb3 = ',prt_params%allom_agb3 write(fates_log(),fmt0) 'allom_agb4 = ',prt_params%allom_agb4 + write(fates_log(),fmt0) 'allom_h2cd1 = ',prt_params%allom_h2cd1 + write(fates_log(),fmt0) 'allom_h2cd2 = ',prt_params%allom_h2cd2 write(fates_log(),fmt0) 'allom_zroot_max_dbh = ',prt_params%allom_zroot_max_dbh write(fates_log(),fmt0) 'allom_zroot_max_z = ',prt_params%allom_zroot_max_z @@ -990,7 +1032,6 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'stoich_phos = ',prt_params%phos_stoich_p1 write(fates_log(),fmt0) 'alloc_organ_priority = ',prt_params%alloc_priority write(fates_log(),fmt0) 'woody = ',prt_params%woody - write(fates_log(),fmt0) 'crown_depth_frac = ',prt_params%crown_depth_frac write(fates_log(),fmt0) 'roota_par = ',prt_params%fnrt_prof_a write(fates_log(),fmt0) 'rootb_par = ',prt_params%fnrt_prof_b write(fates_log(),fmt0) 'fnrt_prof_mode = ',prt_params%fnrt_prof_mode @@ -1060,6 +1101,9 @@ subroutine PRTCheckParams(is_master) logical :: is_season_decid ! Is the PFT cold-deciduous? logical :: is_stress_decid ! Is the PFT drought-deciduous? logical :: is_semi_decid ! Is the PFT drought semi-deciduous? + integer :: nerror ! Count number of errors. If this is not + ! zero by theend of the subroutine, stop + ! the run. npft = size(prt_params%evergreen,1) @@ -1070,14 +1114,21 @@ subroutine PRTCheckParams(is_master) if(.not.is_master) return + ! Initialise nerror with zero. If anything is incorrectly set, nerror will be + ! positive, but we will hold on until all checks are performed before stopping + ! the run. + nerror = 0 if( any(prt_params%organ_id(:)<1) .or. & any(prt_params%organ_id(:)>num_organ_types) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'prt_organ_ids should match the global ids' write(fates_log(),*) 'of organ types found in PRTGenericMod.F90' write(fates_log(),*) 'organ_ids: ',prt_params%organ_id(:) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if ! Check to make sure the organ ids are valid if this is the @@ -1087,22 +1138,28 @@ subroutine PRTCheckParams(is_master) do io = 1,norgans if(prt_params%organ_id(io) == repro_organ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'reproductive tissues are a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' write(fates_log(),*) 'fates_prt_organ_id: ',prt_params%organ_id(:) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if if(prt_params%organ_id(io) == store_organ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'with flexible cnp or c-only alloc hypotheses' write(fates_log(),*) 'storage is a special case' write(fates_log(),*) 'and therefore should not be included in' write(fates_log(),*) 'the parameter file organ list' write(fates_log(),*) 'fates_prt_organ_id: ',prt_params%organ_id(:) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end do @@ -1112,10 +1169,13 @@ subroutine PRTCheckParams(is_master) ! between 0 and 1 if (hlm_parteh_mode .eq. prt_cnp_flex_allom_hyp) then if(any(prt_params%nfix_mresp_scfrac(:)<0._r8) .or. any(prt_params%nfix_mresp_scfrac(:)>1.0_r8)) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'The N fixation surcharge nfix_mresp_sfrac (fates_nfix1) must be between 0-1.' write(fates_log(),*) 'here are the values: ',prt_params%nfix_mresp_scfrac(:) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1136,14 +1196,17 @@ subroutine PRTCheckParams(is_master) ( is_evergreen .and. is_stress_decid ) .or. & ( is_season_decid .and. is_stress_decid ) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'PFT # ',ipft,' must be defined as having one of three' write(fates_log(),*) 'phenology habits, ie, only one of the flags below should' write(fates_log(),*) 'be different than ',ifalse write(fates_log(),*) 'stress_decid: ',prt_params%stress_decid(ipft) write(fates_log(),*) 'season_decid: ',prt_params%season_decid(ipft) write(fates_log(),*) 'evergreen: ',prt_params%evergreen(ipft) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if @@ -1155,6 +1218,7 @@ subroutine PRTCheckParams(is_master) ! In case the product of the lower and upper thresholds is negative, the ! thresholds are inconsistent as both should be defined using the same ! quantity. + write(fates_log(),*) '---~---' write(fates_log(),*) ' When using drought semi-deciduous phenology,' write(fates_log(),*) ' the moist threshold must have the same sign as' write(fates_log(),*) ' the dry threshold. Positive = soil water content [m3/m3],' @@ -1163,10 +1227,13 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' Stress_decid = ',prt_params%stress_decid(ipft) write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) write(fates_log(),*) ' fates_phen_moist_threshold = ',prt_params%phen_moist_threshold (ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 elseif ( prt_params%phen_drought_threshold(ipft) >= prt_params%phen_moist_threshold(ipft) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) ' When using drought semi-deciduous phenology,' write(fates_log(),*) ' the moist threshold must be greater than the dry threshold.' write(fates_log(),*) ' By greater we mean more positive or less negative, and' @@ -1175,8 +1242,10 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) ' Stress_decid = ',prt_params%stress_decid(ipft) write(fates_log(),*) ' fates_phen_drought_threshold = ',prt_params%phen_drought_threshold(ipft) write(fates_log(),*) ' fates_phen_moist_threshold = ',prt_params%phen_moist_threshold (ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1186,13 +1255,16 @@ subroutine PRTCheckParams(is_master) ! is bounded between 0 and 1 (exactly 0 and 1 are acceptable). if ( ( prt_params%phen_fnrt_drop_fraction(ipft) < 0.0_r8 ) .or. & ( prt_params%phen_fnrt_drop_fraction(ipft) > 1.0_r8 ) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) ' Abscission rate for fine roots must be between 0 and 1 for ' write(fates_log(),*) ' deciduous PFTs.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' evergreen flag: (should be 0):',prt_params%evergreen(ipft) write(fates_log(),*) ' phen_fnrt_drop_fraction: ', prt_params%phen_fnrt_drop_fraction(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if @@ -1203,21 +1275,27 @@ subroutine PRTCheckParams(is_master) if ( ( prt_params%woody(ipft) == itrue ) .and. & ( ( prt_params%phen_stem_drop_fraction(ipft) < 0.0_r8 ) .or. & ( prt_params%phen_stem_drop_fraction(ipft) > nearzero ) ) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) ' Non-zero stem-drop fractions are not allowed for woody plants' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' part_params%woody:',prt_params%woody(ipft) write(fates_log(),*) ' phen_stem_drop_fraction: ', prt_params%phen_stem_drop_fraction(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 elseif ( ( prt_params%phen_stem_drop_fraction(ipft) < 0.0_r8 ) .or. & ( prt_params%phen_stem_drop_fraction(ipft) > 1.0_r8 ) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) ' Deciduous non-wood plants must keep 0-100% of their stems' write(fates_log(),*) ' during the deciduous period.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' evergreen flag: (should be 0):',prt_params%evergreen(ipft) write(fates_log(),*) ' phen_stem_drop_fraction: ', prt_params%phen_stem_drop_fraction(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1229,14 +1307,16 @@ subroutine PRTCheckParams(is_master) if ( ( prt_params%seed_alloc(ipft) + & prt_params%seed_alloc_mature(ipft)) > 1.0_r8 ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'The sum of seed allocation from base and mature trees may' write(fates_log(),*) ' not exceed 1.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' seed_alloc: ',prt_params%seed_alloc(ipft) write(fates_log(),*) ' seed_alloc_mature: ',prt_params%seed_alloc_mature(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if ! Check if woody plants have a structural biomass (agb) intercept @@ -1244,21 +1324,46 @@ subroutine PRTCheckParams(is_master) if ( ( prt_params%allom_agb1(ipft) <= tiny(prt_params%allom_agb1(ipft)) ) .and. & ( prt_params%woody(ipft) .eq. 1 ) ) then + write(fates_log(),*) '---~---' write(fates_log(),*) 'Woody plants are expected to have a non-zero intercept' write(fates_log(),*) ' in the diameter to AGB allometry equations' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' allom_agb1: ',prt_params%allom_agb1(ipft) write(fates_log(),*) ' woody: ',prt_params%woody(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) '---~---' + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if + ! Check if parameter 2 for dbh -> height is negative when allom_hmode is 2 + ! (Weibull function / Poorter et al. (2006) + ! ML: FATES definition for parameter 2 is a bit unusual, which is why I added + ! the check. Normally the minus sign is left outside the parameter for + ! Weibull functions. + ! ---------------------------------------------------------------------------------- + if ( ( prt_params%allom_hmode(ipft) == 2 ) .and. & + ( prt_params%allom_d2h2 (ipft) > 0._r8 ) ) then + write(fates_log(),*) "---~---" + write(fates_log(),*) " Incorrect settings for height allometry." + write(fates_log(),*) ' PFT index: ',ipft + write(fates_log(),*) ' allom_hmode: ',prt_params%allom_hmode(ipft) + write(fates_log(),*) ' allom_d2h2: ',prt_params%allom_d2h2 (ipft) + write(fates_log(),*) " Parameter ""allom_d2h2"" must be negative when using" + write(fates_log(),*) " allom_hmode = 2." + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 + end if + ! Check if non-woody plants have structural biomass (agb) intercept ! ---------------------------------------------------------------------------------- ! if ( ( prt_params%allom_agb1(ipft) > tiny(prt_params%allom_agb1(ipft)) ) .and. & ! ( iprt_params%woody(ipft) .ne. 1 ) ) then ! +! write(fates_log(),*) "---~---" ! write(fates_log(),*) 'Non-woody plants are expected to have a zero intercept' ! write(fates_log(),*) ' in the diameter to AGB allometry equations' ! write(fates_log(),*) ' This is because the definition of AGB (as far as allometry)' @@ -1267,8 +1372,10 @@ subroutine PRTCheckParams(is_master) ! write(fates_log(),*) ' PFT#: ',ipft ! write(fates_log(),*) ' allom_agb1: ',prt_params%allom_agb1(ipft) ! write(fates_log(),*) ' woody: ',prt_params%woody(ipft) -! write(fates_log(),*) ' Aborting' -! call endrun(msg=errMsg(sourcefile, __LINE__)) +! write(fates_log(),*) "---~---" +! write(fates_log(),*) '' +! write(fates_log(),*) '' +! nerror = nerror + 1 ! ! end if @@ -1278,13 +1385,16 @@ subroutine PRTCheckParams(is_master) if ( ( prt_params%leaf_stor_priority(ipft) < 0.0_r8 ) .or. & ( prt_params%leaf_stor_priority(ipft) > 1.0_r8 ) ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'Prioritization of carbon allocation to leaf' write(fates_log(),*) ' and root turnover replacement, must be between' write(fates_log(),*) ' 0 and 1' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) 'leaf_stor_priority: ',prt_params%leaf_stor_priority(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if @@ -1293,18 +1403,26 @@ subroutine PRTCheckParams(is_master) ! Make sure nutrient storage fractions are positive if( prt_params%nitr_store_ratio(ipft) < 0._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'With parteh allometric CNP hypothesis' write(fates_log(),*) 'nitr_store_ratio must be > 0' write(fates_log(),*) 'PFT#: ',ipft write(fates_log(),*) 'nitr_store_ratio = ',prt_params%nitr_store_ratio(ipft) - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if if( prt_params%phos_store_ratio(ipft) < 0._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'With parteh allometric CNP hypothesis' write(fates_log(),*) 'phos_store_ratio must be > 0' write(fates_log(),*) 'PFT#: ',ipft write(fates_log(),*) 'nitr_store_ratio = ',prt_params%phos_store_ratio(ipft) - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if do i = 1,norgans @@ -1312,33 +1430,45 @@ subroutine PRTCheckParams(is_master) if(io == sapw_organ) then if ((prt_params%turnover_nitr_retrans(ipft,i) > nearzero)) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' Retranslocation of sapwood tissues should be zero.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' nitrogen retrans: ',prt_params%turnover_nitr_retrans(ipft,i) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if if ((prt_params%turnover_phos_retrans(ipft,i) > nearzero)) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' Retranslocation of sapwood tissues should be zero.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' phosphorus retrans: ',prt_params%turnover_nitr_retrans(ipft,i) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if elseif(io == struct_organ) then if ((prt_params%turnover_nitr_retrans(ipft,i) > nearzero)) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' Retranslocation of structural tissues should be zero.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' carbon retrans: ',prt_params%turnover_nitr_retrans(ipft,i) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if if ((prt_params%turnover_phos_retrans(ipft,i) > nearzero)) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' Retranslocation of structural tissues should be zero.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' phosphorus retrans: ',prt_params%turnover_nitr_retrans(ipft,i) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1347,13 +1477,16 @@ subroutine PRTCheckParams(is_master) (prt_params%turnover_phos_retrans(ipft,i) > 1.0_r8) .or. & (prt_params%turnover_nitr_retrans(ipft,i) < 0.0_r8) .or. & (prt_params%turnover_phos_retrans(ipft,i) < 0.0_r8)) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' Retranslocation should range from 0 to 1.' write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' parameter file organ index: ',i,' global index: ',io write(fates_log(),*) ' nitr: ',prt_params%turnover_nitr_retrans(ipft,i) write(fates_log(),*) ' phos: ',prt_params%turnover_phos_retrans(ipft,i) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end do @@ -1365,18 +1498,24 @@ subroutine PRTCheckParams(is_master) ! if (parteh_mode .eq. prt_carbon_allom_hyp) then if ( ( prt_params%grperc(ipft) < 0.0_r8) .or. & ( prt_params%grperc(ipft) > 1.0_r8 ) ) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' Growth respiration must be between 0 and 1: ',prt_params%grperc(ipft) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if ! elseif(parteh_mode .eq. prt_cnp_flex_allom_hyp) then ! if ( ( any(prt_params%grperc_organ(ipft,:) < 0.0_r8)) .or. & ! ( any(prt_params%grperc_organ(ipft,:) >= 1.0_r8)) ) then +! write(fates_log(),*) "---~---" ! write(fates_log(),*) ' PFT#: ',ipft ! write(fates_log(),*) ' Growth respiration must be between 0 and 1: ',prt_params%grperc_organ(ipft,:) -! write(fates_log(),*) ' Aborting' -! call endrun(msg=errMsg(sourcefile, __LINE__)) +! write(fates_log(),*) "---~---" +! write(fates_log(),*) '' +! write(fates_log(),*) '' +! nerror = nerror + 1 ! end if ! end if @@ -1385,11 +1524,14 @@ subroutine PRTCheckParams(is_master) ! The first nitrogen stoichiometry is used in all cases if ( (any(prt_params%nitr_stoich_p1(ipft,:) < 0.0_r8)) .or. & (any(prt_params%nitr_stoich_p1(ipft,:) >= 1.0_r8))) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' N per C stoichiometry must bet between 0-1' write(fates_log(),*) prt_params%nitr_stoich_p1(ipft,:) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end select @@ -1401,6 +1543,7 @@ subroutine PRTCheckParams(is_master) (prt_params%phos_stoich_p1(ipft,i) < 0._r8) .or. & (prt_params%nitr_stoich_p1(ipft,i) > 1._r8) .or. & (prt_params%phos_stoich_p1(ipft,i) > 1._r8) ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'When the C,N,P allocation hypothesis with flexible' write(fates_log(),*) 'stoichiometry is turned on (prt_cnp_flex_allom_hyp),' write(fates_log(),*) 'all stoichiometries must be greater than or equal to zero,' @@ -1410,18 +1553,23 @@ subroutine PRTCheckParams(is_master) write(fates_log(),*) 'organ index (see head of PRTGenericMod): ',io write(fates_log(),*) 'nitr_stoich: ',prt_params%nitr_stoich_p1(ipft,i) write(fates_log(),*) 'phos_stoich: ',prt_params%phos_stoich_p1(ipft,i) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end do if ( any(prt_params%alloc_priority(ipft,:) < 0) .or. & any(prt_params%alloc_priority(ipft,:) > 6) ) then + write(fates_log(),*) "---~---" write(fates_log(),*) ' PFT#: ',ipft write(fates_log(),*) ' Allocation priorities should be 0-6 for CNP flex hypothesis' write(fates_log(),*) prt_params%alloc_priority(ipft,:) - write(fates_log(),*) ' Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end select @@ -1437,35 +1585,45 @@ subroutine PRTCheckParams(is_master) ! Check that leaf turnover doesn't exeed 1 day if ( (years_per_day / prt_params%leaf_long(ipft,iage)) > 1._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'Leaf turnover time-scale is greater than 1 day!' write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage write(fates_log(),*) 'leaf_long(ipft,iage): ',prt_params%leaf_long(ipft,iage),' [years]' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if ! Check to make sure that all other age-classes for this PFT also ! have non-zero entries, it wouldn't make sense otherwise if ( any(prt_params%leaf_long(ipft,:) <= nearzero) ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'You specified a leaf_long that is zero or' write(fates_log(),*) 'invalid for a particular age class.' write(fates_log(),*) 'Yet, other age classes for this PFT are non-zero.' write(fates_log(),*) 'this doesnt make sense.' write(fates_log(),*) 'ipft = ',ipft write(fates_log(),*) 'leaf_long(ipft,:) = ',prt_params%leaf_long(ipft,:) - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if else if (prt_params%evergreen(ipft) .eq. itrue) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'You specified zero leaf turnover: ' write(fates_log(),*) 'ipft: ',ipft,' iage: ',iage write(fates_log(),*) 'leaf_long(ipft,iage): ',prt_params%leaf_long(ipft,iage) write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' write(fates_log(),*) 'that an evergreen would have leaf maintenance turnover' write(fates_log(),*) 'disable this error if you are ok with this' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1478,23 +1636,30 @@ subroutine PRTCheckParams(is_master) if ( (years_per_day / & (prt_params%leaf_long(ipft,nleafage) * & prt_params%senleaf_long_fdrought(ipft))) > 1._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'Drought-senescent turnover time-scale is greater than 1 day!' write(fates_log(),*) 'ipft: ',ipft write(fates_log(),*) 'leaf_long(ipft,nleafage)*senleaf_long_fdrought: ', & prt_params%leaf_long(ipft,nleafage)*prt_params%senleaf_long_fdrought(ipft),' [years]' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if if ( prt_params%senleaf_long_fdrought(ipft)1._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'senleaf_long_fdrought(ipft) must be greater than 0 ' write(fates_log(),*) 'or less than or equal to 1.' write(fates_log(),*) 'Set this to 1 if you want no accelerated senescence turnover' write(fates_log(),*) 'ipft = ',ipft write(fates_log(),*) 'senleaf_long_fdrought(ipft) = ',prt_params%senleaf_long_fdrought(ipft) - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if @@ -1502,22 +1667,29 @@ subroutine PRTCheckParams(is_master) ! Check that root turnover doesn't exeed 1 day if ( (years_per_day / prt_params%root_long(ipft)) > 1._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'Root turnover time-scale is greater than 1 day!' write(fates_log(),*) 'ipft: ',ipft write(fates_log(),*) 'root_long(ipft): ',prt_params%root_long(ipft),' [years]' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if else if (prt_params%evergreen(ipft) .eq. itrue) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'You specified zero root turnover: ' write(fates_log(),*) 'ipft: ',ipft write(fates_log(),*) 'root_long(ipft): ',prt_params%root_long(ipft) write(fates_log(),*) 'yet this is an evergreen PFT, and it only makes sense' write(fates_log(),*) 'that an evergreen would have root maintenance turnover' write(fates_log(),*) 'disable this error if you are ok with this' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1526,11 +1698,14 @@ subroutine PRTCheckParams(is_master) ! Check that branch turnover doesn't exeed 1 day if ( (years_per_day / prt_params%branch_long(ipft)) > 1._r8 ) then + write(fates_log(),*) "---~---" write(fates_log(),*) 'Branch turnover time-scale is greater than 1 day!' write(fates_log(),*) 'ipft: ',ipft write(fates_log(),*) 'branch_long(ipft): ',prt_params%branch_long(ipft),' [years]' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) + write(fates_log(),*) "---~---" + write(fates_log(),*) '' + write(fates_log(),*) '' + nerror = nerror + 1 end if end if @@ -1538,6 +1713,14 @@ subroutine PRTCheckParams(is_master) end do pftloop + ! If any error was found, abort. We add a single point to abort the run after all + ! checks so users can get all the errors and address them in one go (as opposed to + ! multiple submissions). + if (nerror > 0) then + write(fates_log(),*) 'One or more parameter errors found. Aborting.' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + return end subroutine PRTCheckParams