Updated quadratic solvers in ED2 to use the numerical recipe approach…

…, which avoids severe round-off errors when terms are different by orders of magnitude.
mpaiao · Jan 2, 2024 · ee2bd4f · ee2bd4f
1 parent 4809be5
commit ee2bd4f
Show file tree

Hide file tree

Showing 6 changed files with 812 additions and 186 deletions.
diff --git a/BRAMS/src/lib/numutils.f90 b/BRAMS/src/lib/numutils.f90
diff --git a/ED/build/make/include.mk.macosx b/ED/build/make/include.mk.macosx
@@ -20,8 +20,8 @@ BASE=$(ED_ROOT)/build/
 # libraries compiled with the same compiler you set for F_COMP and C_COMP.  You may still  #
 # be able to compile without HDF5 but it will not run.                                     #
 #------------------------------------------------------------------------------------------#
-ZLIB_PATH=/Users/mlongo/Util/ED2_Libs/zlib-1.2.11
-HDF5_PATH=/Users/mlongo/Util/ED2_Libs/hdf5-1.10.2
+ZLIB_PATH=/usr/local
+HDF5_PATH=/usr/local
 HDF5_INCS=-I$(HDF5_PATH)/include
 HDF5C_INCS=$(HDF5_INCS)
 HDF5_LIBS=-L$(ZLIB_PATH)/lib -lz -L$(HDF5_PATH)/lib -lhdf5_fortran -lhdf5 -lhdf5_hl
@@ -51,10 +51,10 @@ USE_INTERF=1
 #################################### COMPILER SETTINGS #####################################
 CMACH=MAC_OS_X
 FC_TYPE=GNU
-F_COMP=gfortran
-C_COMP=gcc
-LOADER=gfortran
-C_LOADER=gcc
+F_COMP=gfortran-13
+C_COMP=gcc-13
+LOADER=gfortran-13
+C_LOADER=gcc-13
 LIBS=
 MOD_EXT=mod
 ############################################################################################

diff --git a/ED/src/dynamics/farq_katul.f90 b/ED/src/dynamics/farq_katul.f90
@@ -27,6 +27,13 @@
 !==========================================================================================!
 module farq_katul
 
+   !---------------------------------------------------------------------------------------!
+   !     This is a flag used in various sub-routines and functions and denote that we      !
+   ! should ignore the result.                                                             !
+   !---------------------------------------------------------------------------------------!
+   real(kind=8), parameter :: discard = huge(1.d0)
+   !---------------------------------------------------------------------------------------!
+
    contains
    !=======================================================================================!
    !=======================================================================================!
@@ -409,9 +416,6 @@ subroutine optimization_solver8(ib,ipft,lambda,
         real(kind=8)                :: test_dcidg
         real(kind=8)                :: test_dfcdg
         real(kind=8)                :: test_dfedg
-        real(kind=8)                :: test_fc_light
-        real(kind=8)                :: test_fc_rubp
-        real(kind=8)                :: test_fc_3rd
         real(kind=8)                :: opt_ci_light
         real(kind=8)                :: opt_ci_rubp
         real(kind=8)                :: opt_ci_3rd
@@ -851,17 +855,9 @@ subroutine photosynthesis_stomata_solver8(ib,gsc,limit_case,
         real(kind=8)                :: k1,k2        !! Variable used in photosynthesis equation
         real(kind=8)                :: a,b,c        !! Coefficients of the quadratic equation to solve ci
         real(kind=8)                :: rad          !! sqrt(b2-4ac)
+        real(kind=8)                :: ciroot1      !! First root of ci
+        real(kind=8)                :: ciroot2      !! Second root of ci
         real(kind=8)                :: dbdg,dcdg    !! derivatives of b,c wrt. gsc
-        real(kind=8)                :: ci_rubp      !! ci for rubp-limited scenario
-        real(kind=8)                :: dcidg_rubp   !! derivative of ci wrt. gsc for rubp-limited scenario
-        real(kind=8)                :: dfcdg_rubp   !! derivative of fc wrt. gsc for rubp-limited scenario
-        real(kind=8)                :: ci_light     !! ci for light-limited scenario
-        real(kind=8)                :: dcidg_light  !! derivative of ci wrt. gsc for light-limited scenario
-        real(kind=8)                :: dfcdg_light  !! derivative of fc wrt. gsc for light-limited scenario
-        real(kind=8)                :: ci_3rd       !! ci for TPU/CO2-limited scenario
-        real(kind=8)                :: dcidg_3rd    !! derivative of ci wrt. gsc for TPU/CO2-limited scenario
-        real(kind=8)                :: dfcdg_3rd    !! derivative of fc wrt. gsc for TPU/CO2-limited scenario
-
         !------------------------------------------------------------------------------------!
 
 
@@ -889,7 +885,8 @@ subroutine photosynthesis_stomata_solver8(ib,gsc,limit_case,
 
             ! solve the quadratic equation
             rad = sqrt(b ** 2 - 4.d0 * a * c)
-            ci = - (b - rad) / (2.d0 * a)
+            call solve_quadratic8(a,b,c,-discard,ciroot1,ciroot2)
+            ci = max(ciroot1,ciroot2)
             fc = gsbc * (met(ib)%can_co2 - ci)
 
             ! calculate derivatives
@@ -917,7 +914,8 @@ subroutine photosynthesis_stomata_solver8(ib,gsc,limit_case,
             c = (-k1 * aparms(ib)%compp - k2 * aparms(ib)%leaf_resp) / gsbc - k2 * met(ib)%can_co2
 
             rad = sqrt(b ** 2 - 4.d0 * a * c)
-            ci = - (b - rad) / (2.d0 * a)
+            call solve_quadratic8(a,b,c,-discard,ciroot1,ciroot2)
+            ci = max(ciroot1,ciroot2)
             fc = gsbc * (met(ib)%can_co2 - ci)
 
             ! calculate derivatives