text:
!> This module contains methods for various cubic equation of states The cubic
!! eos's are formulated using the m1 and the m2.
!! Suppoerted cubic EOS's
!!   SRK - Soave-Redlich-Kwong
!!   PR    - Peng Robinson
!!   VdW   - Van Der Waals
!!   RK    - Redlich-Kwong
!!   SW    - Schmidt and Wenzel
!!   PT    - Patel-Teja

module cubic
  use thermopack_constants, only: LIQPH, VAPPH, SINGLEPH
  implicit none
  save

  type :: cbBig ! Cubic solver parameters
    sequence
    real :: A,B,C
  end type cbBig

  public ::  &
       cbCalcZfac, cbCalcPressure, cbCalcFug, cbCalcEntropy,&
       cbGres,cbCalcEnthalpy, cbCalcPseudo, cbCalcFreeEnergy
  public :: cbCalcDerivatives_svol, cbCalcFres
  public :: cbDumpEosData
contains

  subroutine cbCalcDerivatives(nc,cbeos,T,p,Zfac)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent(in) :: T,p,Zfac
    real :: v
    ! General part all cubic EOS
    v = Zfac * kRgas*T/p
    call cbCalcDerivatives_svol(nc,cbeos,T,v)
  end subroutine cbCalcDerivatives

  ! Calculate derivatives given temperature [K] and specific volume [L/mol].
  subroutine cbCalcDerivatives_svol(nc,cbeos,T,v)
    use thermopack_constants, only: kRgas
    use eosdata
    use cubic_eos, only: cb_eos
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent(in) :: T, v
    ! Locals
    real :: n0,n1,n2,den,lnn, m, sumn
    real :: D,DD,B,C,R,mBR

    real :: f, fV, fB, fC, fVV, fVB, fVC, fBB,fBC,fCC
    real :: g, gV, gB, gVV, gVB, gBB
    real :: h, hV, hB, hC, hVV, hVB, hVC, hBB, hBC, hCC

    real :: m1B, m2B, mB,m1BB,m2BB,m1BC,m2BC
    real :: m1C, m2C, mC,m1CC,m2CC
    real :: mBB, mBC, mCC
    real :: n1B, n2B, n1V, n2V
    real :: n1C, n2C

    real :: n1BB,n2BB,n1BC,n2BC,n1CC,n2CC

    ! General part all cubic EOS
    sumn = cbeos%sumn
    n0 = v - cbeos%sumb
    n1 = v - cbeos%m1 * cbeos%sumb
    n2 = v - cbeos%m2 * cbeos%sumb
    den = n1*n2
    lnn = n2/n1

    if (abs(n0) <= 1D-20)  then
       print *, "In cbCalcDerivatives_svol: v, b = ", v, cbeos%sumb
       Call StopError &
            ("The (V-m1*B) or (V-m2*B) term has become negative for selected CubicEOS")
    else
       lnn = log(lnn)
       cbeos%pn = kRgas*T/n0
       cbeos%pa = (-1.0)/den
       cbeos%pb = kRgas*T/n0**2 &
            -cbeos%suma*((cbeos%m1 + cbeos%dm1dB*cbeos%sumb)*n2 &
            +(cbeos%m2+cbeos%dm2dB*cbeos%sumb)*n1)/den**2

       cbeos%pt = cbeos%sumn*kRgas/n0  !dp/dt=cbeos%pt+cbeos%pa*cbeos%at + cbeos%pb*cbeos%bt
       cbeos%pv = cbeos%sumn*(-kRgas)*T/n0**2+cbeos%suma*(n1+n2)/den**2  !dp/dv =cbeos%pv
       cbeos%pn = (kRgas*T)/n0 !dP/dn_i = cbeos%pn + cbeos%pa * cbeos%ai(i)
       !    +cbeos%pb * cbeos%bi(i) + cbeos%pc * cbeos%ci(i)

       cbeos%ffn = v/n0
       if (abs(cbeos%ffn) <= 1.0D-20) then
          print *, "In cbCalcDerivatives_svol: inputted T [K], v [L/mol] = ", T, v
          call StopError("Expression v/n <= 0 in selected cubic EOS")
       else

          select case (cbeos%subeosidx)
          case (cbVDW) !< Where m1 = m2 = 0.0
             DD = cbeos%suma/(krGas*T)
             cbeos%ffn = log(cbeos%ffn)

             cbeos%ffnv = -cbeos%sumb/(v*n0)
             cbeos%ffnb = 1.0/n0

             cbeos%ff = sumn*cbeos%ffn-DD/n1

             cbeos%ffa = -1.0/(kRgas*T*n1)
             cbeos%ffb = sumn/n0 - DD * cbeos%m1/(n1*n1)

             cbeos%fft = DD/(T*n1)
             cbeos%fftt = -2.0*DD/(n1*n1)

             cbeos%ffab = -1.0*cbeos%m1/(kRgas*T*n1*n1)
             cbeos%ffat = 1.0/(kRgas*T*T*n1)
             cbeos%ffbb = sumn/n0**2 - 2*DD*cbeos%m1/(n1*n1)
             cbeos%ffbt = DD * cbeos%m1 /(T*n1*n1)

             cbeos%ffv = -sumn*cbeos%sumb/(v*n0) + DD/(n1*n1)
             cbeos%ffvt = -DD/(n1*n1*T)

             cbeos%ffva = 1.0/(n1*n1*kRgas*T)
             cbeos%ffvv = -sumn/(v*v) + 1.0/(n0*n0) - 2.0*DD/(n1*n1*n1)
             cbeos%ffvb = -sumn/(n0*n0) + DD*cbeos%m1/(n1*n1*n1)

             !          else !< Rest of the cubic EOS family
          case (cbSRK,cbPR,cbSW,cbPT,eosLK,cspSRK,cspPR,cpaSRK,cpaPR)  ! 2 and 3 paramter eos
                                                         ! LK added when SRK is used to
                                                         ! generate initial values ....

!>
!! For three param like Patel-Teja, m1 = m1(sumb,sumc) nad m2 = m2(sumb,sumc)
!! For Schmidt and Wenzel, m1 = m1(sumc) and m2 = m2(sumc) - or m2(acfmix)
!!
!! Define reduced Helmholz energy: as F = -g - (amix(T)/T) * f
!!
!!    n = v - cbeos%sumb
!!    n1 = v - cbeos%m1 * cbeos%sumb
!!    n2 = v - cbeos%m2 * cbeos%sumb
!!    den = n1*n2
!!    lnn = log(n2/n1) or log n2 - log n1
!!    h = lnn
!!

!>
!!------------------------------
!! Michelsen | Thermopack
!!  p88-91   |
!!-----------------------------
!!     D     | cbeos%suma,
!!     Di    | cbeos%a(i)
!!     Dj    | cbeos%a(j)
!!     Dt    | cbeos%at
!!     Dtt   | cbeos%att
!!     Dij   | cbeos%aij(i,j)
!!     Dit   | cbeos%ait(i)
!!     B     | cbeos%sumb,
!!     Bi    | cbeos%b(i),
!!     Bij   | cbeos%bij = 0
!!-----------------------------

             m = cbeos%m1 - cbeos%m2

             D = cbeos%suma
             B = cbeos%sumb
             C = cbeos%sumc

             R = kRgas

             mBR = m*B*R

             h = lnn
             f = h/mBR

             g = log(n0) - log(v)

! Derivatives of helper functions

             if (cbeos%subeosidx == cbPT .or. cbeos%subeosidx == cbSW) then
! 1.st derivatives
                m1B = cbeos%dm1dB
                m2B = cbeos%dm2dB

                m1C = cbeos%dm1dC
                m2C = cbeos%dm2dC

                mB = m1B - m2B
                mC = m1C - m2C

                n1B = -m1B*B-cbeos%m1
                n2B = -m2B*B-cbeos%m2

                n1C = -m1C*B
                n2C = -m2C*B

! 2nd derivatives
                m1BB = cbeos%d2m1dB2
                m2BB = cbeos%d2m2dB2

                m1BC = cbeos%d2m1dBdC
                m2BC = cbeos%d2m2dBdC

                m1CC = cbeos%d2m1dC2
                m2CC = cbeos%d2m2dC2

                mBB = m1BB - m2BB
                mBC = m1BC - m2BC
                mCC = m1CC - m2CC

                n1BB = -m1BB * B - 2.0*m1B
                n2BB = -m2BB * B - 2.0*m2B
                n1BC = -m1BC*B - m1C
                n2BC = -m2BC*B - m2C
                n1CC = -m1CC*B
                n2CC = -m2CC*B
             else
                m1B = 0.0D0
                m2B = 0.0D0
                m1C = 0.0D0
                m2C = 0.0D0
                mB = 0.0D00
                mC = 0.0D00
                m1BB = 0.0D0
                m2BB = 0.0D0
                m1C = 0.0D0
                m2C = 0.0D0
                m1CC = 0.0D0
                m2CC = 0.0D0
                mBB = 0.0D00
                mBC = 0.0D00
                mCC = 0.0D00
                n1B = -cbeos%m1
                n2B = -cbeos%m2
                n1C = 0.0D00
                n2C = 0.0D00
                n1BB = 0.0D00
                n2BB = 0.0D00
                n1BC = 0.0D00
                n2BC = 0.0D00
                n1CC = 0.0D00
                n2CC = 0.0D00
             endif

! General for all cubic eos
             n1V = 1.0
             n2V = 1.0

             hV = 1.0/n2 - 1.0/n1
             hB = n2B/n2 - n1B/n1
             hC = n2C/n2 - n1C/n1

             hVV = (-1.0/(n2*n2)) + (1.0/(n1*n1))
             hVB = n1B/(n1*n1) - n2B/(n2*n2)
             hVC = n1C/(n1*n1) - n2C/(n2*n2)

             hBB = n2BB/n2 - n1BB/n1 &
                  - (n2B/n2)*(n2B/n2) &
                  + (n1B/n1)*(n1B/n1)

             hBC = n2BC/n2 - n1BC/n1 &
                  + (n1B*n1C)/(n1*n1) &
                  - (n2B*n2C)/(n2*n2)

             hCC =  n2CC/n2 - n1CC/n1 &
                  - (n2C/n2)*(n2C/n2) &
                  + (n1C/n1)*(n1C/n1)

             ! 1'st
             gV = 1.0/n0 - 1.0/V
             gB = -1.0/n0

             fV = hV/mBR
             fB = (1.0/mBR) * (hB - h*mB/m - h/B)

! Or alternative short form
!             fB = hB/mBR - f/B - f*mB/m
! From michelsen - 2 param with constant m1 and m2
!             fB = -(f+V*fV)/B

             fc = hC/mBR - f*mC/m

             ! 2'nd
             gVV = -1.0/(n0*n0)+1.0/(V*V)
             gVB = 1.0/(n0*n0)
             gBB = -gVB

             fVV = hVV/mBR

! From, mich - 2 param const m1 and m2
             fVB = -(2.0*fV+V*fVV)/B

! General alternative formulation
             fVB = (1.0/mBR) * (hVB - hV*mB/m - hV/B)

             fVC = (1.0/mBR) * (hVC - hV*mC/m)

             fBB = -fB*(1.0/B + mB/m) + (1.0/mBR)*(hBB - hB*(mB/m+1.0/B) &
                  + h*(1.0/B**2 + mB**2/m**2 - mBB/m))

             fBC = -fB*mC/m + (1.0/mBR) * (hBC - (hC*mB + h*mBC)/m &
                  + h*mB*mC/(m*m) - hC/B)

             fCC = (1.0/mBR) * (hCC - (2*hC*mC + h*mCC)/m &
                  + 2*h*mC*mC/(m*m))

             ! The reduced Helmholz function
             cbeos%ff = -g*sumn - (D/T)*f

             ! Now for the Helmholz-derivatives

             ! 1'st deriv
             cbeos%ffn = -g
             cbeos%ffv = -gV*sumn - (D/T)*fV
             cbeos%ffa = -f/T
             cbeos%ffb = -gB*sumn -(D/T)*fB
             cbeos%ffc = -(D/T) * fC
             cbeos%fft = D*f/(T*T)

             ! 2'nd deriv
             cbeos%fftt = -2.0*cbeos%fft/T

             cbeos%ffaa = 0.0D00
             cbeos%ffab = -fB/T
             cbeos%ffac = -fC/T
             cbeos%ffat = f/(T*T)

             cbeos%ffbb = -gBB*sumn - (D/T)*fBB
             cbeos%ffbc = -(D/T)*fBC
             cbeos%ffbt = D*fB/(T*T)

             cbeos%ffcc = -(D/T) * fCC
             cbeos%ffct = (D/(T*T))*fC

             cbeos%ffnv = -gV
             cbeos%ffnb = -gB
             cbeos%ffna = 0.0D0
             cbeos%ffnc = 0.0
             cbeos%ffnn = 0.0D0
             cbeos%ffnt = 0.0D0

             cbeos%ffvt =  D*fV/(T*T)
             cbeos%ffvv = -gVV*sumn - (D/T)* fVV
             cbeos%ffva = -fV/T
             cbeos%ffvb = -gVB*sumn - (D/T)* fVB
             cbeos%ffvc = -(D/T) * fVC

             ! Pressure derivative dp/dC

             cbeos%pc = -R * T * cbeos%ffvc

          end select ! Choice of eos
       endif ! ffn < 0
    endif ! n < 0
  end subroutine cbCalcDerivatives_svol

  !> Explicit calculation of pressure including the volume and temperature derivative
  ! using the current cubic equation of state.
  !
  ! \param T - Temperature [K]
  ! \param v - Specific volume [m3/kmol] (you read correctly, it is not SI)
  !
  !
  subroutine cbCalcPressure (nc,cbeos,T,v,z,p,dpdv,dpdt,d2pdv2,dpdz,recalculate)
    use cubic_eos, only: cb_eos
    use cbmix, only: cbCalcMixtureParams
    use cbHelm, only: cbPress, cbPrst, cbPv, cbPvv, cbPi
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent(in) :: t, v
    real, intent(in), dimension(nc) :: z
    real, intent(out) :: p
    real, optional, intent(out) :: dpdv, dpdt, d2pdv2
    real, dimension(nc), optional, intent(out) :: dpdz
    logical, optional, intent(in) :: recalculate
    ! Locals
    logical :: recalculate_loc
    !
    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    end if

    if (recalculate_loc) then
      call cbCalcMixtureParams(nc,cbeos,t,z)
      call cbCalcDerivatives_svol(nc, cbeos, T, v)
    end if

    p = cbPress(cbeos,T,V)
    if (present(dpdv)) then
      dpdv = cbPv(cbeos)
    endif
    if (present(dpdt)) then
      dpdt = cbPrsT(cbeos)
    endif
    if (present(d2pdv2)) then
      d2pdv2 = cbPvv(cbeos,T,v)
    endif
    if (present(dpdz)) then
      call cbPi(nc,cbeos,dpdz)
    endif
  end subroutine cbCalcPressure

  !-----------------------------------------------------------------------------
  !> Find a root z0 of the cubic polynomial
  !! \f[
  !!   f(z) = z^3 + p2 z^2 + p1 z + p0,
  !! \f]
  !!
  !! to solve for for the compressibility factor.
  !!
  !! \author MAG, 2013-09
  !!
  !! \param p2     - Polynomial coefficient
  !! \param p1     - Polynomial coefficient
  !! \param p0     - Polynomial coefficient
  !!
  !! \return z     - Root found, initial guess as input
  !! \return lconverged - True if solver converged
  !-----------------------------------------------------------------------------
  subroutine cb_solve_cubic_zfac(p2,p1,p0,z,lconverged)
    use numconstants, only: machine_prec
    implicit none
    ! Arguments:
    real,    intent(in)    :: p2,p1,p0
    real,    intent(inout) :: z
    logical, intent(out)   :: lconverged
    ! Locals:
    real               :: zn,dzn,fn,dfn,ddfn,discr,fnl
    integer            :: iiter,irelax
    integer, parameter :: maxiters=50,maxrelaxes=3
    real,    parameter :: dztol=100.0*machine_prec
    !---------------------------------------------------------------------------
    zn = z
    lconverged=.false.
    !
    do iiter = 1,maxiters
      ! Evaluate f and its derivatives.
      fn    = f(zn)
      dfn   = df(zn)
      ddfn  = ddf(zn)
      ! Find the discriminant of the quadratic equation for dz.
      discr = 1.0 - 2.0*fn*ddfn/dfn**2
      if (discr >= 0.0) then
        ! Find dz according to a third-order Newton method.
        dzn = -(dfn/ddfn)*(1.0-sqrt(discr))
      else
        ! The quadratic equation for dz has complex solutions.
        ! Revert to the regular, second-order Newton's method.
        dzn = -fn/dfn
      end if
      zn = zn + dzn
      !
      ! Check for convergence.
      if ( abs(dzn/zn) < dztol ) then
        lconverged = .true.
        exit
      end if
      !
      ! Relax. Shorten the step if it did not decrease f.
      dzn = -dzn
      do irelax = 1,maxrelaxes
        fnl = f(zn)
        if (abs(fnl) >= abs(fn)) then
        !  if (fnl*fn < 0.0) then
            ! Secant line search
        !    if (dzn > 0.0) then
        !      dzn = max(0.05*dzn,-fn*dzn/(fnl-fn))
        !    else
        !      dzn = min(0.05*dzn,-fn*dzn/(fnl-fn))
        !    endif
        !  else
            dzn = 0.5*dzn
        !  endif
          zn  = zn + dzn
        else
          exit
        end if
      end do
      !
    end do
    !
    z = zn
    !
  contains
    !
    pure function f(z)
      real             :: f
      real, intent(in) :: z
      f = p0 + z*( p1 + z*( p2 + z ))
    end function f
    !
    pure function df(z)
      real             :: df
      real, intent(in) :: z
      df = p1 + z*( 2.0*p2 + 3.0*z)
    end function df
    !
    pure function ddf(z)
      real             :: ddf
      real, intent(in) :: z
      ddf = 6.0*z + 2.0*p2
    end function ddf
    !
  end subroutine cb_solve_cubic_zfac

  !-----------------------------------------------------------------------------
  !> Having found the compressibility factor z = z0, that satisfies f(z0) = 0
  !! for the cubic polinomial
  !!
  !! \f[
  !!   f(z) = z^3 + p2 z^2 + p1 z + p0,
  !! \f]
  !!
  !! we can find the others analytically by solving the quadratic
  !! equation
  !!
  !! \f{align*}{
  !!   g(z) &= z^2 + (p2 + z0) z + (p2 z0 + z0^2 +p1 ) \\
  !!   & = z^2 +    q1     z +        q0\\
  !!   & = 0.\\
  !! \f}
  !!
  !! \author MAG, 2013-09
  !!
  !! \param p2     - Polynomial coefficient
  !! \param p1     - Polynomial coefficient
  !! \param z0     - Root of polynomial
  !!
  !! \return zl    - Smallest real root found
  !! \return zg    - Largest real root found
  !! \return iflag - Return code
  !!   ifalg = 1 : Converged to a single, real root
  !!   iflag = 2 : Converged to three, real though maybe degenerate, roots
  !-----------------------------------------------------------------------------
  subroutine cb_cubic_second_zfac(p2,p1,z0,zl,zg,iflag)
    implicit none
    ! Arguments:
    real,    intent(in)    :: p2,p1,z0
    real,    intent(out)   :: zl,zg
    integer, intent(out)   :: iflag
    ! Locals:
    real               :: q0,q1
    real               :: z1,z2,discr
    !---------------------------------------------------------------------------
    !
    ! Solve quadratic equation g(z) = 0 to find the other roots.
    q0    = z0**2 + p2*z0 + p1
    q1    = p2 + z0
    discr = q1**2 - 4*q0
    if (discr < 0.0) then
      ! The two other roots of f are a complex conjugate pair,
      ! so f has only one real root.
      iflag = 1
    else
      ! The two other roots of f are also real, so f has
      ! three real, though maybe degenerate, roots.
      iflag = 2
      discr = sqrt(discr)
      z1 = 0.5*( -q1 + discr )
      ! Refine by a single Newton step.
      z1 = z1-g(z1)/dg(z1)
      z2 = 0.5*( -q1 - discr )
      ! Refine by a single Newton step.
      z2 = z2-g(z2)/dg(z2)
    end if
    !
    select case (iflag)
    case (1)
      ! Found one real root.
      zl = z0
      zg = z0
    case (2)
      ! Found three real, though maybe degenerate, roots.
      zl = min(z0,z1,z2)
      zg = max(z0,z1,z2)
    end select
    !
  contains
    !
    pure function g(z)
      real             :: g
      real, intent(in) :: z
      g = q0 + z*( q1 + z )
    end function g
    !
    pure function dg(z)
      real             :: dg
      real, intent(in) :: z
      dg = q1 + 2.0*z
    end function dg
    !
  end subroutine cb_cubic_second_zfac

  !>     This subroutine solves a cubic polynomial
  !!
  !!     f(z) = z**3 + pp*z**2 + qq*z + rr = 0
  !!
  !!     with a third order method
  !!
  !!     Relaxation is necessarry in the critical region because of
  !!     both the first and second derivative approach zero
  !!
  !!     Original TPlib routine converted to use new data structure
  !!
  !!     \param pp - 1'st polynmominal coefficient
  !!     \param qq - 2'nd polynominal coefficient
  !!     \param rr - 3'rd polynominal coefficient
  !!
  !!     \param  z - Initial and return value for the compressibility factor
  !!
  !!     \return ifail - Return code
  !!       ifail = 0 : Converged, the gradient toward the solution is positive
  !!       ifail = 1 : Not converged after maximum number of iterations
  !!       ifail = 2 : Converged, the gradient toward the solution is negative
  !!       ifail = 3 : Not converged and the maximum number of relaxtions has been reach wi
  !!
  subroutine cbSolveCubicZfac(pp,qq,rr,z,ifail)
    implicit none
    ! Arguments:
    real,    intent (in)    :: pp,qq,rr
    real,    intent (inout) :: z
    integer, intent (out)   :: ifail
    ! Local constants:
    integer, parameter :: maxiter=50
    integer, parameter :: maxrelax=15
    real,    parameter :: eps=10.0*epsilon(0.0) ! 1.0E-18
    ! Local variables
    real :: fz,dfz,ddfz,dz,argum,fznew,znew
    integer ::iter,nrelax
    !
    fz=rr+z*(qq+z*(pp+z))
    iter=0
    ifail=1
    dz=1.0e10*eps
    !
    !do while (iter < maxiter .and. abs(dz) >= eps)
    do while (iter < maxiter .and. abs(dz/z) >= eps) ! Here we check for convergence in dz, but dz may be modified by the relaxing procedure!
      dfz=qq+z*(2.0*pp+3.0*z)
      ddfz=2.0*pp+6.0*z
      argum=1.0-2.0*fz*ddfz/dfz**2
      if (argum<0.0) argum=1.0
      argum=sqrt(argum)
      dz=-(fz/dfz)*2.0/(1.0+argum)
      znew=z+dz
      !if (abs(dz) < eps) then ! finish
      if (abs(dz/z) < eps) then ! finish
        if (dfz < 0.0) then
          ifail = 2
        else
          ifail = 0
        end if
        z = znew
      else ! abs(dz/z) > eps
        fznew = rr+znew*(qq+znew*(pp+znew))
        nrelax = 0
        do while (abs(fznew) >= abs(fz) .and. nrelax < maxrelax)
          nrelax = nrelax + 1
          if (nrelax.ge.maxrelax) then
            ifail = 3
          endif
          dz = dz*0.5
          znew = z + dz
          fznew = rr+znew*(qq+znew*(pp+znew))
        end do
        z = znew
        fz = fznew
        iter = iter + 1
      end if
    end do
    return
  end subroutine cbSolveCubicZfac

  subroutine cbSolveCubicZfacMinimumGibb (cbeos, T, p, pp, qq, rr, big, zfac, ifail, phase)
    use cubic_eos, only: cb_eos
    implicit none
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in) :: T, p
    real, intent (in) :: pp, qq, rr
    type(cbBig), intent(in) :: big
    real, intent (inout) :: zfac
    integer, intent (out) :: ifail, phase

    ! Local variables
    real :: zl, zg, gl, gg
    logical :: lconverged

    call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
    call cb_cubic_second_zfac(pp,qq,zfac,zl,zg,ifail)
    if (zl < big%B) then
      gl = 1.0E20
    else
      call cbGres(cbeos,T,p,zl,gl)
    end if
    if(zg < big%B) then
      gg = 1.0E20
    else
      call cbGres(cbeos,T,p,zg,gg)
    end if
    if (gl < gg) then ! Select the root with lowest Gibbs Free Energy
      zfac = zl
      phase = LIQPH
    else
      zfac = zg
      phase = VAPPH
    end if
    if (ifail == 1) then
      phase = SINGLEPH
    endif

  end subroutine cbSolveCubicZfacMinimumGibb

  !< Routine for calculation of the compressibility factor for a cubic equation of states
  !! This routine is based on the original TPlib routine CB_ZFAC
  !!
  !!  \param T - Temperature [K]
  !!  \param p - Pressure [Pa]
  !!  \param z - Composition vector (1..nc) [mol/mol]
  !!  \param iphase - Phase for which root to return.
  !!            1: Smallest possible (Liquid)
  !!            2: Largest possible (Vapour)
  !!  \param gflag - The TPlib "Gunder" flag
  !!            1: Normal
  !!            2: if the metastable maxima or minima of the z-factor are to be returned
  !!            3: If possibilities of having three roots, return the one having minimum
  !!               Gibbs free energy is to be returned - calls cbGres
  !! \param dZdt - Compressibility factor differential wrpt. temperature
  !! \param dZdp - Compressibility factor differential wrpt. pressure
  !! \param dZdz - Compressibility factor differential wrpt. mole numbers
  !! \param zfac - The compressibility factor
  !!
  !! We are here solving the cubic equation:
  !! z^3 + p*z^2 + q*z + r = 0
  !! Differentiating we get the quadratic equation:
  !! 3z^2 + 2p*z + q = 0
  !! With solution:
  !! z = ((-2p) +- sqrt(argum))/6
  !! argum = (2p)^2-12q
  !! if argum <= 0; The cubic equation will only have one root.
  !!
  subroutine cbCalcZfac(nc,cbeos,T,p,z,iphase,zfac,gflag_opt,dZdt,dZdp,dZdz,minGphase)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbmix, only: cbCalcMixtureParams
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,P,z(:)
    real, intent (out) :: zfac
    integer, intent(in) :: iphase, gflag_opt
    real, intent(out),optional :: dZdt,dZdp
    real, dimension(nc), intent(out),optional :: dZdz
    integer, optional, intent(out) :: minGphase
    ! Local variables
    integer :: iphase_loc
    real :: pp,qq,rr,z1,z2,f1,f2,f3,argum
    type(cbBig) :: big
    integer :: IFAIL1
    integer :: gunder
    real :: zFacGas
    !logical :: lconverged

    ! Calculate first the mixing rules
    call cbCalcMixtureParams (nc,cbeos,T,Z)

    ! Common for all cubic eos ?
    Big%A = cbeos%suma*p/(T*kRgas*t*kRgas)
    Big%B = cbeos%sumb*p/(T*kRgas)
    Big%C = cbeos%sumc*p/(T*kRgas)
    !
    zFacGas = max(2.0,Big%B*1.001)
    !
    !     Initiate ifail
    gunder = gflag_opt
    iphase_loc = iphase
    if (.not. (iphase == VAPPH .or. iphase == LIQPH) ) then
      if (iphase == SINGLEPH) then
        iphase_loc = VAPPH ! shouldn't matter if it's set to VAPPH or LIQPH
      else
        print *, "phase", iphase
        call stoperror('cbCalcZfac: You must give phase , 1 or 2 when' &
             //' calculating the compressibillity factor')
      end if
    endif

    cbeos%extrm(iphase_loc) = 0
    !
    !     Calculate the coefficients for third order equation based on m1 and m2
    pp = big%B*(-cbeos%m1 - cbeos%m2 - 1) - 1.0
    qq = big%B *((cbeos%m1+cbeos%m2)*(big%B+1) &
         + big%B * cbeos%m1 * cbeos%m2) + big%A
    rr = -(big%B**3 + big%B**2)*(cbeos%m1*cbeos%m2) - big%A*big%B
    argum = (2.0*pp)**2-12.0*qq
    f3 = rr+cbeos%single(1)%zcrit*(qq+cbeos%single(1)%zcrit*(pp+cbeos%single(1)%zcrit)) ! Beware for 3 param eq.

    if (cbeos%cubic_verbose) then
      write (*,*) 'zcrit ',cbeos%single%zcrit
      write (*,*) 'm1 ',cbeos%m1
      write (*,*) 'm2 ',cbeos%m2
      write (*,*) 'bigA ',big%A
      write (*,*) 'bigB ',big%B
      write (*,*) 'bigC ',big%C
      write (*,*) 'pp ',pp
      write (*,*) 'qq ',qq
      write (*,*) 'rr ',rr
      write (*,*) 'argum ',argum
      write (*,*) 'f3 ',f3
    endif

    if (argum <= 0.0) then !  Only one root.
      cbeos%nextrm(iphase_loc) = 0
      cbeos%nzfac(iphase_loc) = 1
      if (f3 > 0.0) then
        zfac = big%B ! Initial value from liquid side
      else
        zfac = zFacGas ! Initial value from ideal gas
      end if
      !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
      call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
    else ! Two or three roots
      cbeos%nextrm(iphase_loc) = 2
      z1 = ((-2.0)*pp-sqrt(argum))/6.0
      z2 = ((-2.0)*pp+sqrt(argum))/6.0
      f1 = rr+z1*(qq+z1*(pp+z1))
      f2 = rr+z2*(qq+z2*(pp+z2))
      if (f1*f2 > 0.0) then ! Only one root possible ?
        cbeos%nzfac(iphase_loc) = 1
        if (f3.gt.0.0d0) then
          if (iphase_loc.eq.1) then
            zfac = big%B
            !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
            call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
          else
            if (gunder == 2) then
              zfac = z2
              zfac = max(zfac,big%B*1.001) ! Ensure v > b
              cbeos%extrm(iphase_loc) = 1
            else
              zfac = big%B
              !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
              call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
            end if
          end if
        else ! f3 <= 0
          if (iphase_loc == 1) then
            if (gunder == 2) then
              zfac = Z1
              zfac = max(zfac,big%B*1.001) ! Ensure v > b
              cbeos%extrm(iphase_loc) = 1
            else
              zfac = zFacGas
              !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
              call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
            end if
          else
            zfac = zFacGas
            !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
            call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
          end if
        end if
      else ! f1*f2 <= 0 Three roots possible
        cbeos%nzfac(iphase_loc) = 3
        if (gunder == 3) then
          call cbSolveCubicZfacMinimumGibb (cbeos,T,p,pp,qq,rr,big,zfac,ifail1,minGphase)
        else ! Gunderflag 1 or 2
          if (iphase_loc == 1) then
            zfac = big%B
            !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
            call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
            if(zfac < big%B) then ! Re-initialize from vapour side ....
              zfac = zFacGas
              !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
              call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
            end if
          else
            zfac = zFacGas
            !call cb_solve_cubic_zfac(pp,qq,rr,zfac,lconverged)
            call cbSolveCubicZfac(pp,qq,rr,zfac,ifail1)
          end if
        end if
      end if
    end if

    call cbCalcZfacDiff(nc,cbeos,T,p,zfac,dZdt,dZdp,dZdz)

  end subroutine cbCalczfac

  !< Routine for calculation of the compressibility factor for a cubic equation of states
  !! This routine is based on the original TPlib routine CB_ZFAC
  !!
  !!  \param T - Temperature [K]
  !!  \param p - Pressure [Pa]
  !!  \param z - Composition vector (1..nc) [mol/mol]
  !!  \param iphase - Phase for which root to return.
  !!            1: Smallest possible (Liquid)
  !!            2: Largest possible (Vapour)
  !!  \param gflag - The TPlib "Gunder" flag
  !!            1: Normal
  !!            2: if the metastable maxima or minima of the z-factor are to be returned
  !!            3: If possibilities of having three roots, return the one having minimum
  !!               Gibbs free energy is to be returned - calls cbGres
  !! \param dZdt - Compressibility factor differential wrpt. temperature
  !! \param dZdp - Compressibility factor differential wrpt. pressure
  !! \param dZdz - Compressibility factor differential wrpt. mole numbers
  !! \param zfac - The compressibility factor
  !!
  !!
  !!
  subroutine cbCalcZfacDiff(nc,cbeos,T,p,zfac,dZdt,dZdp,dZdz)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbPi,cbPrst
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,P
    real, intent (in) :: zfac
    real, intent(out),optional :: dZdt,dZdp
    real, dimension(nc), intent(out),optional :: dZdz
    ! Locals
    real ::  dpdt,dvdt,dpdni(nc)

    ! Why is it not possible to disable this call?
    !if (present(dZdt) .or. present(dZdp) .or. present(dZdz)) then
      call cbCalcDerivatives(nc,cbeos,T,p,zfac)
    !endif

    if (present(dzdt)) then
      ! Start development of Z-derivatives here
      dpdt = cbPrst(cbeos)
      dvdt = -dpdt/cbeos%pv
      dzdt = p*dvdt/(kRgas*T) - zfac/T
    endif
    if (present(dzdp)) then
      dzdp = p/(cbeos%pv*kRgas*T) + zfac/p
    endif

    if (present(dZdz)) then
      call cbPi(nc,cbeos,dpdni)
      !dvdni = -dpdni/cbeos%pv
      dZdz =  -p/(kRgas*T)*dpdni/cbeos%pv - zfac
    endif
  end subroutine cbCalcZfacDiff

  !> The function to find the residual Gibbs free energy for the cubic equation of state
  !!
  !! This function is called if the "gflag" is set to 3 in the cbCalcZfac routine.
  !!
  !! \param T - Temperature [K]
  !! \param p - Pressure [Pa]
  !! \param zfac - Compressibility factor
  !! \param gr - The residual Gibbs free energy
  !! \param dgrdt - The residual Gibbs free energy temperature differential
  !! \param dgrdp - The residual Gibbs free energy pressure differential
  !!
  subroutine cbGres(cbeos, T, p, zfac, gr, dgrdt, dgrdp)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbFt
    implicit none
    class(cb_eos), intent(inout) :: cbeos
    real, intent(in) :: t,p,zfac
    real, intent(out) :: gr
    real, optional, intent(out) :: dgrdt, dgrdp

    !     local declarations
    real :: v,n,n1,n2,lnn,m,ares
    real :: ff, ffn, dfdt

    v = zfac*kRgas*T/P
    n = v-cbeos%sumb
    n1 = v-cbeos%m1*cbeos%sumb
    n2 = v-cbeos%m2*cbeos%sumb
    lnn = n2/n1

    if (lnn < 0.0) then
      call StopError ("Parameter zfac or v are out of range")
      gr = 0.0
      return
    end if

    lnn = log(lnn)

    ffn = v/n
    if (ffn <= 0) then
      call StopError("Parameter ffn has become negative in cbGres")
      gr = 0.0
      return
    endif
    ffn = log(ffn)

    if (cbeos%m1 == cbeos%m2) then
      ff = ffn - cbeos%suma /(kRgas*t*n1)
    else
      m = cbeos%m1-cbeos%m2
      ff = ffn - cbeos%suma*lnn/(m*cbeos%sumb*kRgas*t)
    endif

    if (zfac > 0.0) then
      ares = kRgas*t*(ff - log(zfac))
    else
      ares = kRgas*t*(ff - log(-zfac))
    end if

    gr =  ares + kRgas*t*(zfac-1.0)

    if (present(dgrdt)) then
      dfdt = cbFt( cbeos )
      dgrdt = kRgas*(T*dfdt+ff-log(abs(zfac)))
    endif

    if (present(dgrdp)) then
      dgrdp = (zfac-1.0)*kRgas*T/P
    endif
  end subroutine cbGres

  !> The subroutine finds the fugacity coefficients and derivatives
  !!
  !! \param T - Temperature [K]
  !! \param P - Pressure [Pa]
  !! \param Zfac - Compressibility factor
  !! \param dlnfdt - Fugacity coefficients differential wrpt. temperature
  !! \param dlnfdp - Fugacity coefficients differential wrpt. pressure
  !! \param dlnfdz - Fugacity coefficients differential wrpt. mole numbers
  !! \param res_fug - Fugacity coefficients
  !!
  !! \author Oivind Wilhelmsen
  !! \date 2012-06-13

  subroutine cbCalcFug(nc,cbeos,T,p,Zfac,res_fug,dlnfdt,dlnfdp,dlnfdz)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbFt, cbFit, cbPrst, cbPi, cbFij, cbFi
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,P,Zfac   ! Only compressibility without derivatives
    real, dimension(nc), intent (out) :: res_fug
    real, dimension(nc), optional, intent (out) :: dlnfdt, dlnfdp
    real, dimension(nc,nc), optional, intent (out) :: dlnfdz

    ! Local variables
    real :: dfdni(nc), Term, dpdni(nc), d2fdnidt(nc), d2fdnidnj(nc,nc), dpdt
    integer :: i,j

    call cbFi(nc,cbeos,dfdni)
    do i=1,nc
      Term = dfdni(i)-log(Zfac)

      if(Term .GT. 500.0) then
        Term = 500.0
      endif

      res_fug(i)=exp(Term)  ! The fugacity coefficients
    end do

    if(present(dlnfdt)) then  ! Temperature derivatives
      dpdt = cbPrst(cbeos)
      call cbFit(nc, cbeos, d2fdnidt)
      call cbPi(nc, cbeos, dpdni)
      dlnfdt=d2fdnidt+dpdni*dpdt/(kRgas*T*cbeos%pv)+1.0/T
    end if

    if(present(dlnfdp)) then ! Pressure derivative
      call cbPi(nc, cbeos, dpdni)
      dlnfdp=-dpdni/(kRgas*T*cbeos%pv)-1.0/P
    end if

    if(present(dlnfdz)) then ! Composition derivative
      call cbPi(nc, cbeos, dpdni)
      call cbFij(nc,cbeos,d2fdnidnj)
      do i=1,nc
        do j=1,nc
          dlnfdz(j,i)= d2fdnidnj(j,i)+dpdni(i)*dpdni(j)/(kRgas*T*cbeos%pv)+1.0
        end do
      end do
    end if

  end subroutine cbCalcFug

  !> The subroutine finds the entropy residual of cubic Equations of State
  !! optionally the derivatives
  !!
  !! \param T - Temperature [K]
  !! \param P - Pressure [Pa]
  !! \param Z - Composition [-]
  !! \param phase - Phase (1=liquid, 2=vapour)
  !! \param gflag_opt The TPlib "Gunder" flag
  !!            1: Normal
  !!            2: if the metastable maxima or minima of the z-factor are to be returned
  !!            3: If possibilities of having three roots, return the one having minimum
  !!               Gibbs free energy is to be returned - calls cbGres
  !! \param res_entropy - Residual entropy
  !! \param dsdt - Residual entropy differential wrpt. temperature
  !! \param dsdp - Residual entropy differential wrpt. pressure
  !! \param dsdz - Residual entropy differential wrpt. mole numbers
  !!
  !! \author Oivind Wilhelmsen
  !! \date 2012-06-13

  subroutine cbCalcEntropy(nc,cbeos,T,p,Z,phase,res_entropy,gflag_opt,dsdt,dsdp,dsdz)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbFtt, cbFit, cbPrst, cbFt, cbPi, cbFi
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,P
    real, dimension(nc), intent (in)  :: Z
    real, intent (out) :: res_entropy
    integer, intent(in) :: gflag_opt
    integer, intent(in) :: phase
    real, intent(out),optional :: dsdt,dsdp
    real, dimension(nc), intent(out),optional ::dsdz
    ! Local variables
    real :: dfdt, dpdt, d2fdt2, CVres, dpdni(nc), d2fdnidt(nc), dfdni(nc), Term ! Nc+3
    integer :: i
    real :: Zfac

    call cbCalcZfac(nc,cbeos,T,P,Z,phase,Zfac,gflag_opt)
    dpdt=cbPrst(cbeos)
    dfdt=cbFt(cbeos)
    res_entropy=(-kRgas)*(T*dfdt+cbeos%ff-log(abs(Zfac)))

    if(present(dsdt)) then  ! Temperature derivative
      d2fdt2 = cbFtt(cbeos)
      CVres=(-kRgas)*T*(T*d2fdt2+2.0*dfdt)
      dsdt=CVres/T-dpdt**2/cbeos%pv-kRgas/T
    end if

    if(present(dsdp)) then    !Pressure derivative
      dsdp=dpdt/cbeos%pv+kRgas/P
    end if

    if(present(dsdz)) then ! Composition derivative
      call cbFiT(nc, cbeos, d2fdnidt)
      call cbPi(nc, cbeos, dpdni)
      call cbFi(nc, cbeos, dfdni)
      do i=1,nc
        dsdz(i)=d2fdnidt(i)+dpdni(i)*dpdt/(kRgas*T*cbeos%pv)+1.0/T
        dsdz(i)=-dsdz(i)*kRgas*T**2  ! dHresdT

        term = dfdni(i)-log(abs(Zfac)) ! log(fug_i)
        dsdz(i)=dsdz(i)/T-kRgas*term ! dHresdT/T-Rlog(fug_i)
      end do
    end if

  end subroutine cbCalcEntropy

  !> The subroutine finds the enthalpy residual of cubic Equations of State
  !! optionally the derivatives
  !!
  !! \param T - Temperature [K]
  !! \param P - Pressure [Pa]
  !! \param Z - Composition [-]
  !! \param phase - Phase (1=liquid, 2=vapour)
  !! \param gflag_opt The TPlib "Gunder" flag
  !!            1: Normal
  !!            2: if the metastable maxima or minima of the z-factor are to be returned
  !!            3: If possibilities of having three roots, return the one having minimum
  !!               Gibbs free energy is to be returned - calls cbGres
  !! \param res_enthalpy - Residual enthalpy
  !! \param dhdt - Residual enthalpy differential wrpt. temperature
  !! \param dhdp - Residual enthalpy differential wrpt. pressure
  !! \param dhdz - Residual enthalpy differential wrpt. mole numbers
  !!
  !! \author Oivind Wilhelmsen
  !! \date 2012-06-14

  subroutine cbCalcEnthalpy(nc,cbeos,T,p,Z,phase,res_enthalpy,gflag_opt,dhdt,dhdp,dhdz)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbFt, cbFtt, cbFit, cbPrst, cbPi
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,P
    real, dimension(nc), intent (in)  :: Z
    integer, intent(in) :: gflag_opt
    real, intent (out) :: res_enthalpy
    integer, intent(in) :: phase
    real, intent(out), optional :: dhdt, dhdp
    real, dimension(nc), intent(out), optional ::dhdz

    ! Local variables
    real :: Ures, Ares, Sres, Cvres, dfdt, d2fdt2, dpdt, dpdni(nc), &
         d2fdnidt(nc)
    real :: Zfac

    call cbCalcZfac(nc,cbeos,T,P,Z,phase,Zfac,gflag_opt)

    dfdt = cbFt(cbeos)
    Ares = kRgas*T*(cbeos%ff-log(abs(Zfac)))
    Sres = (-kRgas)*(T*dfdt+cbeos%ff-log(abs(Zfac)))
    Ures = (-kRgas)*T**2*(dfdt)
    res_enthalpy = Ures + kRgas*T*(Zfac-1.0)
    if (cbeos%cubic_verbose) then
      write (*,*) 'cbCalcEnthalpy'
      write (*,*) 'dfdt: ',dfdt
      write (*,*) 'Ares: ',Ares
      write (*,*) 'Sres: ',Sres
      write (*,*) 'Ures: ',Ures
      write (*,*) 'Zfac: ',Zfac
    endif
    if (present(dhdt)) then ! Temperature derivative
      dpdt = cbPrst(cbeos)
      d2fdt2 = cbFtt( cbeos)
      CVres = (-kRgas)*T*(T*d2fdt2+2.0*dfdt)
      dhdt = CVres-T*dpdt**2/cbeos%pv-kRgas
    end if

    if (present(dhdp)) then ! Pressure derivative
      dpdt = cbPrst(cbeos)
      dhdp = Zfac*kRgas*T/P+T*dpdt/cbeos%pv
    end if

    if (present(dhdz)) then ! Composition derivative
      dpdt=cbPrst(cbeos)
      call cbFit(nc, cbeos, d2fdnidt)
      call cbPi(nc, cbeos, dpdni)
      dhdz=d2fdnidt+dpdni*dpdt/(kRgas*T*cbeos%pv)+1.0/T
      dhdz=(-dhdz)*kRgas*T**2
    end if

  end subroutine cbCalcEnthalpy

  !----------------------------------------------------------------------
  !> The subroutine finds the innerenergy residual of cubic Equations of State
  !! optionally the derivatives
  !!
  !! \param T The temperature [K]
  !! \param v The specific volume [m3/kmole]
  !! \param Z The mole fraction [-]
  !! \param dudt Temperature derivative [J/kmoleK]
  !! \param dudv Pressure derivative [J/m3]
  !!       0: No derivatives
  !!       1: Analytical derivatives
  !!
  !! \retval energy The internal energy with derivatives [-]
  !!
  !! \author Morten Hammer
  subroutine cbCalcInnerEnergy(nc,cbeos,T,v,Z,u,dudt,dudv,recalculate)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbmix, only: cbCalcMixtureParams
    use cbHelm, only: cbFtt, cbFvt, cbFt
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T,v,Z(1:nc)
    real, intent (out) :: u
    real, optional, intent(out) :: dudt, dudv
    logical, optional, intent(in) :: recalculate

    ! Local variables
    real :: Cvres, dfdt, d2fdt2, d2fdtdv
    logical :: recalculate_loc
    !

    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    end if

    if (recalculate_loc) then
      call cbCalcMixtureParams(nc,cbeos,T,Z)
      call cbCalcDerivatives_svol(nc,cbeos,T,v)
    end if

    dfdt = cbFt(cbeos)
    u = (-kRgas)*T**2*(dfdt)

    if (present(dudt)) then
      d2fdt2 = cbFtt(cbeos)
      CVres = (-kRgas)*T*(T*d2fdt2+2.0*dfdt)
      dudt = CVres
    endif

    if (present(dudv)) then
      d2fdtdv = cbFvt(cbeos)
      dudv = (-kRgas)*T**2*(d2fdtdv)
    endif

  end subroutine cbCalcInnerEnergy

  !> The subroutine calculates the pseudocritical temperature, pressure, acentric factor
  !! compressibility and volume as function of composition.
  !! The pseudo-critical point is the state where with constant Temperature:
  !! (dP/dV)=(d2P/dV^2)=0
  !!
  !! \param Z - Composition [-]
  !! \return Tpc - Pseudo critical temperature [K]
  !! \return Ppc - Pseudo critical pressure [Pa]
  !! \return Zpc - Pseudo critical compressibility [-]
  !! \return Vpc - Pseudo critical volume [m^3/kmole]
  !!
  !! \author Oivind Wilhelmsen
  !! \date 2012-06-14
  subroutine cbCalcPseudo(nc,cbeos,Z,Tpc,Ppc,Zpc,Vpc)
    use thermopack_constants, only: kRgas
    use eosdata
    use cubic_eos, only: cb_eos
    use cbmix, only: cbCalcMixtureParams
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: Z(1:nc)
    real, intent (out) :: Tpc, Ppc, Zpc, Vpc

    ! Local variables
    integer :: i, iter
    real :: ft, dft, ddft, Argum, dTpc, omegaA_mix, OmegaB_mix, zcrit_mix
    real :: zcomb, dzcomb, n, n1, n2, f, fd, fdd, argum1
    real, parameter :: Eps=1.0D-10
    integer, parameter :: Maxiter=20

    ! Initial values

    Tpc=0.0
    Ppc=0.0
    iter=0
    zcomb = 5.0
    dzcomb = 1.0E10

    do i=1,nc
      Tpc=Tpc+Z(i)*cbeos%single(i)%tc
      Ppc=Ppc+Z(i)*cbeos%single(i)%pc
    end do

    ! make sure m1 and m2 are calculated for pseudeo critical t and p
    !call cbCalcParameters (nc, comp, cbeos, Tpc, Z)

    ! First, find the Omega_B and Omega_B of the "mixture" which is such that
    ! (dP/dV)=(d2P/dV2)=0 in a critical point. This is not necessary for three-
    ! parameter cubic EoS because m1 and m2 do not depend on composition

    if (cbeos%subeosidx == cbSW) then !> For the Schmidt and Wenzel EOS
      do while (iter < maxiter .and. abs(dzcomb) > eps)
        n  = zcomb - 1.0
        n1 = zcomb - cbeos%m1
        n2 = zcomb - cbeos%m2
        f = n*n1**2+n*n1*n2+n*n2**2-n1**2*n2-n1*n2**2
        fd = 3.0*(n*n1+n*n2-n1*n2)
        fdd = 6.0*n
        argum1 = 1.0-2.0*f*fdd/fd**2
        if (argum1.gt.0.0) then
          argum1 = sqrt(argum1) ! changed sqrt(argum) to sqrt(argum1) - GS 2012-12-22
        else
          argum1 = 1.0
        end if

        dzcomb = (-2.0)*f/(fd*(1.0+argum1))
        zcomb = zcomb+dzcomb
        iter = iter + 1
      enddo

      ! if (iter == maxiter) call StopError ("Third order method failed to calculate zcrit")
      !  write (*,*) 'Calculation of ohmega/zcrit converged after ',iter-1,' iterations for CB  EOS'
      ! endif

      argum1 = n1*n2/(n*n1+n*n2)
      omegaB_mix = (1.0-argum1)/(zcomb-1.0)
      zcrit_mix = omegaB_mix*zcomb
      omegaA_mix = omegaB_mix*(1.0/n-omegaB_mix)*n1*n2
    elseif (cbeos%subeosidx == cbPT) then
       ! for now
      omegaA_mix=cbeos%single(1)%omegaA
      omegaB_mix=cbeos%single(1)%omegaB
      zcrit_mix=cbeos%single(1)%zcrit
    else
      omegaA_mix=cbeos%single(1)%omegaA
      omegaB_mix=cbeos%single(1)%omegaB
      zcrit_mix=cbeos%single(1)%zcrit
    end if

    ! reinitialize
    iter=0
    dTpc=20.0

    !  Start the third order Newton-Raphsons method to find the correspoding pseudo
    ! critical point

    do while(iter<Maxiter .AND. abs(dTpc)>Eps)
      iter=iter+1
      call cbCalcMixtureParams(nc,cbeos,Tpc,Z)
      ft=cbeos%suma-omegaA_mix/omegaB_mix*kRgas*Tpc*cbeos%sumb
      dft=cbeos%at-omegaA_mix/omegaB_mix*kRgas* ( cbeos%sumb + Tpc * cbeos%bt)
      ddft=cbeos%att
      Argum=(1.0D-2)*ft*ddft/dft**2
      if(Argum < 0.0) Argum=1.0
      Argum=sqrt(Argum)
      dTpc=(-ft)/dft*2.0/(1.0+Argum)
      Tpc=Tpc+dTpc
    end do

    if(iter.GE.Maxiter) then
      call StopError ("Pseudo-critical pressure exceeded Maxiter")
    end if

    Zpc=zcrit_mix
    Vpc=Zpc*cbeos%sumb/omegaB_mix

    Ppc=kRgas*Tpc/(Vpc-cbeos%sumb)-cbeos%suma/((Vpc-cbeos%m1*cbeos%sumb)*&
         (Vpc-cbeos%m2*cbeos%sumb))

  end subroutine cbCalcPseudo

!! ":"-delimited text-output with all data from the current EOS
!!
  subroutine cbDumpEosData(nc,cbeos,T, P, zfac)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbFi, cbFit
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos

    real, intent(in) :: T,P, zfac
    integer :: i,j
    real :: dpdt,dfdni(nc),dpdni,dpdnj,d2fdnidt(nc),d2fdnidnj

    write (*,*) 'Calculated parameters for selected EOS: ',cbeos%eosid
    write (*,*) 'Mixing rule:                            ',cbeos%mruleid

    write (*,*) 'm1  :',cbeos%m1
    write (*,*) 'm2  :',cbeos%m2
    do i=1,nc
       write (*,*) 'alphacorr(',i,') :',cbeos%single(i)%alphaCorrName
       write (*,*) 'omegaA(',i,')    :',cbeos%single(i)%omegaA
       write (*,*) 'omegaB(',i,')    :',cbeos%single(i)%omegaB
       write (*,*) 'omegaC(',i,')    :',cbeos%single(i)%omegaC
       write (*,*) 'zcrit(',i,')     :',cbeos%single(i)%zcrit
    enddo

    write (*,*) 'sumA :',cbeos%suma
    write (*,*) 'sumB :',cbeos%sumb
    write (*,*) 'sumC :',cbeos%sumc

    write (*,*) 'At :',cbeos%At
    write (*,*) 'Att :',cbeos%Att

    write (*,*) 'Bt :',cbeos%Bt
    write (*,*) 'Btt :',cbeos%Btt

! Desireable to keep the values in a single column.
! Use this data to import into for instance EXCEL to directly compare values
! old - new
!
    do i=1,nc
       write (*,*) 'ai(',i,') :',cbeos%ai(i)
       write (*,*) 'bi(',i,') :',cbeos%bi(i)
       write(*,*) 'ait(',i,') :',cbeos%ait(i)
       write (*,*) 'bit(',i,') :',cbeos%bit(i)

       do j=1,nc
          write(*,*) 'aij(',i,',',j,') :',cbeos%aij(i,j)
          write(*,*) 'bij(',i,',',j,') :',cbeos%bij(i,j)
! Could also include mixing rule interactionparameters - but these are really constants
! but could be useful to look at when doing regression ...
!
!          write(*,*) 'bij(',i,',',j,') :',cbeos%kij(i,j)
       end do
    end do

    write (*,*) 'pn :',cbeos%pn
    write (*,*) 'pa :',cbeos%pa
    write (*,*) 'pb :',cbeos%pb
    write (*,*) 'pc :',cbeos%pc
    write (*,*) 'pt :',cbeos%pt
    write (*,*) 'pv :',cbeos%pv
    write (*,*) 'ff :',cbeos%ff
    write (*,*) 'ffn :',cbeos%ffn
    write (*,*) 'ffa :',cbeos%ffa
    write (*,*) 'ffb :',cbeos%ffb
    write (*,*) 'ffc :',cbeos%ffc
    write (*,*) 'fft :',cbeos%fft
    write (*,*) 'ffv :',cbeos%ffv
    write (*,*) 'ffnv :',cbeos%ffnv
    write (*,*) 'ffna :',cbeos%ffna
    write (*,*) 'ffnb :',cbeos%ffnb
    write (*,*) 'ffnc :',cbeos%ffnc
    write (*,*) 'ffaa :',cbeos%ffaa
    write (*,*) 'ffab :',cbeos%ffab
    write (*,*) 'ffac :',cbeos%ffac
    write (*,*) 'ffat :',cbeos%ffat
    write (*,*) 'ffbb :',cbeos%ffbb
    write (*,*) 'ffbc :',cbeos%ffbc
    write (*,*) 'ffbt :',cbeos%ffbt
    write (*,*) 'ffcc :',cbeos%ffcc
    write (*,*) 'ffct :',cbeos%ffct
    write (*,*) 'ffnt :',cbeos%ffnt
    write (*,*) 'fftt :',cbeos%fftt
    write (*,*) 'ffvt :',cbeos%ffvt
    write (*,*) 'ffva :',cbeos%ffva
    write (*,*) 'ffvb :',cbeos%ffvb
    write (*,*) 'ffvc :',cbeos%ffvc
    write (*,*) 'ffvv :',cbeos%ffvv

    write (*,*) 'dm1dB  :',cbeos%dm1dB
    write (*,*) 'dm1dC  :',cbeos%dm1dC
    write (*,*) 'dm2dB  :',cbeos%dm2dB
    write (*,*) 'dm2dC  :',cbeos%dm2dC

    write (*,*) 'd2m1dB2  :',cbeos%d2m1dB2
    write (*,*) 'd2m1dC2  :',cbeos%d2m1dC2
    write (*,*) 'd2m2dB2  :',cbeos%d2m2dB2
    write (*,*) 'd2m2dC2  :',cbeos%d2m2dC2

    ! Other differentials needed in departure functions and fugacities
    call cbFi(nc, cbeos, dfdni)
    do i=1,nc
       write (*,*) 'dfdni(',i,') :',dfdni(i)
       write (*,*) 'ln(fug) :', dfdni(i) - log(zfac)
    enddo

    dpdt=cbeos%pt+cbeos%pa*cbeos%at
    write (*,*) 'dpdt :',dpdt
    call cbFit(nc, cbeos, d2fdnidt)
    do i=1,nc
       dpdni=cbeos%pn+cbeos%pa*cbeos%ai(i)+cbeos%pb*cbeos%bi(i)+cbeos%pc*cbeos%ci(i)

!       d2fdnidt =  &
!            +  cbeos%ffat * cbeos%ai(i) + cbeos%ffa*cbeos%ait(i) &
!            + (cbeos%ffbt + cbeos%ffab*cbeos%at)* cbeos%bi(i) &
!            + (cbeos%ffct + cbeos%ffac*cbeos%at)* cbeos%ci(i)

       write (*,*) 'dpdni(',i,') :',dpdni
       write (*,*) 'd2fdnidt(',i,') :',d2fdnidt(i)
       write (*,*) 'ln(dfug/dt) :', d2fdnidt(i)+dpdni*dpdt/(kRgas*T*cbeos%pv)+1.0/T
    end do

    do i=1,nc
       dpdni=cbeos%pn+cbeos%pa*cbeos%ai(i) &
            + cbeos%pb*cbeos%bi(i) &
            + cbeos%pc*cbeos%ci(i)

       write (*,*) 'dpdni(',i,') :',dpdni
       write (*,*) 'ln(dfug/dp) :', -dpdni/(kRgas*T*cbeos%pv)-1.0/P

    end do

    do i=1,nc
       dpdni=cbeos%pn+cbeos%pa*cbeos%ai(i) &
            +cbeos%pb*cbeos%bi(i) &
            +cbeos%pc*cbeos%ci(i)
       write (*,*) 'dpdni(',i,') :',dpdni
       do j=1,nc
          dpdnj=cbeos%pn+cbeos%pa*cbeos%ai(j) &
               + cbeos%pb*cbeos%bi(j) &
               +cbeos%pc*cbeos%ci(j)
          write (*,*) 'dpdnj(',j,') :',dpdnj
          d2fdnidnj = cbeos%ffnb*cbeos%bi(j) &
               + (cbeos%ffab*cbeos%bi(j) + cbeos%ffac*cbeos%ci(j)) * cbeos%ai(i) + cbeos%ffa*cbeos%aij(i,j) &
               + cbeos%ffb*cbeos%bij(i,j) & !Term used in Wong Sandler
               + (cbeos%ffnb + cbeos%ffab*cbeos%ai(j) + cbeos%ffbb*cbeos%bi(j) + cbeos%ffbc*cbeos%ci(j)) * cbeos%bi(i)  &
               + ( cbeos%ffac*cbeos%ai(j) + cbeos%ffbc*cbeos%bi(j) + cbeos%ffcc*cbeos%ci(j)) * cbeos%ci(i) +&
               cbeos%ffc*cbeos%cij(i,j)
          write (*,*) 'd2fdnidnj(',j,',',i,') :',d2fdnidnj
          write (*,*) 'dlnfdz(',j,',',i,') :', d2fdnidnj+dpdni*dpdnj/(kRgas*T*cbeos%pv)+1.0
       end do
    end do

! New - add also interaction parameters
! Put this at the end for now, (GEIR)

    write(*,*) "kij (i:row, j:col):"
    do i=1,nc
       do j=1,nc
          write(*,"(ES20.10)",advance="no") cbeos%kij(i,j)
       end do
       write(*,*)
    end do
    write(*,*)

  end subroutine cbDumpEosData

  !----------------------------------------------------------------------
  !> Subroutine finds Helmholtz free energy residual of cubic EoS
  !! optionally the derivatives
  !! with respect to T, v
  !!
  !! \param T Temperature [K]
  !! \param v Specific volume [m3/kmol]
  !! \param Z Mole fraction [-]
  !! \param dYdt Temperature derivative [J/kmolK]
  !! \param dYdv Pressure derivative [J/m3]
  !!       0: No derivatives
  !!       1: Analytical derivatives
  !!
  !! \retval Y Free energy with derivatives [-]
  !!
  !! \author GL, 2015-01-22
  subroutine cbCalcFreeEnergy(nc,cbeos,T,v,Z,Y,dYdt,dYdv,recalculate)
    use thermopack_constants, only: kRgas
    use cubic_eos, only: cb_eos
    use cbmix, only: cbCalcMixtureParams
    use cbHelm, only: cbFt
    implicit none
    integer, intent(in) :: nc
    class(cb_eos), intent(inout) :: cbeos
    real, intent (in)  :: T, v, Z(1:nc)
    real, intent (out) :: Y
    real, optional, intent(out) :: dYdt, dYdv
    logical, optional, intent(in) :: recalculate
    ! Local variables
    real :: f, dfdt, dfdv
    logical :: recalculate_loc

    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    end if

    if (recalculate_loc) then
      call cbCalcMixtureParams(nc, cbeos, T, Z)
      call cbCalcDerivatives_svol(nc, cbeos, T, v)
    end if

    f = cbeos%ff
    dfdt = cbFt(cbeos)
    dfdv = cbeos%ffv

    Y = kRgas*T*f

    if (present(dYdt)) then
      dYdt = kRgas*(f + T*dfdt)
    endif

    if (present(dYdv)) then
      dYdv = kRgas*T*dfdv
    endif

  end subroutine cbCalcFreeEnergy

  !> Calculate resudial reduced Helmholtz and differentials
  !!
  !! \param T - Temperature [K]
  !! \param v - Specific volume [m3/kmol]
  !! \param n - Mole numbers [mol]
  !! \param F... - Resudial reduced Helmholtz differentiuals
  !!
  !! \author Morten Hammer
  !! \date 2015-09
  subroutine cbCalcFres(nc,cbeos,T,V,n,F,F_T,F_V,F_n,F_TT,&
       F_TV,F_VV,F_Tn,F_Vn,F_nn,recalculate)
    use cubic_eos, only: cb_eos
    use cbHelm, only: cbF, cbFv, cbFvv, cbFt, cbFtt, cbFvt, &
         cbFi, cbFij, cbFiv, cbFit
    use cbmix, only: cbCalcMixtureParams
    implicit none
    integer, intent(in) :: nc
    real, intent(in) :: T,V,n(nc)
    ! Output.
    class(cb_eos), intent(inout) :: cbeos
    real, optional, intent(out) :: F,F_T,F_V,F_n(nc)
    real, optional, intent(out) :: F_TT,F_TV,F_Tn(nc),F_VV,F_Vn(nc),F_nn(nc,nc)
    logical, optional, intent(in) :: recalculate
    ! Locals
    logical :: recalculate_loc
    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    endif
    if (recalculate_loc) then
      call cbCalcMixtureParams(nc,cbeos,t,n)
      call cbCalcDerivatives_svol(nc, cbeos, T, v)
    endif

    if (present(F)) F = cbF(cbeos)
    if (present(F_T)) F_T = cbFt(cbeos)
    if (present(F_V)) F_V = cbFv(cbeos)
    if (present(F_n)) call cbFi(nc,cbeos,F_n)
    if (present(F_TT)) F_TT = cbFtt(cbeos)
    if (present(F_TV)) F_TV = cbFvt(cbeos)
    if (present(F_Tn)) call cbFit(nc,cbeos,F_Tn)
    if (present(F_VV)) F_VV = cbFvv(cbeos)
    if (present(F_Vn)) call cbFiv(nc,cbeos,F_Vn)
    if (present(F_nn)) call cbFij(nc,cbeos,F_nn)
  end subroutine cbCalcFres

    !> Calculates the contibution to the reduced residual Helmholtz energy F
  !> coming from the cubic eos, along with its derivatives.
  !> Note!!! Input and output in SI units
  subroutine calcCbFder_res_SI(nc,cbeos,T,V,n,F,F_T,F_V,F_n,F_TT,&
       F_TV,F_VV,F_Tn,F_Vn,F_nn,F_VVV,recalculate)
    use cbhelm, only: cbF, cbFv, cbFvv, cbFt, cbFtt, cbFvt, cbFvvv, &
         cbFi, cbFij, cbFiv, cbFit
    use cbmix, only: cbCalcMixtureParams
    use cubic_eos, only: cb_eos
    ! Input
    integer, intent(in) :: nc
    type(cb_eos), intent(inout) :: cbeos
    real, intent(in) :: T,V,n(nc)
    ! Output.
    real, optional, intent(out) :: F,F_T,F_V,F_n(nc),F_VVV
    real, optional, intent(out) :: F_TT,F_TV,F_Tn(nc),F_VV,F_Vn(nc),F_nn(nc,nc)
    logical, optional, intent(in) :: recalculate
    ! Locals
    logical :: recalculate_loc
    real :: sumn
    sumn = sum(n)

    if (present(recalculate)) then
      recalculate_loc = recalculate
    else
      recalculate_loc = .true.
    endif
    if (recalculate_loc) then
      ! Calculate contributions from the non-association part.
      call cbCalcMixtureParams(nc,cbeos,T,n/sumn)
      call cbCalcDerivatives_svol(nc,cbeos,T,1000*V/sumn)
    endif

    if (present(F)) F = sumn*cbF(cbeos)
    if (present(F_T)) F_T = sumn*cbFt(cbeos)
    if (present(F_V)) F_V = 1000*cbFv(cbeos)
    if (present(F_n)) call cbFi(nc,cbeos,F_n)
    if (present(F_TT)) F_TT = cbFtt(cbeos)*sumn
    if (present(F_TV)) F_TV = 1000*cbFvt(cbeos)
    if (present(F_Tn)) call cbFiT(nc,cbeos,F_Tn)
    if (present(F_VV)) F_VV = 1000*1000*cbFvv(cbeos)/sumn
    if (present(F_VVV)) F_VVV = 1.0e9*cbFvvv(cbeos,T,1000*V/sumn)/sumn**2
    if (present(F_Vn)) then
      call cbFiv(nc,cbeos,F_Vn)
      F_Vn = 1000*F_Vn/sumn
    end if
    if (present(F_nn)) then
      call cbFij(nc,cbeos,F_nn)
      F_nn = F_nn/sumn
    end if
  end subroutine calcCbFder_res_SI

end module cubic

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


text:
!
!
!
   INTEGER FUNCTION LJ (J)
   IMPLICIT NONE
   INTEGER J, L, M, K
   REAL*4 Z
!
! Given the harmonic index j=l(l+1)/2+m+1, returns l (degree)
! -GS 10.09.07
!
 If(j<0) then
	  Write(88,*) "LJ: The J-index is out of bounds"
 	  Write(88,*) "The program will STOP ----------"
 	  Write(*, *) "LJ: The J-index is out of bounds"
 	  Write(*, *) "The program will STOP ----------"
          call stop_config
	  Stop
 endif
!
   m=0
   k=j
10 z=(-1.0+sqrt(8.*float(k)-7.))/2.
   lj=aint(z)
   if(z-aint(z)==0.) return
   k=k-1
   m=m+1
   goto 10
   return
   end function LJ
!
!
!

max_stars_repo_name: gassmoeller/selen

max_issues_repo_name: gassmoeller/selen

max_forks_repo_name: gassmoeller/selen


text:
      subroutine  rpn2_cons_update_zero(meqn,maux, idir, iface, q,
     &                                  auxvec_center,
     &                                  auxvec_edge,flux)

      implicit none

      integer meqn,maux,idir, iface
      double precision q(meqn), flux(meqn)
      double precision auxvec_center(maux), auxvec_edge(maux)
      integer m

c     #  f(q) = (n dot u)*q
      do m = 1,meqn
c         # No flux function available for equations in non-conservative form
          flux(m) = 0
      enddo

      end

ext: f

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


subroutine broyden(iter)
  use inputoutput
  use salmon_parallel, only: nproc_group_global,nproc_group_orbital
  use salmon_communication, only: comm_summation
  use salmon_math
  use scf_data
  implicit none
  integer,parameter :: iter_mb=0
  real(8),parameter :: romega0=0.01d0
  integer :: iter_s,iter_e
  integer :: iter,i,j,ix,iy,iz
  integer :: is,is_sta,is_end
  integer :: ibox
  real(8),allocatable :: fb(:,:,:,:)
  real(8),allocatable :: del_fb(:,:,:,:),del_x(:,:,:,:)
  real(8),allocatable :: romega_mb(:)
  real(8) :: s,rho_temp(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),mg_sta(3):mg_end(3))
  real(8),allocatable :: aa(:,:),beta(:,:)
  real(8) :: rbox,rbox2

  if(ilsda==0)then
    is_sta=1
    is_end=1
  else if(ilsda==1)then
    is_sta=1
    is_end=2
  end if

  if (iter <= iter_mb+1) then
    if(ilsda==0)then
      rho(:,:,:)=rho_in(:,:,:,num_rho_stock)+alpha_mb*(rho_out(:,:,:,num_rho_stock)-rho_in(:,:,:,num_rho_stock))
    else
      rho_s(:,:,:,1:2)=rho_s_in(:,:,:,1:2,num_rho_stock)  &
                         +alpha_mb*(rho_s_out(:,:,:,1:2,num_rho_stock)-rho_s_in(:,:,:,1:2,num_rho_stock))
    end if
  else
    iter_s=max(iter_MB+1,iter-nmemory_mb)
    iter_e=iter-1
    allocate(fb(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),mg_sta(3):mg_end(3),iter_s:iter))
    allocate(del_fb(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),mg_sta(3):mg_end(3),iter_s:iter_e))
    allocate(del_x(mg_sta(1):mg_end(1),mg_sta(2):mg_end(2),mg_sta(3):mg_end(3),iter_s:iter_e))
    allocate(romega_mb(iter_s:iter_e))
    allocate(beta(iter_s:iter_e,iter_s:iter_e))
    allocate(aa(1:iter_e-iter_s+1,1:iter_e-iter_s+1))
    do is=is_sta,is_end
      if(ilsda==0)then
        do i=iter_s,iter
          ibox=num_rho_stock-(iter-i)
          fb(:,:,:,i)=rho_out(:,:,:,ibox)-rho_in(:,:,:,ibox)
        end do
      else if(ilsda==1)then
        do i=iter_s,iter
          ibox=num_rho_stock-(iter-i)
          fb(:,:,:,i)=rho_s_out(:,:,:,is,ibox)-rho_s_in(:,:,:,is,ibox)
        end do
      end if
      romega_mb(iter_s:iter_e)=1.d0
  !     &        /sqrt(sum((rho_out_b(:,iter_mb+1:iter-1)-rho_in_b(:,1:iter_mb+1:iter-1))**2))
      if(ilsda==0)then
        do i=iter_s,iter_e
          ibox=num_rho_stock-(iter-i)
          del_x(:,:,:,i)=rho_in(:,:,:,ibox+1)-rho_in(:,:,:,ibox)
          del_fb(:,:,:,i)=fb(:,:,:,i+1)-fb(:,:,:,i)
        end do
      else if(ilsda==1)then
        do i=iter_s,iter_e
          ibox=num_rho_stock-(iter-i)
          del_x(:,:,:,i)=rho_s_in(:,:,:,is,ibox+1)-rho_s_in(:,:,:,is,ibox)
          del_fb(:,:,:,i)=fb(:,:,:,i+1)-fb(:,:,:,i)
        end do
      end if
      do i=iter_s,iter_e
        rbox=0.d0
        do iz=ng_sta(3),ng_end(3)
        do iy=ng_sta(2),ng_end(2)
        do ix=ng_sta(1),ng_end(1)
          rbox=rbox+del_fb(ix,iy,iz,i)**2
        end do
        end do
        end do

        call comm_summation(rbox,rbox2,nproc_group_global)

        s=sqrt(rbox2)

        del_x(:,:,:,i)=del_x(:,:,:,i)/s
        del_fb(:,:,:,i)=del_fb(:,:,:,i)/s
      end do
      do i=1,iter_e-iter_s+1
        do j=1,iter_e-iter_s+1
          rbox=0.d0
          do iz=ng_sta(3),ng_end(3)
          do iy=ng_sta(2),ng_end(2)
          do ix=ng_sta(1),ng_end(1)
            rbox=rbox+del_fb(ix,iy,iz,iter_s-1+i)*del_fb(ix,iy,iz,iter_s-1+j)
          end do
          end do
          end do
          call comm_summation(rbox,rbox2,nproc_group_global)
          aa(i,j)=romega_mb(iter_s-1+i)*romega_mb(iter_s-1+j)*rbox2
          if (i == j) then
            aa(i,j)=aa(i,j)+romega0**2
          end if
        end do
      end do
      call matrix_inverse(aa)
      beta(iter_s:iter_e,iter_s:iter_e)=aa(1:iter_e-iter_s+1,1:iter_e-iter_s+1)
      rho_temp(:,:,:)=0.d0
      do i=iter_s,iter_e
        do j=iter_s,iter_e
          rbox=0.d0
          do iz=ng_sta(3),ng_end(3)
          do iy=ng_sta(2),ng_end(2)
          do ix=ng_sta(1),ng_end(1)
            rbox=rbox+del_fb(ix,iy,iz,i)*fb(ix,iy,iz,iter)
          end do
          end do
          end do
          call comm_summation(rbox,rbox2,nproc_group_orbital)
          rho_temp(:,:,:)=rho_temp(:,:,:)&
              &+romega_mb(i)*romega_mb(j)*beta(i,j)*rbox2*(alpha_mb*del_fb(:,:,:,j)+del_x(:,:,:,j))
        end do
      end do

      if(ilsda==0)then
        call add_in_broyden(rho(:,:,:),rho_in(:,:,:,num_rho_stock),fb(:,:,:,iter),rho_temp(:,:,:))
      else if(ilsda==1)then
        call add_in_broyden(rho_s(:,:,:,is),rho_s_in(:,:,:,is,num_rho_stock),fb(:,:,:,iter),rho_temp(:,:,:))
      end if

    end do

    deallocate(fb,del_fb,del_x,romega_mb,beta,aa)
  end if

  if(ilsda==1)then
    rho(:,:,:)=rho_s(:,:,:,1)+rho_s(:,:,:,2)
  end if

  return
end subroutine broyden


text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X QUIP Module
!X
!% This is a container module. ``use" it if you want to use QUIP stuff
!%>       use QUIP_module
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

module QUIP_internals_module
  use functions_module
  use quip_common_module
  use ewald_module
  use IPModel_lj_module
  use IPModel_sw_module
  use IPModel_tersoff_module
#ifdef HAVE_TB
  use tb_mixing_module
  use tb_common_module
  use tbmodel_bowler_module
  use tbmodel_dftb_module
  use tbmodel_nrl_tb_module
#endif
end module QUIP_internals_module

!% This is a container module. ``use" it if you want to use QUIP stuff
!%>       use QUIP_module
module QUIP_module
  use mpi_context_module
  use scalapack_module
  use matrix_module
  use rs_sparsematrix_module
#ifdef HAVE_TB
  use approxfermi_module
  use tb_kpoints_module
  use tbmodel_module
  use tbmatrix_module
  use tbsystem_module
  use tb_greensfunctions_module
  use tb_module
#endif
  use ip_module
  use potential_module
end module QUIP_module

max_stars_repo_name: albapa/QUIP

max_issues_repo_name: albapa/QUIP

max_forks_repo_name: albapa/QUIP


text:
module mod_monolis_c2d3_shape
  use mod_monolis_prm
  use mod_monolis_shape_util
  implicit none

  private

  real(kdouble), parameter :: gsp(2,3) = reshape([ &
    0.166666666666667d0, 0.166666666666667d0, &
    0.666666666666667d0, 0.166666666666667d0, &
    0.166666666666667d0, 0.666666666666667d0  &
    ], [2,3])

  real(kdouble), parameter :: np(2,3) = reshape([ &
     0.0d0, 0.0d0, &
     1.0d0, 0.0d0, &
     0.0d0, 1.0d0  &
    ], [2,3])

  real(kdouble), parameter :: weight(3) = [ &
     0.166666666666666d0,0.166666666666666d0,0.166666666666666d0 &
    ]

    public :: monolis_C2D3_num_gauss_point
    public :: monolis_C2D3_weight
    public :: monolis_C2D3_integral_point
    public :: monolis_C2D3_node_point
    public :: monolis_C2D3_shapefunc
    public :: monolis_C2D3_shapefunc_deriv
    public :: monolis_C2D3_get_global_deriv

contains

  function monolis_C2D3_num_gauss_point()
    implicit none
    integer(kint) :: monolis_C2D3_num_gauss_point
    monolis_C2D3_num_gauss_point = 3
  end function monolis_C2D3_num_gauss_point

  function monolis_C2D3_weight(i)
    implicit none
    integer(kint) :: i
    real(kdouble) :: monolis_C2D3_weight
    monolis_C2D3_weight = weight(i)
  end function monolis_C2D3_weight

  subroutine monolis_C2D3_integral_point(i, r)
    implicit none
    integer(kint) :: i
    real(kdouble) :: r(2)

    r(1) = gsp(1,i)
    r(2) = gsp(2,i)
  end subroutine monolis_C2D3_integral_point

  subroutine monolis_C2D3_node_point(i, r)
    implicit none
    integer(kint) :: i
    real(kdouble) :: r(2)

    r(1) = np(1,i)
    r(2) = np(2,i)
  end subroutine monolis_C2D3_node_point

  subroutine monolis_C2D3_shapefunc(local, func)
    implicit none
    real(kdouble) :: local(2), func(3)

    func(1) = 1.0d0 - local(1) - local(2)
    func(2) = local(1)
    func(3) = local(2)
  end subroutine monolis_C2D3_shapefunc

  subroutine monolis_C2D3_shapefunc_deriv(local, func)
    implicit none
    real(kdouble) :: local(2), func(3,2)

    func(1,1) = -1.0d0
    func(2,1) =  1.0d0
    func(3,1) =  0.0d0

    func(1,2) = -1.0d0
    func(2,2) =  0.0d0
    func(3,2) =  1.0d0
  end subroutine monolis_C2D3_shapefunc_deriv

  subroutine monolis_C2D3_get_global_deriv(node, r, dndx, det)
    implicit none
    real(kdouble) :: node(2,3), r(2), dndx(3,2), deriv(3,2)
    real(kdouble) :: xj(2,2), inv(2,2), det

    call monolis_C2D3_shapefunc_deriv(r, deriv)
    xj = matmul(node, deriv)
    call monolis_get_inverse_matrix_2d(xj, inv, det)
    dndx = matmul(deriv, inv)
  end subroutine monolis_C2D3_get_global_deriv

end module mod_monolis_c2d3_shape

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
ASUCD Resolution 26 states ASUCDs support for the Coffee House and the UC Davis campus community to become more socially and environmentally responsible purchasers of goods. See also Fair Trade and Fair Trade Coffee.

 Resolution 26
WHEREAS, many students are looking for more socially conscious and environmentally responsible ways to conduct their lives, ways which do not degrade people, animals or the environment; and

WHEREAS, many coffee farmers receive market payments for their coffee which are less than the costs of production, forcing them into a cycle of poverty and debt; and

WHEREAS, intensive coffee farming also leads to environmental problems, such as pesticide pollution, deforestation and extinction of some bird species through habitat destruction; and

WHEREAS, Fair Trade works to correct these imbalances by guaranteeing a minimum wage for the small producers harvest, as well as encouraging the cultivation of coffee without the use of pesticides; and

WHEREAS, with the profit generated from receiving a fair wage, coffee growers are able to invest in such areas as health, education, and environmental protection; and

WHEREAS, more than 500,000 farmers in 20 countries produce and sell more than 32 million pounds of coffee each year through the fair trade network; and

WHEREAS, the United States consumes a fifth of all coffee  produced in the world, and as a whole consumes the greatest  amount of goods and services in the world, invests the most  capital, and thus wields tremendous influence in the global  economy; and

WHEREAS, cooperatively owned companies are more socially responsible; and

WHEREAS, a current example of a cooperatively owned company,  Equal Exchange, reinvests 80% of its profits in its fair trade mission of supporting small farmers and sustainable agriculture;  and

WHEREAS, Equal Exchange is a fullservice company which is  competitive with the current ASUCD supplier and which has experience with larger institutions; and

WHEREAS, ASUCD, as a purchaser of goods and services, has a  responsibility to ensure that funds are spent in a manner  consistent with decent ethical and labor standards, including  assuring that workers are paid a living wage, and purchasing  goods made in the most fair manner when possible; and

WHEREAS, a change towards sustainable coffee will bear positively on the University and its students, by setting an example in the community and furthering education about responsible consumer choices;

THEREFORE BE IT RESOLVED THAT, insofar as is judged practical  and feasible by the Coffee House management, ASUCD will purchase  coffee that is certified as Fair Trade, Organic organic, and shadegrown  by independent certifying agencies and that is from the most socially and environmentally responsible company available; and

THEREFORE BE IT RESOLVED THAT, ASUCD promote Fair Trade coffee  to the greatest extent at the Coffee House and continue to  investigate environmentally and socially responsible  alternatives for its purchases; and

THEREFORE BE IT FURTHER RESOLVED THAT, ASUCD encourage  university and community leaders to adopt a similar policy of  social and environmental responsibility; and

THEREFORE BE IT FINALLY RESOLVED THAT, copies of this  resolution be sent to Chancellor Larry Vanderhoef; Sharon Coulson and Alex Park of the Coffee House; Global Exchange;  The California Aggie; The Davis Enterprise; The Sacramento  Bee; KDVS; and La Palabra, Because People Matter,  Java City, Starbucks, Café Roma, Delta of Venus,  Mishkas, City of Davis Chamber of Commerce, Davis City Council,  Davis City Manager, Retail Food Service Director of Sodexho Marriott Services, CURB Chair, and Oakland City Council Member  Nadel.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Test import-stmt
!      import-stmt  is  IMPORT [ [::] import-name-list ]
!
! Not tested here: import-name-list
!
! import-stmt is only allowed in an interface-body (C1210).
MODULE import_stmt
  INTEGER :: a, b, c

INTERFACE
   SUBROUTINE sub()
     ! No optional parts.
!TODO-F08     IMPORT
   END SUBROUTINE sub

   SUBROUTINE sub1()
     ! Optional import-name-list
!TODO-F08     IMPORT a, b
   END SUBROUTINE sub1

   SUBROUTINE sub2()
     ! Optional import-name-list with :: separater.
!TODO-F08     IMPORT :: a, b, c
   END SUBROUTINE sub2
END INTERFACE

END MODULE

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
!!**.****|****.****|****.****|****.****|****.****|****.****|****.****|****.****|
!! * bremss_spec.f *                               galprop package * 4/14/2000
!!**"****!****"****!****"****!****"****!****"****!****"****!****"****!****"****|

      subroutine bremss_spec(Egam,E0,IZ1,Ne1,dSdK)
c***********************************************************************
c          *** I.Moskalenko (MPE,Garching) *** version of 14.05.1998 ***
c PURPOSE:
c Calculation of the spectrum of e-bremsstrahlung in the fully ionized,
c hydrogen-like, or helium-like gas as well as in neutral H, He, and
c heavier in Born approximation; valid for E0,Egam >~ 0.01 MeV and
c beta0, beta1 >> 2*pi*Af*Z, corresponds to gamma=1.001
c REFERENCES:
c [1] Blumenthal & Gould 1970,Rev.Mod.Phys.42,237 (BG70);
c [2] Gould 1969, Phys.Rev.185,72 (G69);
c [3] Koch & Motz 1959, Rev.Mod.Phys.31,920 (KM59).
c [4] Strong A.W., Moskalenko I.V., Reimer O. 2000, ApJ 537, 763
c INPUT/OUTPUT parameters:
c Egam, E0 [MeV] are the energy of a gamma and total electron energy,
c respectively;
c IZ1 is the nucleus charge;
c Ne1 = 0 for unshielded charge, = 1 for 1-e atoms, = 2 for 2-e atoms.
c dSdK [barn/MeV] is the differential cross section (output).
c Note:
c #1 a one-parameter Hylleraas function is used for helium-like ions,
c    while Hartree-Fock functions are used for heutral He atoms;
c #2 a contribution of free electrons could be estimated by calling
c    the subroutine with Z=1, Ne=0 and multiplying dSdk by ionization level.
c***********************************************************************
      implicit real*8 (a-h,m,o-z)
      real*8 DD1(11)
      common/delta/ delta, IZ, Ne
     2      /Hartree/ DD(11),PH1(11),PH2(11)  ! DD = delta/(2.d0*Af)
      external FPHI1, FPHI2
!Gulli20070821
!$omp threadprivate(/delta/)

      dSdK = 0.d0
      mc2 = 0.511   ! MeV, electron rest mass
      E1 = E0-Egam  ! total energy of an electron after scattering
      if(E0 .le. 1.02*mc2 .or. E1 .le. mc2) return

      IZ = IZ1
      Ne = Ne1
      Pi = 3.14159265359d0
      Af = 7.2974d-3      ! fine structure constant
      R02 = 2.81794d-1**2 ! classical electron radius square (barn)
      EB1 = 0.07 ! MeV,kinetic energy; the upper bound for NonRelat.approx.
      EB2 = 2.  ! MeV, -"-; the lower bound for HE approx. (nonscreened)
      delta = Egam*mc2/(2.d0*E0*E1)
      Ef = 0.
      FS = 0.d0
      FS90 = 0.90
      T0 = E0-mc2   ! Ekin
      P0 = mc2*dsqrt((E0/mc2)**2-1.d0) ! initial electron momentum
      P1 = mc2*dsqrt((E1/mc2)**2-1.d0) ! final -"-
      beta0 = P0/E0
      beta1 = P1/E1

c### Fano-Sauter high-frequency limit (Egam ->E0-mc2, KM59,p.933) ###
      if(Egam/T0 .ge. FS90)
     2   FS =4.*Pi*IZ**3*R02*Af**2/T0*beta0*E0/mc2/(E0/mc2-1.d0)**2
     3      *(4./3. +E0/mc2*(E0/mc2-2.d0)/(E0/mc2+1.d0) *(1.d0
     4      -dlog(1.d0+2.d0*beta0/(1.d0-beta0))*(mc2/E0)**2/2.d0/beta0))

c## Elwert factor (Ef) is valid if   Z*Af*(1./beta1-1./beta0)<<1 ###
      if(T0 .lt. EB2 .and. IZ*Af*(1./beta1-1./beta0) .le. 0.5)
     2   Ef = beta0/beta1 *(1.d0-dexp(-2.d0*Pi*Af*IZ/beta0))
     3      /(1.d0-dexp(-2.d0*Pi*Af*IZ/beta1))

c### NONRELATIVISTIC ENERGIES ( T0 < 0.07 MeV) ###
c# Born approx.(KM59,p.925,3BNa) with Elwert factor (p.931)
c# Born approximation: 2*Pi*Z*Af/beta{0,1}<<1;
c# nonscreened: 2*delta/[Af*Z^(1/3)] >> 1;
      if(T0 .lt. EB1) then
         A = IZ*IZ                    ! no electron contribution
         dSdK = Ef*R02*Af*A*16./3. /Egam
     2      *(mc2/P0)**2*dlog(1.d0+2.d0*P1/(P0-P1))
         if(dSdK .lt. FS) dSdK = (FS+dSdK)/2.   !# 1/2 Fano-Sauter limit
         RETURN
      endif

c### LOW ENEGRY (T0 < 2. MeV) OR delta/(2.*Af*Z) > 4. ###
c# Born approximation (KM59,p.925,3BN): 2*Pi*Z*Af/beta{0,1}<<1;
c# nonscreened: 2*delta/Af >> 1
      if(T0 .lt. EB2 .or. delta/(2.*Af*IZ) .ge. 4.) then
         eps0 = dlog(1.d0+2.d0*P0/(E0-P0))  !=dlog((E0+P0)/(E0-P0))
         eps1 = dlog(1.d0+2.d0*P1/(E1-P1))  !=dlog((E1+P1)/(E1-P1))
         AL = 2.d0*dlog((E0*E1+P0*P1-mc2**2)/(Egam*mc2))
         T1 = 2.d0*E0*E1*((P0/P1)**2+1.d0)/P0**2
         T2 = mc2**2*(eps0*E1/P0**3+eps1*E0/P1**3-eps0*eps1/P0/P1)
         T10 = 8.d0/3.*(E0*E1)/(P0*P1)
         T11 = Egam**2*((E0/P0*E1/P1)**2+1.d0)/(P0*P1)
         T12 = mc2**2 *Egam/(2.d0*P0*P1) *(
     2      eps0*(E0/P0*E1/P0+1.d0)/P0 -eps1*(E0/P1*E1/P1+1.d0)/P1
     3      +2.d0*Egam*E0/P0*E1/P1/(P0*P1) )
         A = IZ*IZ +Ne            ! correction for atomic electrons
         if(T0 .lt. EB2)
     2      A =( IZ*IZ +Ne*(1.d0-dexp(-(T0-EB1)/9./EB1)) ) ! screening
     3         *(1.-0.3*dexp(-Egam/0.033)) *Ef  ! left & right connection
         dSdK = R02*Af*A/Egam*P1/P0*(4.d0/3.-T1+T2+AL*(T10+T11+T12))
         if(dSdK .lt. FS) dSdK = (FS+dSdK)/2.   !# 1/2 Fano-Sauter limit
         RETURN
      endif

c### HIGH  ENERGY (T0 > 2. MeV) AND delta/(2.*Af*Z) < 4. ###

c# Schiff approx. for neutral atoms heavier then He (KM59,p.925,3BSe)
      if(IZ .gt. 2 .and. IZ .eq. Ne) then
         b = IZ**(1.d0/3.d0)/111.d0/delta
         OM = 1.d0/delta**2/(1.d0+b*b)
         dSdK = 2.d0*IZ*IZ*R02*Af/Egam *(
     2      (1.d0+(E1/E0)**2-2.d0/3.*E1/E0)
     3         *(dlog(OM)+1.d0-2.d0/b*datan(b))
     4      +E1/E0*(2.d0/b/b*dlog(1.d0+b*b)
     5         +4./3.*(2.d0-b*b)/b**3*datan(b)-8.d0/3./b/b+2.d0/9.) )
         RETURN
      endif

c# arbitrary screening (G69, BG70; KM59,p.925,3BSb)
      A = IZ*IZ +Ne            !correction for atomic electrons
      phiu = -dlog(delta)-1.d0/2.d0
c# Hartree-Fock approximation for neutral He atoms
      if(IZ .eq. 2 .and. IZ .eq. Ne) then
         DD1(1) = dlog(1.d-3)
         do i=2,11
            DD1(i) = dlog(DD(i))
         enddo
         if(delta/(2.d0*Af) .ge. DD(2)) then
            if(delta/(2.d0*Af) .le. DD(11)) then
               call DINTER(DD1,PH1,11,dlog(delta/(2.d0*Af)),phi1)
               call DINTER(DD1,PH2,11,dlog(delta/(2.d0*Af)),phi2)
            else
               phi1 = 4.*A *phiu  ! asymptotics
               phi2 = 4.*A *phiu
            endif
         else
            phi1 = PH1(1)
            phi2 = PH2(1)
         endif
         phi1 = phi1/4.
         phi2 = phi2/4.
      else
         if(Ne .eq. 0) then  !# UNSHIELDED charge
            phi1 = A*phiu
            phi2 = A*phiu
         else
c# H-like atoms & Hylleraas-1 approximation for He-like atoms
            if(delta/(2.d0*Af) .le. 1.d-4) delta = 1.d-4*(2.d0*Af)
            call SIM1(1.d0,delta,1.d-3,1.d-5,1.d-4,FPHI1,p1)
            call SIM1(1.d0,delta,1.d-3,1.d-5,1.d-4,FPHI2,p2)
            phi1 = 2.*IZ*(-p1       +1.+(Ne-1.)/IZ)  +(IZ-Ne)**2*phiu
            phi2 = 2.*IZ*(-p2+5./6.*(1.+(Ne-1.)/IZ) )+(IZ-Ne)**2*phiu
         endif
      endif
      if(phi1 .gt. A*phiu) phi1 = A *phiu  ! asymptotics
      if(phi2 .gt. A*phiu) phi2 = A *phiu

      dSdK =4.*R02*Af/Egam *((1.d0+(E1/E0)**2)*phi1 -2./3.*E1/E0*phi2)
      if(dSdK .lt. FS) dSdK = (FS+dSdK)/2.      !# 1/2 Fano-Sauter limit
      return
      end

      real*8 function FPHI1(Q)
c***********************************************************************
c used for calculation of PHI1 function (G69,p.74,76)
c          *** I.Moskalenko (MPE,Garching) *** version of 23.09.1997 ***
c***********************************************************************
      implicit real*8 (a-h,o-z)
      common/delta/ delta, IZ, Ne
!Gulli20070821
!$omp threadprivate(/delta/)

      Af = 7.2974d-3          ! fine structure constant
      if(Ne .eq. 1) then      ! one-electron atoms
         AZ = 1.d0/(2.d0*Af*IZ)
         FZ = 1.d0/(1.d0+(AZ*Q)**2)**2 ! form-factor for one-electron atom
         FPHI1 = (Q-delta)**2/Q**3 *(1.d0-FZ)
         return
      endif                   ! two-electron atoms (Hylleraas-1 function)
      AZ = 1.d0/(2.d0*Af*(IZ-5.d0/16.d0))
      FZ = 1.d0/(1.d0+(AZ*Q)**2)**2 ! form-factor for two-electron atom
      FPHI1 = (Q-delta)**2/Q**3 *(2.*(1.d0-FZ) -(1.d0-FZ**2)/IZ)
      return
      end

      real*8 function FPHI2(Q)
c***********************************************************************
c used for calculation of PHI2 function (G69,p.74,76)
c          *** I.Moskalenko (MPE,Garching) *** version of 23.09.1997 ***
c***********************************************************************
      implicit real*8 (a-h,o-z)
      common/delta/ delta, IZ, Ne
!Gulli20070821
!$omp threadprivate(/delta/)

      Af = 7.2974d-3          ! fine structure constant
      if(Ne .eq. 1) then      ! one-electron atoms
         AZ = 1.d0/(2.d0*Af*IZ)
         FZ = 1.d0/(1.d0+(AZ*Q)**2)**2 ! atomic form-factor
         FPHI2=(Q**3+(3.d0-6.d0*dlog(Q/delta))*delta**2*Q
     2      -4.d0*delta**3)/Q**4 *(1.d0-FZ)
         return
      endif                   ! two-electron atoms (Hylleraas-1 function)
      AZ = 1.d0/(2.d0*Af*(IZ-5.d0/16.d0))
      FZ = 1.d0/(1.d0+(AZ*Q)**2)**2 ! form-factor for two-electron atom
      FPHI2=(Q**3+(3.d0-6.d0*dlog(Q/delta))*delta**2*Q
     2   -4.d0*delta**3)/Q**4 *(2.*(1.d0-FZ) -(1.d0-FZ**2)/IZ)
      return
      end


text:
! @expect error

program main
  use smack
  implicit none
  integer, pointer :: box
  integer, target :: x = 3
  box => x
  print *, box == 3
  call assert(box /= 3)
end program main

max_stars_repo_name: soarlab/gandalv

max_issues_repo_name: soarlab/gandalv

max_forks_repo_name: soarlab/gandalv


text:
module type_timer
#include <messenger.h>
    use mod_kinds,      only: rk, ik
    use mod_constants,  only: ZERO
    use mod_chidg_mpi,  only: GLOBAL_MASTER
    use mpi_f08,        only: MPI_WTime
    implicit none

    !>  Timer type for timing procedures
    !!
    !!  @author Nathan A. Wukie
    !!
    !!
    !!
    !-----------------------------------------------------------
    type, public    :: timer_t

        real(rk)    :: start_time           !< Timer starting time
        real(rk)    :: stop_time            !< Timer ending time
        real(rk)    :: elapsed_time = ZERO  !< Timer elapsed time

        logical     :: started = .false.    !< Logical state indicating if the timer was started
        logical     :: stopped = .false.    !< Logical state indicating if the timer was stopped

    contains
        procedure   :: start                !< Start timer
        procedure   :: stop                 !< Stop timer
        procedure   :: elapsed              !< Compute elapsed time
        procedure   :: report               !< Print elapsed time
        procedure   :: reset                !< Reset timer

    end type timer_t
    !-----------------------------------------------------------

contains

    !> Start the timer
    !!
    !!  @author Nathan A. Wukie
    !!
    !---------------------------------------------------------------------------
    subroutine start(self)
        class(timer_t), intent(inout)  :: self

        if (.not. self%started) then
            !call cpu_time(self%start_time)
            self%start_time = MPI_WTime()
        else
            call chidg_signal(WARN,'type_timer: Timer was already started')
        end if

        self%started = .true.

    end subroutine start
    !***************************************************************************

    !> Stop the timer
    !!
    !!  @author Nathan A. Wukie
    !!
    !---------------------------------------------------------------------------
    subroutine stop(self)
        class(timer_t), intent(inout)  :: self

        if (self%started) then
            !call cpu_time(self%stop_time)
            self%stop_time = MPI_Wtime()
        else
            call chidg_signal(WARN,'type_timer: Timer was never started')
        end if

        self%elapsed_time = self%elapsed_time + real(self%stop_time - self%start_time,rk)

        self%stopped = .true.
        self%started = .false.

    end subroutine stop
    !***************************************************************************

    !> Return the elapsed time
    !!
    !!  @author Nathan A. Wukie
    !!
    !----------------------------------------------------------------------------
    function elapsed(self) result(elapsed_time)
        class(timer_t), intent(inout)  :: self

        real(rk)    :: elapsed_time

        if (self%stopped) then
            !elapsed_time = real(self%stop_time - self%start_time,rk)
            elapsed_time = self%elapsed_time
        else
            call chidg_signal(WARN,'type_timer: Timer was not stopped')
            elapsed_time = 123456789._rk
        end if

    end function elapsed
    !****************************************************************************

    !> Procedure for reporting timer information
    !!
    !!  @author Nathan A. Wukie
    !!
    !!  @param[in]  string  Optional character string for reporting contextual information about the timing
    !!
    !------------------------------------------------------------------------------------------------------
    subroutine report(self,messg)
        class(timer_t), intent(inout)           :: self
        character(*),   intent(in),  optional   :: messg

        real(rk)    :: elapsed

        ! Compute the elapsed time
        elapsed =  self%elapsed()

        ! If a message was included, print and append timing. Else, just print the timing
        if (present(messg)) then
            call write_line(messg, ' ', elapsed, delimiter='', io_proc=GLOBAL_MASTER)
        else
            call write_line(elapsed, io_proc=GLOBAL_MASTER)
        end if

    end subroutine report
    !*******************************************************************************************************

    !> Procedure for reseting the clock
    !!
    !!  @author Nathan A. Wukie
    !!
    !!
    !!
    !-------------------------------------------------------------------
    subroutine reset(self)
        class(timer_t), intent(inout)       :: self

        self%start_time   = 0.
        self%stop_time    = 0.
        self%elapsed_time = 0.

        self%started      = .false.
        self%stopped      = .false.

    end subroutine reset
    !********************************************************************************************************

end module type_timer

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
MODULE finish

  USE partlist
  USE laser

  IMPLICIT NONE

  PRIVATE

  PUBLIC :: finalise

CONTAINS

  SUBROUTINE finalise

    CALL deallocate_memory
    CALL MPI_FINALIZE(errcode)

  END SUBROUTINE finalise

  SUBROUTINE deallocate_memory

    INTEGER :: i, stat
    TYPE(particle_species), POINTER :: species, next_species
    TYPE(particle_species), POINTER :: iolist, next_iolist

    DEALLOCATE(x, x_global, xb_global)
    DEALLOCATE(y, y_global, yb_global)
    DEALLOCATE(z, z_global, zb_global)
    DEALLOCATE(ex, ey, ez, bx, by, bz, jx, jy, jz)

    DEALLOCATE(npart_each_rank)
    DEALLOCATE(x_grid_mins, x_grid_maxs, cell_x_min, cell_x_max)
    DEALLOCATE(y_grid_mins, y_grid_maxs, cell_y_min, cell_y_max)
    DEALLOCATE(z_grid_mins, z_grid_maxs, cell_z_min, cell_z_max)

    DEALLOCATE(ex_x_min, ex_x_max, ey_x_min, ey_x_max, ez_x_min, ez_x_max)
    DEALLOCATE(bx_x_min, bx_x_max, by_x_min, by_x_max, bz_x_min, bz_x_max)
    DEALLOCATE(ex_y_min, ex_y_max, ey_y_min, ey_y_max, ez_y_min, ez_y_max)
    DEALLOCATE(bx_y_min, bx_y_max, by_y_min, by_y_max, bz_y_min, bz_y_max)
    DEALLOCATE(ex_z_min, ex_z_max, ey_z_min, ey_z_max, ez_z_min, ez_z_max)
    DEALLOCATE(bx_z_min, bx_z_max, by_z_min, by_z_max, bz_z_min, bz_z_max)

    next_species => species_list
    next_iolist => io_list_data
    DO i = 1, n_species
      species => next_species
      next_species => species%next

      CALL destroy_partlist(species%attached_list)
      DEALLOCATE(species%ext_temp_x_min, STAT=stat)
      DEALLOCATE(species%ext_temp_x_max, STAT=stat)
      DEALLOCATE(species%ext_temp_y_min, STAT=stat)
      DEALLOCATE(species%ext_temp_y_max, STAT=stat)
      DEALLOCATE(species%ext_temp_z_min, STAT=stat)
      DEALLOCATE(species%ext_temp_z_max, STAT=stat)
      DEALLOCATE(species%density, STAT=stat)
      DEALLOCATE(species%temp, STAT=stat)
      DEALLOCATE(species%drift, STAT=stat)
      DEALLOCATE(species, STAT=stat)

      iolist => next_iolist
      next_iolist => iolist%next
      DEALLOCATE(iolist, STAT=stat)
    ENDDO

    IF (ALLOCATED(field_probes)) DEALLOCATE(field_probes)

    CALL deallocate_lasers

    CALL MPI_COMM_FREE(comm, errcode)

  END SUBROUTINE deallocate_memory

END MODULE finish

max_stars_repo_name: TomGoffrey/minepoch

max_issues_repo_name: TomGoffrey/minepoch

max_forks_repo_name: TomGoffrey/minepoch


text:
!
! The Laboratory of Algorithms
!
! The MIT License
!
! Copyright 2011-2015 Andrey Pudov.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the 'Software'), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!

module MGraph

    use MArrays
    use MVertex

    implicit none
    private

    integer, parameter :: DEFAULT_NUMBER_OF_VERTICES = 16

    type, public :: TGraph
        private
        logical, dimension(:,:), pointer     :: adjacencyMatrix
        type(TVertex), dimension(:), pointer :: vertexList

        integer  :: vertexesNumber
    contains
        procedure :: addVertex
        procedure :: addEdge

        procedure :: getAdjacencyMatrix
        procedure :: getVertexList
        procedure :: getVertexNumber

        procedure :: init
        procedure :: destroy
    end type
contains
    subroutine addVertex(instance, value)
        class(TGraph), intent(in out) :: instance
        integer, intent(in)           :: value

        type(TArrays) arrays
        type(TVertex) vertex

        type(TVertex), dimension(:), pointer :: temporary_array1d
        logical, dimension(:,:), pointer     :: temporary_array2d
        integer length
        integer index

        call vertex%init(value, .false.)
        instance%vertexesNumber = instance%vertexesNumber + 1

        ! increase array size when required
        if (instance%vertexesNumber > size(instance%vertexList)) then
            ! increase vertex list
            length = size(instance%vertexList) * 3 / 2
            allocate(temporary_array1d(length))
            temporary_array1d(1:size(instance%vertexList)) = instance%vertexList
            deallocate(instance%vertexList)
            instance%vertexList => temporary_array1d

            ! increase adjacency matrix
            length = size(instance%vertexList)
            allocate(temporary_array2d(length, length))
            temporary_array2d(1:size(instance%adjacencyMatrix, 1), &
                    1:size(instance%adjacencyMatrix, 2)) = instance%adjacencyMatrix
            deallocate(instance%adjacencyMatrix)
            instance%adjacencyMatrix => temporary_array2d
        end if

        instance%vertexList(instance%vertexesNumber) = vertex
    end subroutine

    subroutine addEdge(instance, start, end)
        class(TGraph), intent(in out) :: instance
        integer, intent(in)           :: start
        integer, intent(in)           :: end

        instance%adjacencyMatrix(start, end) = .true.
        instance%adjacencyMatrix(end, start) = .true.
    end subroutine

    function getAdjacencyMatrix(instance) result(value)
        class(TGraph), intent(in)        :: instance
        logical, dimension(:,:), pointer :: value

        value => instance%adjacencyMatrix
    end function

    function getVertexList(instance) result(value)
        class(TGraph), intent(in)            :: instance
        type(TVertex), dimension(:), pointer :: value

        value => instance%vertexList
    end function

    function getVertexNumber(instance) result(value)
        class(TGraph), intent(in) :: instance
        integer :: value

        value = instance%vertexesNumber
    end function

    subroutine init(instance)
        class(TGraph), intent(in out) :: instance

        integer index
        integer jndex

        allocate(instance%adjacencyMatrix(DEFAULT_NUMBER_OF_VERTICES, DEFAULT_NUMBER_OF_VERTICES))
        allocate(instance%vertexList(DEFAULT_NUMBER_OF_VERTICES))

        instance%vertexesNumber  = 0

        do index = 1, size(instance%adjacencyMatrix, 1)
            do jndex = 1, size(instance%adjacencyMatrix, 2)
                instance%adjacencyMatrix(index, jndex) = 0
            end do
        end do
    end subroutine

    subroutine destroy(instance)
        class(TGraph), intent(in out) :: instance

        deallocate(instance%adjacencyMatrix)
        deallocate(instance%vertexList)

        instance%vertexesNumber = 0
    end subroutine
end module

max_stars_repo_name: andreypudov/Algorithms

max_issues_repo_name: jlokimlin/Algorithms

max_forks_repo_name: jlokimlin/Algorithms