LCOV - code coverage report
Current view: | top level - physics/clubb/src/CLUBB_core - calc_roots.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 36 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 3 | 0.0 %

              Line data    Source code
           1             : !---------------------------------------------------------------------------
           2             : ! $Id$
           3             : !===============================================================================
           4             : module calc_roots
           5             :
           6             :   implicit none
           7             :
           8             :   public :: cubic_solve, &
           9             :             quadratic_solve, &
          10             :             cube_root
          11             :
          12             :   private ! Set Default Scope
          13             :
          14             :   contains
          15             :
          16             :   !=============================================================================
          17           0 :   function cubic_solve( nz, a_coef, b_coef, c_coef, d_coef ) &
          18           0 :   result( roots )
          19             :
          20             :     ! Description:
          21             :     ! Solve for the roots of x in a cubic equation.
          22             :     !
          23             :     ! The cubic equation has the form:
          24             :     !
          25             :     ! f(x) = a*x^3 + b*x^2 + c*x + d;
          26             :     !
          27             :     ! where a /= 0.  When f(x) = 0, the cubic formula is used to solve:
          28             :     !
          29             :     ! a*x^3 + b*x^2 + c*x + d = 0.
          30             :     !
          31             :     ! The cubic formula is also called Cardano's Formula.
          32             :     !
          33             :     ! The three solutions for x are:
          34             :     !
          35             :     ! x(1) = -(1/3)*(b/a) + ( S + T );
          36             :     ! x(2) = -(1/3)*(b/a) - (1/2) * ( S + T ) + (1/2)i * sqrt(3) * ( S - T );
          37             :     ! x(3) = -(1/3)*(b/a) - (1/2) * ( S + T ) - (1/2)i * sqrt(3) * ( S - T );
          38             :     !
          39             :     ! where:
          40             :     !
          41             :     ! S = ( R + sqrt( D ) )^(1/3); and
          42             :     ! T = ( R - sqrt( D ) )^(1/3).
          43             :     !
          44             :     ! The determinant, D, is given by:
          45             :     !
          46             :     ! D = R^2 + Q^3.
          47             :     !
          48             :     ! The values of R and Q relate back to the a, b, c, and d coefficients:
          49             :     !
          50             :     ! Q = ( 3*(c/a) - (b/a)^2 ) / 9; and
          51             :     ! R = ( 9*(b/a)*(c/a) - 27*(d/a) - 2*(b/a)^3 ) / 54.
          52             :     !
          53             :     ! When D < 0, there are three unique, real-valued roots.  When D = 0, there
          54             :     ! are three real-valued roots, but one root is a double root or a triple
          55             :     ! root.  When D > 0, there is one real-valued root and there are two roots
          56             :     ! that are complex conjugates.
          57             :
          58             :     ! References:
          59             :     ! http://mathworld.wolfram.com/CubicFormula.html
          60             :     !-----------------------------------------------------------------------
          61             :
          62             :     use constants_clubb, only: &
          63             :         three,     & ! Constant(s)
          64             :         two,       &
          65             :         one_half,  &
          66             :         one_third
          67             :
          68             :     use clubb_precision, only: &
          69             :         core_rknd    ! Variable(s)
          70             :
          71             :     implicit none
          72             :
          73             :     ! Input Variables
          74             :     integer, intent(in) :: &
          75             :       nz    ! Number of vertical levels
          76             :
          77             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
          78             :       a_coef, & ! Coefficient a (of x^3) in a*x^3 + b*x^2 + c^x + d = 0    [-]
          79             :       b_coef, & ! Coefficient b (of x^2) in a*x^3 + b*x^2 + c^x + d = 0    [-]
          80             :       c_coef, & ! Coefficient c (of x) in a*x^3 + b*x^2 + c^x + d = 0      [-]
          81             :       d_coef    ! Coefficient d in a*x^3 + b*x^2 + c^x + d = 0             [-]
          82             :
          83             :     ! Return Variables
          84             :     complex( kind = core_rknd ), dimension(nz,3) :: &
          85             :       roots    ! Roots of x that satisfy a*x^3 + b*x^2 + c*x + d = 0       [-]
          86             :
          87             :     ! Local Variables
          88             :     real( kind = core_rknd ), dimension(nz) :: &
          89           0 :       cap_Q_coef,  & ! Coefficient Q in cubic formula     [-]
          90           0 :       cap_R_coef,  & ! Coefficient R in cubic formula     [-]
          91           0 :       determinant    ! Determinant D in cubic formula     [-]
          92             :
          93             :     complex( kind = core_rknd ), dimension(nz) :: &
          94           0 :       sqrt_det,   & ! Square root of determinant D in cubic formula     [-]
          95           0 :       cap_S_coef, & ! Coefficient S in cubic formula                    [-]
          96           0 :       cap_T_coef    ! Coefficient T in cubic formula                    [-]
          97             :
          98             :     complex( kind = core_rknd ), parameter :: &
          99             :       i_cmplx = ( 0.0_core_rknd, 1.0_core_rknd )  ! i = sqrt(-1)
         100             :
         101             :     complex( kind = core_rknd ) :: &
         102             :       sqrt_3,          & ! Sqrt 3 (complex data type)
         103             :       one_half_cmplx,  & ! 1/2 (complex data type)
         104             :       one_third_cmplx    ! 1/3 (complex data type)
         105             :
         106             :
         107             :     ! Declare some constants as complex data types in order to prevent
         108             :     ! data-type conversion warning messages.
         109           0 :     sqrt_3 = cmplx( sqrt( three ), kind = core_rknd )
         110           0 :     one_half_cmplx = cmplx( one_half, kind = core_rknd )
         111           0 :     one_third_cmplx = cmplx( one_third, kind = core_rknd )
         112             :
         113             :     ! Find the value of the coefficient Q; where
         114             :     ! Q = ( 3*(c/a) - (b/a)^2 ) / 9.
         115             :     cap_Q_coef = ( three * (c_coef/a_coef) - (b_coef/a_coef)**2 ) &
         116           0 :                  / 9.0_core_rknd
         117             :
         118             :     ! Find the value of the coefficient R; where
         119             :     ! R = ( 9*(b/a)*(c/a) - 27*(d/a) - 2*(b/a)^3 ) / 54.
         120             :     cap_R_coef = ( 9.0_core_rknd * (b_coef/a_coef) * (c_coef/a_coef) &
         121             :                    - 27.0_core_rknd * (d_coef/a_coef) &
         122           0 :                    - two * (b_coef/a_coef)**3 ) / 54.0_core_rknd
         123             :
         124             :     ! Find the value of the determinant D; where
         125             :     ! D = R^2 + Q^3.
         126           0 :     determinant = cap_Q_coef**3 + cap_R_coef**2
         127             :
         128             :     ! Calculate the square root of the determinant.  This will be a complex
         129             :     ! number.
         130           0 :     sqrt_det = sqrt( cmplx( determinant, kind = core_rknd ) )
         131             :
         132             :     ! Find the value of the coefficient S; where
         133             :     ! S = ( R + sqrt( D ) )^(1/3).
         134             :     cap_S_coef &
         135           0 :     = ( cmplx( cap_R_coef, kind = core_rknd ) + sqrt_det )**one_third_cmplx
         136             :
         137             :     ! Find the value of the coefficient T; where
         138             :     ! T = ( R - sqrt( D ) )^(1/3).
         139             :     cap_T_coef &
         140           0 :     = ( cmplx( cap_R_coef, kind = core_rknd ) - sqrt_det )**one_third_cmplx
         141             :
         142             :     ! Find the values of the roots.
         143             :     ! This root is always real-valued.
         144             :     ! x(1) = -(1/3)*(b/a) + ( S + T ).
         145             :     roots(:,1) &
         146             :     = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
         147           0 :       + ( cap_S_coef + cap_T_coef )
         148             :
         149             :     ! This root is real-valued when D < 0 (even though the square root of the
         150             :     ! determinant is a complex number), as well as when D = 0 (when it is part
         151             :     ! of a double or triple root).  When D > 0, this root is a complex number.
         152             :     ! It is the complex conjugate of roots(3).
         153             :     ! x(2) = -(1/3)*(b/a) - (1/2) * ( S + T ) + (1/2)i * sqrt(3) * ( S - T ).
         154             :     roots(:,2) &
         155             :     = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
         156             :       - one_half_cmplx * ( cap_S_coef + cap_T_coef ) &
         157           0 :       + one_half_cmplx * i_cmplx * sqrt_3 * ( cap_S_coef - cap_T_coef )
         158             :
         159             :     ! This root is real-valued when D < 0 (even though the square root of the
         160             :     ! determinant is a complex number), as well as when D = 0 (when it is part
         161             :     ! of a double or triple root).  When D > 0, this root is a complex number.
         162             :     ! It is the complex conjugate of roots(2).
         163             :     ! x(3) = -(1/3)*(b/a) - (1/2) * ( S + T ) - (1/2)i * sqrt(3) * ( S - T ).
         164             :     roots(:,3) &
         165             :     = - one_third_cmplx * cmplx( b_coef/a_coef, kind = core_rknd ) &
         166             :       - one_half_cmplx * ( cap_S_coef + cap_T_coef ) &
         167           0 :       - one_half_cmplx * i_cmplx * sqrt_3 * ( cap_S_coef - cap_T_coef )
         168             :
         169             :
         170           0 :     return
         171             :
         172           0 :   end function cubic_solve
         173             :
         174             :   !=============================================================================
         175           0 :   function quadratic_solve( nz, a_coef, b_coef, c_coef ) &
         176           0 :   result( roots )
         177             :
         178             :     ! Description:
         179             :     ! Solve for the roots of x in a quadratic equation.
         180             :     !
         181             :     ! The equation has the form:
         182             :     !
         183             :     ! f(x) = a*x^2 + b*x + c;
         184             :     !
         185             :     ! where a /= 0.  When f(x) = 0, the quadratic formula is used to solve:
         186             :     !
         187             :     ! a*x^2 + b*x + c = 0.
         188             :     !
         189             :     ! The two solutions for x are:
         190             :     !
         191             :     ! x(1) = ( -b + sqrt( b^2 - 4*a*c ) ) / (2*a); and
         192             :     ! x(2) = ( -b - sqrt( b^2 - 4*a*c ) ) / (2*a).
         193             :     !
         194             :     ! The determinant, D, is given by:
         195             :     !
         196             :     ! D = b^2 - 4*a*c.
         197             :     !
         198             :     ! When D > 0, there are two unique, real-valued roots.  When D = 0, there
         199             :     ! are two real-valued roots, but they are a double root.  When D < 0, there
         200             :     ! there are two roots that are complex conjugates.
         201             :
         202             :     ! References:
         203             :     !-----------------------------------------------------------------------
         204             :
         205             :     use constants_clubb, only: &
         206             :         four, & ! Constant(s)
         207             :         two
         208             :
         209             :     use clubb_precision, only: &
         210             :         core_rknd    ! Variable(s)
         211             :
         212             :     implicit none
         213             :
         214             :     ! Input Variables
         215             :     integer, intent(in) :: &
         216             :       nz    ! Number of vertical levels
         217             :
         218             :     real( kind = core_rknd ), dimension(nz), intent(in) :: &
         219             :       a_coef, & ! Coefficient a (of x^2) in a*x^2 + b*x + c = 0    [-]
         220             :       b_coef, & ! Coefficient b (of x) in a*x^2 + b*x + c = 0      [-]
         221             :       c_coef    ! Coefficient c in a*x^2 + b*x + c = 0             [-]
         222             :
         223             :     ! Return Variables
         224             :     complex( kind = core_rknd ), dimension(nz,2) :: &
         225             :       roots    ! Roots of x that satisfy a*x^2 + b*x + c = 0       [-]
         226             :
         227             :     ! Local Variables
         228             :     real( kind = core_rknd ), dimension(nz) :: &
         229           0 :       determinant    ! Determinant D in quadratic formula     [-]
         230             :
         231             :     complex( kind = core_rknd ), dimension(nz) :: &
         232           0 :       sqrt_det    ! Square root of determinant D in quadratic formula     [-]
         233             :
         234             :
         235             :     ! Find the value of the determinant D; where
         236             :     ! D = b^2 - 4*a*c.
         237           0 :     determinant = b_coef**2 - four * a_coef * c_coef
         238             :
         239             :     ! Calculate the square root of the determinant.  This will be a complex
         240             :     ! number.
         241           0 :     sqrt_det = sqrt( cmplx( determinant, kind = core_rknd ) )
         242             :
         243             :     ! Find the values of the roots.
         244             :     ! This root is real-valued when D > 0, as well as when D = 0 (when it is
         245             :     ! part of a double root).  When D < 0, this root is a complex number.  It is
         246             :     ! the complex conjugate of roots(2).
         247             :     ! x(1) = ( -b + sqrt( b^2 - 4*a*c ) ) / (2*a); and
         248             :     roots(:,1) = ( -cmplx( b_coef, kind = core_rknd ) + sqrt_det ) &
         249           0 :                  / cmplx( two * a_coef, kind = core_rknd )
         250             :
         251             :     ! This root is real-valued when D > 0, as well as when D = 0 (when it is
         252             :     ! part of a double root).  When D < 0, this root is a complex number.  It is
         253             :     ! the complex conjugate of roots(1).
         254             :     ! x(2) = ( -b - sqrt( b^2 - 4*a*c ) ) / (2*a).
         255             :     roots(:,2) = ( -cmplx( b_coef, kind = core_rknd ) - sqrt_det ) &
         256           0 :                  / cmplx( two * a_coef, kind = core_rknd )
         257             :
         258             :
         259           0 :     return
         260             :
         261           0 :   end function quadratic_solve
         262             :
         263             :   !=============================================================================
         264           0 :   function cube_root( x )
         265             :
         266             :     ! Description:
         267             :     ! Calculates the cube root of x.
         268             :     !
         269             :     ! When x >= 0, this code simply calculates x^(1/3).  When x < 0, this code
         270             :     ! uses x^(1/3) = -|x|^(1/3).  This eliminates numerical errors when the
         271             :     ! exponent of 1/3 is not treated as exactly 1/3, which would sometimes
         272             :     ! result in values of NaN.
         273             :     !
         274             :     ! References:
         275             :     !-----------------------------------------------------------------------
         276             :
         277             :     use constants_clubb, only: &
         278             :         one_third, & ! Constant(s)
         279             :         zero
         280             :
         281             :     use clubb_precision, only: &
         282             :         core_rknd    ! Variable(s)
         283             :
         284             :     implicit none
         285             :
         286             :     ! Input Variables
         287             :     real( kind = core_rknd ), intent(in) :: &
         288             :       x  ! Variable x
         289             :
         290             :     ! Return Variables
         291             :     real( kind = core_rknd ) :: &
         292             :       cube_root  ! Cube root of x
         293             :
         294             :
         295           0 :     if ( x >= zero ) then
         296           0 :        cube_root = x**one_third
         297             :     else ! x < 0
         298           0 :        cube_root = -abs(x)**one_third
         299             :     endif ! x >= 0
         300             :
         301             :
         302             :     return
         303             :
         304             :   end function cube_root
         305             :
         306             : !===============================================================================
         307             :
         308             : end module calc_roots

Generated by: LCOV version 1.14