Skip to content
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

Revised numerical methods for a few functions #409

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Revised numerical methods for a few functions #409

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ee2bd4f
Updated quadratic solvers in ED2 to use the numerical recipe approach…
mpaiao Jan 2, 2024
cd53bad
Revised some functions and sub-routines in numutils.f90 for improved …
mpaiao Dec 31, 2024
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
754 changes: 682 additions & 72 deletions BRAMS/src/lib/numutils.f90
View file
Open in desktop

Large diffs are not rendered by default.

15 changes: 13 additions & 2 deletions ED/src/dynamics/farq_katul.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
!=======================================================================================!
!=======================================================================================!
Expand Down Expand Up @@ -848,6 +855,8 @@ 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 !! derivatives of b wrt. gsc
real(kind=8) :: dcdg !! derivatives of c wrt. gsc
!------------------------------------------------------------------------------------!
Expand Down Expand Up @@ -877,7 +886,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
Expand Down Expand Up @@ -905,7 +915,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
Expand Down
193 changes: 42 additions & 151 deletions ED/src/dynamics/farq_leuning.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -403,7 +403,6 @@ subroutine comp_photo_tempfun(ib,leaf_aging_factor,green_leaf_factor)
real(kind=8) :: aterm !> 'a' term for quadratic equation [ ---]
real(kind=8) :: bterm !> 'b' term for quadratic equation [ ---]
real(kind=8) :: cterm !> 'c' term for quadratic equation [ ---]
real(kind=8) :: discr !> discriminant for quadratic eqn [ ---]
real(kind=8) :: jroot1 !> root for quadratic equation [ ---]
real(kind=8) :: jroot2 !> root for quadratic equation [ ---]
!------------------------------------------------------------------------------------!
Expand Down Expand Up @@ -659,54 +658,41 @@ subroutine comp_photo_tempfun(ib,leaf_aging_factor,green_leaf_factor)



!----- Solve the quadratic equation. ---------------------------------------------!
call solve_quadratic8(aterm,bterm,cterm,discard,jroot1,jroot2)
!---------------------------------------------------------------------------------!
if (aterm == 0.d0) then
!----- Equation was reduced to linear equation. -------------------------------!
aparms(ib)%jact = - cterm / bterm
!------------------------------------------------------------------------------!



!---------------------------------------------------------------------------------!
! Make sure at least one root is valid. !
!---------------------------------------------------------------------------------!
if ( ( jroot1 /= discard ) .or. ( jroot2 /= discard ) ) then
aparms(ib)%jact = min( jroot1, jroot2 )
else
!----- Check whether discriminant is positive. --------------------------------!
discr = bterm * bterm - 4.d0 * aterm * cterm
if (discr == 0.d0) then
aparms(ib)%jact = - bterm / (2.d0 * aterm)
elseif (discr > 0.d0) then
!----- Always pick the smallest of the roots. ------------------------------!
jroot1 = ( - bterm - sqrt(discr) ) / ( 2.d0 * aterm)
jroot2 = ( - bterm + sqrt(discr) ) / ( 2.d0 * aterm)
if (jroot1 <= jroot2) then
aparms(ib)%jact = jroot1
else
aparms(ib)%jact = jroot2
end if
!---------------------------------------------------------------------------!
else
!----- Broadcast the bad news. ---------------------------------------------!
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) '========================================'
write(unit=*,fmt='(a)' ) '========================================'
write(unit=*,fmt='(a)' ) ' Discriminant is negative.'
write(unit=*,fmt='(a)' ) '----------------------------------------'
write(unit=*,fmt='(a,1x,f12.5)' ) 'PAR = ', met(ib)%par * mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'QYIELD = ', thispft(ib)%phi_psII
write(unit=*,fmt='(a,1x,f12.5)' ) 'CURVATURE = ', thispft(ib)%curvpar
write(unit=*,fmt='(a,1x,f12.5)' ) 'JMAX = ', aparms(ib)%jm &
* mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'JM0 = ', thispft(ib)%jm0 &
* mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'LEAF_TEMP = ', met(ib)%leaf_temp + t008
write(unit=*,fmt='(a,1x,es12.5)') 'A = ', aterm
write(unit=*,fmt='(a,1x,es12.5)') 'B = ', bterm
write(unit=*,fmt='(a,1x,es12.5)') 'C = ', cterm
write(unit=*,fmt='(a,1x,es12.5)') 'DISCR = ', discr
write(unit=*,fmt='(a)' ) '========================================'
write(unit=*,fmt='(a)' ) '========================================'
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) ''
call fatal_error(' Negative discriminant when seeking the J rate.' &
,'comp_photo_tempfun','farq_leuning.f90')
!---------------------------------------------------------------------------!
end if
!----- Broadcast the bad news. ------------------------------------------------!
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) '==========================================='
write(unit=*,fmt='(a)' ) '==========================================='
write(unit=*,fmt='(a)' ) ' Discriminant is negative.'
write(unit=*,fmt='(a)' ) '-------------------------------------------'
write(unit=*,fmt='(a,1x,f12.5)' ) 'PAR = ', met(ib)%par * mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'QYIELD = ', thispft(ib)%phi_psII
write(unit=*,fmt='(a,1x,f12.5)' ) 'CURVATURE = ', thispft(ib)%curvpar
write(unit=*,fmt='(a,1x,f12.5)' ) 'JMAX = ', aparms(ib)%jm * mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'JM0 = ', thispft(ib)%jm0 * mol_2_umol8
write(unit=*,fmt='(a,1x,f12.5)' ) 'LEAF_TEMP = ', met(ib)%leaf_temp + t008
write(unit=*,fmt='(a,1x,es12.5)') 'A = ', aterm
write(unit=*,fmt='(a,1x,es12.5)') 'B = ', bterm
write(unit=*,fmt='(a,1x,es12.5)') 'C = ', cterm
write(unit=*,fmt='(a,1x,es12.5)') 'DISCR = ', bterm*bterm-4.d0*aterm*cterm
write(unit=*,fmt='(a)' ) '==========================================='
write(unit=*,fmt='(a)' ) '==========================================='
write(unit=*,fmt='(a)' ) ''
write(unit=*,fmt='(a)' ) ''
call fatal_error(' Negative discriminant when seeking the J rate.' &
,'comp_photo_tempfun','farq_leuning.f90')
!------------------------------------------------------------------------------!
end if
!---------------------------------------------------------------------------------!
Expand Down Expand Up @@ -1004,7 +990,6 @@ subroutine solve_aofixed_case(ib,answer)
real(kind=8) :: aquad
real(kind=8) :: bquad
real(kind=8) :: cquad
real(kind=8) :: discr
real(kind=8) :: gswroot1
real(kind=8) :: gswroot2
real(kind=8) :: ciroot1
Expand Down Expand Up @@ -1071,31 +1056,7 @@ subroutine solve_aofixed_case(ib,answer)
bquad = qterm1 * qterm2 - aquad * thispft(ib)%b - qterm3
cquad = - qterm1 * qterm2 * thispft(ib)%b - qterm3 * met(ib)%blyr_cond_h2o
!----- Solve the quadratic equation for gsw. -------------------------------------!
if (aquad == 0.d0) then
!----- Not really a quadratic equation. ---------------------------------------!
gswroot1 = -cquad / bquad
gswroot2 = discard
else
!----- A quadratic equation, find the discriminant. ---------------------------!
discr = bquad * bquad - 4.d0 * aquad * cquad
!----- Decide what to do based on the discriminant. ---------------------------!
if (discr == 0.d0) then
!----- Double root. --------------------------------------------------------!
gswroot1 = - bquad / (2.d0 * aquad)
gswroot2 = discard
!---------------------------------------------------------------------------!
elseif (discr > 0.d0) then
!----- Two distinct roots. -------------------------------------------------!
gswroot1 = (- bquad - sqrt(discr)) / (2.d0 * aquad)
gswroot2 = (- bquad + sqrt(discr)) / (2.d0 * aquad)
!---------------------------------------------------------------------------!
else
!----- None of the roots are real, this solution failed. -------------------!
return
!---------------------------------------------------------------------------!
end if
!------------------------------------------------------------------------------!
end if
call solve_quadratic8(aquad,bquad,cquad,discard,gswroot1,gswroot2)
!---------------------------------------------------------------------------------!


Expand Down Expand Up @@ -1868,7 +1829,6 @@ subroutine find_lint_co2_bounds(ib,cimin,cimax,bounded)
real(kind=8) :: ytmp ! variable for ciQ
real(kind=8) :: ztmp ! variable for ciQ
real(kind=8) :: wtmp ! variable for ciQ
real(kind=8) :: discr ! The discriminant of the quad. eq. [mol2/mol2]
real(kind=8) :: ciroot1 ! 1st root for the quadratic eqn. [ mol/mol]
real(kind=8) :: ciroot2 ! 2nd root for the quadratic eqn. [ mol/mol]
!------------------------------------------------------------------------------------!
Expand Down Expand Up @@ -1898,46 +1858,9 @@ subroutine find_lint_co2_bounds(ib,cimin,cimax,bounded)
+ met(ib)%blyr_cond_co2 * aparms(ib)%tau + aparms(ib)%rho
cquad = aparms(ib)%tau * (aparms(ib)%nu - met(ib)%blyr_cond_co2 * met(ib)%can_co2 ) &
+ aparms(ib)%sigma
!----- 2. Decide whether this is a true quadratic case or not. ----------------------!
if (aquad /= 0.d0) then
!---------------------------------------------------------------------------------!
! This is a true quadratic case, the first step is to find the discriminant. !
!---------------------------------------------------------------------------------!
discr = bquad * bquad - 4.d0 * aquad * cquad
if (discr == 0.d0) then
!------------------------------------------------------------------------------!
! Discriminant is zero, both roots are the same. We save only one, and !
! make the other negative, which will make the guess discarded. !
!------------------------------------------------------------------------------!
ciroot1 = - bquad / (2.d0 * aquad)
ciroot2 = - discard
elseif (discr > 0.d0) then
!----- Positive discriminant, two solutions are possible. ---------------------!
ciroot1 = (- bquad + sqrt(discr)) / (2.d0 * aquad)
ciroot2 = (- bquad - sqrt(discr)) / (2.d0 * aquad)
!------------------------------------------------------------------------------!
else
!----- Discriminant is negative. Impossible to solve. ------------------------!
ciroot1 = - discard
ciroot2 = - discard
!------------------------------------------------------------------------------!
end if
!---------------------------------------------------------------------------------!
else
!---------------------------------------------------------------------------------!
! This is a linear case, the xi term is zero. There is only one number that !
! works for this case. !
!---------------------------------------------------------------------------------!
ciroot1 = - cquad / bquad
ciroot2 = - discard
!----- Not used, just for the debugging process. ---------------------------------!
discr = bquad * bquad
!---------------------------------------------------------------------------------!
end if
!------------------------------------------------------------------------------------!
! Save the largest of the values. In case both were discarded, we switch it to !
! the positive discard so this will never be chosen. !
!------------------------------------------------------------------------------------!
!----- 2. Solve the quadratic equation. ---------------------------------------------!
call solve_quadratic8(aquad,bquad,cquad,-discard,ciroot1,ciroot2)
!----- 3. Save the largest value (or flag them as invalid when both are invalid). ---!
cigsw=max(ciroot1, ciroot2)
if (cigsw == -discard) cigsw = discard
!------------------------------------------------------------------------------------!
Expand All @@ -1962,44 +1885,12 @@ subroutine find_lint_co2_bounds(ib,cimin,cimax,bounded)
bquad = ytmp * aparms(ib)%rho - aparms(ib)%xi * wtmp + ztmp * aparms(ib)%tau
cquad = - aparms(ib)%tau * wtmp + ytmp * aparms(ib)%sigma

!----- 3. Decide whether this is a true quadratic case or not. ----------------------!
if (aquad /= 0.d0) then
!---------------------------------------------------------------------------------!
! This is a true quadratic case, the first step is to find the discriminant. !
!---------------------------------------------------------------------------------!
discr = bquad * bquad - 4.d0 * aquad * cquad
if (discr == 0.d0) then
!------------------------------------------------------------------------------!
! Discriminant is zero, both roots are the same. We save only one, and !
! make the other negative, which will make the guess discarded. !
!------------------------------------------------------------------------------!
ciroot1 = - bquad / (2.d0 * aquad)
ciroot2 = -discard
!------------------------------------------------------------------------------!
elseif (discr > 0.d0) then
!----- Positive discriminant, two solutions are possible. ---------------------!
ciroot1 = (- bquad + sqrt(discr)) / (2.d0 * aquad)
ciroot2 = (- bquad - sqrt(discr)) / (2.d0 * aquad)
!------------------------------------------------------------------------------!
else
!----- Discriminant is negative. Impossible to solve. ------------------------!
ciroot1 = -discard
ciroot2 = -discard
!------------------------------------------------------------------------------!
end if
else
!---------------------------------------------------------------------------------!
! This is a linear case, the xi term is zero. There is only one number !
! that works for this case. !
!---------------------------------------------------------------------------------!
ciroot1 = - cquad / bquad
ciroot2 = -discard
!----- Not used, just for the debugging process. ---------------------------------!
discr = bquad * bquad
end if
!----- 3. Solve the quadratic equation. ---------------------------------------------!
call solve_quadratic8(aquad,bquad,cquad,-discard,ciroot1,ciroot2)

!------------------------------------------------------------------------------------!
! Save the largest of the values. In case both were discarded, we switch it to !
! the positive discard so this will never be chosen. !
! 4. Save the largest of the values. In case both were discarded, we switch it !
! to the positive discard so this will never be chosen. !
!------------------------------------------------------------------------------------!
ciQ=max(ciroot1, ciroot2)
if (ciQ == -discard) ciQ = discard
Expand Down
Loading
Loading