Skip to content

Change quadratic code solver #4255

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ekluzek merged 13 commits into from
Oct 15, 2026
+119 −100
Merged

Change quadratic code solver #4255

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
7501171
Change name of quadratic subroutine with a _roots suffix to make it's…
ekluzek Oct 9, 2026
3540c55
Remove mention of reference for quadratic_roots, since it's in it's o…
ekluzek Oct 9, 2026
16847de
Change names of temporary variables for clarity, show working in PF u…
ekluzek Oct 9, 2026
00d84d8
Send in line and file for quadratic_roots to give a message about the…
ekluzek Oct 9, 2026
b1e8f30
Merge remote-tracking branch 'escomp/b4b-dev' into change_code_recipes
ekluzek Oct 9, 2026
b39995a
Make the reference to the solution more clear, fix description of ban…
ekluzek Oct 9, 2026
151e5ca
Update src/biogeophys/TridiagonalMod.F90
ekluzek Oct 12, 2026
04281b3
Update src/biogeophys/TridiagonalMod.F90
ekluzek Oct 12, 2026
4c05c03
Update src/utils/quadraticMod.F90
ekluzek Oct 12, 2026
94ad439
Update src/utils/quadraticMod.F90
ekluzek Oct 12, 2026
c130f70
Add back the revision history in
ekluzek Oct 12, 2026
c5a5493
Merge branch 'b4b-dev' into change_code_recipes
ekluzek Oct 13, 2026
2aacc8e
Merge branch 'b4b-dev' into change_code_recipes
ekluzek Oct 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/biogeophys/BandDiagonalMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ module BandDiagonalMod
subroutine BandDiagonal(bounds, lbj, ubj, jtop, jbot, numf, filter, nband, b, r, u)
!
! !DESCRIPTION:
! Tridiagonal matrix solution
! Band-diagonal matrix solution from LAPACK
!
! !ARGUMENTS:
implicit none
Expand Down
42 changes: 5 additions & 37 deletions src/biogeophys/LunaMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,6 @@ module LunaMod
private :: VcmxTKattge !Calculate the temperature response for Vcmax, based on Kattge and Knorr 2007
private :: RespTBernacchi !Calculate the temperature response for respiration, following Bernacchi PCE 2001
private :: Photosynthesis_luna !calculate the photosynthetic rate for nitrogen allocation
private :: Quadratic !Calculate the soultion using the quadratic formula

!------------------------------------------------------------------------------
!Constants
Expand Down Expand Up @@ -214,7 +213,6 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
use clm_varpar , only : nlevsoi
use perf_mod , only : t_startf, t_stopf
use clm_varctl , only : use_cn
use quadraticMod , only : quadratic
use CNSharedParamsMod , only : CNParamsShareInst
use shr_infnan_mod, only : isnan => shr_infnan_isnan
use OzoneBaseMod, only : ozone_base_type
Expand Down Expand Up @@ -1091,6 +1089,7 @@ end subroutine Nitrogen_investments
! two phases. First phase is where Rubisco is limiting (Wc <= Wj) and second phase is where light is limiting (Wj > Wc).

subroutine Photosynthesis_luna(forc_pbot, tleafd, relh, CO2a,O2a, rb, Vcmax, JmeanL, ci, Kc, Kj, A)
use quadraticMod, only : quadratic_roots
implicit none
real(r8), intent (in) :: forc_pbot !air presure (Pa)
real(r8), intent (in) :: tleafd !daytime leaf temperature (oC)
Expand Down Expand Up @@ -1165,7 +1164,8 @@ subroutine Photosynthesis_luna(forc_pbot, tleafd, relh, CO2a,O2a, rb, Vcmax, Jme
phi = forc_pbot * (1.37_r8 * gs_mol + 1.6_r8 * gb_mol) / (gb_mol * gs_mol)
bquad = awc - CO2c + phi * Vcmax
cquad = -(c_p * phi * Vcmax + awc * CO2c)
call Quadratic(aquad, bquad, cquad, r1, r2)
call quadratic_roots(aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ci = max(r1, r2)
if (ci < 0.0_r8) ci = c_p + 0.5_r8 * ciold
end do
Expand All @@ -1183,7 +1183,8 @@ subroutine Photosynthesis_luna(forc_pbot, tleafd, relh, CO2a,O2a, rb, Vcmax, Jme
phi = forc_pbot * (1.37_r8 * gs_mol + 1.6_r8 * gb_mol) / (gb_mol * gs_mol)
bquad = 2.0_r8 * c_p - CO2c + phi * JmeanL / 4.0_r8
cquad = -(c_p * phi * JmeanL / 4.0_r8 + 2.0_r8 * c_p * CO2c)
call Quadratic(aquad, bquad, cquad, r1, r2)
call quadratic_roots(aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ci = max(r1, r2)
if (ci < 0.0_r8) ci = c_p + 0.5_r8 * ciold
Kj = max(ci - c_p, 0.0_r8) / (4.0_r8 * ci + 8.0_r8 * c_p)
Expand Down Expand Up @@ -1372,38 +1373,5 @@ real(r8) function RespTBernacchi(tleaf)
end function RespTBernacchi


!********************************************************************************************************************
!Calculate the soultion using the quadratic formula

subroutine Quadratic(a,b,c,r1,r2)
implicit none
real(r8), intent(in) :: a !coefficient a
real(r8), intent(in) :: b !coefficient b
real(r8), intent(in) :: c !coefficient c
real(r8), intent(out) :: r1 !root one
real(r8), intent(out) :: r2 !root one
real(r8) :: q ! temporary term for quadratic solution

r1 = 1.0e36_r8
r2 = 1.0e36_r8

if (a == 0.0_r8) return

if (b .GE. 0.0_r8) then
q = -0.5_r8 * (b + sqrt(b*b - 4.0_r8*a*c))
else
q = -0.5_r8 * (b - sqrt(b*b - 4.0_r8*a*c))
end if

r1 = q / a

if (q .NE. 0.0_r8)then
r2 = c / q
else
r2 = 1.0e36_r8
end if

end subroutine Quadratic

end module LunaMod

50 changes: 33 additions & 17 deletions src/biogeophys/PhotosynthesisMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ module PhotosynthesisMod
use clm_varpar , only : nlevcan, nvegwcs, mxpft
use clm_varcon , only : c14ratio, spval, isecspday
use decompMod , only : bounds_type, subgrid_level_patch
use QuadraticMod , only : quadratic
use QuadraticMod , only : quadratic_roots
use pftconMod , only : pftcon
use atm2lndType , only : atm2lnd_type
use CanopyStateType , only : canopystate_type
Expand Down Expand Up @@ -1898,7 +1898,8 @@ subroutine Photosynthesis ( bounds, fn, filterp, &
aquad = params_inst%theta_psii
bquad = -(qabs + jmax_z(p,iv))
cquad = qabs * jmax_z(p,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
je = min(r1,r2)

! Iterative loop for ci beginning with initial guess
Expand Down Expand Up @@ -2375,8 +2376,9 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,&
!!DESCRIPTION:
!Use Brent's method to find the root of a single variable function ci_func, which is known to exist between x1 and x2.
!The found root will be updated until its accuracy is tol.

!
!!REVISION HISTORY:
Comment thread
ekluzek marked this conversation as resolved.
!
!Dec 14/2012: Jinyun Tang, modified from numerical recipes in F90 by press et al. 1188-1189
!
!!ARGUMENTS:
Expand Down Expand Up @@ -2650,13 +2652,15 @@ subroutine ci_func(ci, fval, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,&
aquad = params_inst%theta_cj(ivt(p))
bquad = -(ac(p,iv) + aj(p,iv))
cquad = ac(p,iv) * aj(p,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ai = min(r1,r2)

aquad = params_inst%theta_ip
bquad = -(ai + ap(p,iv))
cquad = ai * ap(p,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ag(p,iv) = max(0._r8,min(r1,r2))

! Net photosynthesis. Exit iteration if an < 0
Expand All @@ -2679,13 +2683,15 @@ subroutine ci_func(ci, fval, p, iv, c, gb_mol, je, cair, oair, lmr_z, par_z,&
(2.0*medlynintercept(patch%itype(p))*1.e-06_r8 + term * &
(1.0 - medlynslope(patch%itype(p))* medlynslope(patch%itype(p)) / rh_can)) * term

call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol = max(r1,r2) * 1.e06_r8
else if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 )then
aquad = cs
bquad = cs*(gb_mol - bbb(p)) - mbb(p)*an(p,iv)*forc_pbot(c)
cquad = -gb_mol*(cs*bbb(p) + mbb(p)*an(p,iv)*forc_pbot(c)*rh_can)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol = max(r1,r2)
end if

Expand Down Expand Up @@ -3566,15 +3572,17 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, &
aquad = params_inst%theta_psii
bquad = -(qabs + jmax_z(p,sun,iv))
cquad = qabs * jmax_z(p,sun,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
je_sun = min(r1,r2)

! sha
qabs = 0.5_r8 * (1._r8 - params_inst%fnps) * par_z_sha(p,iv) * 4.6_r8
aquad = params_inst%theta_psii
bquad = -(qabs + jmax_z(p,sha,iv))
cquad = qabs * jmax_z(p,sha,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
je_sha = min(r1,r2)

! Iterative loop for ci beginning with initial guess
Expand Down Expand Up @@ -4344,26 +4352,30 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
aquad = params_inst%theta_cj(ivt(p))
bquad = -(ac(p,sun,iv) + aj(p,sun,iv))
cquad = ac(p,sun,iv) * aj(p,sun,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ai = min(r1,r2)

aquad = params_inst%theta_ip
bquad = -(ai + ap(p,sun,iv))
cquad = ai * ap(p,sun,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ag(p,sun,iv) = max(0._r8,min(r1,r2))

! Shaded
aquad = params_inst%theta_cj(ivt(p))
bquad = -(ac(p,sha,iv) + aj(p,sha,iv))
cquad = ac(p,sha,iv) * aj(p,sha,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ai = min(r1,r2)

aquad = params_inst%theta_ip
bquad = -(ai + ap(p,sha,iv))
cquad = ai * ap(p,sha,iv)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
ag(p,sha,iv) = max(0._r8,min(r1,r2))

! Net photosynthesis. Exit iteration if an < 0
Expand Down Expand Up @@ -4421,7 +4433,8 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
(2.0_r8*medlynintercept(patch%itype(p))*1.e-06_r8 + term * &
(1.0_r8 - medlynslope(patch%itype(p))* medlynslope(patch%itype(p)) / rh_can)) * term

call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol_sun = max(r1,r2) * 1.e06_r8
end if

Expand All @@ -4438,7 +4451,8 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
(2.0_r8*medlynintercept(patch%itype(p))*1.e-06_r8 + term * (1.0 - medlynslope(patch%itype(p))* &
medlynslope(patch%itype(p)) / rh_can)) * term

call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol_sha = max(r1,r2)* 1.e06_r8
end if
else if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 )then
Expand All @@ -4450,7 +4464,8 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
aquad = cs_sun
bquad = cs_sun*(gb_mol - max(bsun*bbb(p),1._r8)) - mbb(p)*an_sun(p,iv)*forc_pbot(c)
cquad = -gb_mol*(cs_sun*max(bsun*bbb(p),1._r8) + mbb(p)*an_sun(p,iv)*forc_pbot(c)*rh_can)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol_sun = max(r1,r2)
end if

Expand All @@ -4462,7 +4477,8 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b
aquad = cs_sha
bquad = cs_sha*(gb_mol - max(bsha*bbb(p),1._r8)) - mbb(p)*an_sha(p,iv)*forc_pbot(c)
cquad = -gb_mol*(cs_sha*max(bsha*bbb(p),1._r8) + mbb(p)*an_sha(p,iv)*forc_pbot(c)*rh_can)
call quadratic (aquad, bquad, cquad, r1, r2)
call quadratic_roots (aquad, bquad, cquad, r1, r2, &
file=sourcefile, line=__LINE__)
gs_mol_sha = max(r1,r2)
end if
end if
Expand Down
7 changes: 5 additions & 2 deletions src/biogeophys/TridiagonalMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,16 +24,19 @@ subroutine Tridiagonal (bounds, lbj, ubj, jtop, numf, filter, a, b, c, r, u)
!
! !DESCRIPTION:
! Tridiagonal matrix solution
! This implements the Crank-Nicholson method for solving a tridiagonal matrix system of equations, customized to CLM Fortran structures.
! The method is based on Press, et. al. 1992. Numerical Recipes in FORTRAN:
! The Art of Scientific Computing. Cambridge University Press, New York.
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use clm_varctl , only : iulog
use decompMod , only : bounds_type
use ColumnType , only : col
use ColumnType , only : col
!
! !ARGUMENTS:
implicit none
type(bounds_type), intent(in) :: bounds
type(bounds_type), intent(in) :: bounds
integer , intent(in) :: lbj, ubj ! lbinning and ubing level indices
integer , intent(in) :: jtop( bounds%begc: ) ! top level for each column [col]
integer , intent(in) :: numf ! filter dimension (should not include hydrologically inactive points)
Expand Down
1 change: 0 additions & 1 deletion src/utils/clmfates_interfaceMod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2646,7 +2646,6 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, &
use clm_varcon , only : tfrz, namep
use clm_varctl , only : iulog
use PatchType , only : patch
use quadraticMod , only : quadratic
use EDtypesMod , only : ed_site_type
use FatesPatchMod, only : fates_patch_type
use FatesCohortMod , only : fates_cohort_type
Expand Down
Loading
Loading