text:
!> The module eosdata contains the definitions of the equation of state, mixing
!> rule and the interaction parameters.

module cubic_eos
  use eosdata
  use thermopack_constants, only: eosid_len, bibref_len, ref_len, uid_len, mix_len
  use thermopack_var, only: base_eos_param, base_eos_dealloc
  use unifac, only: unifacdb
  implicit none

  type :: kijdatadb
    character(len=eosid_len) :: eosid
    character(len=eosid_len) :: mruleid
    character(len=ref_len) :: ref
    character(len=bibref_len) :: bib_ref
    character(len=uid_len) :: uid1
    character(len=uid_len) :: uid2
    real :: kijvalue
  end type kijdatadb

  type :: interGEdatadb
    character(len=eosid_len) :: eosid
    character(len=eosid_len) :: mruleid
    character(len=ref_len) :: ref
    character(len=bibref_len) :: bib_ref
    character(len=uid_len) :: uid1
    character(len=uid_len) :: uid2
    real :: kijvalue
    integer :: correlation
    real :: alphaijvalue(2)
    real :: polyij(3)
    real :: polyji(3)
  end type interGEdatadb

  type :: mixExcessGibbs
    integer :: mGE ! Variant: MHV1 (cbMixHuronVidal),
    !          MHV2 (cbMixHuronVidal2),
                   !          NRTL (cbMixNRTL)
    integer, allocatable, dimension(:,:) :: correlation ! 0: Inactive (Use kij)
    ! 1: default, 2: Maribo-Mogensen (2015) - 10.1002/aic.14829
    real, allocatable, dimension(:,:) :: alpha
    real, allocatable, dimension(:,:) :: aGE, bGE, cGE
    ! tau*T = dg/R = aGE + bGE*T + cGE*T^2
    ! tau is dimensionless
    ! TODO: General NRTL form: tau = a + b/T + c/T^2 + d ln(T) + e T^f
  contains
    procedure, public :: dealloc => excess_gibbs_deallocate
    procedure, public :: excess_gibbs_allocate_and_init
    ! Assignment operator
    procedure  :: assign_excess_gibbs_mix
    generic,   public   :: assignment(=) => assign_excess_gibbs_mix
  end type mixExcessGibbs

  integer, parameter :: nDegreePoly=2
  type :: Fraction
    real, dimension(nDegreePoly+1):: pNum !<Numerator polynom
    real, dimension(nDegreePoly+1):: pDen !<Denominator polynom
  end type Fraction
  !< y = (pNum(1)+pNum(2)*x+...pNum(n-1)*x**n)/(pDen(1)+pDen(2)*x+...pDen(n-1)*x**n)
  !< where n = nDegreePoly.  Generally pDen(1) = 1.0

  type :: mixWongSandler
    real, allocatable, dimension(:,:) :: alphaij
    type(Fraction), allocatable, dimension(:,:) :: f_kij, f_tauij
  contains
    procedure, public :: dealloc => WS_deallocate
    procedure, public :: WS_allocate_and_init
    ! Assignment operator
    procedure  :: assign_WS_mix
    generic,   public   :: assignment(=) => assign_WS_mix
  end type mixWongSandler

  type :: singleData
     !> Cubic EoS quantities for individual components
     real :: omegaA, omegaB, omegaC, zcrit ! Allow for 3-paramter EOS
     real :: a !< Energy constant in front of alpha. Units: Pa*L^2/mol^2.
     real :: b !< Covolume parameter. Units: L/mol.
     real :: c !< Parameter in Patel-Teja EoS. Units: L/mol.

     real :: acf !< Acentric factor to use for the cubic EoS [-]
     real :: tc  !< Critical temperature to use in the cubic EoS [K]
     real :: pc  !< Critical pressure to use in the cubic EoS [Pa]

     ! Alpha correlations for a factor
     integer :: alphaMethod=-1
     character (len=short_label_len) :: alphaCorrName
     real :: alpha, dalphadt, d2alphadt2 !< alpha is dimensionless
     real :: alphaParams(3)=-1e10 ! Fitted coeffs in correlation. Max 3.

     ! Beta correlations for b factor
     integer :: betaMethod=1
     character (len=12) :: betaCorrName = "Classic     "
     real :: beta, dbetadt, d2betadt2 !< beta is dimensionless
     real :: betaParams(3)=-1e10 ! Fitted coeffs in correlation. Max 3.
  end type singleData

  type, extends(base_eos_param) :: cb_eos
    character (len=mix_len) :: mruleid
    character (len=20) :: name
    integer :: mruleidx
    logical  :: cubic_verbose = .false.

    ! Parameters m1 and m2 appearing in the denominator of the attractive term
    real :: m1, m2
    real :: dm1dB, dm1dC, dm2dB, dm2dC
    real :: d2m1dB2, d2m1dC2, d2m2dB2, d2m2dC2,d2m2dBdC,d2m1dBdC
    real :: delta ! m1 and m2 are usually functions of this EoS-dependent constant

    real :: a, b, c !< Dependent of component
    real :: sumN
    real :: sumA !< [Pa*L^2/mol^2]
    real :: sumB !< Molfraction average of single-component b_i [L/mol]
    real :: sumC !< [L/mol]

    ! Pressure differentials
    real :: pn,pa,pb,pc,pt,pv

    ! Reduced residual Helmholtz free energy ff, and its derivatives wrt
    ! explicit variables T,v,n and explicit parameters a,b,c. These are
    ! combined to yield total derivatives wrt T,v,n in the cbhelm module.
    real :: ff,fft,fftt
    real :: ffn,ffnv,ffna,ffnb,ffnc,ffnn,ffnt
    real :: ffa,ffaa,ffab,ffac,ffat
    real :: ffb,ffbb,ffbc,ffbt
    real :: ffc,ffcc,ffct
    real :: ffv,ffvt,ffvv,ffva,ffvb,ffvc

    ! Parameter differentials that will be combined with the ff differentials
    real :: at, att, bt, btt ! derivatives wrt temperature
    real, allocatable  :: ai(:),ait(:), aij(:,:) ! derivatives wrt mole numbers and temperature
    real, allocatable  :: bi(:),bit(:),bitt(:), ci(:), bij(:,:), cij(:,:) ! derivatives wrt mole numbers and temperature

    ! Parameters for numerical solution of cubic EoS
    integer :: extrm(2), nextrm(2), nzfac(2)

    ! Cubic EoS quantities for individual components
    type(singleData), allocatable :: single(:)

    ! Mixing parameters
    real, allocatable, dimension(:,:) :: kij ! interaction parameter for energy
    real, allocatable, dimension(:,:) :: lij ! interaction parameter for size
    real, allocatable, dimension(:,:) :: lowcase_bij
    logical :: simple_covolmixing
    type(mixExcessGibbs) :: mixGE
    type(mixWongSandler) :: mixWS
    type(unifacdb), pointer :: unifdb => NULL()

  contains
    procedure, public :: dealloc => cubic_eos_dealloc
    procedure, public  :: allocate_and_init => allocate_and_init_cubic_eos
    !procedure, public :: de_allocate_cubic_eos
    ! Assignment operator
    procedure, pass(This), public :: assign_eos => assign_cubic_eos
  end type cb_eos

  type, extends(cb_eos) :: lk_eos

  !   private

  ! contains
  end type lk_eos

  type, extends(cb_eos) :: cpa_eos

  !   private

  ! contains
  end type cpa_eos

  integer, parameter :: cbMixClassicGroup = 1 !< Classic kij type mixing
  integer, parameter :: cbMixVdW = 11 !< Classic vdW mixing rule for am and bm - using k_ij == k_ji
  integer, parameter :: cbMixVdWCPA = 12 !< CPA mixing rule (same as cbMixVdW, but kij from another db)
  integer, parameter :: cbMixReid = 13    !< Unsymmertic mixing rule where k_ij can be different than k_ji
  integer, parameter :: cbMixGEGroup = 2
  integer, parameter :: cbMixHuronVidal = 21 !< Huron vidal mixing rule
  integer, parameter :: cbMixHuronVidal2 = 22 !< Huron vidal mixing rule
  integer, parameter :: cbMixNRTL = 23 !< NRTL mixing rule
  integer, parameter :: cbMixUNIFAC = 24 !< UNIFAC mixing rule
  integer, parameter :: cbMixHVCPA = 25 !< Huron Vidal mixing rule (classic, but kij from another db)
  integer, parameter :: cbMixHVCPA2 = 26 !< Huron Vidal mixing rule (classic, but kij from another db)
  integer, parameter :: cbMixWongSandler = 3 !< Wong Sandler mixing rule

  type mix_label_mapping
    integer :: mix_idx_group
    integer :: mix_idx
    character(len=short_label_len) :: short_label
    character(len=label_len) :: label
    character(len=label_len) :: alias
  end type mix_label_mapping

  integer, parameter :: n_mix_rules = 10
  type(mix_label_mapping), dimension(n_mix_rules), parameter :: mix_label_db = (/&
       mix_label_mapping(mix_idx_group=cbMixClassicGroup,&
       mix_idx=cbMixVdW, short_label="VDW", label="Classic",&
       alias = "CLASSIC"), &
       mix_label_mapping(mix_idx_group=cbMixClassicGroup,&
       mix_idx=cbMixVdWCPA, short_label="Classic(CPA)", label="Classic(CPA)",&
       alias = ""), &
       mix_label_mapping(mix_idx_group=cbMixClassicGroup,&
       mix_idx=cbMixReid, short_label="Reid", label="Reid",&
       alias = ""), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixHuronVidal, short_label="HV", label="Huron-Vidal",&
       alias = "HV1/HV0"), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixHuronVidal2, short_label="HV2", label="Huron-Vidal2",&
       alias = ""), &
       mix_label_mapping(mix_idx_group=cbMixWongSandler,&
       mix_idx=cbMixWongSandler, short_label="WongSandler", label="Wong-Sandler",&
       alias = "WS"), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixNRTL, short_label="NRTL", label="NRTL",&
       alias = ""), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixUNIFAC, short_label="UNIFAC", label="UNIFAC",&
       alias = "UMR/VTPR"), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixHVCPA, short_label="HVCPA", label="HVCPA",&
       alias = ""), &
       mix_label_mapping(mix_idx_group=cbMixGEGroup,&
       mix_idx=cbMixHVCPA2, short_label="HVCPA2", label="HVCPA2",&
       alias = "") &
       /)

  integer, parameter :: nHVCorrs = 4
  integer, parameter, dimension(nHVCorrs) :: HVCorrIndices = (/&
       cbMixHuronVidal,&
       cbMixHuronVidal2,&
       cbMixHVCPA,&
       cbMixHVCPA2/)

  integer, parameter :: nGECorrs = 6
  integer, parameter, dimension(nGECorrs) :: GECorrIndices = (/&
       HVCorrIndices,&
       cbMixNRTL,&
       cbMixGEGroup/)

  ! Alphacorrelations. To add new correlations, add it below and update module cbAlpha
  integer, parameter :: cbAlphaClassicIdx = 1 !< Classic alpha-correlation VdW, RK, SRK and PR
  integer, parameter :: cbAlphaTwuIdx = 2     !< Twu-Coon-Bluck-Cunninghan exponential formulation
  integer, parameter :: cbAlphaMcIdx = 3     !< Mathias-Copeman expression for polar substances
  integer, parameter :: cbAlphaGergIdx = 4     !< Gerg-Water(PR) expression for polar substances Q1,Q2 and Q3
  integer, parameter :: cbAlphaClassicFitIdx = 5 !< Classic alpha-corr where with replaced by a fitted parameter. Used in CPA.
  integer, parameter :: cbAlphaUMRIdx = 6 !< Classic alpha-corr with another function m(acf)

  integer, parameter :: cbAlphaGBIdx = 7 !< Alpha-corr of Graboski and Daubert (10.1021/i260068a009)
  integer, parameter :: cbAlphaRKIdx = 8 !< Alpha-corr of Redlich-Kwong (10.1021/cr60137a013)
  integer, parameter :: cbAlphaSoaveIdx = 9 !< Alpha-corr of Soave ()
  integer, parameter :: cbAlphaPRIdx = 10 !< Alpha-corr of Peng-Robinson
  integer, parameter :: cbAlphaPTIdx = 11 !< Alpha-corr of Patel-Teja
  integer, parameter :: cbAlphaSWIdx = 12 !< Alpha-corr of Schmidt-Wensel
  integer, parameter :: cbAlphaVDWIdx = 13 !< Alpha-corr of van der Waals

  type alpha_label_mapping
    integer :: alpha_idx
    integer :: n_param
    character(len=short_label_len) :: short_label
    character(len=label_len) :: description
    integer :: classic_for_eos_idx
  end type alpha_label_mapping

  integer, parameter :: n_alpha_corrs = 13
  type(alpha_label_mapping), dimension(n_alpha_corrs), parameter :: alpha_corr_db = (/&
       alpha_label_mapping(alpha_idx=cbAlphaClassicIdx, n_param=1, short_label="CLASSIC", &
       description="Classic alpha-correlation VdW, SRK, PR, PT and SW", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaTwuIdx, n_param=3, short_label="TWU", &
       description="Twu-Coon-Bluck-Cunninghan exponential formulation", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaMcIdx, n_param=3, short_label="MC", &
       description="Mathias-Copeman expression for polar substances", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaGergIdx, n_param=3, short_label="GERG", &
       description="Gerg-Water(PR) expression for polar substances", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaClassicFitIdx, n_param=1, &
       short_label="CLASSICFIT", &
       description="Classic alpha-corr with fitted parameter. Used in CPA.", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaUMRIdx, n_param=1, short_label="UMR", &
       description="Classic alpha-corr with another function m(acf)", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaGBIdx, n_param=1, short_label="GD", &
       description="Alpha-corr of Graboski and Daubert (10.1021/i260068a009)", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaRKIdx, n_param=1, short_label="RK", &
       description="Alpha-corr of Redlich-Kwong (10.1021/cr60137a013)", &
       classic_for_eos_idx=-1), &
       alpha_label_mapping(alpha_idx=cbAlphaSoaveIdx, n_param=1, short_label="SOAVE", &
       description="Alpha-corr of Soave", &
       classic_for_eos_idx=cbSRK), &
       alpha_label_mapping(alpha_idx=cbAlphaPRIdx, n_param=1, short_label="PR", &
       description="Alpha-corr of Peng-Robinson", &
       classic_for_eos_idx=cbPR), &
       alpha_label_mapping(alpha_idx=cbAlphaPTIdx, n_param=1, short_label="PT", &
       description="Alpha-corr of Patel-Teja", &
       classic_for_eos_idx=cbPT), &
       alpha_label_mapping(alpha_idx=cbAlphaSWIdx, n_param=1, short_label="SW", &
       description="Alpha-corr of cbAlphaSWIdx", &
       classic_for_eos_idx=cbSW), &
       alpha_label_mapping(alpha_idx=cbAlphaVDWIdx, n_param=1, short_label="VDW", &
       description="Alpha-corr of cbAlphaVDWIdx", &
       classic_for_eos_idx=cbVDW) &
       /)

  ! Beta correlations for cubic EoS "b" parameter. To add new correlations, add
  ! it below and update module cbBeta.
  integer, parameter :: nBetaCorrs = 2
  integer, parameter :: cbBetaClassicIdx = 1 !< Classic beta-correlation yielding the classic, temperature-independent covolume
  integer, parameter :: cbBetaQuantumIdx = 2 !< Beta-correlation for quantum fluids He, Ne, H2, D2
  character(len=11), dimension(nBetaCorrs), parameter :: betaCorrNames = (/ &
       "CLASSIC    ", "QUANTUM    " /)
  integer, parameter, dimension(nBetaCorrs) :: betaCorrNumParams = (/0,2/)

contains

  subroutine WS_deallocate(mixWS)
    class(mixWongSandler), intent(inout) :: mixWS
    ! Locals
    integer :: err
    !
    err = 0
    if (allocated(mixWS%f_kij)) deallocate (mixWS%f_kij, stat=err)
    if (err /= 0) call stoperror('WS_deallocate: could not deallocate array: mixWS%f_kij')
    !
    if (allocated(mixWS%f_tauij)) deallocate (mixWS%f_tauij, stat=err)
    if (err /= 0) call stoperror('WS_deallocate: could not deallocate array: mixWS%f_tauij')
    if (allocated(mixWS%alphaij)) deallocate (mixWS%alphaij, stat=err)
    if (err /= 0) call stoperror('WS_deallocate: could not allocate array: mixWS%alphaij')

  end subroutine WS_deallocate

  subroutine WS_allocate_and_init(mixWS,nc)
    use utilities, only: allocate_nc_x_nc
    class(mixWongSandler), intent(inout) :: mixWS
    integer, intent(in) :: nc
    ! Locals
    type(Fraction):: ZeroFraction !<= 0.0/1.0
    integer :: i, j, err
    ZeroFraction%pNum = 0.0
    ZeroFraction%pDen = 0.0
    ZeroFraction%pDen(1) = 1.0
    !
    err = 0
    call mixWS%dealloc()
    allocate (mixWS%f_kij(nc,nc), stat=err)
    if (err /= 0) call stoperror('WS_allocate_and_init: could not allocate array: mixWS%f_kij')
    allocate (mixWS%f_tauij(nc,nc), stat=err)
    if (err /= 0) call stoperror('WS_allocate_and_init: could not allocate array: mixWS%f_tauij')
    call allocate_nc_x_nc(mixWS%Alphaij,nc,"mixWS%Alphaij")
    do i=1,nc
      do j=1,nc
        mixWS%f_kij(i,j) = ZeroFraction
        mixWS%f_Tauij(i,j) = ZeroFraction
        mixWS%Alphaij(i,j) = 0.0
      enddo
    enddo
  end subroutine WS_allocate_and_init

  subroutine assign_WS_mix(mixWS1,mixWS2)
    class(mixWongSandler), intent(inout) :: mixWS1
    class(mixWongSandler), intent(in) :: mixWS2
    if (allocated(mixWS2%f_kij)) then
      call mixWS1%WS_allocate_and_init(size(mixWS2%f_kij,dim=1))
      mixWS1%f_kij = mixWS2%f_kij
      mixWS1%f_Tauij = mixWS2%f_Tauij
      mixWS1%Alphaij = mixWS2%Alphaij
    endif
  end subroutine assign_WS_mix

  subroutine excess_gibbs_deallocate(mixGE)
    class(mixExcessGibbs), intent(inout) :: mixGE
    ! Locals
    integer :: err
    if (allocated (mixGE%alpha)) deallocate (mixGE%alpha, STAT=err)
    if (allocated (mixGE%aGE)) deallocate (mixGE%aGE, STAT=err)
    if (allocated (mixGE%bGE)) deallocate (mixGE%bGE, STAT=err)
    if (allocated (mixGE%cGE)) deallocate (mixGE%cGE, STAT=err)
    if (allocated (mixGE%correlation)) deallocate (mixGE%correlation, STAT=err)
  end subroutine excess_gibbs_deallocate

  subroutine excess_gibbs_allocate_and_init(mixGE,nc)
    class(mixExcessGibbs), intent(inout) :: mixGE
    integer, intent(in) :: nc
    ! Locals
    integer :: i, j
    call mixGE%dealloc()
    allocate (mixGE%alpha(nc,nc))
    allocate (mixGE%aGE(nc,nc))
    allocate (mixGE%bGE(nc,nc))
    allocate (mixGE%cGE(nc,nc))
    allocate (mixGE%correlation(nc,nc))
    do i=1,nc
      do j=1,nc
        mixGE%alpha(i,j) = 0.0
        mixGE%aGE(i,j) = 0.0
        mixGE%bGE(i,j) = 0.0
        mixGE%cGE(i,j) = 0.0
        mixGE%correlation(i,j) = 0
      enddo
   enddo
  end subroutine excess_gibbs_allocate_and_init

  subroutine assign_excess_gibbs_mix(mixGE1,mixGE2)
    class(mixExcessGibbs), intent(inout) :: mixGE1
    class(mixExcessGibbs), intent(in) :: mixGE2
    if (allocated(mixGE2%alpha)) then
      call mixGE1%excess_gibbs_allocate_and_init(size(mixGE2%alpha,dim=1))
      mixGE1%alpha = mixGE2%alpha
      mixGE1%aGE = mixGE2%aGE
      mixGE1%bGE = mixGE2%bGE
      mixGE1%cGE = mixGE2%cGE
      mixGE1%correlation = mixGE2%correlation
    endif
  end subroutine assign_excess_gibbs_mix

  !> Allocate memory for cubic eos
  function cubic_eos_constructor(nc,eos_label) result(cb)
    ! Input:
    integer, intent(in) :: nc
    character(len=*), intent(in) :: eos_label
    ! Created object:
    type(cb_eos) :: cb
    ! Locals

    call cb%allocate_and_init(nc,eos_label)

  end function cubic_eos_constructor

  !> Allocate memory for LK eos
  function lk_eos_constructor(nc,eos_label) result(lk)
    ! Input:
    integer, intent(in) :: nc
    character(len=*), intent(in) :: eos_label
    ! Created object:
    type(lk_eos) :: lk
    !
    call lk%cb_eos%allocate_and_init(nc,eos_label)

  end function lk_eos_constructor

  !> Allocate memory for CPA eos
  function cpa_eos_constructor(nc,eos_label) result(cpa)
    ! Input:
    integer, intent(in) :: nc
    character(len=*), intent(in) :: eos_label
    ! Created object:
    type(cpa_eos) :: cpa
    ! Locals

    call cpa%cb_eos%allocate_and_init(nc,eos_label)

  end function cpa_eos_constructor

  subroutine assign_cubic_eos(this,other)
    use utilities, only: allocate_nc, allocate_nc_x_nc
    class(cb_eos), intent(inout) :: this
    class(*), intent(in) :: other
    ! Locals
    integer :: istat
    select type (other)
    class is (cb_eos)
      call this%assign_base_eos_param(other)
      this%eosid = other%eosid
      this%mruleid = other%mruleid
      this%name = other%name
      this%eosidx = other%eosidx
      this%mruleidx = other%mruleidx
      this%subeosidx = other%subeosidx
      this%cubic_verbose = other%cubic_verbose

      this%m1 = other%m1
      this%m2 = other%m2
      this%dm1dB = other%dm1dB
      this%dm1dC = other%dm1dC
      this%dm2dB = other%dm2dB
      this%dm2dC = other%dm2dC
      this%d2m1dB2 = other%d2m1dB2
      this%d2m1dC2 = other%d2m1dC2
      this%d2m2dB2 = other%d2m2dB2
      this%d2m2dC2 = other%d2m2dC2
      this%d2m2dBdC = other%d2m2dBdC
      this%d2m1dBdC = other%d2m1dBdC
      this%delta = other%delta

      this%a = other%a
      this%b = other%b
      this%c = other%c
      this%sumN = other%sumN
      this%sumA = other%sumA
      this%sumB = other%sumB
      this%sumC = other%sumC

      this%pn = other%pn
      this%pa = other%pa
      this%pb = other%pb
      this%pc = other%pc
      this%pt = other%pt
      this%pv = other%pv

      this%ff = other%ff
      this%fft = other%fft
      this%fftt = other%fftt
      this%ffn = other%ffn
      this%ffnv = other%ffnv
      this%ffna = other%ffna
      this%ffnb = other%ffnb
      this%ffnc = other%ffnc
      this%ffnn = other%ffnn
      this%ffnt = other%ffnt
      this%ffa = other%ffa
      this%ffaa = other%ffaa
      this%ffab = other%ffab
      this%ffac = other%ffac
      this%ffat = other%ffat
      this%ffb = other%ffb
      this%ffbb = other%ffbb
      this%ffbc = other%ffbc
      this%ffbt = other%ffbt
      this%ffc = other%ffc
      this%ffcc = other%ffcc
      this%ffct = other%ffct
      this%ffv = other%ffv
      this%ffvt = other%ffvt
      this%ffvv = other%ffvv
      this%ffva = other%ffva
      this%ffvb = other%ffvb
      this%ffvc = other%ffvc

      this%at = other%at
      this%att = other%att
      this%bt = other%bt
      this%btt = other%btt

      if (allocated(other%aij)) then
        call allocate_nc_x_nc(this%aij,size(other%aij,dim=1),"aij")
        this%aij = other%aij
      endif
      if (allocated(other%bij)) then
        call allocate_nc_x_nc(this%bij,size(other%bij,dim=1),"bij")
        this%bij = other%bij
      endif
      if (allocated(other%cij)) then
        call allocate_nc_x_nc(this%cij,size(other%cij,dim=1),"cij")
        this%cij = other%cij
      endif
      if (allocated(other%ai)) then
        call allocate_nc(this%ai,size(other%ai,dim=1),"ai")
        this%ai = other%ai
      endif
      if (allocated(other%ait)) then
        call allocate_nc(this%ait,size(other%ait,dim=1),"ait")
        this%ait = other%ait
      endif
      if (allocated(other%bi)) then
        call allocate_nc(this%bi,size(other%bi,dim=1),"bi")
        this%bi = other%bi
      endif
      if (allocated(other%bit)) then
        call allocate_nc(this%bit,size(other%bit,dim=1),"bit")
        this%bit = other%bit
      endif
      if (allocated(other%bitt)) then
        call allocate_nc(this%bitt,size(other%bitt,dim=1),"bitt")
        this%bitt = other%bitt
      endif
      if (allocated(other%ci)) then
        call allocate_nc(this%ci,size(other%ci,dim=1),"ci")
        this%ci = other%ci
      endif

      this%extrm = other%extrm
      this%nextrm = other%nextrm
      this%nzfac = other%nzfac

      if (allocated(other%single)) then
        this%single = other%single
      else if (allocated(this%single)) then
        deallocate(this%single, stat=istat)
        if (istat /= 0) call stoperror("cubic_eos: Not able to deallocate single")
      endif

      if (allocated(other%kij)) then
        call allocate_nc_x_nc(this%kij,size(other%kij,dim=1),"kij")
        this%kij = other%kij
      endif
      if (allocated(other%lij)) then
         call allocate_nc_x_nc(this%lij,size(other%lij,dim=1),"lij")
         this%lij = other%lij
      endif
      if (allocated(other%lowcase_bij)) then
        call allocate_nc_x_nc(this%lowcase_bij,size(other%lowcase_bij,dim=1),"lowcase_bij")
        this%lowcase_bij = other%lowcase_bij
      endif
      this%simple_covolmixing = other%simple_covolmixing
      this%mixGE = other%mixGE
      this%mixWS = other%mixWS
      if (associated(other%unifdb)) then
        if (.not. associated(this%unifdb)) then
          allocate(unifacdb :: this%unifdb, stat=istat)
          if (istat /= 0) call stoperror("cubic_eos: Not able to allocate unifacdb")
        endif
        this%unifdb = other%unifdb
      endif
    class default
      print *,"assign_cubic_eos: Should not be here"
    end select

  end subroutine assign_cubic_eos

  subroutine allocate_and_init_cubic_eos(eos,nc,eos_label)
    use utilities, only: allocate_nc, allocate_nc_x_nc
    class(cb_eos), intent(inout) :: eos
    integer, intent(in) :: nc
    character(len=*), intent(in) :: eos_label !< EOS label
    ! Locals
    integer :: err
    err = 0 ! needs initialization due to if clauses below

    call allocate_nc(eos%ai,nc,"eos%ai")
    call allocate_nc(eos%ait,nc,"eos%ait")
    call allocate_nc(eos%bi,nc,"eos%bi")
    call allocate_nc(eos%bit,nc,"eos%bit")
    call allocate_nc(eos%bitt,nc,"eos%bitt")
    call allocate_nc(eos%ci,nc,"eos%ci")

    call allocate_nc_x_nc(eos%aij,nc,"eos%aij")
    call allocate_nc_x_nc(eos%bij,nc,"eos%bij")
    call allocate_nc_x_nc(eos%cij,nc,"eos%cij")

    call allocate_nc_x_nc(eos%kij,nc,"eos%kij")
    call allocate_nc_x_nc(eos%lij,nc,"eos%lij")
    call allocate_nc_x_nc(eos%lowcase_bij,nc,"eos%lowcase_bij")

    if (allocated (eos%single)) deallocate(eos%single, STAT=err)
    if (err /= 0) call stoperror('allocate_cbeos: could not deallocate single')
    allocate (eos%single(nc), STAT=err)
    if (err /= 0) call stoperror('allocate_cbeos: could not allocate single')

    eos%eosid = ""
    eos%mruleid = ""
    eos%name = ""
    eos%eosidx = 0
    eos%mruleidx = 0
    eos%subeosidx = 0
    eos%cubic_verbose = .false.

    eos%m1 = 0
    eos%m2 = 0
    eos%dm1dB = 0
    eos%dm1dC = 0
    eos%dm2dB = 0
    eos%dm2dC = 0
    eos%d2m1dB2 = 0
    eos%d2m1dC2 = 0
    eos%d2m2dB2 = 0
    eos%d2m2dC2 = 0
    eos%d2m2dBdC = 0
    eos%d2m1dBdC = 0
    eos%delta = 0

    eos%a = 0
    eos%b = 0
    eos%c = 0
    eos%sumN = 0
    eos%sumA = 0
    eos%sumB = 0
    eos%sumC = 0

    eos%pn = 0
    eos%pa = 0
    eos%pb = 0
    eos%pc = 0
    eos%pt = 0
    eos%pv = 0

    eos%ff = 0
    eos%fft = 0
    eos%fftt = 0
    eos%ffn = 0
    eos%ffnv = 0
    eos%ffna = 0
    eos%ffnb = 0
    eos%ffnc = 0
    eos%ffnn = 0
    eos%ffnt = 0
    eos%ffa = 0
    eos%ffaa = 0
    eos%ffab = 0
    eos%ffac = 0
    eos%ffat = 0
    eos%ffb = 0
    eos%ffbb = 0
    eos%ffbc = 0
    eos%ffbt = 0
    eos%ffc = 0
    eos%ffcc = 0
    eos%ffct = 0
    eos%ffv = 0
    eos%ffvt = 0
    eos%ffvv = 0
    eos%ffva = 0
    eos%ffvb = 0
    eos%ffvc = 0

    eos%at  = 0
    eos%att = 0
    eos%bt  = 0
    eos%btt = 0

    eos%aij  = 0
    eos%bij  = 0
    eos%cij  = 0
    eos%ai   = 0
    eos%ait  = 0
    eos%bi   = 0
    eos%bit  = 0
    eos%bitt = 0
    eos%ci   = 0

    eos%extrm  = 0
    eos%nextrm = 0
    eos%nzfac  = 0

    eos%kij = 0
    eos%lij = 0
    eos%lowcase_bij = 0
    eos%simple_covolmixing = .true.
  end subroutine allocate_and_init_cubic_eos

  !! \author Morten H
  subroutine cubic_eos_dealloc(eos)
    use utilities, only: deallocate_real, deallocate_real_2
    ! Passed object:
    class(cb_eos), intent(inout) :: eos
    ! Loclas
    integer :: stat
    call base_eos_dealloc(eos)

    call deallocate_real(eos%ai,"eos%ai")
    call deallocate_real(eos%ait,"eos%ait")
    call deallocate_real(eos%bi,"eos%bi")
    call deallocate_real(eos%bit,"eos%bit")
    call deallocate_real(eos%bitt,"eos%bitt")
    call deallocate_real(eos%ci,"eos%ci")

    call deallocate_real_2(eos%aij,"eos%aij")
    call deallocate_real_2(eos%bij,"eos%bij")
    call deallocate_real_2(eos%cij,"eos%cij")

    call deallocate_real_2(eos%kij,"eos%kij")
    call deallocate_real_2(eos%lij,"eos%lij")
    call deallocate_real_2(eos%lowcase_bij,"eos%lowcase_bij")

    stat = 0
    if (allocated(eos%single)) deallocate(eos%single, STAT=stat)
    if (stat /= 0) write (*,*) 'Error deallocating single'

    call eos%mixGE%dealloc()
    call eos%mixWS%dealloc()
    if (associated(eos%unifdb)) then
      call eos%unifdb%dealloc()
      stat = 0
      deallocate(eos%unifdb,STAT=stat)
      if (stat /= 0) write (*,*) 'Error deallocating unifdb'
    endif

  end subroutine cubic_eos_dealloc

  function isHVmixModel(mix_idx) result(isHV)
    integer, intent(in) :: mix_idx
    logical :: isHV
    ! Locals
    integer :: i
    isHV = .false.
    do i=1,nHVCorrs
      if (HVCorrIndices(i) == mix_idx) then
        isHV = .true.
        return
      endif
    enddo
  end function isHVmixModel

  function isGEmixModel(mix_idx) result(isGE)
    integer, intent(in) :: mix_idx
    logical :: isGE
    ! Locals
    integer :: i
    isGE = .false.
    do i=1,nGECorrs
      if (GECorrIndices(i) == mix_idx) then
        isGE = .true.
        return
      endif
    enddo
  end function isGEmixModel

  function get_mix_db_idx(short_label) result(idx)
    character(len=*), intent(in) :: short_label
    integer :: idx
    ! Locals
    integer :: i
    idx = -1
    do i=1,n_mix_rules
      if (str_eq(short_label,mix_label_db(i)%short_label) .or. &
           string_match(short_label,mix_label_db(i)%alias)) then
        idx = i
        return
      endif
    enddo
  end function get_mix_db_idx

  function get_alpha_db_idx(short_label) result(idx)
    character(len=*), intent(in) :: short_label
    integer :: idx
    ! Locals
    integer :: i
    idx = -1
    do i=1,n_alpha_corrs
      if ( str_eq(short_label,alpha_corr_db(i)%short_label)) then
        idx = i
        return
      endif
    enddo
    print *, "Possible alphacorrs:"
    do i=1,n_alpha_corrs
      print *,trim(alpha_corr_db(i)%short_label)
    enddo
    call StopError ('nameToIdx::unknown alpha corr name ')
  end function get_alpha_db_idx

  function get_alpha_db_idx_from_alpha_idx(alpha_idx) result(idx)
    integer, intent(in) :: alpha_idx
    integer :: idx
    ! Locals
    integer :: i
    idx = -1
    do i=1,n_alpha_corrs
      if (alpha_idx == alpha_corr_db(i)%alpha_idx) then
        idx = i
        return
      endif
    enddo
  end function get_alpha_db_idx_from_alpha_idx

  function eos_to_classic_alpha_db_idx(eosidx) result(alpha_idx)
    integer, intent(in) :: eosidx
    integer :: alpha_idx
    ! Locals
    integer :: i
    alpha_idx = -1
    do i=1,n_alpha_corrs
      if (eosidx == alpha_corr_db(i)%classic_for_eos_idx) then
        alpha_idx = i
        exit
      endif
    enddo
  end function eos_to_classic_alpha_db_idx

  function is_classic_alpha(alphaidx) result(is_classic)
    integer, intent(in) :: alphaidx
    logical :: is_classic
    !
    is_classic = (&
         alphaidx == cbAlphaClassicFitIdx .OR. &
         alphaidx == cbAlphaSoaveIdx .OR. &
         alphaidx == cbAlphaPRIdx .OR. &
         alphaidx == cbAlphaPTIdx .OR. &
         alphaidx == cbAlphaSWIdx .OR. &
         alphaidx == cbAlphaVDWIdx)

  end function is_classic_alpha

end module cubic_eos

max_stars_repo_name: jabirali/thermopack

max_issues_repo_name: jabirali/thermopack

max_forks_repo_name: jabirali/thermopack


text:
************************************************************************
* This routine checks the singularity cancellation, given the two      *
* appropriate routines - real_g, sub_gs - and the momentum array p.    *
* A given sub-process (j,k) is determined to be singular if            *
*  (npart=4) msq(j,k) > 1d3 * (smallest value of msq(j,k))             *
*  (npart=5) msq(j,k) > 1d4 * (smallest value of msq(j,k))             *
* This seems to be a reasonable criterion for this set of singular     *
* points. The singularity is considered uncancelled if                 *
*  msq(j,k)/(sum of dipoles) is more than 2% different from 1          *
************************************************************************
      subroutine singcheck(real_g,sub_gs,p)
      implicit none
      include 'constants.f'
      include 'nflav.f'
      include 'npart.f'
      include 'ptilde.f'
      include 'process.f'
      include 'jetlabel.f'
      include 'noglue.f'
      include 'nqcdjets.f'
      include 'incldip.f'
      include 'notag.f'
      integer j,k,jj,kk,jmax,kmax,nd,icount
      double precision p(mxpart,4),q(mxpart,4),pjet(mxpart,4),
     . msq(-nf:nf,-nf:nf),msqs(-nf:nf,-nf:nf),msqc(maxd,-nf:nf,-nf:nf),
     . s(mxpart,mxpart),rcut,debugsmall,debuglarge,xtoler,pttwo,
     . premsq(-nf:nf,-nf:nf),premsqs(-nf:nf,-nf:nf),smin,dot
      character*32 debugmsg
      logical includedipole
      external real_g,sub_gs
      common/rcut/rcut
      data icount/0/
      save icount,premsq,premsqs

      if     (npart .eq. 3) then
        jmax=5
c        if (case .eq. 'dirgam') then
c          kmax=5
c        elseif (case .eq. 'twojet') then
c          kmax=9
        if ((case .eq. 'tt_tot') .or. (case .eq. 'cc_tot')
     .   .or. (case .eq. 'bb_tot')) then
          kmax=4
        else
        jmax=1
        kmax=1
        xtoler=1d3
c        smin=min(-dot(p,1,5),-dot(p,2,5))
c        if (smin. gt. 1d-3) return
        endif
        xtoler=1d3
      elseif (npart .eq. 4) then
        jmax=5
        kmax=5
        xtoler=1d3
        if ( (case .eq. 'qg_tbq') .or .(case .eq. 'qq_tbg')
     .   .or.(case .eq. 'epem3j') ) then
          jmax=1
          kmax=1
        endif
        jmax=1
        kmax=1
      elseif (npart .eq. 5) then
        if ((case .eq. 't_bbar') .or. (case .eq. 'bq_tpq')) then
        jmax=5
        kmax=6
        xtoler=1d3
        elseif (case .eq. 'W_bjet') then
          jmax=1
          kmax=1
          xtoler=1d3
        else
        jmax=1
        kmax=1
        xtoler=1d3
        smin=min(-dot(p,1,5),-dot(p,2,5))
        smin=min(smin,-dot(p,1,6))
        smin=min(smin,-dot(p,2,6))
        smin=min(smin,-dot(p,1,7))
        smin=min(smin,-dot(p,2,7))
        smin=min(smin,+dot(p,5,6))
        smin=min(smin,+dot(p,5,7))
        smin=min(smin,+dot(p,6,7))
c        if (smin. gt. 1d-1) return
        endif
      elseif (npart .eq. 7) then
        jmax=5
        kmax=5
        xtoler=1d3
        if ((case .eq. 'tt_bbl') .or .(case .eq. 'tt_bbh')) then
          jmax=1
          kmax=1
        endif
      else
        write(6,*) 'No singularity check implemented yet'
        stop
      endif

      jmax=1
      kmax=1

      do j=1,jmax
      do k=1,kmax
         if     (npart .eq. 3) then
           if ((case .eq. 'tt_tot') .or. (case .eq. 'cc_tot')
     .   .or. (case .eq. 'bb_tot')) then
c             call coll3m(p,k,j)
           elseif (case .eq. 'dirgam') then
c             call coll3gam(p,k,j)
           elseif (case .eq. 'twojet') then
c             call coll3jet(p,k,j)
           else
c             call coll3(p,k,j)
           endif
         elseif (npart .eq. 4) then
           if     (case .eq. 'W_cjet') then
c             call coll4c(p,k,j)
           elseif (case .eq. 'W_tndk') then
             if (k .eq. 3) goto 68
c             call coll4t(p,k,j)
           elseif ( (case .eq. 'qg_tbq') .or. (case .eq. 'qq_tbg')
     .         .or. (case .eq. 'epem3j') ) then
             write(6,*) 'p(5,4)=',p(5,4)
             write(6,*) 'p(6,4)=',p(6,4)
           else
c             call coll4a(p,k,j)
           endif
         elseif (npart .eq. 5) then
           if ((case .eq. 't_bbar') .or. (case .eq. 'bq_tpq')) then
c             call coll_stop(p,k,j)
           elseif (case .eq. 'W_bjet') then
             write(6,*) 'p(5,4)=',p(5,4)
             write(6,*) 'p(6,4)=',p(6,4)
             write(6,*) 'p(7,4)=',p(7,4)
           else
c             call coll5(p,k,j)
           endif
         elseif (npart .eq. 7) then
           if ((case .eq. 'tt_bbl') .or. (case .eq. 'tt_bbh')) then
             write(6,*) 'p(9,4)=',p(9,4)
           else
c             call coll5(p,k,j)
           endif
         endif
         call real_g(p,msq)
         call sub_gs(p,msqc)
         do jj=-nf,nf
         do kk=-nf,nf
            msqs(jj,kk)=0d0
         enddo
         enddo
         call dotem(9,p,s)
         call genclust2(p,rcut,pjet,0)
c         write(6,*) 'rl,pt34',0,pttwo(3,4,p)
c         write(6,*) 'nd,jets, msq(5,0)',0,nd,jets,msq(5,0)
         if (jets .eq. -1) then
           write(6,*) 'This point does not have a final state b-jet'
           do jj=-nf,nf
           do kk=-nf,nf
              msq(jj,kk)=0d0
           enddo
           enddo
           goto 68
         endif
         if (jets .lt. nqcdjets-notag) then
c           write(6,*) 'This point does not have enough jets'
           do jj=-nf,nf
           do kk=-nf,nf
              msq(jj,kk)=0d0
           enddo
           enddo
           goto 68
         endif

         write(*,*) 'Point ',j,':'
         write(6,*) '2d0*p1.p4',s(1,4)
         write(6,*) '2d0*p2.p4',s(2,4)
         if (npart .gt. 2) then
         write(6,*) '2d0*p1.p5',s(1,5)
         write(6,*) '2d0*p2.p5',s(2,5)
         write(6,*) '2d0*p3.p4',s(3,4)
         write(6,*) '2d0*p3.p5',s(3,5)
         write(6,*) '2d0*p4.p5',s(4,5)
         endif
         if (npart .gt. 3) then
         write(6,*) '2d0*p1.p6',s(1,6)
         write(6,*) '2d0*p2.p6',s(2,6)
         write(6,*) '2d0*p3.p6',s(3,6)
         write(6,*) '2d0*p4.p6',s(4,6)
         write(6,*) '2d0*p5.p6',s(5,6)
         endif
         if (npart .gt. 4) then
         write(6,*) '2d0*p1.p7',s(1,7)
         write(6,*) '2d0*p2.p7',s(2,7)
         write(6,*) '2d0*p5.p7',s(5,7)
         write(6,*) '2d0*p6.p7',s(6,7)
         endif

c         write(6,*) '0',jets,msq(2,0),(pjet(5,4)+pjet(6,4))**2
c     .                      -(pjet(5,1)+pjet(6,1))**2
c     .                      -(pjet(5,2)+pjet(6,2))**2
c     .                      -(pjet(5,3)+pjet(6,3))**2
         do nd=1,ndmax
            do jj=1,9
            do kk=1,4
               q(jj,kk)=ptilde(nd,jj,kk)
            enddo
            enddo
            if (incldip(nd) .eqv. .false.) goto 72
c            write(6,*) 'nd,pt34',nd,pttwo(3,4,q)
c            if ((nd .eq. 1) .or. (nd .eq. 8)) then
c            write(6,*) 'KINEMATICS FOR DIPOLE ',nd
c            call writeout(q)
c            endif
            call dotem(9,q,s)
c            call genclust2(q,rcut,pjet,1)
c            if (jets .ge. nqcdjets) then

            if (includedipole(nd,q)) then
               write(6,*) nd
c            write(6,*) nd,jets,msqc(nd,5,0)
               do jj=-nflav,nflav
               do kk=-nflav,nflav
                  msqs(jj,kk)=msqs(jj,kk)+msqc(nd,jj,kk)
               enddo
               enddo
            endif
c            write(6,*) 'nd,jets,msqc(nd,2,0)',nd,jets,msqc(nd,2,0)
   72    continue
         enddo

c--- find smallest value of msq
         debugsmall=msq(0,0)
         debuglarge=msq(0,0)
         do jj=-nflav,nflav
         do kk=-nflav,nflav
           if ((msq(jj,kk) .lt. debugsmall)
     .     .and. (msq(jj,kk) .gt. 0d0)) debugsmall=msq(jj,kk)
           if ((msq(jj,kk) .gt. debuglarge)
     .     .and. (msq(jj,kk) .gt. 0d0)) debuglarge=msq(jj,kk)
         enddo
         enddo

         if ((debuglarge/debugsmall .lt. xtoler)
     .  .or. (debuglarge .lt. 1d-10)) then
c           write(6,*) ' OK - no singular configurations'
c           goto 68
         endif

         if (debugsmall .gt. 1d-3) debugsmall=debugsmall/xtoler/2d0
c         if (debuglarge/debugsmall .lt. xtoler*1d2)
c     .       debugsmall=debuglarge/xtoler

         do jj=-nflav,nflav
         do kk=-nflav,nflav

c         if ((gqonly) .and. (jj.ne.0) .and. (kk.ne.0)) goto 69
c         if ((jj .gt. 0) .or. (kk .gt. 0)) goto 69 ! DEBUG: QBG & GQB ONLY
c         if ((jj .le. 0) .or. (kk .le. 0)) goto 69 ! DEBUG: QQ only

         if (msq(jj,kk) .eq. 0d0) then
            if (msqs(jj,kk) .eq. 0d0) then
               goto 69
               debugmsg='   OK   zero msq and subtraction'
            else
               debugmsg=' FAILED subtraction with msq=0'
            endif
            write(*,22) jj,kk,'    n/a   ',msq(jj,kk),msqs(jj,kk),
     .      debugmsg
            goto 69
            endif

         if ((msq(jj,kk)/debuglarge) .lt. 1d0/xtoler) then
            if (abs(abs(msq(jj,kk)/msqs(jj,kk))-1d0) .lt. 0.02d0)then
              debugmsg='   OK   cancelled'
            else
              debugmsg='   OK   not singular'
            endif
            write(*,21) jj,kk,msq(jj,kk)/msqs(jj,kk),
     .                    msq(jj,kk),msqs(jj,kk),debugmsg
         else
             if (msqs(jj,kk) .eq. 0d0) then
                debugmsg=' FAILED singular, no subtraction'
                write(*,22) jj,kk,'    n/a   ',msq(jj,kk),msqs(jj,kk),
     . debugmsg
                goto 69
             endif
             if (abs(abs(msq(jj,kk)/msqs(jj,kk))-1d0) .gt. 0.02d0)then
                debugmsg=' FAILED singularity uncancelled'
             else
                debugmsg='   OK   cancelled'
             endif
             write(*,21) jj,kk,msq(jj,kk)/msqs(jj,kk),
     .                  msq(jj,kk),msqs(jj,kk),debugmsg
         endif

   69    continue
         enddo
         enddo

   68 continue
      enddo
      pause
      enddo

c---c--- This block of code is useful for checking the symmetry of
c---c--- the real matrix elements and subtractions
c---      icount=icount+1
c---      if (icount .eq. 1) then
c---c--- save contents of arrays for reals and subtractions
c---        do j=-nf,nf
c---       do k=-nf,nf
c---       premsq(j,k)=msq(j,k)
c---       premsqs(j,k)=msqs(j,k)
c---       enddo
c---       enddo
c---      else
c---c--- compare reals and subtractions with previous iteration
c---        do j=-nf,nf
c---        do k=-nf,nf
c---c        if (j*k .gt. 0) goto 67 ! NO QQ OR QBQB (for Z+2jet 1<>2 symm)
c---c        if (j*k .le. 0) goto 67 ! ONLY QQ/QBQB (for Z+2jet 1<>2,5<>6 symm)
c---c       if (j*k .lt. 0) goto 67 ! DEBUG: NO QQB OR QBQ
c---c       if (j*k .ge. 0) goto 67 ! DEBUG: QQB & QBQ only
c---c      if (j*k .le. 0) goto 67 ! DEBUG: QQ & QBQB only
c---c      if ((j .le. 0) .or. (k .le. 0)) goto 67 ! DEBUG: QQ only
c---c       if ((j*k .gt. 0).or.(abs(j).eq.5).or.(abs(k).eq.5)) goto 67 ! NO QQ OR QBQB, OR B QUARK
c---        if ((abs(msq(j,k)) .lt. 1d-15) .and. (abs(msqs(j,k)).lt.1d-15))
c---     .                    goto 67 ! skip zero contributions
c---        write(6,23) j,k,premsq(j,k),msq(k,j),premsqs(j,k),msqs(k,j),
c---     .                  premsq(j,k)/msq(k,j),premsqs(j,k)/msqs(k,j)
c---   67        continue
c---        enddo
c---        enddo
c---        icount=0
c---        pause
c---      endif

c      write(6,*) 'Singularity check completed'
c      write(6,*) 'Do ''mcfm | grep -v OK'' to look for failures'
c      stop

      return

   21 format(1x,2i3,f10.6,2e14.6,1x,a32)
   22 format(1x,2i3,a10,2e14.6,1x,a32)
   23 format(2i4,4e12.5,2f10.6)

      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!
! wavy - A spectral ocean wave modeling and development framework
! Copyright (c) 2017, Wavebit Scientific LLC
! All rights reserved.
!
! Licensed under the BSD 3-clause license. See LICENSE for details.

module mod_time_integration

use mod_precision,only:rk => realkind
use mod_spectrum,only:spectrum_type
use mod_domain,only:domain_type

implicit none

private

public :: backward_euler
public :: exact_exponential
public :: forward_euler
public :: integrate

interface integrate
  module procedure :: integrate_spectrum
  module procedure :: integrate_domain
endinterface integrate

interface exact_exponential
  module procedure :: exact_exponential_spectrum
  module procedure :: exact_exponential_domain
endinterface exact_exponential

interface backward_euler
  module procedure :: backward_euler_spectrum
  module procedure :: backward_euler_domain
endinterface backward_euler

interface forward_euler
  module procedure :: forward_euler_spectrum
  module procedure :: forward_euler_domain
endinterface forward_euler

contains

!-------------------------------------------------------------------------------
pure type(spectrum_type) function integrate_spectrum(func,initial,tendency,dt)
  !! Integrates spectrum forward in time using a time integration method
  !! provided as the argument `func`.
  interface
    pure type(spectrum_type) function func(initial,tendency,dt)
      import :: spectrum_type,rk
      type(spectrum_type),intent(in) :: initial
      type(spectrum_type),intent(in) :: tendency
      real(kind=rk),intent(in) :: dt
    endfunction func
  endinterface
  type(spectrum_type),intent(in) :: initial !! Initial spectrum instance
  type(spectrum_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  integrate_spectrum = func(initial,tendency,dt)
endfunction integrate_spectrum
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure type(domain_type) function integrate_domain(func,initial,tendency,dt)
  !! Integrates domain forward in time using a time integration method
  !! provided as the argument `func`.
  interface
    pure type(domain_type) function func(initial,tendency,dt)
      import :: domain_type,rk
      type(domain_type),intent(in) :: initial
      type(domain_type),intent(in) :: tendency
      real(kind=rk),intent(in) :: dt
    endfunction func
  endinterface
  type(domain_type),intent(in) :: initial !! Initial domain instance
  type(domain_type),intent(in) :: tendency !! Spectrum domain instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  integrate_domain = func(initial,tendency,dt)
endfunction integrate_domain
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(spectrum_type) function backward_euler_spectrum(initial,&
  tendency,dt) result(spec)
  !! Integrates a spectrum forward in time using a 1st order implicit backward
  !! Euler integration scheme.
  type(spectrum_type),intent(in) :: initial !! Initial spectrum instance
  type(spectrum_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  spec = initial / (- tendency * dt + 1._rk)
endfunction backward_euler_spectrum
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(domain_type) function backward_euler_domain(initial,&
  tendency,dt) result(domain)
  !! Integrates a domain instance forward in time using a 1st order implicit
  !! backward Euler integration scheme.
  type(domain_type),intent(in) :: initial !! Initial spectrum instance
  type(domain_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  domain = initial
  call domain % setSpectrum(backward_euler(initial % getSpectrum(),&
    tendency % getSpectrum(),dt))
endfunction backward_euler_domain
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(spectrum_type) function exact_exponential_spectrum(initial,&
  tendency,dt) result(spec)
  !! Integrates a spectrum instance forward in time using the exact exponential.
  type(spectrum_type),intent(in) :: initial !! Initial spectrum instance
  type(spectrum_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  spec = initial
  call spec % setSpectrum(exp(tendency % getSpectrum() * dt))
  !spec = (initial * exp(tendency % getSpectrum() * dt))
  spec = initial * spec
endfunction exact_exponential_spectrum
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(domain_type) function exact_exponential_domain(initial,&
  tendency,dt) result(domain)
  !! Integrates a domain instance forward in time using the exact exponential.
  type(domain_type),intent(in) :: initial !! Initial spectrum instance
  type(domain_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  domain = initial
  domain = exact_exponential(initial % getSpectrum(),&
    tendency % getSpectrum(),dt)
endfunction exact_exponential_domain
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(spectrum_type) function forward_euler_spectrum(initial,&
  tendency,dt) result(spec)
  !! Integrates a spectrum forward in time using a 1st order Euler integration
  !! scheme.
  type(spectrum_type),intent(in) :: initial !! Initial spectrum instance
  type(spectrum_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  spec = initial + tendency * dt
endfunction forward_euler_spectrum
!-------------------------------------------------------------------------------

!-------------------------------------------------------------------------------
pure elemental type(domain_type) function forward_euler_domain(initial,&
  tendency,dt) result(domain)
  !! Integrates a domain instance forward in time using a 1st order Euler
  !! integration scheme.
  type(domain_type),intent(in) :: initial !! Initial spectrum instance
  type(domain_type),intent(in) :: tendency !! Spectrum tendency instance
  real(kind=rk),intent(in) :: dt !! Time step [s]
  domain = initial
  domain = forward_euler(initial % getSpectrum(),tendency % getSpectrum(),dt)
endfunction forward_euler_domain
!-------------------------------------------------------------------------------
endmodule mod_time_integration

max_stars_repo_name: dksasaki/wavy

max_issues_repo_name: dksasaki/wavy

max_forks_repo_name: dksasaki/wavy