From 7501171612409419a6172a0f622d2f8b69f68bd3 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 9 Oct 2026 11:56:39 -0600 Subject: [PATCH 01/10] Change name of quadratic subroutine with a _roots suffix to make it's purpose more clear, add some comments, add option to send in file and line number on error, get the unit test working with the changes --- src/biogeophys/LunaMod.F90 | 40 +--------- src/biogeophys/PhotosynthesisMod.F90 | 32 ++++---- src/utils/clmfates_interfaceMod.F90 | 1 - src/utils/quadraticMod.F90 | 75 +++++++++++++------ .../test/quadratic_test/test_quadratic.pf | 21 +++--- 5 files changed, 83 insertions(+), 86 deletions(-) diff --git a/src/biogeophys/LunaMod.F90 b/src/biogeophys/LunaMod.F90 index 950fefbc32..3f662371af 100644 --- a/src/biogeophys/LunaMod.F90 +++ b/src/biogeophys/LunaMod.F90 @@ -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 @@ -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 @@ -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) @@ -1165,7 +1164,7 @@ 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) ci = max(r1, r2) if (ci < 0.0_r8) ci = c_p + 0.5_r8 * ciold end do @@ -1183,7 +1182,7 @@ 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) 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) @@ -1372,38 +1371,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 diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 2f96a20d5d..b2a287b2d7 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -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 @@ -1898,7 +1898,7 @@ 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) je = min(r1,r2) ! Iterative loop for ci beginning with initial guess @@ -2650,13 +2650,13 @@ 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) 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) ag(p,iv) = max(0._r8,min(r1,r2)) ! Net photosynthesis. Exit iteration if an < 0 @@ -2679,13 +2679,13 @@ 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) 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) gs_mol = max(r1,r2) end if @@ -3566,7 +3566,7 @@ 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) je_sun = min(r1,r2) ! sha @@ -3574,7 +3574,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & 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) je_sha = min(r1,r2) ! Iterative loop for ci beginning with initial guess @@ -4344,26 +4344,26 @@ 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) 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) 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) 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) ag(p,sha,iv) = max(0._r8,min(r1,r2)) ! Net photosynthesis. Exit iteration if an < 0 @@ -4421,7 +4421,7 @@ 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) gs_mol_sun = max(r1,r2) * 1.e06_r8 end if @@ -4438,7 +4438,7 @@ 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) gs_mol_sha = max(r1,r2)* 1.e06_r8 end if else if ( stomatalcond_mtd == stomatalcond_mtd_bb1987 )then @@ -4450,7 +4450,7 @@ 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) gs_mol_sun = max(r1,r2) end if @@ -4462,7 +4462,7 @@ 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) gs_mol_sha = max(r1,r2) end if end if diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index 3f9a662c02..35c4bb02a3 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -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 diff --git a/src/utils/quadraticMod.F90 b/src/utils/quadraticMod.F90 index 87bb250945..ebcfb7c1b5 100644 --- a/src/utils/quadraticMod.F90 +++ b/src/utils/quadraticMod.F90 @@ -1,76 +1,105 @@ module quadraticMod use abortutils , only: endrun - use shr_kind_mod, only: r8 => shr_kind_r8 - use shr_log_mod , only: errMsg => shr_log_errMsg + use shr_kind_mod, only: r8 => shr_kind_r8, CX => shr_kind_cx use clm_varctl , only: iulog implicit none - public :: quadratic + private + + public :: quadratic_roots ! Solve for the two roots of a quadratic equation character(len=*), parameter, private :: sourcefile = & __FILE__ contains - subroutine quadratic (a, b, c, r1, r2) + subroutine quadratic_roots (a, b, c, root1, root2, file, line) ! ! !DESCRIPTION: !==============================================================================! !----------------- Solve quadratic equation for its two roots -----------------! !==============================================================================! - ! Solution from Press et al (1986) Numerical Recipes: The Art of Scientific - ! Computing (Cambridge University Press, Cambridge), pp. 145. + ! Implements the numerically stable formulation from Press et al (1986) + ! Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge) + ! + ! NOTE: Special handling for these cases... + ! Will truncate the square root term to zero if it is very small and negative, otherwise will error out + ! Root2 will be set to spval if it would be undefined by a division by zero ! ! !REVISION HISTORY: - ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson + ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson + ! 5/14/18: Modify endrun handling and handle truncation to zero of small negative square root term EBK + ! 10/9/26: Clarify variable names, add comments and pass in file and line number of calling routine EBK ! ! !USES: + use clm_varcon, only: spval + implicit none ! ! !ARGUMENTS: - real(r8), intent(in) :: a,b,c ! Terms for quadratic equation - real(r8), intent(out) :: r1,r2 ! Roots of quadratic equation - ! + real(r8), intent(in) :: a, b, c ! Coefficients of the quadratic equation x = a*x^2 + b*x + c + real(r8), intent(out) :: root1, root2 ! The two roots of the quadratic equation + character(len=*), intent(in), optional :: file ! File name of calling routine + integer, intent(in), optional :: line ! Line number of calling routine ! !LOCAL VARIABLES: - real(r8) :: q ! Temporary term for quadratic solution - real(r8) :: root ! Term that will have a square root taken - character(len=*), parameter :: subname = 'quadratic' + real(r8) :: q ! Temporary term for the solution + real(r8) :: root ! Temporary Term that will have a square root taken + integer :: pline ! Line number to use in endrun calls + character(len=CX) :: pfile ! File name to use in endrun calls !------------------------------------------------------------------------------ + + if ( present(file) )then + pfile = file + else + pfile = sourcefile + end if + if ( present(line) )then + pline = line + else + pline = __LINE__ + end if + + ! Initialize the roots to spval in case it exits early due to an error + root1 = spval + root2 = spval + ! If the "a" coefficient is zero, then this is linear rather than quadratic which we assume is a mistake if (a == 0._r8) then - write (iulog,*) subname//' ERROR: Quadratic solution error: a = ',a - write (iulog,*) errmsg(sourcefile, __LINE__) - call endrun(msg=subname//' ERROR: Quadratic solution error' ) + write (iulog,*) ' ERROR: Quadratic solution error: the "a" coefficient is zero = ',a + call endrun(msg='Quadratic solution error', file=pfile, line=pline) return end if + ! Compute the term that will have a square root taken root = b*b - 4._r8*a*c if ( root < 0.0 )then + ! Truncate to zero if the term is very small and negative (relative to b), otherwise error out if ( -root < 3.0_r8*epsilon(b) )then root = 0.0_r8 else - write (iulog,*) subname//' ERROR: Quadratic solution error: b^2 - 4ac is negative = ', root - write (iulog,*) errmsg(sourcefile, __LINE__) - call endrun( msg=subname//' ERROR: Quadratic solution error: b^2 - 4ac is negative' ) + write (iulog,*) ' ERROR: Quadratic solution error: root would be complex b^2 - 4*a*c = ', root + call endrun( msg=' ERROR: Quadratic solution error: root would be complex', file=pfile, line=pline ) return end if end if + ! Compute the term that will be the first root multiplied by the "a" coefficient if (b >= 0._r8) then q = -0.5_r8 * (b + sqrt(root)) else q = -0.5_r8 * (b - sqrt(root)) end if - r1 = q / a + ! Solve for the two roots to return + root1 = q / a if (q /= 0._r8) then - r2 = c / q + root2 = c / q else - r2 = 1.e36_r8 + root2 = spval end if - end subroutine quadratic + end subroutine quadratic_roots end module quadraticMod diff --git a/src/utils/test/quadratic_test/test_quadratic.pf b/src/utils/test/quadratic_test/test_quadratic.pf index 7acfbae6e9..651749e686 100644 --- a/src/utils/test/quadratic_test/test_quadratic.pf +++ b/src/utils/test/quadratic_test/test_quadratic.pf @@ -58,7 +58,7 @@ contains b = 0.0_r8 c = 0.0_r8 do i = 1, size(a) - call quadratic (a(i), b, c, r1, r2) + call quadratic_roots (a(i), b, c, r1, r2) call check_root(a(i), b, c, r1) @assertEqual(0.0_r8, r1, tolerance=tol) @assertEqual(1.e+36_r8, r2, tolerance=tol) @@ -84,7 +84,7 @@ contains expectedr1 = (/ -2.0_r8, -2.0_r8 /) expectedr2 = (/ -1.0_r8, +2.0_r8 /) do i = 1, size(a) - call quadratic (a(i), b(i), c(i), r1, r2) + call quadratic_roots (a(i), b(i), c(i), r1, r2) call check_root(a(i), b(i), c(i), r1) call check_root(a(i), b(i), c(i), r2) @assertEqual(expectedr1(i), r1, tolerance=tol) @@ -117,7 +117,7 @@ contains do j = 1, size(pert) c_coef = c(i)+pert(j) print *, 'a, b, c_ = ', a(i), b(i), c_coef - call quadratic (a(i), b(i), c_coef, r1, r2) + call quadratic_roots (a(i), b(i), c_coef, r1, r2) call check_root(a(i), b(i), c_coef, r1) call check_root(a(i), b(i), c_coef, r2) @assertEqual(expectedr1(i), r1, tolerance=tol_big) @@ -142,9 +142,10 @@ contains a = 0.0_r8 b = a c = a - call quadratic (a, b, c, r1, r2) + call quadratic_roots (a, b, c, r1, r2, & + file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - 'quadratic ERROR: Quadratic solution error') + 'Quadratic solution error') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsazero @@ -159,9 +160,10 @@ contains a = 1.0_r8 b = 2.0_r8 c = 5.0_r8 - call quadratic (a, b, c, r1, r2) + call quadratic_roots (a, b, c, r1, r2, & + file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - 'quadratic ERROR: Quadratic solution error: b^2 - 4ac is negative') + ' ERROR: Quadratic solution error: root would be complex') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsimaginary @@ -175,9 +177,10 @@ contains a = 1.0_r8 b = 4.0_r8 c = 4.0_r8 + 100.0_r8*epsilon(b) - call quadratic (a, b, c, r1, r2) + call quadratic_roots (a, b, c, r1, r2, & + file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - 'quadratic ERROR: Quadratic solution error: b^2 - 4ac is negative') + ' ERROR: Quadratic solution error: root would be complex') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsbarelyimaginary From 3540c55998e8d64827ccf1d356a1c151a0d1f551 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 9 Oct 2026 12:00:41 -0600 Subject: [PATCH 02/10] Remove mention of reference for quadratic_roots, since it's in it's own module now --- src/biogeophys/PhotosynthesisMod.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index b2a287b2d7..e783b004a4 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -2376,9 +2376,6 @@ subroutine brent(x, x1,x2,f1, f2, tol, ip, iv, ic, gb_mol, je, cair, oair,& !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: - !Dec 14/2012: Jinyun Tang, modified from numerical recipes in F90 by press et al. 1188-1189 - ! !!ARGUMENTS: real(r8), intent(out) :: x ! indepedent variable of the single value function ci_func(x) real(r8), intent(in) :: x1, x2, f1, f2 ! minimum and maximum of the variable domain to search for the solution ci_func(x1) = f1, ci_func(x2)=f2 From 16847deed2a675107b4bba2b0b346542d32a5c00 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 9 Oct 2026 12:31:27 -0600 Subject: [PATCH 03/10] Change names of temporary variables for clarity, show working in PF unit tests --- src/utils/quadraticMod.F90 | 38 +++++++++---------- .../test/quadratic_test/test_quadratic.pf | 6 +-- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/src/utils/quadraticMod.F90 b/src/utils/quadraticMod.F90 index ebcfb7c1b5..4be725f3c6 100644 --- a/src/utils/quadraticMod.F90 +++ b/src/utils/quadraticMod.F90 @@ -21,7 +21,7 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) !==============================================================================! !----------------- Solve quadratic equation for its two roots -----------------! !==============================================================================! - ! Implements the numerically stable formulation from Press et al (1986) + ! Implements the numerically stable formulation from Press et al (1986) ! Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge) ! ! NOTE: Special handling for these cases... @@ -44,8 +44,8 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) character(len=*), intent(in), optional :: file ! File name of calling routine integer, intent(in), optional :: line ! Line number of calling routine ! !LOCAL VARIABLES: - real(r8) :: q ! Temporary term for the solution - real(r8) :: root ! Temporary Term that will have a square root taken + real(r8) :: root1_x_a ! Temporary term for the solution first root multiplied by the "a" coefficient + real(r8) :: discriminant ! Temporary Term that will have a square root taken integer :: pline ! Line number to use in endrun calls character(len=CX) :: pfile ! File name to use in endrun calls !------------------------------------------------------------------------------ @@ -64,38 +64,38 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) ! Initialize the roots to spval in case it exits early due to an error root1 = spval root2 = spval - + ! If the "a" coefficient is zero, then this is linear rather than quadratic which we assume is a mistake if (a == 0._r8) then - write (iulog,*) ' ERROR: Quadratic solution error: the "a" coefficient is zero = ',a - call endrun(msg='Quadratic solution error', file=pfile, line=pline) + write (iulog,*) 'ERROR: Quadratic solution error: the "a" coefficient is zero = ',a + call endrun(msg='Quadratic solution error: input equation is linear not quadratic', file=pfile, line=pline) return end if ! Compute the term that will have a square root taken - root = b*b - 4._r8*a*c - if ( root < 0.0 )then + discriminant = b*b - 4._r8*a*c + if ( discriminant < 0.0 )then ! Truncate to zero if the term is very small and negative (relative to b), otherwise error out - if ( -root < 3.0_r8*epsilon(b) )then - root = 0.0_r8 + if ( -discriminant < 3.0_r8*epsilon(b) )then + discriminant = 0.0_r8 else - write (iulog,*) ' ERROR: Quadratic solution error: root would be complex b^2 - 4*a*c = ', root - call endrun( msg=' ERROR: Quadratic solution error: root would be complex', file=pfile, line=pline ) + write (iulog,*) 'ERROR: Quadratic solution error: root would be complex b^2 - 4*a*c = ', discriminant + call endrun( msg='Quadratic solution error: root would be complex', file=pfile, line=pline ) return end if end if - + ! Compute the term that will be the first root multiplied by the "a" coefficient if (b >= 0._r8) then - q = -0.5_r8 * (b + sqrt(root)) + root1_x_a = -0.5_r8 * (b + sqrt(discriminant)) else - q = -0.5_r8 * (b - sqrt(root)) + root1_x_a = -0.5_r8 * (b - sqrt(discriminant)) end if - + ! Solve for the two roots to return - root1 = q / a - if (q /= 0._r8) then - root2 = c / q + root1 = root1_x_a / a + if (root1_x_a /= 0._r8) then + root2 = c / root1_x_a else root2 = spval end if diff --git a/src/utils/test/quadratic_test/test_quadratic.pf b/src/utils/test/quadratic_test/test_quadratic.pf index 651749e686..1042aae3c7 100644 --- a/src/utils/test/quadratic_test/test_quadratic.pf +++ b/src/utils/test/quadratic_test/test_quadratic.pf @@ -145,7 +145,7 @@ contains call quadratic_roots (a, b, c, r1, r2, & file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - 'Quadratic solution error') + 'Quadratic solution error: input equation is linear not quadratic') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsazero @@ -163,7 +163,7 @@ contains call quadratic_roots (a, b, c, r1, r2, & file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - ' ERROR: Quadratic solution error: root would be complex') + 'Quadratic solution error: root would be complex') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsimaginary @@ -180,7 +180,7 @@ contains call quadratic_roots (a, b, c, r1, r2, & file=__FILE__, line=__LINE__) expected_msg = endrun_msg( & - ' ERROR: Quadratic solution error: root would be complex') + 'Quadratic solution error: root would be complex') @assertExceptionRaised(expected_msg) end subroutine check_errorcondsbarelyimaginary From 00d84d87d216976323cc1ee328beabace66d6411 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 9 Oct 2026 12:50:07 -0600 Subject: [PATCH 04/10] Send in line and file for quadratic_roots to give a message about the caller if endrun happens --- src/biogeophys/LunaMod.F90 | 6 ++-- src/biogeophys/PhotosynthesisMod.F90 | 45 ++++++++++++++++++---------- src/utils/quadraticMod.F90 | 1 + 3 files changed, 35 insertions(+), 17 deletions(-) diff --git a/src/biogeophys/LunaMod.F90 b/src/biogeophys/LunaMod.F90 index 3f662371af..755025149a 100644 --- a/src/biogeophys/LunaMod.F90 +++ b/src/biogeophys/LunaMod.F90 @@ -1164,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_roots(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 @@ -1182,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_roots(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) diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index e783b004a4..564dd145b2 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -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_roots (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 @@ -2647,13 +2648,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_roots (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_roots (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 @@ -2676,13 +2679,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_roots (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_roots (aquad, bquad, cquad, r1, r2) + call quadratic_roots (aquad, bquad, cquad, r1, r2, & + file=sourcefile, line=__LINE__) gs_mol = max(r1,r2) end if @@ -3563,7 +3568,8 @@ 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_roots (aquad, bquad, cquad, r1, r2) + call quadratic_roots (aquad, bquad, cquad, r1, r2, & + file=sourcefile, line=__LINE__) je_sun = min(r1,r2) ! sha @@ -3571,7 +3577,8 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & aquad = params_inst%theta_psii bquad = -(qabs + jmax_z(p,sha,iv)) cquad = qabs * jmax_z(p,sha,iv) - call quadratic_roots (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 @@ -4341,26 +4348,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_roots (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_roots (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_roots (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_roots (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 @@ -4418,7 +4429,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_roots (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 @@ -4435,7 +4447,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_roots (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 @@ -4447,7 +4460,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_roots (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 @@ -4459,7 +4473,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_roots (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 diff --git a/src/utils/quadraticMod.F90 b/src/utils/quadraticMod.F90 index 4be725f3c6..ec8f838cfb 100644 --- a/src/utils/quadraticMod.F90 +++ b/src/utils/quadraticMod.F90 @@ -23,6 +23,7 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) !==============================================================================! ! Implements the numerically stable formulation from Press et al (1986) ! Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge) + ! Allows for roots that are technically complex if they are close to rounding to zero. ! ! NOTE: Special handling for these cases... ! Will truncate the square root term to zero if it is very small and negative, otherwise will error out From b39995a44f447729119f4d676d92861679a088e1 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Fri, 9 Oct 2026 14:46:49 -0600 Subject: [PATCH 05/10] Make the reference to the solution more clear, fix description of banddiagonal, and remove whitespace after a few lines --- src/biogeophys/BandDiagonalMod.F90 | 2 +- src/biogeophys/TridiagonalMod.F90 | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/biogeophys/BandDiagonalMod.F90 b/src/biogeophys/BandDiagonalMod.F90 index 13529a21ef..7acaf7ad0e 100644 --- a/src/biogeophys/BandDiagonalMod.F90 +++ b/src/biogeophys/BandDiagonalMod.F90 @@ -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 diff --git a/src/biogeophys/TridiagonalMod.F90 b/src/biogeophys/TridiagonalMod.F90 index 46532b0d8b..40b2326ade 100644 --- a/src/biogeophys/TridiagonalMod.F90 +++ b/src/biogeophys/TridiagonalMod.F90 @@ -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. + ! The method is from 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) From 151e5ca2d1ebf3d1c9689e3ce5a7c49e4e5ffadd Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 12 Oct 2026 13:48:34 -0600 Subject: [PATCH 06/10] Update src/biogeophys/TridiagonalMod.F90 Co-authored-by: Samuel Levis --- src/biogeophys/TridiagonalMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biogeophys/TridiagonalMod.F90 b/src/biogeophys/TridiagonalMod.F90 index 40b2326ade..af2ed54104 100644 --- a/src/biogeophys/TridiagonalMod.F90 +++ b/src/biogeophys/TridiagonalMod.F90 @@ -24,7 +24,7 @@ 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. + ! This implements the Crank-Nicholson method for solving a tridiagonal matrix system of equations, customized to CLM Fortran structures. ! The method is from Press, et. al. 1992. Numerical Recipes in FORTRAN: ! The Art of Scientific Computing. Cambridge University Press, New York. ! From 04281b396af74a1f724bd3dbd9dafe8abc4e7fd4 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 12 Oct 2026 13:48:48 -0600 Subject: [PATCH 07/10] Update src/biogeophys/TridiagonalMod.F90 Co-authored-by: Samuel Levis --- src/biogeophys/TridiagonalMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/biogeophys/TridiagonalMod.F90 b/src/biogeophys/TridiagonalMod.F90 index af2ed54104..4d1eb1dbb9 100644 --- a/src/biogeophys/TridiagonalMod.F90 +++ b/src/biogeophys/TridiagonalMod.F90 @@ -25,7 +25,7 @@ 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 from Press, et. al. 1992. Numerical Recipes in FORTRAN: + ! The method is based on Press, et. al. 1992. Numerical Recipes in FORTRAN: ! The Art of Scientific Computing. Cambridge University Press, New York. ! ! !USES: From 4c05c03667efa61d29f783c065ba40eefa45bd7e Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 12 Oct 2026 13:50:14 -0600 Subject: [PATCH 08/10] Update src/utils/quadraticMod.F90 I might revert this for going on master if it changes answers. Co-authored-by: Samuel Levis --- src/utils/quadraticMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/quadraticMod.F90 b/src/utils/quadraticMod.F90 index ec8f838cfb..1eb2cd43c8 100644 --- a/src/utils/quadraticMod.F90 +++ b/src/utils/quadraticMod.F90 @@ -75,7 +75,7 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) ! Compute the term that will have a square root taken discriminant = b*b - 4._r8*a*c - if ( discriminant < 0.0 )then + if ( discriminant < 0.0_r8 )then ! Truncate to zero if the term is very small and negative (relative to b), otherwise error out if ( -discriminant < 3.0_r8*epsilon(b) )then discriminant = 0.0_r8 From 94ad43928bf50ba3548f5d9223b7a948d90aacf4 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 12 Oct 2026 13:50:25 -0600 Subject: [PATCH 09/10] Update src/utils/quadraticMod.F90 Co-authored-by: Samuel Levis --- src/utils/quadraticMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/quadraticMod.F90 b/src/utils/quadraticMod.F90 index 1eb2cd43c8..107adaab22 100644 --- a/src/utils/quadraticMod.F90 +++ b/src/utils/quadraticMod.F90 @@ -21,7 +21,7 @@ subroutine quadratic_roots (a, b, c, root1, root2, file, line) !==============================================================================! !----------------- Solve quadratic equation for its two roots -----------------! !==============================================================================! - ! Implements the numerically stable formulation from Press et al (1986) + ! Implements the numerically stable formulation as discussed in Press et al (1986) ! Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge) ! Allows for roots that are technically complex if they are close to rounding to zero. ! From c130f70c3faa07b1e8d80f1095d22a13841274f5 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 12 Oct 2026 13:54:09 -0600 Subject: [PATCH 10/10] Add back the revision history in --- src/biogeophys/PhotosynthesisMod.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 564dd145b2..5ca4f0c6bd 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -2376,7 +2376,11 @@ 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: + ! + !Dec 14/2012: Jinyun Tang, modified from numerical recipes in F90 by press et al. 1188-1189 + ! !!ARGUMENTS: real(r8), intent(out) :: x ! indepedent variable of the single value function ci_func(x) real(r8), intent(in) :: x1, x2, f1, f2 ! minimum and maximum of the variable domain to search for the solution ci_func(x1) = f1, ci_func(x2)=f2