eckstein_yields.F90

Go to the documentation of this file.

module eckstein_yields

  use itm_types
        
  real(R8), parameter :: mathpi = 3.1415926535897932384626433832795
  real(R8), parameter :: degtorad = 0.01745329251994329576923690768489

contains
        
  function etf(M1, M2, Z1, Z2)
    real(R8) :: etf, m1, m2, z1, z2
    etf = (30.74*(m1 + m2)*z1*sqrt(z1**0.6666666666666666 + z2**0.6666666666666666)*z2)/m2
  end function etf
        
  function omegal(epsl)
    real(R8) :: omegal, epsl
    omegal = 0.008*(epsl**0.1504) + 0.1728*sqrt(epsl) + epsl
  end function omegal
        
  function sn(epsl)
    real(R8) :: sn, epsl
    sn = (0.5*log(1 + 1.2288*epsl))/omegal(epsl)
  end function sn
        
  function seyield(E0, M1, M2, Z1, Z2, q, lambda, u, ETh)
    real(R8) :: seyield, e0, m1, m2, z1, z2, q, lambda, u, eth, epsl, snucl, omega, etoeth
    epsl = e0/etf(m1, m2, z1, z2)
    snucl = sn(epsl)
    omega = omegal(epsl)
    etoeth =  ((e0/eth) - 1.0)**u
    seyield = (etoeth*q*snucl)/(etoeth + (lambda/omega))
  end function seyield
        
  function sayield(E0, theta, f, b, c, Esp)
    real(R8) :: sayield, e0, theta, thetarad, f, b, c, esp
    real(R8) :: thetastar, cosfact
    thetarad = theta * degtorad
    thetastar = mathpi - acos(sqrt(1/(1 + e0/esp)))
    cosfact = cos(((mathpi/2.)**c)*((thetarad/thetastar)**c))
    sayield = exp(b*(1 - 1/cosfact))/(cosfact**f)
  end function sayield

  function reyieldlight(E0, M1, M2, Z1, Z2, A1, A2, A3, A4)
    real(R8) :: reyield, e0, m1, m2, z1, z2, a1, a2, a3, a4, eps

    eps = e0/etf(m1, m2, z1, z2)

    reyieldlight = (a1 * eps**a2) / (1.0 + a3 * eps**a4)

  end function reyieldlight

  function reyieldself(E0, M1, M2, Z1, Z2, A1, A2, A3, A4)
    real(R8) :: reyield, e0, m1, m2, z1, z2, a1, a2, a3, a4, eps

    eps = e0/etf(m1, m2, z1, z2)

    if(a3 * eps**a4 .lt. 700) then  ! prevent overflow DPC
       reyieldself = exp(a1*eps**a2)/(1.0 + exp(a3 * eps**a4))
    else
       reyieldself = exp(a1*eps**a2 - a3 * eps**a4)
    endif

  end function reyieldself

  function reyield(E0, M1, M2, Z1, Z2, A1, A2, A3, A4)
    real(R8) :: reyield, e0, m1, m2, z1, z2, a1, a2, a3, a4, eps
    if (m1 .ge. m2) then
!       write(*,*) 'M1 = ', M1, ' is >= M2 = ',M2,' Using  ryieldself'
       reyield = reyieldself(e0, m1, m2, z1, z2, a1, a2, a3, a4)
    else
!       write(*,*) 'M1 = ', M1, ' is <  M2 = ',M2,' Using  ryieldlight'
       reyield = reyieldlight(e0, m1, m2, z1, z2, a1, a2, a3, a4)
    end if
  end function reyield

  function rayield(angledeg, C1, C2, C3, C4)
    real(R8) :: rayield, angledeg, c1, c2, c3, c4

    rayield = c1 + c2 * tanh(c3 * angledeg * degtorad + c4)

  end function rayield

end module eckstein_yields