Skip to content

Commit

Permalink
Merge pull request #11 from rgknox/mpaiao-pr-allom-modes-moreparams
Browse files Browse the repository at this point in the history 
added parameters and merged up
  • Loading branch information
@mpaiao
mpaiao authored Mar 27, 2024
2 parents 37db010 + 97dc27a commit 6573572
Show file tree
Hide file tree
Showing 19 changed files with 9,612 additions and 6,775 deletions.
11 changes: 6 additions & 5 deletions biogeophys/FatesPlantHydraulicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -937,6 +937,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
real(r8), parameter :: min_trim = 0.1_r8 ! The lower cap on trimming function used
! to estimate maximum leaf carbon


ccohort_hydr => ccohort%co_hydr
ft = ccohort%pft
nlevrhiz = csite_hydr%nlevrhiz
Expand Down Expand Up @@ -980,13 +981,13 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
! Get the target, or rather, maximum leaf carrying capacity of plant
! Lets also avoid super-low targets that have very low trimming functions

! efleaf_coh hard-coded to 1 in the call below to avoid zero leaf volume
call bleaf(ccohort%dbh,ccohort%pft,ccohort%crowndamage, &
max(ccohort%canopy_trim,min_trim),ccohort%efleaf_coh, leaf_c_target)
max(ccohort%canopy_trim,min_trim),1.0_r8, leaf_c_target)

ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * &
prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8)

if( (ccohort%status_coh == leaves_on) .or. ccohort_hydr%is_newly_recruited ) then
ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * &
prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8)
end if

! Step sapwood volume
! -----------------------------------------------------------------------------------
Expand Down
152 changes: 7 additions & 145 deletions biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ module FATESPlantRespPhotosynthMod
use FatesInterfaceTypesMod, only : hlm_parteh_mode
use FatesInterfaceTypesMod, only : numpft
use FatesInterfaceTypesMod, only : nleafage
use FatesUtilsMod, only : QuadraticRoots => QuadraticRootsSridharachary
use EDParamsMod, only : maxpft
use EDParamsMod, only : nlevleaf
use EDParamsMod, only : nclmax
Expand Down Expand Up @@ -1379,7 +1380,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
aquad = theta_psii
bquad = -(qabs + jmax)
cquad = qabs * jmax
call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
je = min(r1,r2)

! Initialize intercellular co2
Expand Down Expand Up @@ -1409,7 +1410,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
aquad = theta_cj_c3
bquad = -(ac + aj)
cquad = ac * aj
call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
agross = min(r1,r2)

else
Expand Down Expand Up @@ -1440,13 +1441,13 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
aquad = theta_cj_c4
bquad = -(ac + aj)
cquad = ac * aj
call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
ai = min(r1,r2)

aquad = theta_ip
bquad = -(ai + ap)
cquad = ai * ap
call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
agross = min(r1,r2)

end if
Expand Down Expand Up @@ -1485,7 +1486,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
(2.0*stomatal_intercept_btran + term * &
(1.0 - medlyn_slope(ft)* medlyn_slope(ft) / vpd)) * term

call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)

else if ( stomatal_model == ballberry_model ) then !stomatal conductance calculated from Ball et al. (1987)
Expand All @@ -1494,7 +1495,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in
cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + &
bb_slope(ft)*anet*can_press * ceair/ veg_esat )

call quadratic_f (aquad, bquad, cquad, r1, r2)
call QuadraticRoots(aquad, bquad, cquad, r1, r2)
gs_mol = max(r1,r2)
end if

Expand Down Expand Up @@ -1922,145 +1923,6 @@ end function fth25_f

! =====================================================================================

subroutine quadratic_f (a, b, c, r1, r2)
!
! !DESCRIPTION:
!==============================================================================!
!----------------- Solve quadratic equation for its two roots -----------------!
!==============================================================================!
! This solution is mostly derived from:
! Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 1992. Numerical Recipes
! in Fortran77: The Art of Scientific Computing. 2nd edn. Cambridge
! University Press, Cambridge UK, ISBN 0-521-43064-X.
! Available at: http://numerical.recipes/oldverswitcher.html, section 5.6.
!
! !REVISION HISTORY:
! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson
! 7/23/16: Copied over from CLM by Ryan Knox
! 12/30/23: Instead of issuing errors when a=0, solve the trivial cases too.
! Check determinant sign, and stop the run when it is negative.
!
! !USES:
!
! !ARGUMENTS:
real(r8), intent(in) :: a,b,c ! Terms for quadratic equation
real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation
!
! !LOCAL VARIABLES:
real(r8) :: discriminant ! Discriminant
real(r8) :: q ! Temporary term for quadratic solution
logical :: a_offzero ! Is a close to zero?
logical :: b_offzero ! Is b close to zero?
logical :: c_offzero ! Is c close to zero?
! ! Local constants:
real(r8), parameter :: discard = 1.e36_r8 ! Large number for discarding answer
!------------------------------------------------------------------------------

! Save logical tests.
a_offzero = abs(a) > nearzero
b_offzero = abs(b) > nearzero
c_offzero = abs(c) > nearzero

if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then
! Quadratic equation with two non-zero solutions (but may be complex solutions)
discriminant = b*b - 4._r8 * a * c

! Proceed only when the discriminant is non-negative or only tiny negative
if (discriminant >= - nearzero) then
! Coerce discriminant to non-negative
discriminant = max(0._r8,discriminant)

! Find q as in the numerical recipes. If b or c are non-zero, q cannot
! be zero, no need for additional checks.
q = - 0.5_r8 * (b + sign(sqrt(discriminant),b))
r1 = q / a
r2 = c / q
else
! Negative discriminant, stop the run.
write (fates_log(),'(a)') '---~---'
write (fates_log(),'(a)') ' Fatal error!'
write (fates_log(),'(a)') ' Quadratic equation discriminant is negative.'
write (fates_log(),'(a)') '---~---'
write (fates_log(),'(a,1x,es12.5)') ' a = ',a
write (fates_log(),'(a,1x,es12.5)') ' b = ',b
write (fates_log(),'(a,1x,es12.5)') ' c = ',c
write (fates_log(),'(a,1x,es12.5)') ' discriminant = ',discriminant
call endrun(msg=errMsg(sourcefile, __LINE__))
end if
else if (a_offzero) then
! b and c are nearly zero. Both roots must be zero.
r1 = 0._r8
r2 = 0._r8
else if (b_offzero) then
! "a" is not zero, not a true quadratic equation. Single root.
r1 = - c / b
r2 = discard
else
! Both a and b are zero, this really doesn't make any sense and should never
! happen. If it does, issue an error and stop the run.
write (fates_log(),'(a)') '---~---'
write (fates_log(),'(a)') ' Fatal error!'
write (fates_log(),'(a)') ' This solver expects ''a'' and/or ''b'' to be non-zero.'
write (fates_log(),'(a)') '---~---'
write (fates_log(),'(a,1x,es12.5)') ' a = ',a
write (fates_log(),'(a,1x,es12.5)') ' b = ',b
write (fates_log(),'(a,1x,es12.5)') ' c = ',c
write (fates_log(),'(a)') '---~---'
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

return
end subroutine quadratic_f

! ====================================================================================

subroutine quadratic_fast (a, b, c, r1, r2)
!
! !DESCRIPTION:
!==============================================================================!
!----------------- Solve quadratic equation for its two roots -----------------!
! THIS METHOD SIMPLY REMOVES THE DIV0 CHECK AND ERROR REPORTING !
!==============================================================================!
! Solution from Press et al (1986) Numerical Recipes: The Art of Scientific
! Computing (Cambridge University Press, Cambridge), pp. 145.
!
! !REVISION HISTORY:
! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson
! 7/23/16: Copied over from CLM by Ryan Knox
!
! !USES:
!
! !ARGUMENTS:
real(r8), intent(in) :: a,b,c ! Terms for quadratic equation
real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation
!
! !LOCAL VARIABLES:
real(r8) :: q ! Temporary term for quadratic solution
!------------------------------------------------------------------------------

! if (a == 0._r8) then
! write (fates_log(),*) 'Quadratic solution error: a = ',a
! call endrun(msg=errMsg(sourcefile, __LINE__))
! end if

if (b >= 0._r8) then
q = -0.5_r8 * (b + sqrt(b*b - 4._r8*a*c))
else
q = -0.5_r8 * (b - sqrt(b*b - 4._r8*a*c))
end if

r1 = q / a
! if (q /= 0._r8) then
r2 = c / q
! else
! r2 = 1.e36_r8
! end if

end subroutine quadratic_fast


! ====================================================================================

subroutine UpdateCanopyNCanNRadPresent(currentPatch)

! ---------------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 6573572

Please sign in to comment.