diff --git a/BRAMS/src/lib/numutils.f90 b/BRAMS/src/lib/numutils.f90 index 8d2f47e81..aadaf4a6c 100644 --- a/BRAMS/src/lib/numutils.f90 +++ b/BRAMS/src/lib/numutils.f90 @@ -1811,7 +1811,7 @@ end function expected !==========================================================================================! !==========================================================================================! -! This function simply computes f(x) = exp (-x²). Just to make errorfun neater. ! +! This function simply computes f(x) = exp (-x^2). Just to make errorfun neater. ! !------------------------------------------------------------------------------------------! real function expmsq(x) implicit none @@ -2219,7 +2219,6 @@ end subroutine diagon - !==========================================================================================! !==========================================================================================! ! EIFUN8 -- This function computes the exponential integral function, defined by ! @@ -2228,35 +2227,49 @@ end subroutine diagon ! | exp(t) ! ! Ei(x) = | -------- dt ! ! _| t ! -! 0 ! +! -Inf ! ! ! -! This function is based on: ! +! This function checks for two approaches: series expansion, which typically works ! +! best when x is small, and the asymptotic expansion, which typically works best when x is ! +! large. The approach selects the smallest result (in absolute numbers) as the most ! +! accurate method. Both the series expansion and the asymptotic expansion are provided in ! +! AS72. This approach also checks for some other edge cases, and ignores the results when ! +! the value is very negative. ! ! ! -! Press, W. H., S. A. Teukolsky, W. T. Vetterling, B. P. Flannery: 1992. Numerical recipes ! -! in Fortran 77. Cambridge University Press, section 6.3 p. 215-219. ! +! Reference: ! +! ! +! Abramowitz, M., and I. A. Stegun, Eds., 1972: Handbook of mathematical functions with ! +! formulas, graphs, and mathematical tables. 10th ed., No. 55, Applied Mathematics ! +! Series, National Bureau of Standards, Washington, DC, USA (AS72). ! ! ! -! with the difference that we also solve for negative numbers. Zero cannot be solved, so ! -! if this happens, or if the sought number would lead to infinity, we stop the model. ! !------------------------------------------------------------------------------------------! real(kind=8) function eifun8(x) - use rconstants, only : euler_gam8 & ! intent(in) - , lnexp_min8 & ! intent(in) - , lnexp_max8 & ! intent(in) - , tiny_num8 ! ! intent(in) + use rconstants, only : euler_gam8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 & ! intent(in) + , tiny_num8 & ! intent(in) + , almost_zero8 ! ! intent(in) implicit none !----- Arguments. ----------------------------------------------------------------------! real(kind=8), intent(in) :: x !----- Local variables. ----------------------------------------------------------------! - real(kind=8) :: sum - real(kind=8) :: term - real(kind=8) :: fact - real(kind=8) :: prev - real(kind=8) :: diter - integer :: iter + real(kind=8) :: usum + real(kind=8) :: uxk + real(kind=8) :: uxkm1 + real(kind=8) :: vsum + real(kind=8) :: vxkm1 + real(kind=8) :: vxk + real(kind=8) :: ei_series + real(kind=8) :: ei_asymptote + integer :: k !----- Local constants. ----------------------------------------------------------------! - real(kind=8), parameter :: powerlim = 1.5d+01 - real(kind=8), parameter :: converge = 1.0d-7 - integer , parameter :: maxiter = 100 + real(kind=8), parameter :: discard8 = 1.0d+36 + integer , parameter :: maxiter = 100 + !----- Polynomial coefficients. --------------------------------------------------------! + real(kind=8), dimension(4), parameter :: apoly = (/ 8.5733287401d+00, 1.8059015973d+01 & + , 8.6347608925d+00, 2.6777373430d-01 /) + real(kind=8), dimension(4), parameter :: bpoly = (/ 9.5733223454d+00, 2.5632956149d+01 & + , 2.1099653083d+01, 3.9584969228d+00 /) !---------------------------------------------------------------------------------------! @@ -2266,63 +2279,157 @@ real(kind=8) function eifun8(x) if (x == 0.d0) then !------------------------------------------------------------------------------------! ! Zero. This is a singularity and the user should never call it in this case. ! - ! That's sad, but we ought to quit this run and tell the user why the run crashed. ! - !------------------------------------------------------------------------------------! - call abort_run('Exponential integral cannot be solved for x = 0.' & - ,'eifun8','numutils.f90') - elseif (x >= lnexp_max8) then - !----- Huge value, crash because this is iminent over-flow. -------------------------! - write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x - write(unit=*,fmt='(a,1x,es12.5)') 'Maximum acceptable X =',lnexp_max8 - call abort_run('Exponential integral cannot be solved for x = 0.' & - ,'eifun8','numutils.f90') - elseif (abs(x) <= lnexp_min8) then - !----- Huge negative number, the result can be rounded to zero. ---------------------! + !------------------------------------------------------------------------------------! + stop 'Exponential integral cannot be solved for x = 0.' + !------------------------------------------------------------------------------------! + elseif (x <= lnexp_min8) then + !------------------------------------------------------------------------------------! + ! Huge negative value, the result can be set to zero. ! + !------------------------------------------------------------------------------------! eifun8 = 0.d0 - elseif (abs(x) <= tiny_num8) then - !----- The number is too close to zero, bypass iterative methods. -------------------! - eifun8 = euler_gam8 + log(abs(x)) - elseif (abs(x) <= powerlim) then - !------------------------------------------------------------------------------------! - ! Input x is small, so we use the power method. ! - !------------------------------------------------------------------------------------! - fact = 1.d0 - sum = 0.d0 - powerloop: do iter=1,maxiter - diter = dble(iter) - fact = fact * x / diter - term = fact / diter - sum = sum + term - !----- If the term is tiny, we have reached convergence, quit the loop. ----------! - if (abs(term) < converge * abs(sum)) exit powerloop - end do powerloop - eifun8 = euler_gam8 + log(abs(x)) + sum + !------------------------------------------------------------------------------------! + elseif (x <= -1.d0) then + !------------------------------------------------------------------------------------! + ! For negative values less than -1.0, we use the polynomial approximation ! + ! (Equation 5.1.56 of AS72), by taking that Ei(x) = - E1(-x). ! + !------------------------------------------------------------------------------------! + eifun8 = exp(x)/x & + * ( x * ( x * ( x * ( x - apoly(1) ) + apoly(2) ) - apoly(3) ) + apoly(4) ) & + / ( x * ( x * ( x * ( x - bpoly(1) ) + bpoly(2) ) - bpoly(3) ) + bpoly(4) ) + !------------------------------------------------------------------------------------! else !------------------------------------------------------------------------------------! - ! Input x is large, so we use the asymptotic approximation. ! - !------------------------------------------------------------------------------------! - sum = 0.d0 - term = 1.d0 - asymploop: do iter=1,maxiter - diter = dble(iter) - prev = term - term = term * diter / x - if (abs(term) < converge) then - !----- The term is tiny, we have reached convergence, quit the loop. ----------! - exit asymploop - elseif (abs(term) >= abs(prev)) then - !------------------------------------------------------------------------------! - ! Series is diverging, we are probably reaching round-off errors, we better ! - ! stop now. ! - !------------------------------------------------------------------------------! - sum = sum - prev - exit asymploop + ! Find both the series expansion and the asymptotic expansion, and pick the one ! + ! with the lowest absolute value. ! + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Series expansion: Equation 5.1.11 of AS72, by taking that Ei(x) = - E1(-x). ! + ! ! + ! Inf ! + ! Ei(x) = gamma + ln(x) + SUM u(x,k), ! + ! k=1 ! + ! ! + ! where u(x,k) = [ (-1)^k * (-x)^k / (k * k!) ] ! + ! ! + ! To efficiently compute the terms inside the summation, we use that: ! + ! ! + ! u(x,k) = x * (k -1) / k^2 * u(x,k-1), for k >= 2. ! + !------------------------------------------------------------------------------------! + uxk = x + usum = uxk + do_expansion: do k = 2, maxiter + !----- Update the current summation term. ----------------------------------------! + uxkm1 = uxk + uxk = x * dble( k - 1 ) / dble( k * k ) * uxkm1 + !----- Check for degenerate or very large estimate. ------------------------------! + if ( abs(uxk) > discard8 .or. abs(usum) > discard8) then + usum = sign(discard8,usum) + exit do_expansion + end if + !----- Check for convergence. ----------------------------------------------------! + if ( any(abs(uxk) <= [ almost_zero8 * abs(usum), tiny_num8] ) ) exit do_expansion + !----- Update summation. ---------------------------------------------------------! + usum = usum + uxk + !---------------------------------------------------------------------------------! + end do do_expansion + !----- Find the series solution. ----------------------------------------------------! + if ( abs(usum) == discard8) then + ei_series = sign(discard8,usum) + else + ei_series = euler_gam8 + log(abs(x)) + usum + end if + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Asymptote expansion: Equation 5.1.51 of AS72 by taking AS72's n=1 and ! + ! Ei(x) = -E1(-x)). ! + ! ! + ! Inf ! + ! Ei(x) = exp(x) / x * SUM v(x,k), ! + ! k=0 ! + ! ! + ! where v(x,k) = k! / x^k ! + ! ! + ! To efficiently compute the terms inside the summation, we use that: ! + ! ! + ! v(x,k) = k / x * v(x,n -1), for k >= 1. ! + !------------------------------------------------------------------------------------! + vxk = 1.d0 + vsum = vxk + do_asymptote: do k=1,maxiter + !----- Update the current summation term. ----------------------------------------! + vxkm1 = vxk + vxk = vxkm1 * dble(k) / x + !---------------------------------------------------------------------------------! + ! This method can become degenerate for low x or lead to exceedinly large ! + ! values, in these cases, halt evaluation. ! + !---------------------------------------------------------------------------------! + if ( abs(vxkm1) < abs(vxk) .or. abs(vsum) > discard8) then + vsum = sign(discard8,vsum) + exit do_asymptote + end if + !----- Check for convergence. ----------------------------------------------------! + if ( any(abs(vxk) <= [ almost_zero8 * abs(vsum), tiny_num8] ) ) exit do_asymptote + !----- Update summation. ---------------------------------------------------------! + vsum = vsum + vxk + !---------------------------------------------------------------------------------! + end do do_asymptote + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! If the solution became degenerate, skip value. ! + !------------------------------------------------------------------------------------! + if (abs(vsum) == discard8) then + ei_asymptote = sign(discard8,vsum) + else + ei_asymptote = exp(x) * vsum / x + end if + !------------------------------------------------------------------------------------! + + + + + + !------------------------------------------------------------------------------------! + ! Pick the lowest absolute value as long as the sign is reasonable. ! + !------------------------------------------------------------------------------------! + if (all(abs([ei_series,ei_asymptote]) == discard8)) then + !----- Huge value, crash because this is iminent over-flow. ----------------------! + write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x + stop 'Exponential integral cannot be solved for large absolute x.' + !---------------------------------------------------------------------------------! + elseif (x < 0.d0) then + !---------------------------------------------------------------------------------! + ! Exponential integral is negative when x is negative, however, for some ! + ! values between -15 < x < -14, the solutions become numerically unstable. Check ! + ! for the most reasonable estimate. ! + !---------------------------------------------------------------------------------! + if (ei_series > 0.d0 .and. ei_asymptote > 0.d0) then + write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x + write(unit=*,fmt='(a,1x,es12.5)') 'Series expansion estimate = ',ei_series + write(unit=*,fmt='(a,1x,es12.5)') 'Asymptote expansion estimate = ',ei_asymptote + stop 'Exponential integral failed solving, another method might be needed.' + elseif (ei_series > 0.d0) then + eifun8 = ei_asymptote + elseif (ei_asymptote > 0.d0) then + eifun8 = ei_series + elseif (abs(ei_series) < abs(ei_asymptote)) then + eifun8 = ei_series else - sum = sum + term + eifun8 = ei_asymptote end if - end do asymploop - eifun8 = exp(x) * (1.d0 + sum) / x + !---------------------------------------------------------------------------------! + elseif (abs(ei_series) < abs(ei_asymptote)) then + eifun8 = ei_series + else + eifun8 = ei_asymptote + end if + !------------------------------------------------------------------------------------! end if + !---------------------------------------------------------------------------------------! return end function eifun8 @@ -2333,6 +2440,509 @@ end function eifun8 +!==========================================================================================! +!==========================================================================================! +! Heapsort is a robust and efficient sorting algorithm introduced by W64. The algorithm ! +! implemented here is built from the Wikipedia heapsort pseudocode, which is in turn based ! +! on K97. ! +! ! +! Williams, JWJ (1964). Algorithm 232 - Heapsort, Commun. ACM 7, 347-348. ! +! doi:10.1145/512274.512284 (W64). ! +! ! +! Knuth, D (1997). The Art of Computer Programming - volume 3: sort and searching. ! +! section 5.2.3. Sorting by selection (p. 144-155). ISBN 978-0-201-89685-5 (K97). ! +! ! +! Wikipedia link: https://en.wikipedia.org/wiki/Heapsort ! +!------------------------------------------------------------------------------------------! +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 (inner loop) + integer :: ilwr ! Lower index (inner loop) + integer :: iupr ! Upper index (inner loop) + integer :: olwr ! Lower index (outer loop) + integer :: oupr ! Upper index (outer loop) + real :: aux ! Placeholder for element swapping + !---------------------------------------------------------------------------------------! + + + !----- 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 + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Set initial guess of lower index ilwr to half the size of the vector and iupr to ! + ! the size of the vector. During the heap setting stage, ilwr will be reduced until it ! + ! becomes 0, and then we start decreasing iupr until it becomes 1, at which point the ! + ! vector becomes sorted. ! + !---------------------------------------------------------------------------------------! + olwr = nx/2 + 1 + oupr = nx+1 + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Main loop. ! + !---------------------------------------------------------------------------------------! + outer_loop: do + !------------------------------------------------------------------------------------! + ! Exit outer loop if we reach the upper bound has already reached 1. ! + !------------------------------------------------------------------------------------! + if (oupr == 2) exit outer_loop + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check whether we are in the heap setting phase or in the retirement-and-promotion ! + ! phase. ! + !------------------------------------------------------------------------------------! + if (olwr > 1) then + !----- Heap construction. --------------------------------------------------------! + olwr = olwr - 1 + !---------------------------------------------------------------------------------! + else + !---------------------------------------------------------------------------------! + ! Heap extraction. ! + !---------------------------------------------------------------------------------! + !----- Shift upper side down one step. -------------------------------------------! + oupr = oupr -1 + !----- Swap indices. -------------------------------------------------------------! + aux = xo(oupr) + xo(oupr) = xo(1) + xo(1) = aux + !---------------------------------------------------------------------------------! + end if + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Sift down step. ! + !------------------------------------------------------------------------------------! + i = olwr + inner_loop: do + !----- Find the lower and right elements. ----------------------------------------! + ilwr = 2 * i + iupr = ilwr + 1 + !---------------------------------------------------------------------------------! + + !----- Make sure we do not exceed the heap size. ---------------------------------! + if (iupr > oupr) exit inner_loop + !---------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------! + ! Test whether there is an upper element that is larger, and swap the order. ! + ! Make sure that the elements are bounded before testing vector elements, to ! + ! avoid segmentation violation. ! + !---------------------------------------------------------------------------------! + if (iupr < oupr) then + if (xo(ilwr) < xo(ilwr+1)) ilwr = ilwr + 1 + end if + !---------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------! + ! Test whether or not to swap elements. ! + !---------------------------------------------------------------------------------! + if (xo(i) < xo(ilwr)) then + aux = xo(i) + xo(i) = xo(ilwr) + xo(ilwr) = aux + i = ilwr + else + exit inner_loop + end if + !---------------------------------------------------------------------------------! + end do inner_loop + !------------------------------------------------------------------------------------! + 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). ! +! We test whether or not this is a trivial case that does not require solving the full ! +! equation. For the full equation, we use the approach by H02 to avoid floating point ! +! issues when solving roots. We further check whether or not 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. ! +! ! +! Higham, N. J., 2002: Accuracy and Stability of Numerical Algorithms. 2nd ed., Society ! +! for Industrial and Applied Mathematics, Philadelphia, PA, United States, ! +! doi:10.1137/1.9780898718027 (H02). ! +!------------------------------------------------------------------------------------------! +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 + 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) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Follow H02's approach to find the largest root (absolute value) from the ! + ! traditional quadratic equation, then derive the second root from the first one. ! + ! This is safe whenever b or c are non-zero. ! + !---------------------------------------------------------------------------------! + root1 = - (bquad + sign(sqrt(discr),bquad)) / ( 2. * aquad ) + root2 = cquad / ( aquad * root1 ) + !---------------------------------------------------------------------------------! + 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). ! +! We test whether or not this is a trivial case that does not require solving the full ! +! equation. For the full equation, we use the approach by H02 to avoid floating point ! +! issues when solving roots. We further check whether or not 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. ! +! ! +! Higham, N. J., 2002: Accuracy and Stability of Numerical Algorithms. 2nd ed., Society ! +! for Industrial and Applied Mathematics, Philadelphia, PA, United States, ! +! doi:10.1137/1.9780898718027 (H02). ! +!------------------------------------------------------------------------------------------! +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 + 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) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Follow H02's approach to find the largest root (absolute value) from the ! + ! traditional quadratic equation, then derive the second root from the first one. ! + ! This is safe whenever b or c are non-zero. ! + !---------------------------------------------------------------------------------! + root1 = - (bquad + sign(sqrt(discr),bquad)) / ( 2.d0 * aquad ) + root2 = cquad / ( aquad * root1 ) + !---------------------------------------------------------------------------------! + 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 +!==========================================================================================! +!==========================================================================================! + + return +end subroutine solve_quadratic8 +!==========================================================================================! +!==========================================================================================! + + + + + !==========================================================================================! !==========================================================================================! diff --git a/ED/src/dynamics/farq_katul.f90 b/ED/src/dynamics/farq_katul.f90 index 0e184f0b1..90e769971 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 !=======================================================================================! !=======================================================================================! @@ -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 !------------------------------------------------------------------------------------! @@ -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 @@ -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 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 0a5504f58..eddc7588e 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -4658,11 +4658,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). ----------------------------------! @@ -4681,6 +4680,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. !---------------------------------------------------------------------------------------! @@ -4937,23 +4938,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/memory/consts_coms.F90 b/ED/src/memory/consts_coms.F90 index 2345913b0..7fe069102 100644 --- a/ED/src/memory/consts_coms.F90 +++ b/ED/src/memory/consts_coms.F90 @@ -381,7 +381,7 @@ Module consts_coms !---------------------------------------------------------------------------------------! real, parameter :: srtwo = 1.414213562373095 ! Square root of 2. [ ---] real, parameter :: srthree = 1.732050807568877 ! Square root of 3. [ ---] - real, parameter :: sqrt2o2 = 0.5 * srtwo ! ½ Square root of 2. [ ---] + real, parameter :: sqrt2o2 = 0.5 * srtwo ! 1/2 * Square root of 2. [ ---] real, parameter :: srtwoi = 1./srtwo ! 1./ Square root of 2. [ ---] real, parameter :: srthreei = 1./srthree ! 1./ Square root of 3. [ ---] real, parameter :: onethird = 1./3. ! 1/3 [ ---] @@ -394,12 +394,12 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Universal constants ! !---------------------------------------------------------------------------------------! - real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m²/K^4] - real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m²kg/s²/K] - real, parameter :: t00 = 273.15 ! 0°C [ °C] - real, parameter :: rmol = 8.314510 ! Molar gas constant [ J/mol/K] - real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m³] - real, parameter :: volmoll = volmol*1e3 ! Molar volume at STP [ L] + real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m2/K4] + real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m2 kg/s2/K] + real, parameter :: t00 = 273.15 ! 0 degC [ degC] + real, parameter :: rmol = 8.314510 ! Molar gas constant [ J/mol/K] + real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m3] + real, parameter :: volmoll = volmol*1e3 ! Molar volume at STP [ L] !---------------------------------------------------------------------------------------! @@ -438,11 +438,11 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! General Earth properties ! !---------------------------------------------------------------------------------------! - real, parameter :: vonk = 0.40 ! Von Kármán constant [ ---] - real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s²] + real, parameter :: vonk = 0.40 ! Von Karman constant [ ---] + real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s2] real, parameter :: erad = 6370997. ! Earth radius [ m] real, parameter :: erad2 = 2.*erad ! Earth diameter [ m] - real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m²] + real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m2] real, parameter :: p00 = 1.e5 ! Reference pressure [ Pa] real, parameter :: prefsea = 101325. ! Reference sea level pressure [ Pa] real, parameter :: p00i = 1. / p00 ! 1/p00 [ 1/Pa] @@ -458,7 +458,7 @@ Module consts_coms ! third edition, Academic Press, Amsterdam, 418pp. (Chapters 3 and 10). ! ! ! ! Air diffusion properties. These properties are temperature-dependent in reality, ! - ! but for simplicity we assume them constants, using the value at 20°C. ! + ! but for simplicity we assume them constants, using the value at 20 degC. ! ! ! ! Thermal diffusivity - Computed from equation on page 32 of MU08; ! ! Kinematic viscosity - Computed from equation on page 32 of MU08; ! @@ -466,9 +466,9 @@ Module consts_coms ! 10.11 (MU08). ! ! These terms could be easily made function of temperature in the future if needed be. ! !---------------------------------------------------------------------------------------! - real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m²/s] + real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m2/s] real, parameter :: dth_diff = 0.007 ! Temperature dependency slope [ 1/K] - real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m²/s] + real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m2/s] real, parameter :: dkin_visc = 0.007 ! Temperature dependency slope [ 1/K] !---------------------------------------------------------------------------------------! @@ -516,8 +516,8 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Liquid water properties ! !---------------------------------------------------------------------------------------! - real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m³] - real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m³/kg] + real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m3] + real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cliq = 4.186e3 ! Liquid water specific heat (Cl) [ J/kg/K] real, parameter :: cliqi = 1./cliq ! Inverse of water heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! @@ -527,12 +527,12 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Ice properties ! !---------------------------------------------------------------------------------------! - real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m³] - real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m³/kg] - real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m³] - real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m³/kg] - real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m³] - real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m³/kg] + real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m3] + real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m3/kg] + real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m3] + real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m3/kg] + real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m3] + real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cice = 2.093e3 ! Ice specific heat (Ci) [ J/kg/K] real, parameter :: cicei = 1. / cice ! Inverse of ice heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! @@ -547,8 +547,8 @@ Module consts_coms real, parameter :: t3plei = 1./t3ple ! 1./T3 [ 1/K] real, parameter :: es3ple = 611.65685464 ! Vapour pressure at T3 (es3) [ Pa] real, parameter :: es3plei = 1./es3ple ! 1./es3 [ 1/Pa] - real, parameter :: epes3ple = ep * es3ple ! epsilon × es3 [ Pa kg/kg] - real, parameter :: rh2ot3ple = rh2o * t3ple ! Rv × T3 [ J/kg] + real, parameter :: epes3ple = ep * es3ple ! epsilon * es3 [ Pa kg/kg] + real, parameter :: rh2ot3ple = rh2o * t3ple ! Rv * T3 [ J/kg] real, parameter :: alli = 3.34e5 ! Lat. heat - fusion (Lf)[ J/kg] real, parameter :: alvl3 = 2.50e6 ! Lat. heat - vaporisation (Lv)[ J/kg] real, parameter :: alvi3 = alli + alvl3 ! Lat. heat - sublimation (Ls)[ J/kg] @@ -558,8 +558,8 @@ Module consts_coms real, parameter :: lvordry = alvl3 / rdry ! Lv/Ra [ K] real, parameter :: lvorvap = alvl3 / rh2o ! Lv/Rv [ K] real, parameter :: lsorvap = alvi3 / rh2o ! Ls/Rv [ K] - real, parameter :: lvt3ple = alvl3 * t3ple ! Lv × T3 [ K J/kg] - real, parameter :: lst3ple = alvi3 * t3ple ! Ls × T3 [ K J/kg] + real, parameter :: lvt3ple = alvl3 * t3ple ! Lv * T3 [ K J/kg] + real, parameter :: lst3ple = alvi3 * t3ple ! Ls * T3 [ K J/kg] real, parameter :: uiicet3 = cice * t3ple ! u at triple point, only ice [ J/kg] real, parameter :: uiliqt3 = uiicet3 + alli ! u at triple point, only liq. [ J/kg] real, parameter :: dcpvl = cph2o - cliq ! difference of sp. heat [ J/kg/K] @@ -671,16 +671,16 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Carbon-related unit conversions. ! !---------------------------------------------------------------------------------------! - real, parameter :: mol_2_umol = 1.e6 ! mol => µmol - real, parameter :: umol_2_mol = 1.e-6 ! µmol => mol - real, parameter :: umol_2_kgC = 1.20107e-8 ! µmol(CO2) => kg(C) - real, parameter :: Watts_2_Ein = 4.6e-6 ! W/m2 => mol/m²/s - real, parameter :: Ein_2_Watts = 1./Watts_2_Ein ! mol/m²/s => W/m2 - real, parameter :: kgC_2_umol = 1. / umol_2_kgC ! kg(C) => µmol(CO2) - real, parameter :: kgom2_2_tonoha = 10. ! kg(C)/m² => ton(C)/ha - real, parameter :: tonoha_2_kgom2 = 0.1 ! ton(C)/ha => kg(C)/m² - real, parameter :: umols_2_kgCyr = umol_2_kgC * yr_sec ! µmol(CO2)/s => kg(C)/yr - real, parameter :: kgCday_2_umols = kgC_2_umol / day_sec ! kg(C)/day => µmol(CO2)/s + real, parameter :: mol_2_umol = 1.e6 ! mol => umol + real, parameter :: umol_2_mol = 1.e-6 ! umol => mol + real, parameter :: umol_2_kgC = 1.20107e-8 ! umol(CO2) => kg(C) + real, parameter :: Watts_2_Ein = 4.6e-6 ! W/m2 => mol/m2/s + real, parameter :: Ein_2_Watts = 1./Watts_2_Ein ! mol/m2/s => W/m2 + real, parameter :: kgC_2_umol = 1. / umol_2_kgC ! kg(C) => umol(CO2) + real, parameter :: kgom2_2_tonoha = 10. ! kg(C)/m2 => ton(C)/ha + real, parameter :: tonoha_2_kgom2 = 0.1 ! ton(C)/ha => kg(C)/m2 + real, parameter :: umols_2_kgCyr = umol_2_kgC * yr_sec ! umol(CO2)/s => kg(C)/yr + real, parameter :: kgCday_2_umols = kgC_2_umol / day_sec ! kg(C)/day => umol(CO2)/s !---------------------------------------------------------------------------------------! #endif diff --git a/ED/src/utils/numutils.f90 b/ED/src/utils/numutils.f90 index 22c830115..e8930d092 100644 --- a/ED/src/utils/numutils.f90 +++ b/ED/src/utils/numutils.f90 @@ -431,11 +431,12 @@ end subroutine cumsum !==========================================================================================! ! This subroutine is the double precision version of the linear system solver above. ! ! It will solve the linear system AA . X = Y for given AA and Y, using the Gaussian ! -! elimination method with partial pivoting and back-substitution. This subroutine is ! -! based on: ! +! elimination method with partial pivoting and back-substitution. This subroutine builds ! +! on the algorithm described in P92, but adopting a vector-based approach compatible with ! +! modern Fortran. ! ! ! ! Press, W. H., S. A. Teukolsky, W. T. Vetterling, B. P. Flannery: 1992. Numerical recipes ! -! in Fortran 77. Cambridge University Press. ! +! in Fortran 77. Cambridge University Press (P92). ! !------------------------------------------------------------------------------------------! subroutine lisys_solver8(nsiz,AA,Y,X,sing) implicit none @@ -536,6 +537,8 @@ end subroutine lisys_solver8 + + !==========================================================================================! !==========================================================================================! ! EIFUN8 -- This function computes the exponential integral function, defined by ! @@ -544,35 +547,49 @@ end subroutine lisys_solver8 ! | exp(t) ! ! Ei(x) = | -------- dt ! ! _| t ! -! 0 ! +! -Inf ! ! ! -! This function is based on: ! +! This function checks for two approaches: series expansion, which typically works ! +! best when x is small, and the asymptotic expansion, which typically works best when x is ! +! large. The approach selects the smallest result (in absolute numbers) as the most ! +! accurate method. Both the series expansion and the asymptotic expansion are provided in ! +! AS72. This approach also checks for some other edge cases, and ignores the results when ! +! the value is very negative. ! ! ! -! Press, W. H., S. A. Teukolsky, W. T. Vetterling, B. P. Flannery: 1992. Numerical recipes ! -! in Fortran 77. Cambridge University Press, section 6.3 p. 215-219. ! +! Reference: ! +! ! +! Abramowitz, M., and I. A. Stegun, Eds., 1972: Handbook of mathematical functions with ! +! formulas, graphs, and mathematical tables. 10th ed., No. 55, Applied Mathematics ! +! Series, National Bureau of Standards, Washington, DC, USA (AS72). ! ! ! -! with the difference that we also solve for negative numbers. Zero cannot be solved, so ! -! if this happens, or if the sought number would lead to infinity, we stop the model. ! !------------------------------------------------------------------------------------------! real(kind=8) function eifun8(x) - use consts_coms, only : euler_gam8 & ! intent(in) - , lnexp_min8 & ! intent(in) - , lnexp_max8 & ! intent(in) - , tiny_num8 ! ! intent(in) + use consts_coms, only : euler_gam8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 & ! intent(in) + , tiny_num8 & ! intent(in) + , almost_zero8 ! ! intent(in) implicit none !----- Arguments. ----------------------------------------------------------------------! real(kind=8), intent(in) :: x !----- Local variables. ----------------------------------------------------------------! - real(kind=8) :: sum - real(kind=8) :: term - real(kind=8) :: fact - real(kind=8) :: prev - real(kind=8) :: diter - integer :: iter + real(kind=8) :: usum + real(kind=8) :: uxk + real(kind=8) :: uxkm1 + real(kind=8) :: vsum + real(kind=8) :: vxkm1 + real(kind=8) :: vxk + real(kind=8) :: ei_series + real(kind=8) :: ei_asymptote + integer :: k !----- Local constants. ----------------------------------------------------------------! - real(kind=8), parameter :: powerlim = 1.5d+01 - real(kind=8), parameter :: converge = 1.0d-7 - integer , parameter :: maxiter = 100 + real(kind=8), parameter :: discard8 = 1.0d+36 + integer , parameter :: maxiter = 100 + !----- Polynomial coefficients. --------------------------------------------------------! + real(kind=8), dimension(4), parameter :: apoly = (/ 8.5733287401d+00, 1.8059015973d+01 & + , 8.6347608925d+00, 2.6777373430d-01 /) + real(kind=8), dimension(4), parameter :: bpoly = (/ 9.5733223454d+00, 2.5632956149d+01 & + , 2.1099653083d+01, 3.9584969228d+00 /) !---------------------------------------------------------------------------------------! @@ -582,63 +599,157 @@ real(kind=8) function eifun8(x) if (x == 0.d0) then !------------------------------------------------------------------------------------! ! Zero. This is a singularity and the user should never call it in this case. ! - ! That's sad, but we ought to quit this run and tell the user why the run crashed. ! - !------------------------------------------------------------------------------------! - call fatal_error('Exponential integral cannot be solved for x = 0.' & - ,'eifun8','numutils.f90') - elseif (x >= lnexp_max8) then - !----- Huge value, crash because this is iminent over-flow. -------------------------! - write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x - write(unit=*,fmt='(a,1x,es12.5)') 'Maximum acceptable X =',lnexp_max8 - call fatal_error('Exponential integral cannot be solved for x = 0.' & - ,'eifun8','numutils.f90') - elseif (abs(x) <= lnexp_min8) then - !----- Huge negative number, the result can be rounded to zero. ---------------------! + !------------------------------------------------------------------------------------! + stop 'Exponential integral cannot be solved for x = 0.' + !------------------------------------------------------------------------------------! + elseif (x <= lnexp_min8) then + !------------------------------------------------------------------------------------! + ! Huge negative value, the result can be set to zero. ! + !------------------------------------------------------------------------------------! eifun8 = 0.d0 - elseif (abs(x) <= tiny_num8) then - !----- The number is too close to zero, bypass iterative methods. -------------------! - eifun8 = euler_gam8 + log(abs(x)) - elseif (abs(x) <= powerlim) then - !------------------------------------------------------------------------------------! - ! Input x is small, so we use the power method. ! - !------------------------------------------------------------------------------------! - fact = 1.d0 - sum = 0.d0 - powerloop: do iter=1,maxiter - diter = dble(iter) - fact = fact * x / diter - term = fact / diter - sum = sum + term - !----- If the term is tiny, we have reached convergence, quit the loop. ----------! - if (abs(term) < converge * abs(sum)) exit powerloop - end do powerloop - eifun8 = euler_gam8 + log(abs(x)) + sum + !------------------------------------------------------------------------------------! + elseif (x <= -1.d0) then + !------------------------------------------------------------------------------------! + ! For negative values less than -1.0, we use the polynomial approximation ! + ! (Equation 5.1.56 of AS72), by taking that Ei(x) = - E1(-x). ! + !------------------------------------------------------------------------------------! + eifun8 = exp(x)/x & + * ( x * ( x * ( x * ( x - apoly(1) ) + apoly(2) ) - apoly(3) ) + apoly(4) ) & + / ( x * ( x * ( x * ( x - bpoly(1) ) + bpoly(2) ) - bpoly(3) ) + bpoly(4) ) + !------------------------------------------------------------------------------------! else !------------------------------------------------------------------------------------! - ! Input x is large, so we use the asymptotic approximation. ! - !------------------------------------------------------------------------------------! - sum = 0.d0 - term = 1.d0 - asymploop: do iter=1,maxiter - diter = dble(iter) - prev = term - term = term * diter / x - if (abs(term) < converge) then - !----- The term is tiny, we have reached convergence, quit the loop. ----------! - exit asymploop - elseif (abs(term) >= abs(prev)) then - !------------------------------------------------------------------------------! - ! Series is diverging, we are probably reaching round-off errors, we better ! - ! stop now. ! - !------------------------------------------------------------------------------! - sum = sum - prev - exit asymploop + ! Find both the series expansion and the asymptotic expansion, and pick the one ! + ! with the lowest absolute value. ! + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Series expansion: Equation 5.1.11 of AS72, by taking that Ei(x) = - E1(-x). ! + ! ! + ! Inf ! + ! Ei(x) = gamma + ln(x) + SUM u(x,k), ! + ! k=1 ! + ! ! + ! where u(x,k) = [ (-1)^k * (-x)^k / (k * k!) ] ! + ! ! + ! To efficiently compute the terms inside the summation, we use that: ! + ! ! + ! u(x,k) = x * (k -1) / k^2 * u(x,k-1), for k >= 2. ! + !------------------------------------------------------------------------------------! + uxk = x + usum = uxk + do_expansion: do k = 2, maxiter + !----- Update the current summation term. ----------------------------------------! + uxkm1 = uxk + uxk = x * dble( k - 1 ) / dble( k * k ) * uxkm1 + !----- Check for degenerate or very large estimate. ------------------------------! + if ( abs(uxk) > discard8 .or. abs(usum) > discard8) then + usum = sign(discard8,usum) + exit do_expansion + end if + !----- Check for convergence. ----------------------------------------------------! + if ( any(abs(uxk) <= [ almost_zero8 * abs(usum), tiny_num8] ) ) exit do_expansion + !----- Update summation. ---------------------------------------------------------! + usum = usum + uxk + !---------------------------------------------------------------------------------! + end do do_expansion + !----- Find the series solution. ----------------------------------------------------! + if ( abs(usum) == discard8) then + ei_series = sign(discard8,usum) + else + ei_series = euler_gam8 + log(abs(x)) + usum + end if + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Asymptote expansion: Equation 5.1.51 of AS72 by taking AS72's n=1 and ! + ! Ei(x) = -E1(-x)). ! + ! ! + ! Inf ! + ! Ei(x) = exp(x) / x * SUM v(x,k), ! + ! k=0 ! + ! ! + ! where v(x,k) = k! / x^k ! + ! ! + ! To efficiently compute the terms inside the summation, we use that: ! + ! ! + ! v(x,k) = k / x * v(x,n -1), for k >= 1. ! + !------------------------------------------------------------------------------------! + vxk = 1.d0 + vsum = vxk + do_asymptote: do k=1,maxiter + !----- Update the current summation term. ----------------------------------------! + vxkm1 = vxk + vxk = vxkm1 * dble(k) / x + !---------------------------------------------------------------------------------! + ! This method can become degenerate for low x or lead to exceedinly large ! + ! values, in these cases, halt evaluation. ! + !---------------------------------------------------------------------------------! + if ( abs(vxkm1) < abs(vxk) .or. abs(vsum) > discard8) then + vsum = sign(discard8,vsum) + exit do_asymptote + end if + !----- Check for convergence. ----------------------------------------------------! + if ( any(abs(vxk) <= [ almost_zero8 * abs(vsum), tiny_num8] ) ) exit do_asymptote + !----- Update summation. ---------------------------------------------------------! + vsum = vsum + vxk + !---------------------------------------------------------------------------------! + end do do_asymptote + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! If the solution became degenerate, skip value. ! + !------------------------------------------------------------------------------------! + if (abs(vsum) == discard8) then + ei_asymptote = sign(discard8,vsum) + else + ei_asymptote = exp(x) * vsum / x + end if + !------------------------------------------------------------------------------------! + + + + + + !------------------------------------------------------------------------------------! + ! Pick the lowest absolute value as long as the sign is reasonable. ! + !------------------------------------------------------------------------------------! + if (all(abs([ei_series,ei_asymptote]) == discard8)) then + !----- Huge value, crash because this is iminent over-flow. ----------------------! + write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x + stop 'Exponential integral cannot be solved for large absolute x.' + !---------------------------------------------------------------------------------! + elseif (x < 0.d0) then + !---------------------------------------------------------------------------------! + ! Exponential integral is negative when x is negative, however, for some ! + ! values between -15 < x < -14, the solutions become numerically unstable. Check ! + ! for the most reasonable estimate. ! + !---------------------------------------------------------------------------------! + if (ei_series > 0.d0 .and. ei_asymptote > 0.d0) then + write(unit=*,fmt='(a,1x,es12.5)') 'Attempted X = ',x + write(unit=*,fmt='(a,1x,es12.5)') 'Series expansion estimate = ',ei_series + write(unit=*,fmt='(a,1x,es12.5)') 'Asymptote expansion estimate = ',ei_asymptote + stop 'Exponential integral failed solving, another method might be needed.' + elseif (ei_series > 0.d0) then + eifun8 = ei_asymptote + elseif (ei_asymptote > 0.d0) then + eifun8 = ei_series + elseif (abs(ei_series) < abs(ei_asymptote)) then + eifun8 = ei_series else - sum = sum + term + eifun8 = ei_asymptote end if - end do asymploop - eifun8 = exp(x) * (1.d0 + sum) / x + !---------------------------------------------------------------------------------! + elseif (abs(ei_series) < abs(ei_asymptote)) then + eifun8 = ei_series + else + eifun8 = ei_asymptote + end if + !------------------------------------------------------------------------------------! end if + !---------------------------------------------------------------------------------------! return end function eifun8 @@ -651,8 +762,17 @@ end function eifun8 !==========================================================================================! !==========================================================================================! -! Heapsort is a robust and efficient sorting algorithm. For more details, check ! -! The Numerical Recipes Book (chapter 8). ! +! Heapsort is a robust and efficient sorting algorithm introduced by W64. The algorithm ! +! implemented here is built from the Wikipedia heapsort pseudocode, which is in turn based ! +! on K97. ! +! ! +! Williams, JWJ (1964). Algorithm 232 - Heapsort, Commun. ACM 7, 347-348. ! +! doi:10.1145/512274.512284 (W64). ! +! ! +! Knuth, D (1997). The Art of Computer Programming - volume 3: sort and searching. ! +! section 5.2.3. Sorting by selection (p. 144-155). ISBN 978-0-201-89685-5 (K97). ! +! ! +! Wikipedia link: https://en.wikipedia.org/wiki/Heapsort ! !------------------------------------------------------------------------------------------! subroutine heapsort(nx,xi,increase,xo) implicit none @@ -662,11 +782,12 @@ subroutine heapsort(nx,xi,increase,xo) 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 + integer :: i ! Counter (inner loop) + integer :: ilwr ! Lower index (inner loop) + integer :: iupr ! Upper index (inner loop) + integer :: olwr ! Lower index (outer loop) + integer :: oupr ! Upper index (outer loop) + real :: aux ! Placeholder for element swapping !---------------------------------------------------------------------------------------! @@ -685,13 +806,13 @@ subroutine heapsort(nx,xi,increase,xo) !---------------------------------------------------------------------------------------! - ! 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. ! + ! Set initial guess of lower index ilwr to half the size of the vector and iupr to ! + ! the size of the vector. During the heap setting stage, ilwr will be reduced until it ! + ! becomes 0, and then we start decreasing iupr until it becomes 1, at which point the ! + ! vector becomes sorted. ! !---------------------------------------------------------------------------------------! - l = nx/2 + 1 - ir = nx + olwr = nx/2 + 1 + oupr = nx+1 !---------------------------------------------------------------------------------------! @@ -700,79 +821,74 @@ subroutine heapsort(nx,xi,increase,xo) !---------------------------------------------------------------------------------------! outer_loop: do !------------------------------------------------------------------------------------! - ! Check whether we are in the hiring phase or in the retirement-and-promotion ! + ! Exit outer loop if we reach the upper bound has already reached 1. ! + !------------------------------------------------------------------------------------! + if (oupr == 2) exit outer_loop + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Check whether we are in the heap setting phase or in the retirement-and-promotion ! ! phase. ! !------------------------------------------------------------------------------------! - if (l > 1) then - !---------------------------------------------------------------------------------! - ! Still in hiring phase. ! - !---------------------------------------------------------------------------------! - l = l - 1 - aux = xo(l) + if (olwr > 1) then + !----- Heap construction. --------------------------------------------------------! + olwr = olwr - 1 !---------------------------------------------------------------------------------! else !---------------------------------------------------------------------------------! - ! In the retirement-and-promotion phase. ! + ! Heap extraction. ! !---------------------------------------------------------------------------------! - !----- 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 + !----- Shift upper side down one step. -------------------------------------------! + oupr = oupr -1 + !----- Swap indices. -------------------------------------------------------------! + aux = xo(oupr) + xo(oupr) = xo(1) + xo(1) = aux !---------------------------------------------------------------------------------! end if !------------------------------------------------------------------------------------! + !------------------------------------------------------------------------------------! - ! Whether in hiring phase or promotion phase, we here set up to sift down element ! - ! aux to its proper level. ! + ! Sift down step. ! !------------------------------------------------------------------------------------! - i = l - j = l+1 + i = olwr inner_loop: do - if (j > ir) exit inner_loop + !----- Find the lower and right elements. ----------------------------------------! + ilwr = 2 * i + iupr = ilwr + 1 + !---------------------------------------------------------------------------------! + + !----- Make sure we do not exceed the heap size. ---------------------------------! + if (iupr > oupr) exit inner_loop + !---------------------------------------------------------------------------------! - !----- Compare to the better underling. ------------------------------------------! - if (j < ir) then - if(xo(j) < xo(j+1)) j = j + 1 + + !---------------------------------------------------------------------------------! + ! Test whether there is an upper element that is larger, and swap the order. ! + ! Make sure that the elements are bounded before testing vector elements, to ! + ! avoid segmentation violation. ! + !---------------------------------------------------------------------------------! + if (iupr < oupr) then + if (xo(ilwr) < xo(ilwr+1)) ilwr = ilwr + 1 end if !---------------------------------------------------------------------------------! !---------------------------------------------------------------------------------! - ! Check whether to demote aux or not. ! + ! Test whether or not to swap elements. ! !---------------------------------------------------------------------------------! - if (aux < xo(j)) then - !----- Demote aux. ------------------------------------------------------------! - xo(i) = xo(j) - i = j - j = j + j - !------------------------------------------------------------------------------! + if (xo(i) < xo(ilwr)) then + aux = xo(i) + xo(i) = xo(ilwr) + xo(ilwr) = aux + i = ilwr else - !----- This is aux's level. Set j to terminate the sift-down. ----------------! - j = ir + 1 - !------------------------------------------------------------------------------! + exit inner_loop end if !---------------------------------------------------------------------------------! end do inner_loop !------------------------------------------------------------------------------------! - - !----- Put aux into its slot. -------------------------------------------------------! - xo(i) = aux - !------------------------------------------------------------------------------------! end do outer_loop !---------------------------------------------------------------------------------------! @@ -909,70 +1025,231 @@ end function fquant_mask + !==========================================================================================! !==========================================================================================! -! FUNCTION bpow01 -!\brief Safe power estimate to avoid floating point exceptions -!\author Marcos Longo 3 March 2021 -!\details This function to calculate power functions for numbers bounded between 0 and 1 -!! safely. It uses that y = x ** a = exp(a * ln(x)), and use the safe log limits -!! to avoid FPE errors. This function "should" work for values greater than 1 too, -!! but don't use if for negative numbers. +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! We test whether or not this is a trivial case that does not require solving the full ! +! equation. For the full equation, we use the approach by H02 to avoid floating point ! +! issues when solving roots. We further check whether or not 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. ! +! ! +! Higham, N. J., 2002: Accuracy and Stability of Numerical Algorithms. 2nd ed., Society ! +! for Industrial and Applied Mathematics, Philadelphia, PA, United States, ! +! doi:10.1137/1.9780898718027 (H02). ! !------------------------------------------------------------------------------------------! -real(kind=4) function bpow01(x,a) - use consts_coms, only : lnexp_min & ! intent(in) - , lnexp_max ! ! intent(in) +subroutine solve_quadratic(aquad,bquad,cquad,undef,root1,root2) + use consts_coms, only : tiny_num implicit none !----- Arguments. ----------------------------------------------------------------------! - real(kind=4), intent(in) :: x - real(kind=4), intent(in) :: a + 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) :: lnexp + real(kind=4) :: discr + 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 if x is greater than zero (if it is exactly zero, we cannot use the log ! - ! approach). ! - !---------------------------------------------------------------------------------------! - if (x > 0.) then - !----- Find the bounded term inside the exponential. --------------------------------! - lnexp = max(lnexp_min,min(lnexp_max,a * log(x))) + !---------------------------------------------------------------------------------------! + ! 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) + !---------------------------------------------------------------------------------! - !----- Report result. ---------------------------------------------------------------! - bpow01 = exp(lnexp) + !---------------------------------------------------------------------------------! + ! Follow H02's approach to find the largest root (absolute value) from the ! + ! traditional quadratic equation, then derive the second root from the first one. ! + ! This is safe whenever b or c are non-zero. ! + !---------------------------------------------------------------------------------! + root1 = - (bquad + sign(sqrt(discr),bquad)) / ( 2. * aquad ) + root2 = cquad / ( aquad * root1 ) + !---------------------------------------------------------------------------------! + 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 (x == 0. .and. a > 0.) then - !----- By definition 0^a = 0 as long as a > 0. --------------------------------------! - bpow01 = 0. + else if (b_offzero) then + !------------------------------------------------------------------------------------! + ! "aquad" is not zero, not a true quadratic equation. Single root. ! + !------------------------------------------------------------------------------------! + root1 = - cquad / bquad + root2 = undef !------------------------------------------------------------------------------------! else !------------------------------------------------------------------------------------! - ! Invalid input data, stop everything. ! + ! 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)' ) ' Invalid variables for bpow01!' - write(unit=*,fmt='(a)' ) '-----------------------------------------------' - write(unit=*,fmt='(a,1x,es12.5)') ' x = ',x - write(unit=*,fmt='(a,1x,es12.5)') ' a = ',a - write(unit=*,fmt='(a)' ) '-----------------------------------------------' - write(unit=*,fmt='(a)' ) ' If x >= 0., then make sure a >= 0. ' - write(unit=*,fmt='(a)' ) ' If x = 0., then make sure a > 0. ' - write(unit=*,fmt='(a)' ) ' Negative x values are not allowed. ' - write(unit=*,fmt='(a)' ) '-----------------------------------------------' - call fatal_error('Invalid values for x and/or a.','bpow01','numutils.f90') + 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 function bpow01 +end subroutine solve_quadratic +!==========================================================================================! +!==========================================================================================! + + + + + !==========================================================================================! !==========================================================================================! +! Sub-routine that solves the quadratic equation ( a * x**2 + b * x + c = 0). ! +! We test whether or not this is a trivial case that does not require solving the full ! +! equation. For the full equation, we use the approach by H02 to avoid floating point ! +! issues when solving roots. We further check whether or not 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. ! +! ! +! Higham, N. J., 2002: Accuracy and Stability of Numerical Algorithms. 2nd ed., Society ! +! for Industrial and Applied Mathematics, Philadelphia, PA, United States, ! +! doi:10.1137/1.9780898718027 (H02). ! +!------------------------------------------------------------------------------------------! +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 + 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) + !---------------------------------------------------------------------------------! + + !---------------------------------------------------------------------------------! + ! Follow H02's approach to find the largest root (absolute value) from the ! + ! traditional quadratic equation, then derive the second root from the first one. ! + ! This is safe whenever b or c are non-zero. ! + !---------------------------------------------------------------------------------! + root1 = - (bquad + sign(sqrt(discr),bquad)) / ( 2.d0 * aquad ) + root2 = cquad / ( aquad * root1 ) + !---------------------------------------------------------------------------------! + 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 +!==========================================================================================! +!==========================================================================================!