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