From ee2bd4f3d7a01c71596fecc8a46719853321b38e Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 2 Jan 2024 10:47:10 -0300 Subject: [PATCH] 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. --- BRAMS/src/lib/numutils.f90 | 483 +++++++++++++++++++++++++++++++ ED/build/make/include.mk.macosx | 12 +- ED/src/dynamics/farq_katul.f90 | 28 +- ED/src/dynamics/farq_leuning.f90 | 193 +++--------- ED/src/init/ed_params.f90 | 61 +++- ED/src/utils/numutils.f90 | 221 ++++++++++++++ 6 files changed, 812 insertions(+), 186 deletions(-) diff --git a/BRAMS/src/lib/numutils.f90 b/BRAMS/src/lib/numutils.f90 index 8d2f47e81..3feb04109 100644 --- a/BRAMS/src/lib/numutils.f90 +++ b/BRAMS/src/lib/numutils.f90 @@ -2333,6 +2333,489 @@ end function eifun8 +!==========================================================================================! +!==========================================================================================! +! Heapsort is a robust and efficient sorting algorithm. For more details, check ! +! The Numerical Recipes Book (chapter 8). ! +!------------------------------------------------------------------------------------------! +subroutine heapsort(nx,xi,increase,xo) + implicit none + !----- Arguments. ----------------------------------------------------------------------! + integer , intent(in) :: nx ! Size of input/output vectors + real , dimension(nx), intent(in) :: xi ! Input vector + logical , intent(in) :: increase ! Sort from small to large? + real , dimension(nx), intent(out) :: xo ! Output vector + !----- Local variables. ----------------------------------------------------------------! + integer :: i ! Counter + integer :: ir ! Index of selected data + integer :: j ! Index of selected data + integer :: l ! Index of selected data + real :: aux ! Placeholder + !---------------------------------------------------------------------------------------! + + + !----- Skip routine in case this has only one element. ---------------------------------! + if (nx < 2) then + xo(:) = xi(:) + return + else if (.not. increase) then + !----- Cheat by making numbers negative. We switch values back in before leaving. --! + xo(:) = -xi(:) + !------------------------------------------------------------------------------------! + else + xo(:) = xi(:) + end if + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! The index l will be decremented from its initial value down to 1 during the ! + ! "hiring" (heap creation) phase. Once it reaches 1, the index ir will be decremented ! + ! from its initial value down to 1 during the "retirement-and-promotion" (heap ! + ! selection) phase. ! + !---------------------------------------------------------------------------------------! + l = nx/2 + 1 + ir = nx + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Main loop. ! + !---------------------------------------------------------------------------------------! + outer_loop: do + !------------------------------------------------------------------------------------! + ! Check whether we are in the hiring phase or in the retirement-and-promotion ! + ! phase. ! + !------------------------------------------------------------------------------------! + if (l > 1) then + !---------------------------------------------------------------------------------! + ! Still in hiring phase. ! + !---------------------------------------------------------------------------------! + l = l - 1 + aux = xo(l) + !---------------------------------------------------------------------------------! + else + !---------------------------------------------------------------------------------! + ! In the retirement-and-promotion phase. ! + !---------------------------------------------------------------------------------! + !----- Clear a space at end of array. --------------------------------------------! + aux = xo(ir) + !----- Retire the top of the heap into it. ---------------------------------------! + xo(ir) = xo(1) + !----- Decrease the size of the corporation. -------------------------------------! + ir = ir -1 + !----- Check how we are doing with promotions. -----------------------------------! + if (ir == 1) then + !----- Done with the last promotion. The least competent worker of all! ------! + xo(1) = aux + !------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------! + ! Time to leave the loop (and the sub-routine). ! + !------------------------------------------------------------------------------! + exit outer_loop + !------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------! + end if + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Whether in hiring phase or promotion phase, we here set up to sift down element ! + ! aux to its proper level. ! + !------------------------------------------------------------------------------------! + i = l + j = l+1 + inner_loop: do + if (j > ir) exit inner_loop + + !----- Compare to the better underling. ------------------------------------------! + if (j < ir) then + if(xo(j) < xo(j+1)) j = j + 1 + end if + !---------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------! + ! Check whether to demote aux or not. ! + !---------------------------------------------------------------------------------! + if (aux < xo(j)) then + !----- Demote aux. ------------------------------------------------------------! + xo(i) = xo(j) + i = j + j = j + j + !------------------------------------------------------------------------------! + else + !----- This is aux's level. Set j to terminate the sift-down. ----------------! + j = ir + 1 + !------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------! + end do inner_loop + !------------------------------------------------------------------------------------! + + !----- Put aux into its slot. -------------------------------------------------------! + xo(i) = aux + !------------------------------------------------------------------------------------! + end do outer_loop + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Before we leave, check whether this should be a high-to-low sorting. In case so, ! + ! switch the sign again. ! + !---------------------------------------------------------------------------------------! + if (.not. increase) xo(:) = -xo(:) + !---------------------------------------------------------------------------------------! + + return +end subroutine heapsort +!==========================================================================================! +!==========================================================================================! + + + + +!==========================================================================================! +!==========================================================================================! +! Function that defines the quantile given a vector. This is a rather simple ! +! estimator, it should work reasonably well as long as x is sufficiently large and not a ! +! crazy distribution. ! +!------------------------------------------------------------------------------------------! +real function fquant(nx,x,prob) + implicit none + !----- Arguments. ----------------------------------------------------------------------! + integer , intent(in) :: nx + real , dimension(nx), intent(in) :: x + real , intent(in) :: prob + !----- Internal variables. -------------------------------------------------------------! + real , dimension(nx) :: xsort + integer :: il + integer :: ih + real :: wl + real :: wh + real :: ridx + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Sanity check: prob must be between 0 and 1. If not, crash! ! + !---------------------------------------------------------------------------------------! + if (prob < 0. .or. prob > 1.) then + write(unit=*,fmt='(a)' ) ' ' + write(unit=*,fmt='(a)' ) ' ' + write(unit=*,fmt='(a)' ) '=================================================' + write(unit=*,fmt='(a)' ) '=================================================' + write(unit=*,fmt='(a)' ) ' In function fquant: Invalid PROB!' + write(unit=*,fmt='(a)' ) '-------------------------------------------------' + write(unit=*,fmt='(a,1x,es12.5)') ' -> Provided PROB: ',prob + write(unit=*,fmt='(a)' ) '-------------------------------------------------' + write(unit=*,fmt='(a)' ) ' ' + write(unit=*,fmt='(a)' ) ' ' + call fatal_error('Invalid prob setting','fquant','numutils.f90') + end if + !---------------------------------------------------------------------------------------! + + + !----- Sort output vector. -------------------------------------------------------------! + call heapsort(nx,x,.true.,xsort) + !---------------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------------! + ! Find the quantile position in terms of indices. ! + !---------------------------------------------------------------------------------------! + !----- Position without interpolation. -------------------------------------------------! + ridx = 1. + prob * real(nx-1) + !----- Index just before ridx. ---------------------------------------------------------! + il = max(1,floor(ridx)) + !----- Index just after ridx. ----------------------------------------------------------! + ih = min(nx,ceiling(ridx)) + !----- Quantile is the interpolated value. ---------------------------------------------! + if (il == ih) then + fquant = xsort(il) + else + !----- Weight factors. --------------------------------------------------------------! + wl = ridx - real(il) + wh = real(ih) - ridx + !------------------------------------------------------------------------------------! + + !----- Quantile is the weighted average. --------------------------------------------! + fquant = (wl * xsort(il) + wh * xsort(ih)) / (wl + wh) + !------------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------------! + + + return +end function fquant +!==========================================================================================! +!==========================================================================================! + + + + +!==========================================================================================! +!==========================================================================================! +! Wrapper for the function above, with an additional mask vector to select only some ! +! elements. ! +!------------------------------------------------------------------------------------------! +real function fquant_mask(nx,x,mask,prob) + implicit none + !----- Arguments. ----------------------------------------------------------------------! + integer , intent(in) :: nx + real , dimension(nx), intent(in) :: x + logical, dimension(nx), intent(in) :: mask + real , intent(in) :: prob + !----- Internal variables. -------------------------------------------------------------! + real , dimension(:) , allocatable :: xuse + integer :: nuse + !----- External functions. -------------------------------------------------------------! + real , external :: fquant + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Count elements to be used, then allocate xuse and xsort. ! + !---------------------------------------------------------------------------------------! + nuse = count(mask) + allocate (xuse(nuse)) + xuse = pack(x,mask) + fquant_mask = fquant(nuse,xuse,prob) + deallocate(xuse) + !---------------------------------------------------------------------------------------! + + return +end function fquant_mask +!==========================================================================================! +!==========================================================================================! + + + + + +!==========================================================================================! +!==========================================================================================! +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! This is an extension of the Numeric Recipes in Fortran 90 to account for the trivial ! +! cases and for checking when the discriminant is negative. ! +! The subroutine also requires a "undef" flag to be passed, which will flag cases ! +! in which one or both solutions are not valid. This is an argument so the solver can be ! +! used when either the largest or the smallest root is sought. ! +!------------------------------------------------------------------------------------------! +subroutine solve_quadratic(aquad,bquad,cquad,undef,root1,root2) + use rconstants, only : tiny_num + implicit none + !----- Arguments. ----------------------------------------------------------------------! + real(kind=4), intent(in) :: aquad + real(kind=4), intent(in) :: bquad + real(kind=4), intent(in) :: cquad + real(kind=4), intent(in) :: undef + real(kind=4), intent(out) :: root1 + real(kind=4), intent(out) :: root2 + !----- Internal variables. -------------------------------------------------------------! + real(kind=4) :: discr + real(kind=4) :: qfact + logical :: a_offzero + logical :: b_offzero + logical :: c_offzero + !---------------------------------------------------------------------------------------! + + + + !----- Save logical tests. -------------------------------------------------------------! + a_offzero = abs(aquad) >= tiny_num + b_offzero = abs(bquad) >= tiny_num + c_offzero = abs(cquad) >= tiny_num + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Check for cases to solve. ! + !---------------------------------------------------------------------------------------! + if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then + !------------------------------------------------------------------------------------! + ! Quadratic equation with two non-zero solutions. Find the discriminant to find ! + ! out whether the solutions are real (if negative, then the roots are complex). ! + !------------------------------------------------------------------------------------! + discr = bquad*bquad - 4.0 * aquad * cquad + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check discriminant sign (but allow for round-off errors). ! + !------------------------------------------------------------------------------------! + if (discr >= - tiny_num) then + !----- Coerce discriminant to non-negative. --------------------------------------! + discr = max(0.0,discr) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Find the q factor as in the numerical recipes, which allows for a more ! + ! robust solution. This is safe whenever b or c are non-zero, as q cannot be ! + ! zero in these cases. + !---------------------------------------------------------------------------------! + qfact = - 0.5 * (bquad + sign(sqrt(discr),bquad)) + root1 = qfact / aquad + root2 = cquad / qfact + else + !----- Negative discriminant, return invalid roots. ------------------------------! + root1 = undef + root2 = undef + !---------------------------------------------------------------------------------! + end if + else if (a_offzero) then + !------------------------------------------------------------------------------------! + ! Both bquad and cquad are nearly zero. Double root, and both have to be zero. ! + !------------------------------------------------------------------------------------! + root1 = 0.0 + root2 = 0.0 + !------------------------------------------------------------------------------------! + else if (b_offzero) then + !------------------------------------------------------------------------------------! + ! "aquad" is not zero, not a true quadratic equation. Single root. ! + !------------------------------------------------------------------------------------! + root1 = - cquad / bquad + root2 = undef + !------------------------------------------------------------------------------------! + else + !------------------------------------------------------------------------------------! + ! Both aquad and bquad are zero, this really doesn't make any sense and should ! + ! never happen. If it does, issue an error and stop the run. ! + !------------------------------------------------------------------------------------! + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a)') ' Quadratic equation cannot be solved!' + write (unit=*,fmt='(a)') ' ''aquad'' and/or ''bquad'' must be non-zero.' + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a,1x,es12.5)') ' aquad = ',aquad + write (unit=*,fmt='(a,1x,es12.5)') ' bquad = ',bquad + write (unit=*,fmt='(a,1x,es12.5)') ' cquad = ',cquad + write (unit=*,fmt='(a)') '------------------------------------------------' + call fatal_error(' Invalid coefficients for quadratic equation' & + ,'solve_quadratic','numutils.f90') + !------------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------------! + + return +end subroutine solve_quadratic +!==========================================================================================! +!==========================================================================================! + + + + + +!==========================================================================================! +!==========================================================================================! +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! This is an extension of the Numeric Recipes in Fortran 90 to account for the trivial ! +! cases and for checking when the discriminant is negative. ! +! The subroutine also requires a "undef" flag to be passed, which will flag cases ! +! in which one or both solutions are not valid. This is an argument so the solver can be ! +! used when either the largest or the smallest root is sought. ! +!------------------------------------------------------------------------------------------! +subroutine solve_quadratic8(aquad,bquad,cquad,undef,root1,root2) + use rconstants, only : tiny_num8 + implicit none + !----- Arguments. ----------------------------------------------------------------------! + real(kind=8), intent(in) :: aquad + real(kind=8), intent(in) :: bquad + real(kind=8), intent(in) :: cquad + real(kind=8), intent(in) :: undef + real(kind=8), intent(out) :: root1 + real(kind=8), intent(out) :: root2 + !----- Internal variables. -------------------------------------------------------------! + real(kind=8) :: discr + real(kind=8) :: qfact + logical :: a_offzero + logical :: b_offzero + logical :: c_offzero + !---------------------------------------------------------------------------------------! + + + + !----- Save logical tests. -------------------------------------------------------------! + a_offzero = abs(aquad) >= tiny_num8 + b_offzero = abs(bquad) >= tiny_num8 + c_offzero = abs(cquad) >= tiny_num8 + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Check for cases to solve. ! + !---------------------------------------------------------------------------------------! + if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then + !------------------------------------------------------------------------------------! + ! Quadratic equation with two non-zero solutions. Find the discriminant to find ! + ! out whether the solutions are real (if negative, then the roots are complex). ! + !------------------------------------------------------------------------------------! + discr = bquad*bquad - 4.d0 * aquad * cquad + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check discriminant sign (but allow for round-off errors). ! + !------------------------------------------------------------------------------------! + if (discr >= - tiny_num8) then + !----- Coerce discriminant to non-negative. --------------------------------------! + discr = max(0.d0,discr) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Find the q factor as in the numerical recipes, which allows for a more ! + ! robust solution. This is safe whenever b or c are non-zero, as q cannot be ! + ! zero in these cases. + !---------------------------------------------------------------------------------! + qfact = - 5.d-1 * (bquad + sign(sqrt(discr),bquad)) + root1 = qfact / aquad + root2 = cquad / qfact + else + !----- Negative discriminant, return invalid roots. ------------------------------! + root1 = undef + root2 = undef + !---------------------------------------------------------------------------------! + end if + else if (a_offzero) then + !------------------------------------------------------------------------------------! + ! Both bquad and cquad are nearly zero. Double root, and both have to be zero. ! + !------------------------------------------------------------------------------------! + root1 = 0.d0 + root2 = 0.d0 + !------------------------------------------------------------------------------------! + else if (b_offzero) then + !------------------------------------------------------------------------------------! + ! "aquad" is not zero, not a true quadratic equation. Single root. ! + !------------------------------------------------------------------------------------! + root1 = - cquad / bquad + root2 = undef + !------------------------------------------------------------------------------------! + else + !------------------------------------------------------------------------------------! + ! Both aquad and bquad are zero, this really doesn't make any sense and should ! + ! never happen. If it does, issue an error and stop the run. ! + !------------------------------------------------------------------------------------! + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a)') ' Quadratic equation cannot be solved!' + write (unit=*,fmt='(a)') ' ''aquad'' and/or ''bquad'' must be non-zero.' + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a,1x,es12.5)') ' aquad = ',aquad + write (unit=*,fmt='(a,1x,es12.5)') ' bquad = ',bquad + write (unit=*,fmt='(a,1x,es12.5)') ' cquad = ',cquad + write (unit=*,fmt='(a)') '------------------------------------------------' + call fatal_error(' Invalid coefficients for quadratic equation' & + ,'solve_quadratic8','numutils.f90') + !------------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------------! + + return +end subroutine solve_quadratic8 +!==========================================================================================! +!==========================================================================================! + + + + + !==========================================================================================! !==========================================================================================! diff --git a/ED/build/make/include.mk.macosx b/ED/build/make/include.mk.macosx index 69ed3365c..9d0d94a2d 100644 --- 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 index 31ac9b186..877910f0e 100644 --- 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 diff --git a/ED/src/dynamics/farq_leuning.f90 b/ED/src/dynamics/farq_leuning.f90 index 8b77da71d..329afbc8f 100644 --- a/ED/src/dynamics/farq_leuning.f90 +++ b/ED/src/dynamics/farq_leuning.f90 @@ -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 [ ---] !------------------------------------------------------------------------------------! @@ -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 !---------------------------------------------------------------------------------! @@ -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 @@ -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) !---------------------------------------------------------------------------------! @@ -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] !------------------------------------------------------------------------------------! @@ -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 !------------------------------------------------------------------------------------! @@ -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 diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index 98b14d90b..3edbea93f 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -5368,11 +5368,10 @@ subroutine init_pft_mort_params() real :: aquad real :: bquad real :: cquad - real :: discr real :: lambda_ref real :: lambda_eff - real :: leff_neg - real :: leff_pos + real :: leff_one + real :: leff_two real, dimension(n_pft) :: rho_use integer :: ipft !----- Local constants for C18 mortality (see below). ----------------------------------! @@ -5391,6 +5390,8 @@ subroutine init_pft_mort_params() ! near-extinction but not so ! high that it would be ! difficult to display in output + real, parameter :: discard = huge(1.0) ! Meaningless value for discarding + ! a root of the quadratic solver. !---------------------------------------------------------------------------------------! @@ -5647,23 +5648,55 @@ subroutine init_pft_mort_params() ! instead, we will wait until the patch age is older than time2canopy. We want, ! ! however, to make the mean patch age to be 1/treefall_disturbance_rate. The ! ! equation below can be retrieved by integrating the steady-state probability ! - ! distribution function. The equation is quadratic and the discriminant will ! - ! never be zero and the treefall_disturbance_rate will be always positive because ! - ! the values of time2canopy and treefall_disturbance_rate have already been ! - ! tested in ed_opspec.F90. ! + ! distribution function. ! !---------------------------------------------------------------------------------! aquad = time2canopy * time2canopy * lambda_ref - 2. * time2canopy bquad = 2. * time2canopy * lambda_ref - 2. cquad = 2. * lambda_ref - !------ Find the discriminant. ---------------------------------------------------! - discr = bquad * bquad - 4. * aquad * cquad - leff_neg = - 0.5 * (bquad - sqrt(discr)) / aquad - leff_pos = - 0.5 * (bquad + sqrt(discr)) / aquad !---------------------------------------------------------------------------------! - ! Use the maximum value, but don't let the value to be too large otherwise ! - ! the negative exponential will cause underflow. ! + + + !------ Solve the quadratic equation. --------------------------------------------! + call solve_quadratic(aquad,bquad,cquad,discard,leff_one,leff_two) + !---------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------! + ! The discriminant should be always positive and the ! + ! treefall_disturbance_rate will be always positive because the values of ! + ! time2canopy and treefall_disturbance_rate have already been tested in ! + ! ed_opspec.F90. In any case, we add a check in here to ensure at least one of ! + ! the solutions is valid. ! !---------------------------------------------------------------------------------! - lambda_eff = min(lnexp_max,max(leff_neg,leff_pos)) + if ( ( leff_one /= discard ) .or. ( leff_two /= discard ) ) then + !------------------------------------------------------------------------------! + ! Use the maximum value, but don't let the value to be too large other- ! + ! wise the negative exponential will cause underflow. ! + !------------------------------------------------------------------------------! + lambda_eff = min(lnexp_max,max(leff_one,leff_two)) + !---------------------------------------------------------------------------------! + 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)' ) 'TIME2CANOPY = ', time2canopy + write(unit=*,fmt='(a,1x,f12.5)' ) 'LAMBA_REF = ', lambda_ref + write(unit=*,fmt='(a,1x,es12.5)') 'A = ', aquad + write(unit=*,fmt='(a,1x,es12.5)') 'B = ', bquad + write(unit=*,fmt='(a,1x,es12.5)') 'C = ', cquad + write(unit=*,fmt='(a,1x,es12.5)') 'DISCR = ', bquad*bquad-4.*aquad*cquad + write(unit=*,fmt='(a)' ) '===========================================' + write(unit=*,fmt='(a)' ) '===========================================' + write(unit=*,fmt='(a)' ) '' + write(unit=*,fmt='(a)' ) '' + call fatal_error(' Negative discriminant when seeking effective lambda.' & + ,'init_pft_mort_params','ed_params.f90') + !------------------------------------------------------------------------------! + end if !---------------------------------------------------------------------------------! else lambda_eff = lambda_ref diff --git a/ED/src/utils/numutils.f90 b/ED/src/utils/numutils.f90 index de2c50065..5207f9b50 100644 --- a/ED/src/utils/numutils.f90 +++ b/ED/src/utils/numutils.f90 @@ -906,3 +906,224 @@ end function fquant_mask !==========================================================================================! !==========================================================================================! + + + + +!==========================================================================================! +!==========================================================================================! +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! This is an extension of the Numeric Recipes in Fortran 90 to account for the trivial ! +! cases and for checking when the discriminant is negative. ! +! The subroutine also requires a "undef" flag to be passed, which will flag cases ! +! in which one or both solutions are not valid. This is an argument so the solver can be ! +! used when either the largest or the smallest root is sought. ! +!------------------------------------------------------------------------------------------! +subroutine solve_quadratic(aquad,bquad,cquad,undef,root1,root2) + use consts_coms, only : tiny_num + implicit none + !----- Arguments. ----------------------------------------------------------------------! + real(kind=4), intent(in) :: aquad + real(kind=4), intent(in) :: bquad + real(kind=4), intent(in) :: cquad + real(kind=4), intent(in) :: undef + real(kind=4), intent(out) :: root1 + real(kind=4), intent(out) :: root2 + !----- Internal variables. -------------------------------------------------------------! + real(kind=4) :: discr + real(kind=4) :: qfact + logical :: a_offzero + logical :: b_offzero + logical :: c_offzero + !---------------------------------------------------------------------------------------! + + + + !----- Save logical tests. -------------------------------------------------------------! + a_offzero = abs(aquad) >= tiny_num + b_offzero = abs(bquad) >= tiny_num + c_offzero = abs(cquad) >= tiny_num + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Check for cases to solve. ! + !---------------------------------------------------------------------------------------! + if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then + !------------------------------------------------------------------------------------! + ! Quadratic equation with two non-zero solutions. Find the discriminant to find ! + ! out whether the solutions are real (if negative, then the roots are complex). ! + !------------------------------------------------------------------------------------! + discr = bquad*bquad - 4.0 * aquad * cquad + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check discriminant sign (but allow for round-off errors). ! + !------------------------------------------------------------------------------------! + if (discr >= - tiny_num) then + !----- Coerce discriminant to non-negative. --------------------------------------! + discr = max(0.0,discr) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Find the q factor as in the numerical recipes, which allows for a more ! + ! robust solution. This is safe whenever b or c are non-zero, as q cannot be ! + ! zero in these cases. + !---------------------------------------------------------------------------------! + qfact = - 0.5 * (bquad + sign(sqrt(discr),bquad)) + root1 = qfact / aquad + root2 = cquad / qfact + else + !----- Negative discriminant, return invalid roots. ------------------------------! + root1 = undef + root2 = undef + !---------------------------------------------------------------------------------! + end if + else if (a_offzero) then + !------------------------------------------------------------------------------------! + ! Both bquad and cquad are nearly zero. Double root, and both have to be zero. ! + !------------------------------------------------------------------------------------! + root1 = 0.0 + root2 = 0.0 + !------------------------------------------------------------------------------------! + else if (b_offzero) then + !------------------------------------------------------------------------------------! + ! "aquad" is not zero, not a true quadratic equation. Single root. ! + !------------------------------------------------------------------------------------! + root1 = - cquad / bquad + root2 = undef + !------------------------------------------------------------------------------------! + else + !------------------------------------------------------------------------------------! + ! Both aquad and bquad are zero, this really doesn't make any sense and should ! + ! never happen. If it does, issue an error and stop the run. ! + !------------------------------------------------------------------------------------! + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a)') ' Quadratic equation cannot be solved!' + write (unit=*,fmt='(a)') ' ''aquad'' and/or ''bquad'' must be non-zero.' + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a,1x,es12.5)') ' aquad = ',aquad + write (unit=*,fmt='(a,1x,es12.5)') ' bquad = ',bquad + write (unit=*,fmt='(a,1x,es12.5)') ' cquad = ',cquad + write (unit=*,fmt='(a)') '------------------------------------------------' + call fatal_error(' Invalid coefficients for quadratic equation' & + ,'solve_quadratic','numutils.f90') + !------------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------------! + + return +end subroutine solve_quadratic +!==========================================================================================! +!==========================================================================================! + + + + + +!==========================================================================================! +!==========================================================================================! +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! This is an extension of the Numeric Recipes in Fortran 90 to account for the trivial ! +! cases and for checking when the discriminant is negative. ! +! The subroutine also requires a "undef" flag to be passed, which will flag cases ! +! in which one or both solutions are not valid. This is an argument so the solver can be ! +! used when either the largest or the smallest root is sought. ! +!------------------------------------------------------------------------------------------! +subroutine solve_quadratic8(aquad,bquad,cquad,undef,root1,root2) + use consts_coms, only : tiny_num8 + implicit none + !----- Arguments. ----------------------------------------------------------------------! + real(kind=8), intent(in) :: aquad + real(kind=8), intent(in) :: bquad + real(kind=8), intent(in) :: cquad + real(kind=8), intent(in) :: undef + real(kind=8), intent(out) :: root1 + real(kind=8), intent(out) :: root2 + !----- Internal variables. -------------------------------------------------------------! + real(kind=8) :: discr + real(kind=8) :: qfact + logical :: a_offzero + logical :: b_offzero + logical :: c_offzero + !---------------------------------------------------------------------------------------! + + + + !----- Save logical tests. -------------------------------------------------------------! + a_offzero = abs(aquad) >= tiny_num8 + b_offzero = abs(bquad) >= tiny_num8 + c_offzero = abs(cquad) >= tiny_num8 + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Check for cases to solve. ! + !---------------------------------------------------------------------------------------! + if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then + !------------------------------------------------------------------------------------! + ! Quadratic equation with two non-zero solutions. Find the discriminant to find ! + ! out whether the solutions are real (if negative, then the roots are complex). ! + !------------------------------------------------------------------------------------! + discr = bquad*bquad - 4.d0 * aquad * cquad + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check discriminant sign (but allow for round-off errors). ! + !------------------------------------------------------------------------------------! + if (discr >= - tiny_num8) then + !----- Coerce discriminant to non-negative. --------------------------------------! + discr = max(0.d0,discr) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Find the q factor as in the numerical recipes, which allows for a more ! + ! robust solution. This is safe whenever b or c are non-zero, as q cannot be ! + ! zero in these cases. + !---------------------------------------------------------------------------------! + qfact = - 5.d-1 * (bquad + sign(sqrt(discr),bquad)) + root1 = qfact / aquad + root2 = cquad / qfact + else + !----- Negative discriminant, return invalid roots. ------------------------------! + root1 = undef + root2 = undef + !---------------------------------------------------------------------------------! + end if + else if (a_offzero) then + !------------------------------------------------------------------------------------! + ! Both bquad and cquad are nearly zero. Double root, and both have to be zero. ! + !------------------------------------------------------------------------------------! + root1 = 0.d0 + root2 = 0.d0 + !------------------------------------------------------------------------------------! + else if (b_offzero) then + !------------------------------------------------------------------------------------! + ! "aquad" is not zero, not a true quadratic equation. Single root. ! + !------------------------------------------------------------------------------------! + root1 = - cquad / bquad + root2 = undef + !------------------------------------------------------------------------------------! + else + !------------------------------------------------------------------------------------! + ! Both aquad and bquad are zero, this really doesn't make any sense and should ! + ! never happen. If it does, issue an error and stop the run. ! + !------------------------------------------------------------------------------------! + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a)') ' Quadratic equation cannot be solved!' + write (unit=*,fmt='(a)') ' ''aquad'' and/or ''bquad'' must be non-zero.' + write (unit=*,fmt='(a)') '------------------------------------------------' + write (unit=*,fmt='(a,1x,es12.5)') ' aquad = ',aquad + write (unit=*,fmt='(a,1x,es12.5)') ' bquad = ',bquad + write (unit=*,fmt='(a,1x,es12.5)') ' cquad = ',cquad + write (unit=*,fmt='(a)') '------------------------------------------------' + call fatal_error(' Invalid coefficients for quadratic equation' & + ,'solve_quadratic8','numutils.f90') + !------------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------------! + + return +end subroutine solve_quadratic8 +!==========================================================================================! +!==========================================================================================!