text:
! PsychroLib version 2.0.0 (https://github.com/psychrometrics/psychrolib)
! Copyright (c) 2018 D. Thevenard and D. Meyer for the current library implementation
! Copyright (c) 2017 ASHRAE Handbook — Fundamentals for ASHRAE equations and coefficients
! Licensed under the MIT License.

module psychrolib
  !+ Module overview
  !+  Contains functions for calculating thermodynamic properties of gas-vapor mixtures
  !+  and standard atmosphere suitable for most engineering, physical, and meteorological
  !+  applications.
  !+
  !+  Most of the functions are an implementation of the formulae found in the
  !+  2017 ASHRAE Handbook - Fundamentals, in both International System (SI),
  !+  and Imperial (IP) units. Please refer to the information included in
  !+  each function for their respective reference.
  !+
  !+ Example
  !+  use psychrolib, only: GetTDewPointFromRelHum, SI
  !+  ! Set the unit system, for example to SI (can be either 'SI' or 'IP')
  !+  SetUnitSystem(SI)
  !+  ! Calculate the dew point temperature for a dry bulb temperature of 25 C and a relative humidity of 80%
  !+  print *, GetTDewPointFromRelHum(25.0, 0.80)
  !+ 21.3094
  !+
  !+ Copyright
  !+  - For the current library implementation
  !+     Copyright (c) 2018 D. Thevenard and D. Meyer.
  !+  - For equations and coefficients published ASHRAE Handbook — Fundamentals, Chapter 1
  !+     Copyright (c) 2017 ASHRAE Handbook — Fundamentals (https://www.ashrae.org)
  !+
  !+ License
  !+  MIT (https://github.com/psychrometrics/psychrolib/LICENSE.txt)
  !+
  !+ Note from the Authors
  !+  We have made every effort to ensure that the code is adequate, however, we make no
  !+  representation with respect to its accuracy. Use at your own risk. Should you notice
  !+  an error, or if you have a suggestion, please notify us through GitHub at
  !+  https://github.com/psychrometrics/psychrolib/issues.

  implicit none

  private
  public :: IP
  public :: SI
  public :: SetUnitSystem
  public :: GetUnitSystem
  public :: isIP
  public :: GetTRankineFromTFahrenheit
  public :: GetTKelvinFromTCelsius
  public :: GetTWetBulbFromTDewPoint
  public :: GetTWetBulbFromRelHum
  public :: GetRelHumFromTDewPoint
  public :: GetRelHumFromTWetBulb
  public :: GetTDewPointFromRelHum
  public :: GetTDewPointFromTWetBulb
  public :: GetVapPresFromRelHum
  public :: GetRelHumFromVapPres
  public :: GetTDewPointFromVapPres
  public :: GetVapPresFromTDewPoint
  public :: GetTWetBulbFromHumRatio
  public :: GetHumRatioFromTWetBulb
  public :: GetHumRatioFromRelHum
  public :: GetRelHumFromHumRatio
  public :: GetHumRatioFromTDewPoint
  public :: GetTDewPointFromHumRatio
  public :: GetHumRatioFromVapPres
  public :: GetVapPresFromHumRatio
  public :: GetDryAirEnthalpy
  public :: GetDryAirDensity
  public :: GetDryAirVolume
  public :: GetSatVapPres
  public :: GetSatHumRatio
  public :: GetSatAirEnthalpy
  public :: GetVaporPressureDeficit
  public :: GetDegreeOfSaturation
  public :: GetMoistAirEnthalpy
  public :: GetMoistAirVolume
  public :: GetMoistAirDensity
  public :: GetStandardAtmPressure
  public :: GetStandardAtmTemperature
  public :: GetSeaLevelPressure
  public :: GetStationPressure
  public :: CalcPsychrometricsFromTWetBulb
  public :: CalcPsychrometricsFromTDewPoint
  public :: CalcPsychrometricsFromRelHum

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Global constants
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  real, parameter ::  R_DA_IP = 53.350
    !+ Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

  real, parameter ::  R_DA_SI = 287.042
    !+ Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

  integer, parameter :: IP = 1
  integer, parameter :: SI = 2

  integer  :: PSYCHROLIB_UNITS = 0 ! 0 = undefined
    !+ Unit system to use.

  real ::  PSYCHROLIB_TOLERANCE = 1.0
    !+ Tolerance of temperature calculations

  contains

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Helper functions
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine SetUnitSystem(UnitSystem)
  !+ Set the system of units to use (SI or IP).
  !+ Notes: this function *HAS TO BE CALLED* before the library can be used

    integer, intent(in)    :: UnitSystem
    !+ Units: string indicating the system of units chosen (SI or IP)

    if (.not. (UnitSystem == SI .or. UnitSystem == IP)) then
      error stop "The system of units has to be either SI or IP."
    end if

    PSYCHROLIB_UNITS = UnitSystem

    ! Define tolerance on temperature calculations
    ! The tolerance is the same in IP and SI
    if (UnitSystem == IP) then
      PSYCHROLIB_TOLERANCE = 0.001 * 9.0 / 5.0
    else
      PSYCHROLIB_TOLERANCE = 0.001
    end if
  end subroutine SetUnitSystem

  function GetUnitSystem() result(UnitSystem)
    !+ Return the system of units in use.
    integer :: UnitSystem
    UnitSystem = PSYCHROLIB_UNITS
  end function GetUnitSystem

  function isIP()
    !+ Check whether the system in use is IP or SI
    logical :: isIP
    if (PSYCHROLIB_UNITS == IP) then
      isIP = .true.
    else if (PSYCHROLIB_UNITS == SI) then
      isIP = .false.
    else
      error stop "The system of units has not been defined."
    end if
  end function isIP

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Conversion between temperature units
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetTRankineFromTFahrenheit(TFahrenheit) result(TRankine)
    real, intent(in)  :: TFahrenheit
      !+ Temperature in degree Fahrenheit
    real              :: TRankine
      !+ Temperature in degree Rankine
    real              :: ZERO_FAHRENHEIT_AS_RANKINE
      !+ Zero degree Fahrenheit (°F) expressed as degree Rankine (°R)

    ZERO_FAHRENHEIT_AS_RANKINE = 459.67

    TRankine = TFahrenheit + ZERO_FAHRENHEIT_AS_RANKINE
  end function GetTRankineFromTFahrenheit

  function GetTKelvinFromTCelsius(TCelsius) result(TKelvin)
    real, intent(in)  :: TCelsius
      !+ Temperature in degree Celsius
    real              :: TKelvin
      !+ Tempearatyre in Kelvin
    real              :: ZERO_CELSIUS_AS_KELVIN
      ! Zero degree Fahrenheit (°F) expressed as Kelvin (K)

    ZERO_CELSIUS_AS_KELVIN = 273.15

    TKelvin = TCelsius + ZERO_CELSIUS_AS_KELVIN
  end function GetTKelvinFromTCelsius

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Conversions between dew point, wet bulb, and relative humidity
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetTWetBulbFromTDewPoint(TDryBulb, TDewPoint, Pressure) result(TWetBulb)
    !+ Return wet-bulb temperature given dry-bulb temperature, dew-point temperature, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  :: TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  :: TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real, intent(in)  :: Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              :: TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real              :: HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (TDewPoint > TDryBulb) then
      error stop "Error: dew point temperature is above dry bulb temperature"
    end if

    HumRatio = GetHumRatioFromTDewPoint(TDewPoint, Pressure)
    TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure)
  end function GetTWetBulbFromTDewPoint

  function GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) result(TWetBulb)
    !+ Return wet-bulb temperature given dry-bulb temperature, relative humidity, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  RelHum
      !+ Relative humidity in range [0, 1]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (RelHum < 0.0 .or. RelHum > 1.0) then
      error stop "Error: relative humidity is outside range [0,1]"
    end if

    HumRatio = GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure)
    TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure)
  end function GetTWetBulbFromRelHum

  function GetRelHumFromTDewPoint(TDryBulb, TDewPoint) result(RelHum)
    !+ Return relative humidity given dry-bulb temperature and dew-point temperature.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 22

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real              ::  RelHum
      !+ Relative humidity in range [0, 1]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real              ::  SatVapPres
      !+ Vapor pressure of saturated air in Psi [IP] or Pa [SI]

    if (TDewPoint > TDryBulb) then
      error stop "Error: dew point temperature is above dry bulb temperature"
    end if

    VapPres     = GetSatVapPres(TDewPoint)
    SatVapPres  = GetSatVapPres(TDryBulb)
    RelHum      = VapPres / SatVapPres
  end function GetRelHumFromTDewPoint

  function GetRelHumFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(RelHum)
    !+ Return relative humidity given dry-bulb temperature, wet bulb temperature and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  RelHum
      !+ Relative humidity in range [0, 1]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (TWetBulb > TDryBulb) then
      error stop "Error: wet bulb temperature is above dry bulb temperature"
    end if

    HumRatio = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)
    RelHum   = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure)
  end function GetRelHumFromTWetBulb

  function GetTDewPointFromRelHum(TDryBulb, RelHum) result(TDewPoint)
    !+ Return dew-point temperature given dry-bulb temperature and relative humidity.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  RelHum
      !+ Relative humidity in range [0, 1]
    real              ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (RelHum < 0.0 .or. RelHum > 1.0) then
      error stop "Error: relative humidity is outside range [0,1]"
    end if

    VapPres   = GetVapPresFromRelHum(TDryBulb, RelHum)
    TDewPoint = GetTDewPointFromVapPres(TDryBulb, VapPres)
  end function GetTDewPointFromRelHum

  function GetTDewPointFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(TDewPoint)
    !+ Return dew-point temperature given dry-bulb temperature, wet-bulb temperature, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (TWetBulb > TDryBulb) then
      error stop "Error: wet bulb temperature is above dry bulb temperature"
    end if

    HumRatio  = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)
    TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure)
  end function GetTDewPointFromTWetBulb

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Conversions between dew point, or relative humidity and vapor pressure
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetVapPresFromRelHum(TDryBulb, RelHum) result(VapPres)
    !+ Return partial pressure of water vapor as a function of relative humidity and temperature.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 12, 22

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  RelHum
      !+ Relative humidity in range [0, 1]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (RelHum < 0.0 .or. RelHum > 1.0) then
      error stop "Error: relative humidity is outside range [0,1]"
    end if

    VapPres = RelHum * GetSatVapPres(TDryBulb)
  end function GetVapPresFromRelHum

  function GetRelHumFromVapPres(TDryBulb, VapPres) result(RelHum)
    !+ Return relative humidity given dry-bulb temperature and vapor pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 12, 22

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real              ::  RelHum
      !+ Relative humidity in range [0, 1]

    if (VapPres < 0.0) then
      error stop "Error: partial pressure of water vapor in moist air cannot be negative"
    end if

    RelHum = VapPres / GetSatVapPres(TDryBulb)
  end function GetRelHumFromVapPres

  function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint)
    !+ Return dew-point temperature given dry-bulb temperature and vapor pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 and 6
    !+ Notes:
    !+ The dew point temperature is solved by inverting the equation giving water vapor pressure
    !+ at saturation from temperature rather than using the regressions provided
    !+ by ASHRAE (eqn. 37 and 38) which are much less accurate and have a
    !+ narrower range of validity.
    !+ The Newton-Raphson (NR) method is used on the logarithm of water vapour
    !+ pressure as a function of temperature, which is a very smooth function
    !+ Convergence is usually achieved in 3 to 5 iterations.
    !+ TDryBulb is not really needed here, just used for convenience.

    real, intent(in)    ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real                ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real                ::  Tdp
      !+ Calculated value of dew point temperatures, solved for iteratively in °F [IP] or °C [SI]
    real                ::  lnVP
      !+ Natural logarithm of partial pressure of water vapor pressure in moist air
    real                ::  d_lnVP
      !+ Derivative of function, calculated numerically
    real                ::  lnVP_c
      !+ Value of log of vapor water pressure used in NR calculation
    real                ::  Tdp_c
      !+ Value of Tdp used in NR calculation
    real                ::  d_Tdp
      !+ Value of temperature step used in NR calculation
    real                ::  STEPSIZE
      !+ Size of timestep (dimensionless)
    real                :: TMidPoint
      !+ Midpoint of domain of validity
    real, dimension(2)  ::  BOUNDS
      !+ Valid temperature range in °F [IP] or °C [SI]
    integer             :: index
      !+ Index used in the calculation

    ! Bounds and step size as a function of the system of units
    if (isIP()) then
        BOUNDS(1) = -148.0
        BOUNDS(2) =  392.0
        STEPSIZE  =  0.01 * 9.0 / 5.0
    else
        BOUNDS(1) = -100.0
        BOUNDS(2) =  200.0
        STEPSIZE  =  0.01
    end if

    TMidPoint = (BOUNDS(2) + BOUNDS(2)) / 2.0

    ! Bounds outside which a solution cannot be found
    if (VapPres <= GetSatVapPres(BOUNDS(1)) .or. VapPres >= GetSatVapPres(BOUNDS(2))) then
      error stop "Error: partial pressure of water vapor is outside range of validity of equations"
    end if

    ! First guess
    Tdp = TDryBulb
    lnVP = log(VapPres)
    Tdp_c = Tdp
    index = 1

    do while ((abs(Tdp - Tdp_c) > PSYCHROLIB_TOLERANCE) .or. (index < 2))
      ! Current point
      Tdp_c = Tdp
      lnVP_c = log(GetSatVapPres(Tdp_c))
      ! Step - negative in the right part of the curve, positive in the left part
      ! to avoid going past the domain of validity of eqn. 5 and 6
      ! when Tdp_c is close to its bounds
      if (Tdp_c > TMidPoint) then
        d_Tdp = -STEPSIZE
      else
        d_Tdp = STEPSIZE
      end if
      ! Derivative of function, calculated numerically
      d_lnVP = (log(GetSatVapPres(Tdp_c + d_Tdp)) - lnVP_c) / d_Tdp
      ! New estimate, bounded by domain of validity of eqn. 5 and 6
      Tdp = Tdp_c - (lnVP_c - lnVP) / d_lnVP
      Tdp = max(Tdp, BOUNDS(1))
      Tdp = min(Tdp, BOUNDS(2))
      index = index + 1
    end do

  TDewPoint = min(Tdp, TDryBulb)
  end function GetTDewPointFromVapPres

  function GetVapPresFromTDewPoint(TDewPoint) result(VapPres)
    !+ Return vapor pressure given dew point temperature.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 36

    real, intent(in)  ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    VapPres = GetSatVapPres(TDewPoint)
  end function GetVapPresFromTDewPoint

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Conversions from wet-bulb temperature, dew-point temperature, or relative humidity to humidity ratio
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) result(TWetBulb)
    !+ Return wet-bulb temperature given dry-bulb temperature, humidity ratio, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35 solved for Tstar

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real              ::  TDewPoint
      !+ TDewPoint : Dew-point temperature in °F [IP] or °C [SI]
    real              ::  TWetBulbSup
      !+ Upper value of wet bulb temperature in bissection method (initial guess is from dry bulb temperature) in °F [IP] or °C [SI]
    real              ::  TWetBulbInf
      !+ Lower value of wet bulb temperature in bissection method (initial guess is from dew point temperature) in °F [IP] or °C [SI]
    real              ::  Wstar
      !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio cannot be negative"
    end if

    TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure)

    ! Initial guesses
    TWetBulbSup = TDryBulb
    TWetBulbInf = TDewPoint
    TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.0

    ! Bisection loop
    do while(TWetBulbSup - TWetBulbInf > PSYCHROLIB_TOLERANCE)

    ! Compute humidity ratio at temperature Tstar
    Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)

    ! Get new bounds
    if (Wstar > HumRatio) then
      TWetBulbSup = TWetBulb
    else
      TWetBulbInf = TWetBulb
    end if

    ! New guess of wet bulb temperature
    TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.0
    end do
  end function GetTWetBulbFromHumRatio

  function GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(HumRatio)
    !+ Return humidity ratio given dry-bulb temperature, wet-bulb temperature, and pressure.
    !+ References:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  Wsstar
      !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (TWetBulb > TDryBulb) then
      error stop "Error: wet bulb temperature is above dry bulb temperature"
    end if

    Wsstar = GetSatHumRatio(TWetBulb, Pressure)

    if (isIP()) then
      if (TWetBulb >= 32.0) then
        HumRatio = ((1093.0 - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb))   &
                   / (1093.0 + 0.444 * TDryBulb - TWetBulb)
      else
        HumRatio = ((1220.0 - 0.04 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb))    &
                   / (1220.0 + 0.444 * TDryBulb - 0.48 * TWetBulb)
      end if
    else
      if (TWetBulb >= 0.0) then
        HumRatio = ((2501.0 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb))   &
                   / (2501.0 + 1.86 * TDryBulb - 4.186 * TWetBulb)
      else
           HumRatio = ((2830.0 - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) &
                      / (2830.0 + 1.86 * TDryBulb - 2.1 * TWetBulb)
      end if
    end if
  end function GetHumRatioFromTWetBulb

  function GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure) result(HumRatio)
    !+ Return humidity ratio given dry-bulb temperature, relative humidity, and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  RelHum
      !+ Relative humidity in range [0, 1]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (RelHum < 0.0 .or. RelHum > 1.0) then
      error stop "Error: relative humidity is outside range [0,1]"
    end if

    VapPres   = GetVapPresFromRelHum(TDryBulb, RelHum)
    HumRatio  = GetHumRatioFromVapPres(VapPres, Pressure)
  end function GetHumRatioFromRelHum

  function GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) result(RelHum)
  !+ Return relative humidity given dry-bulb temperature, humidity ratio, and pressure.
  !+ Reference:
  !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  RelHum
      !+ Relative humidity in range [0, 1]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio cannot be negative"
    end if

    VapPres = GetVapPresFromHumRatio(HumRatio, Pressure)
    RelHum  = GetRelHumFromVapPres(TDryBulb, VapPres)
  end function GetRelHumFromHumRatio

  function GetHumRatioFromTDewPoint(TDewPoint, Pressure) result(HumRatio)
    !+ Return humidity ratio given dew-point temperature and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    VapPres   = GetSatVapPres(TDewPoint)
    HumRatio  = GetHumRatioFromVapPres(VapPres, Pressure)
  end function GetHumRatioFromTDewPoint

  function GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) result(TDewPoint)
    !+ Return dew-point temperature given dry-bulb temperature, humidity ratio, and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio cannot be negative"
    end if

    VapPres   = GetVapPresFromHumRatio(HumRatio, Pressure)
    TDewPoint = GetTDewPointFromVapPres(TDryBulb, VapPres)
  end function GetTDewPointFromHumRatio

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Conversions between humidity ratio and vapor pressure
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetHumRatioFromVapPres(VapPres, Pressure) result(HumRatio)
    !+ Return humidity ratio given water vapor pressure and atmospheric pressure.
    !+ Reference:
    !+ ASHRAE Fundamentals (2005) ch. 6 eqn. 22;
    !+ ASHRAE Fundamentals (2009) ch. 1 eqn. 22.

    real, intent(in)  ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    if (VapPres < 0.0) then
      error stop "Error: partial pressure of water vapor in moist air cannot be negative"
    end if

    HumRatio = 0.621945 * VapPres / (Pressure-VapPres)
  end function GetHumRatioFromVapPres

  function GetVapPresFromHumRatio(HumRatio, Pressure) result(VapPres)
    !+ Return vapor pressure given humidity ratio and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 20 solved for pw

    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    VapPres = Pressure * HumRatio / (0.621945 + HumRatio)
  end function GetVapPresFromHumRatio

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Dry Air Calculations
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetDryAirEnthalpy(TDryBulb) result(DryAirEnthalpy)
    !+ Return dry-air enthalpy given dry-bulb temperature.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 28

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real              ::  DryAirEnthalpy
      !+ Dry air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI]

    if (isIP()) then
      DryAirEnthalpy = 0.240 * TDryBulb
    else
      DryAirEnthalpy = 1006 * TDryBulb
    end if
  end function GetDryAirEnthalpy

  function GetDryAirDensity(TDryBulb, Pressure) result(DryAirDensity)
    !+ Return dry-air density given dry-bulb temperature and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1
    !+ Notes:
    !+ Eqn 14 for the perfect gas relationship for dry air.
    !+ Eqn 1 for the universal gas constant.
    !+ The factor 144 in IP is for the conversion of Psi = lb in⁻² to lb ft⁻².

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  DryAirDensity
      !+ Dry air density in lb ft⁻³ [IP] or kg m⁻³ [SI]

    if (isIP()) then
      DryAirDensity = (144 * Pressure) / R_DA_IP / GetTRankineFromTFahrenheit(TDryBulb)
    else
      DryAirDensity = Pressure / R_DA_SI / GetTKelvinFromTCelsius(TDryBulb)
    end if
  end function GetDryAirDensity

  function GetDryAirVolume(TDryBulb, Pressure) result(DryAirVolume)
    !+ Return dry-air volume given dry-bulb temperature and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1
    !+ Notes:
    !+ Eqn 14 for the perfect gas relationship for dry air.
    !+ Eqn 1 for the universal gas constant.
    !+ The factor 144 in IP is for the conversion of Psi = lb in⁻² to lb ft⁻².

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  DryAirVolume
      !+ Dry air volume in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI]

    if (isIP()) then
      DryAirVolume = GetTRankineFromTFahrenheit(TDryBulb) * R_DA_IP / (144 * Pressure)
    else
      DryAirVolume = GetTKelvinFromTCelsius(TDryBulb) * R_DA_SI / Pressure
    end if
  end function GetDryAirVolume

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Saturated Air Calculations
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetSatVapPres(TDryBulb) result(SatVapPres)
    !+ Return saturation vapor pressure given dry-bulb temperature.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1  eqn 5

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real              ::  SatVapPres
      !+ Vapor pressure of saturated air in Psi [IP] or Pa [SI]
    real              ::  LnPws
      !+ Log of Vapor Pressure of saturated air (dimensionless)
    real              ::  T
      !+ Dry bulb temperature in R [IP] or K [SI]

    if (isIP()) then
      if (TDryBulb < -148.0 .or. TDryBulb > 392.0) then
        error stop "Error: dry bulb temperature must be in range [-148, 392]°F"
      end if

      T = GetTRankineFromTFahrenheit(TDryBulb)

      if (TDryBulb <= 32.) then
        LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T**2    &
                + 3.5575832E-10 * T**3 - 9.0344688E-14 * T**4 + 4.1635019 * log(T))
      else
        LnPws = -1.0440397E+04 / T - 1.1294650E+01 - 2.7022355E-02* T + 1.2890360E-05 * T**2  &
                - 2.4780681E-09 * T**3 + 6.5459673 * log(T)
      end if

      else
        if (TDryBulb < -100.0 .or. TDryBulb > 200.0) then
          error stop "Error: dry bulb temperature must be in range [-100, 200]°C"
        end if

        T = GetTKelvinFromTCelsius(TDryBulb)

        if (TDryBulb <= 0) then
          LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T**2    &
                  + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T)
        else
          LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T**2   &
                  - 1.4452093E-08 * T**3 + 6.5459673 * log(T)
        end if
      end if

    SatVapPres = exp(LnPws)
  end function GetSatVapPres

  function GetSatHumRatio(TDryBulb, Pressure) result(SatHumRatio)
    !+ Return humidity ratio of saturated air given dry-bulb temperature and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 36, solved for W

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  SatHumRatio
      !+ Humidity ratio of saturated air in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  SatVaporPres
      !+ Vapor pressure of saturated air in in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]

    SatVaporPres  = GetSatVapPres(TDryBulb)
    SatHumRatio   = 0.621945 * SatVaporPres / (Pressure-SatVaporPres)
  end function GetSatHumRatio

  function GetSatAirEnthalpy(TDryBulb, Pressure) result(SatAirEnthalpy)
    !+ Return saturated air enthalpy given dry-bulb temperature and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  SatAirEnthalpy
      !+ Saturated air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI]

    SatAirEnthalpy = GetMoistAirEnthalpy(TDryBulb, GetSatHumRatio(TDryBulb, Pressure))
  end function GetSatAirEnthalpy

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Moist Air Calculations
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetVaporPressureDeficit(TDryBulb, HumRatio, Pressure) result(VaporPressureDeficit)
    !+ Return Vapor pressure deficit given dry-bulb temperature, humidity ratio, and pressure.
    !+ Reference:
    !+ Oke (1987) eqn 2.13a

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  VaporPressureDeficit
      !+ Vapor pressure deficit in Psi [IP] or Pa [SI]
    real              ::  RelHum
      !+ Relative humidity in range [0, 1]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    RelHum = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure)
    VaporPressureDeficit = GetSatVapPres(TDryBulb) * (1.0 - RelHum)
  end function GetVaporPressureDeficit

  function GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) result(DegreeOfSaturation)
    !+ Return the degree of saturation (i.e humidity ratio of the air / humidity ratio of the air at saturation
    !+ at the same temperature and pressure) given dry-bulb temperature, humidity ratio, and atmospheric pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2009) ch. 1 eqn 12
    !+ Notes:
    !+ This definition is absent from the 2017 Handbook. Using 2009 version instead.

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  DegreeOfSaturation
      !+ Degree of saturation in arbitrary unit

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    DegreeOfSaturation = HumRatio/GetSatHumRatio(TDryBulb, Pressure)
  end function GetDegreeOfSaturation

  function GetMoistAirEnthalpy(TDryBulb, HumRatio) result(MoistAirEnthalpy)
    !+ Return moist air enthalpy given dry-bulb temperature and humidity ratio.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 30

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real              ::  MoistAirEnthalpy
      !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    if (isIP()) then
        MoistAirEnthalpy = 0.240 * TDryBulb + HumRatio * (1061.0 + 0.444 * TDryBulb)
    else
        MoistAirEnthalpy = (1.006 * TDryBulb + HumRatio * (2501.0 + 1.86 * TDryBulb)) * 1000.0
    end if
  end function GetMoistAirEnthalpy

  function GetMoistAirVolume(TDryBulb, HumRatio, Pressure) result(MoistAirVolume)
    !+ Return moist air specific volume given dry-bulb temperature, humidity ratio, and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 26
    !+ Notes:
    !+ In IP units, R_DA_IP / 144 equals 0.370486 which is the coefficient appearing in eqn 26
    !+ The factor 144 is for the conversion of Psi = lb in⁻² to lb ft⁻².

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  MoistAirVolume
      !+ Specific volume of moist air in ft³ lb⁻¹ of dry air [IP] or in m³ kg⁻¹ of dry air [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    if (isIP()) then
        MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1.0 + 1.607858 * HumRatio) / (144.0 * Pressure)
    else
        MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1.0 + 1.607858 * HumRatio) / Pressure
    end if
  end function GetMoistAirVolume

  function GetMoistAirDensity(TDryBulb, HumRatio, Pressure) result(MoistAirDensity)
    !+ Return moist air density given humidity ratio, dry bulb temperature, and pressure.
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 11

    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)  ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(in)  ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real              ::  MoistAirDensity
      !+ Moist air density in lb ft⁻³ [IP] or kg m⁻³ [SI]

    if (HumRatio < 0.0) then
      error stop "Error: humidity ratio is negative"
    end if

    MoistAirDensity = (1.0 + HumRatio) / GetMoistAirVolume(TDryBulb, HumRatio, Pressure)
  end function GetMoistAirDensity

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Standard atmosphere
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function GetStandardAtmPressure(Altitude) result(StandardAtmPressure)
    !+ Return standard atmosphere barometric pressure, given the elevation (altitude).
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 3

    real, intent(in)  ::  Altitude
      !+ Altitude in ft [IP] or m [SI]
    real              ::  StandardAtmPressure
      !+ Standard atmosphere barometric pressure in Psi [IP] or Pa [SI]

    if (isIP()) then
        StandardAtmPressure = 14.696 * (1.0 - 6.8754e-06 * Altitude)**5.2559
    else
        StandardAtmPressure = 101325 * (1 - 2.25577e-05 * Altitude)**5.2559
    end if
  end function GetStandardAtmPressure

  function GetStandardAtmTemperature(Altitude) result(StandardAtmTemperature)
    !+ Return standard atmosphere temperature, given the elevation (altitude).
    !+ Reference:
    !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 4

    real, intent(in)  ::  Altitude
      !+ Altitude in ft [IP] or m [SI]
    real              ::  StandardAtmTemperature
      !+ Standard atmosphere dry-bulb temperature in °F [IP] or °C [SI]

    if (isIP()) then
        StandardAtmTemperature = 59.0 - 0.00356620 * Altitude
    else
        StandardAtmTemperature = 15.0 - 0.0065 * Altitude
    end if
  end function GetStandardAtmTemperature

  function GetSeaLevelPressure(StnPressure, Altitude, TDryBulb) result(SeaLevelPressure)
    !+ Return sea level pressure given dry-bulb temperature, altitude above sea level and pressure.
    !+ Reference:
    !+ Hess SL, Introduction to theoretical meteorology, Holt Rinehart and Winston, NY 1959,
    !+ ch. 6.5; Stull RB, Meteorology for scientists and engineers, 2nd edition,
    !+ Brooks/Cole 2000, ch. 1.
    !+ Notes:
    !+ The standard procedure for the US is to use for TDryBulb the average
    !+ of the current station temperature and the station temperature from 12 hours ago.

    real, intent(in)  ::  StnPressure
      !+ Observed station pressure in Psi [IP] or Pa [SI]
    real, intent(in)  ::  Altitude
      !+ Altitude in ft [IP] or m [SI]
    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real              ::  SeaLevelPressure
      !+ Sea level barometric pressure in Psi [IP] or Pa [SI]
    real              ::  TColumn
      !+ Average temperature in column of air in R [IP] or K [SI]
    real              ::  H
      !+ scale height (dimensionless)

    if (isIP()) then
      ! Calculate average temperature in column of air, assuming a lapse rate
      ! of 3.6 °F/1000ft
      TColumn = TDryBulb + 0.0036 * Altitude / 2.0

      ! Determine the scale height
      H = 53.351 * GetTRankineFromTFahrenheit(TColumn)
    else
      ! Calculate average temperature in column of air, assuming a lapse rate
      ! of 6.5 °C/km
      TColumn = TDryBulb + 0.0065 * Altitude / 2.0

      ! Determine the scale height
      H = 287.055 * GetTKelvinFromTCelsius(TColumn) / 9.807
    end if

    ! Calculate the sea level pressure
    SeaLevelPressure = StnPressure * exp(Altitude / H)
  end function GetSeaLevelPressure

  function GetStationPressure(SeaLevelPressure, Altitude, TDryBulb) result(StationPressure)
    !+ Return station pressure from sea level pressure.
    !+ Reference:
    !+ See 'GetSeaLevelPressure'
    !+ Notes:
    !+ This function is just the inverse of 'GetSeaLevelPressure'.

    real, intent(in)  ::  SeaLevelPressure
      !+ Sea level barometric pressure in Psi [IP] or Pa [SI]
    real, intent(in)  ::  Altitude
      !+ Altitude in ft [IP] or m [SI]
    real, intent(in)  ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real              ::  StationPressure
      !+ Station pressure in Psi [IP] or Pa [SI]

    StationPressure = SeaLevelPressure / GetSeaLevelPressure(1.0, Altitude, TDryBulb)
  end function GetStationPressure

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Functions to set all psychrometric values
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine CalcPsychrometricsFromTWetBulb(TDryBulb,           &
                                            TWetBulb,           &
                                            Pressure,           &
                                            HumRatio,           &
                                            TDewPoint,          &
                                            RelHum,             &
                                            VapPres,            &
                                            MoistAirEnthalpy,   &
                                            MoistAirVolume,     &
                                            DegreeOfSaturation)

    !+ Utility function to calculate humidity ratio, dew-point temperature, relative humidity,
    !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given
    !+ dry-bulb temperature, wet-bulb temperature, and pressure.

    real, intent(in)    ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real, intent(out)   ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(out)   ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real, intent(out)   ::  RelHum
      !+ Relative humidity in range [0, 1]
    real, intent(out)   ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real, intent(out)   ::  MoistAirEnthalpy
      !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI]
    real, intent(out)   ::  MoistAirVolume
      !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI]
    real, intent(out)   ::  DegreeOfSaturation
      !+ Degree of saturation [unitless]

    HumRatio            = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure)
    TDewPoint           = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure)
    RelHum              = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure)
    VapPres             = GetVapPresFromHumRatio(HumRatio, Pressure)
    MoistAirEnthalpy    = GetMoistAirEnthalpy(TDryBulb, HumRatio)
    MoistAirVolume      = GetMoistAirVolume(TDryBulb, HumRatio, Pressure)
    DegreeOfSaturation  = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure)
  end subroutine CalcPsychrometricsFromTWetBulb

  subroutine CalcPsychrometricsFromTDewPoint(TDryBulb,           &
                                             TDewPoint,          &
                                             Pressure,           &
                                             HumRatio,           &
                                             TWetBulb,           &
                                             RelHum,             &
                                             VapPres,            &
                                             MoistAirEnthalpy,   &
                                             MoistAirVolume,     &
                                             DegreeOfSaturation)

    !+ Utility function to calculate humidity ratio, wet-bulb temperature, relative humidity,
    !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given
    !+ dry-bulb temperature, dew-point temperature, and pressure.

    real, intent(in)    ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real, intent(out)   ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(out)   ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(out)   ::  RelHum
      !+ Relative humidity in range [0, 1]
    real, intent(out)   ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real, intent(out)   ::  MoistAirEnthalpy
      !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI]
    real, intent(out)   ::  MoistAirVolume
      !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI]
    real, intent(out)   ::  DegreeOfSaturation
      !+ Degree of saturation [unitless]

    HumRatio            = GetHumRatioFromTDewPoint(TDewPoint, Pressure)
    TWetBulb            = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure)
    RelHum              = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure)
    VapPres             = GetVapPresFromHumRatio(HumRatio, Pressure)
    MoistAirEnthalpy    = GetMoistAirEnthalpy(TDryBulb, HumRatio)
    MoistAirVolume      = GetMoistAirVolume(TDryBulb, HumRatio, Pressure)
    DegreeOfSaturation  = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure)
  end subroutine CalcPsychrometricsFromTDewPoint

  subroutine CalcPsychrometricsFromRelHum(TDryBulb,           &
                                          RelHum,             &
                                          Pressure,           &
                                          HumRatio,           &
                                          TWetBulb,           &
                                          TDewPoint,          &
                                          VapPres,            &
                                          MoistAirEnthalpy,   &
                                          MoistAirVolume,     &
                                          DegreeOfSaturation)

    !+ Utility function to calculate humidity ratio, wet-bulb temperature, dew-point temperature,
    !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given
    !+ dry-bulb temperature, relative humidity and pressure.

    real, intent(in)    ::  TDryBulb
      !+ Dry-bulb temperature in °F [IP] or °C [SI]
    real, intent(in)    ::  RelHum
      !+ Relative humidity in range [0, 1]
    real, intent(in)    ::  Pressure
      !+ Atmospheric pressure in Psi [IP] or Pa [SI]
    real, intent(out)   ::  HumRatio
      !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI]
    real, intent(out)   ::  TWetBulb
      !+ Wet-bulb temperature in °F [IP] or °C [SI]
    real, intent(out)   ::  TDewPoint
      !+ Dew-point temperature in °F [IP] or °C [SI]
    real, intent(out)   ::  VapPres
      !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI]
    real, intent(out)   ::  MoistAirEnthalpy
      !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI]
    real, intent(out)   ::  MoistAirVolume
      !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI]
    real, intent(out)   ::  DegreeOfSaturation
      !+ Degree of saturation [unitless]

    HumRatio            = GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure)
    TWetBulb            = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure)
    TDewPoint           = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure)
    VapPres             = GetVapPresFromHumRatio(HumRatio, Pressure)
    MoistAirEnthalpy    = GetMoistAirEnthalpy(TDryBulb, HumRatio)
    MoistAirVolume      = GetMoistAirVolume(TDryBulb, HumRatio, Pressure)
    DegreeOfSaturation  = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure)
  end subroutine CalcPsychrometricsFromRelHum

end module psychrolib

max_stars_repo_name: semiamaurotic/psychrolib

max_issues_repo_name: semiamaurotic/psychrolib

max_forks_repo_name: semiamaurotic/psychrolib


! { dg-do compile }
! PR fortran/67525
! Code contributed by Gerhard Steinmetz
!
real function f(x)
   select type (x)         ! { dg-error "shall be polymorphic" }
   end select
end function f

real function g(x)
   select type (x=>null()) ! { dg-error "shall be polymorphic" }
   end select
end function g

subroutine a(x)
   select type (x)         ! { dg-error "shall be polymorphic" }
   end select
end subroutine a


text:
C--------------------------------------------------------------------------
        subroutine ftmrec(ounit,nkey,record,status)

C       modify the nth keyword in the CHU, by replacing it with the
C       input 80 character string.
C
C       ounit   i  fortran output unit number
C       nkey    i  sequence number (starting with 1) of the keyword to read
C       record  c  80-character string to replace the record with
C       OUTPUT PARAMETERS:
C       status  i  output error status (0 = ok)
C
C       written by Wm Pence, HEASARC/GSFC, June 1991

        integer ounit,nkey,status
        character*(*) record
        character rec*80

C       find the old keyword; just use REC as a temporary variable
        call ftgrec(ounit,nkey,rec,status)

        rec=record
C       overwrite the keyword with the new record
        call ftmodr(ounit,rec,status)
        end

ext: f

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
        function rn(idummy)
        real*8 rn,ran
        common/seed/ij,kl
        save init
        data init /1/

        if (init.eq.1) then
        init=0
        if (ij .eq. 0) then
        ij=1802
        kl=9373
        endif
        call rmarin(ij,kl)
        end if
*
  10    call ranmar(ran)
        if (ran.lt.1d-16) goto 10
        rn=ran
*
        end
*
      subroutine ranmar(rvec)
*     -----------------
* universal random number generator proposed by marsaglia and zaman
* in report fsu-scri-87-50
* in this version rvec is a double precision variable.
      implicit real*8(a-h,o-z)
      common/ raset1 / ranu(97),ranc,rancd,rancm
      common/ raset2 / iranmr,jranmr
      save /raset1/,/raset2/
      uni = ranu(iranmr) - ranu(jranmr)
      if(uni .lt. 0d0) uni = uni + 1d0
      ranu(iranmr) = uni
      iranmr = iranmr - 1
      jranmr = jranmr - 1
      if(iranmr .eq. 0) iranmr = 97
      if(jranmr .eq. 0) jranmr = 97
      ranc = ranc - rancd
      if(ranc .lt. 0d0) ranc = ranc + rancm
      uni = uni - ranc
      if(uni .lt. 0d0) uni = uni + 1d0
      rvec = uni
      end

      subroutine rmarin(ij,kl)
*     -----------------
* initializing routine for ranmar, must be called before generating
* any pseudorandom numbers with ranmar. the input values should be in
* the ranges 0<=ij<=31328 ; 0<=kl<=30081
      implicit real*8(a-h,o-z)
      common/ raset1 / ranu(97),ranc,rancd,rancm
      common/ raset2 / iranmr,jranmr
      save /raset1/,/raset2/
* this shows correspondence between the simplified input seeds ij, kl
* and the original marsaglia-zaman seeds i,j,k,l.
* to get the standard values in the marsaglia-zaman paper (i=12,j=34
* k=56,l=78) put ij=1802, kl=9373
      i = mod( ij/177 , 177 ) + 2
      j = mod( ij     , 177 ) + 2
      k = mod( kl/169 , 178 ) + 1
      l = mod( kl     , 169 )
      do 300 ii = 1 , 97
        s =  0d0
        t = .5d0
        do 200 jj = 1 , 24
          m = mod( mod(i*j,179)*k , 179 )
          i = j
          j = k
          k = m
          l = mod( 53*l+1 , 169 )
          if(mod(l*m,64) .ge. 32) s = s + t
          t = .5d0*t
  200   continue
        ranu(ii) = s
  300 continue
      ranc  =   362436d0 / 16777216d0
      rancd =  7654321d0 / 16777216d0
      rancm = 16777213d0 / 16777216d0
      iranmr = 97
      jranmr = 33
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
/* Examples for testing */

true;
if false then true else false;

x/;
x;

x = true;
x;
if x then false else x;

lambda x. x;
(lambda x. x) (lambda x. x x);
lambda _. true;

{x=lambda x.x, y=(lambda x.x)(lambda x.x)};
{x=lambda x.x, y=(lambda x.x)(lambda x.x)}.x;

"hello";

timesfloat (timesfloat 2.0 3.0) (timesfloat 4.0 5.0);

0;
succ (pred 0);
iszero (pred (succ (succ 0)));

let x=true in x;
let _=true in false;

max_stars_repo_name: erszcz/tapl-erlang

max_issues_repo_name: erszcz/tapl-erlang

max_forks_repo_name: erszcz/tapl-erlang


! chk_location.f90 --
!     Check the Fortran 2008 extensions to the localisation functions
!
program chk_location
    implicit none

    real, dimension(5) :: array = (/ 1.0, 5.0, 3.0, 1.0, 5.0 /)

    write( *, '(a)' ) 'Location of maximum (using BACK=.true.): ', maxloc(array, back= .true.)
    write( *, '(a)' ) 'Location of minimum (using BACK=.true.): ', minloc(array, back= .true.)

    write( *, '(a)' ) 'Location of 3.0:                         ', findloc(array, 3.0)
end program chk_location


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:32 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_ADD_SOURCES_ACOUSTIC_ADJOINT__genmod
          INTERFACE
            SUBROUTINE COMPUTE_ADD_SOURCES_ACOUSTIC_ADJOINT
            END SUBROUTINE COMPUTE_ADD_SOURCES_ACOUSTIC_ADJOINT
          END INTERFACE
        END MODULE COMPUTE_ADD_SOURCES_ACOUSTIC_ADJOINT__genmod


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE FFONOF (IFLD, INTYP, CFIELD, ISON, *)
C=======================================================================

C   --*** FFONOF *** (FFLIB) Parse free-field ON/OFF
C   --   Written by Amy Gilkey - revised 02/24/86
C   --
C   --FFONOF parses an on/off option from an input field.  No field is
C   --assumed 'ON'.
C   --
C   --Parameters:
C   --   IFLD - IN/OUT - the index of the current field number, incremented
C   --   INTYP - IN - the input type from the free-field reader
C   --   CFIELD - IN - the input option string
C   --   ISON - OUT - true iff the option is ON, set only if no error
C   --   * - return statement if the field is invalid; message is printed

      INTEGER IFLD
      INTEGER INTYP(*)
      CHARACTER*(*) CFIELD(*)
      LOGICAL ISON

      CHARACTER*4 OPT

      IF (INTYP(IFLD) .EQ. 0) THEN
         OPT = CFIELD(IFLD)
      ELSE IF (INTYP(IFLD) .LE. -1) THEN
         OPT = 'ON'
      ELSE
         OPT = ' '
      END IF
      IF ((OPT(:2) .NE. 'ON') .AND. (OPT(:3) .NE. 'OFF')) THEN
         CALL PRTERR ('CMDERR', 'Expected "ON" or "OFF"')
         GOTO 100
      END IF

      ISON = (OPT(:2) .EQ. 'ON')

      IF (INTYP(IFLD) .GE. -1) IFLD = IFLD + 1
      RETURN

  100 CONTINUE
      IF (INTYP(IFLD) .GE. -1) IFLD = IFLD + 1
      RETURN 1
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
! (C) Copyright 2019-2021 NOAA/NWS/NCEP/EMC.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

module magic_geom_mod

use iso_c_binding
use kinds
use fckit_mpi_module, only: fckit_mpi_comm
use atlas_module

implicit none

private

!> Geometry data structure
type, public :: magic_geom
  integer :: npx, npy, npz
  real(kind=kind_real), allocatable, dimension(:) :: ak, bk
  real(kind=kind_real), allocatable, dimension(:) :: grid_xt, grid_yt
  real(kind=kind_real), allocatable, dimension(:) :: pfull, phalf
  real(kind=kind_real), allocatable, dimension(:,:) :: grid_lon, grid_lat
contains
  procedure :: init => geom_init
  procedure :: delete => geom_delete
  procedure :: clone => geom_clone
  procedure :: info => geom_info
end type magic_geom

! ------------------------------------------------------------------------------
contains
! ------------------------------------------------------------------------------

! ------------------------------------------------------------------------------
!> Initialize and allocate memory for geometry object
subroutine geom_init(self, c_conf)
  class(magic_geom), target, intent(inout) :: self
  type(c_ptr),               intent(in)    :: c_conf

  type(fckit_mpi_comm) :: f_comm

  f_comm = fckit_mpi_comm()

end subroutine geom_init
! ------------------------------------------------------------------------------

! ------------------------------------------------------------------------------
!> Geometry destructor
subroutine geom_delete(self)
  class(magic_geom), intent(inout) :: self

end subroutine geom_delete
! ------------------------------------------------------------------------------

! ------------------------------------------------------------------------------
!> Clone, self = other
subroutine geom_clone(self, other)
  class(magic_geom), intent(inout) :: self
  class(magic_geom), intent(inout) :: other

end subroutine geom_clone
! ------------------------------------------------------------------------------

! ------------------------------------------------------------------------------
!> Print geometry info to stdout
subroutine geom_info(self)
  class(magic_geom), intent(in) :: self

end subroutine geom_info
! ------------------------------------------------------------------------------

end module magic_geom_mod

max_stars_repo_name: JCSDA-internal/magic

max_issues_repo_name: JCSDA-internal/magic

max_forks_repo_name: JCSDA-internal/magic


text:
C    @(#)lastch.f	20.5 1/15/98
        function lastch( text )
        character text *(*)
C
C       FIND the last non_blank character in the RECORD.
C
        if ( len( text ) .eq. 0 ) then
           lastch = 0       ! null record
        else if ( ichar( text(1:1) ) .eq. 0 ) then
           lastch = 0       ! "C" style null string
        else if ( text .eq. ' ' ) then
           lastch = 1
        else
           l = len( text )
           do i = l,2,-1
           if ( text(i:i) .ne. ' ' ) go to 110
           enddo
           i = 1
110        lastch = i
        endif
        return
        end

        function lastch0 ( text )
        character text *(*)
C
C       FIND the last non_blank character in the RECORD.
C
        character * 1 null

        null = char(0)
        if ( len(text) .eq. 0 ) then
           lastch0 = 0       ! null record
        else if ( text(1:1) .eq. null ) then
           lastch0 = 0       ! "C" style null string
        else if ( text .eq. ' ' ) then
           lastch0 = 1
        else
           lc = index( text, null ) - 1
           if ( lc .lt. 1 ) then
             lc = len (text)
           else if (text(lc:lc) .eq. char(10)) then
             lc = lc - 1
           endif
           do i = lc, 2, -1
             if ( text(i:i) .ne. ' ' ) go to 110
           enddo
           i = 1
110        lastch0 = i
        endif
        return
        end

        function lastchln( text )
        character *(*) text

        character * 4  ctst

        lc = len( text )
        lastchln = 0
        if ( lc .eq. 0 ) return   ! null record
        if ( text .eq. ' ' ) then
           lastchln = 1
           return
        endif
c***           null        LF           FF          CR
        ctst = char(0) // char(10) // char(12) // char(13)
        i = 1
        do while ( i .lt. lc+1 )
           if ( index( ctst, text(i:i) ) .ne. 0 ) then
              lc = i - 1
              goto 100
           endif
           i = i + 1
        enddo
  100   continue
        do while ( lc .gt. 1  .and.  text(lc:lc) .eq. ' ' )
           lc = lc - 1
        enddo
        lastchln = lc
        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
Heres a compilation of Gene Shalitlike oneline impressions of places Ive eaten at, or slightly longer lines where warranted. Most of these comments are not worthy of spreading throughout the wiki. Others are summaries of comments or edits Ive previously made on individual pages. If you agree with my general comments about restaurants, you might find this page useful.

 American
 Americana
   Applebees A few years back, I said mediocrity at its worst. The food has gotten better, so Ill revise that to consistently mediocre for most of its food, and good marks for the burgers.
   Bakers Square Decent pies.
   Cindys Greasier than Jerry Lewiss hair
   The Graduate Decent place for dancing which somehow manages to avoid getting a noise violation
   The Grill at G Street Pub Attempts and fails to deliver food classier than bar food and Shadiness Factor shadily takes cover charge right before announcing last call.
   Hot Dogger With the exception of the chicken linguica, I like their stuff, particularly the nice selection of mustards.
   Seasons Good overall, but the rabbit pasta was generous with the bunny and light on the interesting.
   Sudwerk Theyve got a western burger thats pretty tasty, but their appetizer selection was better before the ownership change.

 Sandwiches / Delis
   Noahs Bagels Okay bagels. Despite faux NY atmosphere, these bastards replaced Dr. Browns Kosher soda with more regional Thomas Kemper Soda.
   Plutos I enjoy their salads and the occasional seasonal ingredients. The stuffing is very filling if youre looking for an inexpensive and filling side dish
   Quiznos Pricey but decent for a chain sandwich shop.
   Subway Lowquality, overprocessed, and devoid of taste.
   Sunrise Restaurant Vegetarian Vietnamese sandwich is inexpensive and delicious
   Togos Their blatant overuse of avocado is inexcusable.
   Zias Delicatessen The meatball sandwich w/ extra sauce is one my my favorite lunches. Try also the pepper salame sandwich and the gelato. Yum!

 California / Fusion food
   Bistro 33 The fairly priced panini make this the best restaurant to eat at after 10 PM. The service has improved much since their first six months.
   Cafe Bernardo Overpriced and overrushed for a restaurant where you have to order at the counter.
   Fuzio Mediocre YuppieChain food.
   Mustard Seed I avoid b/c the owner gave my iron chef team a low rating and the wait is too long. Otherwise fine if someone else if paying or you really like mayo.
   Natural Food Works Rose Anns Farmers Kitchen Cafe is refreshingly unlike any other Davis restaurant. Youll enjoy your ultrafresh meal here, but they typically serve food that takes a while to prepare. As such, this is a great place to chill out, but if youre in a rush, I recommend you give them a call ahead of time its worth it!

 Asian
 Afghan
   Kabul Afghan Cuisine Meat mantoo dumplings are my favorite dumpling in Davis, and the Lamb Kebob is some of the tastiest meat Ive had in Davis.

 Chinese
   Davis Noodle City One of my favorite Davis restaurants. The homemade noodles well complement their generously portioned and tasty dishes. They are probably the Panda Express vs. Davis Noodle City fastest madetoorder Chinese place in Davis.
   Ding How Bland and uninteresting. Everyone else seems to like their stuff but Ive never ordered anything here that was all that good.
   Great Wall Generic.
   House of Chang Good Hong Kong Style pan fried noodles.
   Hunan Generic.
   Jade Garden Sizzling oyster dish is distinct and satisfying.
   Yummy Guide Cafe As a Westerner, my thoughts of this place invoke references to the wiki:Wikipedia:Uncanny_Valley Uncanny Valley.
   Noodle Express Red bean noodles are cheap and a welcome departure from standard americanchinese fare
   Orange Hut The potstickers are really good, but someone must have left Captain Hazelwood near the lunch offerings.
   Old Teahouse Gets points for authenticity, but nothing has ever warranted repeat visits. Suggestions are welcome.
   Shanghai Town Best (only ...) Dim Sum in Davis, and hands down the best Chinese in Davis. The deep fried and Americanized dishes come out overly oily, but numerous steamed dishes, noodle dishes, and the Lions Head Casserole are plentiful, wellprepared, and unique to the Davis Chinese scene. Anything dish preceded by Shanghai is well worth getting.
   Silver Dragon Good quantity lunch specials, but generic.
   Wok n Roll Really generic.
   Wok of Flame The house special panfried noodles come with a gazillion different vegetables, which isnt to my taste, but might be to yours.

 Thai
   2k Thai Food Most improved restaurant 2005.
   Sophias Thai Kitchen Fritter appetizers are really good and the Tom Ka Gai (sp?) is probably the best Thai dish in town. They get mad props for including rice into the price of their dishes. On weekend nights they get pretty busy, but the service and atmosphere of the bar is always great, so its no big whoop to wait there for a table.
   Thai Bistro Fresher feel than other Thai places, and the duck is always good.

 Japanese / Sushi
   Fuji Chef Long wait makes it unappealing, but the high variety buffet is worth it if you avoid the sushi boats and get directly from chefs. Get there early the crowds are small Ive gotten nasty fish later in the day.
   Jusco Had bad fish here once and am loathe to return. The fried calamari is good, but the sushi has far too much rice, and not enough seafood. The sushi is about the same quality youll get at a Chinese buffet.
   Osaka Sushi Decent quality, but pricey.
   Nobu Hiro Generally has the best quality sushi for a lunch buffet, but limited offerings and quality can be variable.
   Zen Toro Good atmosphere, good spam sushi.

 Korean
   Manna Ive always enjoyed my food here and should go more frequently. I enjoy the bibimbap.

 Vietnamese
   Hoa Viet Their garlic chicken trumps that of all chinese restaurants.
   Sunrise Restaurant Vegetarian Vietnamese sandwich is cheap and delicious.

 Indian
   Kathmandu Kitchen The vegetarian momo plate is a really good TIbetan(?) dish. The buffet has gone downhill in recent times.
   New Delhi Chaat Cafe I frequent the Panir wrap and that Bhel Puri, but make sure to tell them to spice them up lest they be too bland. The wraps by themselves always leave me hungry, so its best to get something else too.
   Rajas I enjoy the rice pudding which is the strongest element of the buffet. They servers are really nice.

 European
 British
   London Fishn Chips Decent and plentiful fish and chips, really good fried shrimp.

 Czech
   Little Prague Venison with pistachio nuts is yummy as are their breakfast dumplings. The potato pancakes are quite good, but only worth the money during happy hour. Desserts are top notch. Enough of my friends have complained about the owner and his wife that I no longer go here.

 French
   Crepeville Winner of a call to the Dept. of Health after an argument with the unhelpful owner. In 2003, I got a soup and sandwich to go. I finished the soup soon after, but didnt get to the tuna sandwich till the next day. Upon eating it, i discovered something plasticky. I figured it was packaging, but it turned out to be a bandaid minus the gauze. Sick to my stomach, I stormed into Crepeville and demanded to see the owner. I explained my discovery to him, telling him that I wanted to know from whom it came, if they had any diseases, etc. His reaction was completely unsatisfactory. He only spoke in reaction to my specific complaints, and he showed no effort to take control of the situation or to calm me down. He made me feel as if I were somehow inconveniencing him just by being there. His damage control tactics were limited to telling me of their stellar health code record, how he was sorry for my experience but that these things happen, and that I should call the dept. of health if I was concerned. In response to him unhelpfully saying, what do you want me to do about it?, I gave him my phone number and instructed him to contact me when they found out the origin of the BandAid.  I then went outside, called the Dept. of Health, and filed a complaint. After the call, I realized something the manager didnt even have the decency to offer me my money back! I went in and said, arent you forgetting something? He looked at me blankly, and then said of course after I requested a refund. BR
  Im not resentful about the Bandaid itself and thats why Ive not explicated this on the Crepeville page. It was a horrible experience that no one should have to go through, but restaurants do have accidents and Im sympathetic to that. I spoke with the Dept. of Health recently, and it appears that theyve not had any complaints since mine, convincing me Crepeville is a perfectly sanitary restaurant. However, I do harbor resentment with the management and refuse to eat there. The customer might not always be right, but when they find a medical waste product in their food, the least the manager can do is immediately offer them their money back. BTW, Crepeville never called me to tell me from whom the Bandaid came.

 Italian  Primarily Pizza
   Mountain Mikes Pizza Great portions for the price, but nothing special tastewise
   Papa Murphys Um, since when have the Irish been renowned for their pizza making skills?
   Papa Johns No pizzeria should depend on a side of garlic butter for flavor.
      Dont you mean artificiallyflavored, partiallyhydrogenated soy oil?  Users/ArlenAbraham aa
        Ya, that stuff is the worst!  Users/StevenDaubert
   Pizza Guys Okay, with good value when they have specials.
   Steves Place Pizza Blech. The breadsticks are far better than the pizza and their subs are subpar there should not be cold veggies on a hot sub.
   Symposium Spicey Greek Pizza is best pizza in town, and the flaming cheese is hair singeingly good..
   Woodstocks Pizza Dessert pizza impressed me, but financially exploits peoples misunderstanding of Pizza: FSI Calculations geometry
   Village Bakery Really variable, but still consistently the 2nd best pizza in Davis, notably the Eggplant and Kalamata pizza

 Italian  Traditional
   Caffe Italia They cant make any dish right.
   Pasta Pasta? Really good appetizers, but their use of canned pesto and inability to make Pasta Arrabiata is inexcusable
   Strings Italian Cafe Generic and exercises the Law of the Diminishing Plate Size on their allyoucaneat night. I do enjoy the bread though.

 Spanish
   Aioli Bodega Espanola Really friendly service, and they dont seem to mind when you ask for bread 10,000 times to sop up the sauces.

 Taqueria
   El Mariachi Best salsa in town, notably the pumpkin one.
   Taqueria Guadalajara Best Burrito in town. The Wet Chile Relleno Super Burrito is amazing.

 Latininfused
   Chipotle Good by NY standards but equivalent to 711 by Bay Area standards. The only Taqueria in town that doesnt give complementary chips.

 Mediterranean
 MiddleEastern
   Ali Baba Fessenjen is very sweet and awesome. Bacon Mushroom burger is best burger in town.
   Cafe Mediterranee Good muhumarah and the potato side salad is surprisingly worth getting.
   Sams Mediterranean I was reluctant to get the chicken sandwich, but I swear by it now.
   Pita Pit Too sterile for 2AM hunger pangs. Nearly all of their items would be better on a bread other than pita.

 Fast Food
 Burgers
   Burger King In specifically the Davis Burger King, Pallid and flat Whoppers are better dubbed Whimpers
   InNOut Extra syrup shake, animal style double double, and twice fried fries are the perfect fast food meal. And they dont put mustard on burgers, thank goodness!
   Jack in the Box This Crap in the Box is incredible slow.
   McDonalds Atrocious soda and food, and best used as an American Embassy to be used only for shelter, utensils, or a place to take a McShit.
   Redrum Burger The onion ring batter on the calamari is worth trying once, but prepare to wait a long time or call ahead. The ZOOM is really good, but if you sneak it into the movie theater, its presence will be known.

 Chicken
   KFC Kentucky Fried Cramps is what youll get from eating here. I like KFC, but its not worth the intestinal distress.

 Mexican
   Del Taco Makes Taco Bell look authentic.
   Taco Bell I typically get the Nacho Cheese Chalupa.
   Jimboys Tacos Makes Del Taco look authentic. Cheese Quesadilla gave me the first stomach ache Id had since I was 7.

 Chinese
   Panda Express Pandas to Americans sweet tooths and their inability to identify quality Chinese food. I really panda Davis location in the Panda Express vs. Davis Noodle City experiment.

 Dessert
   Baskin Robbins Dollar scoop night makes this place awfully tempting
   Ben & Jerrys Generous with portions and fun to watch them toss the scoop onto the cone.
   The Candy House of Davis Scrumptious, generous with the free samples, and costs only $1 a truffle. Owing to location, its the most underappreciated business in Davis.
   Ciocolat Salmon Croque Monsieur perplexed me at first, but wound up being very good.
   Konditorei Cassia dessert is heavenly. The hours are odd on weekends, but I always make an effort to bring out of towners here.

 Coffee / Cafes
   The Posh Bagel Salt bagel with Lox Trim is affordable and passable for a California bagel place.
see also Doughnuts

 Outside of Davis
 Oakland
   Zacharys Pizza Best Chicago style pizza

 Sacramento
   Buca di Beppo Dont be fooled by the huge family style servings this place should serve Rolaids as after dinner mints.
   Kru The Poki, served with an assortment of three flavored salts, was unbelievably good. Unbelievably good.
   New Canton Best Dim Sum Ive had in California.
   Rice Bowl Good selection of Dim Sum, notably the Fried Meat Vegetable Bun and CharShiu Bow. If you like snails, be warned that crunchy baby snail shells cling to snail meat.
   SF Market Good Har Gow dumplings and Sesame Balls, but the duck was almost indistinguishable from chicken.
   Streets of London Pub Delicious meat pie. Try it on Sunday trivia night.
   Tapa the World Good overall, and I seem to remember the mushrooms being fantastic.
   Thai Basil Better Thai than most places in Davis.

 Winters
   Buckhorn Steak and Roadhouse Best steak house in the area. Try the cajun steak bites appetizer if you want to try the steak, but would prefer a different main course.
   Putah Creek Cafe My favorite breakfast place in the area, Ive never gotten something I didnt enjoy here.

20050910 11:34:09 nbsp I, for one, demand that Craig dresses up like Gene Shalit when writing these reviews.  You can buy surplus mock facial hair at Stagecraft Studios in Berkeley. Users/LeightonHinkley

20051107 00:39:21 nbsp thx Mr. CB Users/MichaelGiardina

20051108 14:42:23 nbsp Given your Chipotle review, its obvious what your Posh Bagel review should be... Users/AlexanderWoo

20051119 14:49:29 nbsp You should check out Steady Eddys in Winters.  Good mochas, good sandwiches, great place to spend time. Users/GrahamFreeman

20051207 22:37:16 nbsp I would like to tag along to dim sum sometime.  Have you been to the Scottish Meat Pie Shop in Dixon? Users/JudithTruman

20060104 22:18:18 nbsp I think Taco Bell makes Del Taco look authentic.  I mean, sure D.T. has fries, but their mainstays are not a product of food science endeavors like taco bell is. Users/JaimeRaba

20060307 19:45:44 nbsp Interesting story on Crepeville. I dont eat there because I consider it Crapville. Crepes are the easiest thing to make at home anyway, so I save a lot of money. Ever gone to Sogas and tried their salads? Theyre a pile of iceberg lettuce with something worse than Hidden Valley Ranch... for $10. Users/KarlMogel

20060308 13:53:09 nbsp In spite of the bland atmosphere, Ding How deserves another drivebyCraiging.  Order off their vegetarian menu (its a different color), and marvel at the delicious vegetarian beef and green beans (thinlysliced, savory seitan).  Lightyears better than the regular beef and peas.  Also, it might be worth adding Farmers Kitchen Table to your Natural Food Works entry, because that is the name of Rose Anns restaurant there.  xox Users/ElHarvey

20061013 17:00:28 nbsp Dude you have to understand this inst The City or New York.  Davis might not have the best food around, but it isnt too bad.  Also on Del Taco, the only reason it exists is to make $.79 Chicken Soft Tacos for kids who are either drunk or stoned late at night.  Quality is not to be expected. Users/RussellAnderson
  Yeah, I kind of wonder why anyone would go to Del Taco sober or during the day Taco Bell is far tastier for cheap greasemeat. Generally, Id say that the food in Davis has consistently improved and become more diverse since I moved here in 2001. Notably, Asian cuisine has really flourished with the addition of an Afghan, 2 Vietnamese sandwich, an Indian, a Filipino, a Taiwanese noodle, a Shanghai style chinese, a Korean, and multiple Thai restaurants. In terms of quality, one need not complain about restaurants by holding them to a big city yardstick. When I go to a restaurant, I judge it based on whether I can cook and easily assembly a superior dish at home. Im neither Italian by heritage nor that great a cook, but I can easily make better Italian food than any place in town. If I, armed with only a cookbook and Safeway ingredients, can do better than people who cook for a living, theres definitely something wrong. That rule is why I praise AliBaba, Farmers Kitchen Cafe, and Shanghai Town and insult Great Wall, Pasta?, Papa Murphys and Crepeville. Well, I insult Crepeville because I found a bandaid in my food and the manager reacted like a douchebag, but you get the point. Users/CraigBrozinsky

20070303 20:03:18 nbsp You need a one liner for the educated eatery! wanna go one of these days? Users/StevenDaubert

20070628 14:30:57 nbsp could you elaborate on Yummy guide as being related to uncanny valley.  Is it that as the food gets more authentic we enjoy it until it reaches an almost authentic status, and then we can only see the differences and become repulsed?  Im a little confused by the comparison, but it sounds interesting. Users/MattHh
  Thats kind of right. If they were completely European, that would be good. If it were completely Chinese, that would be good. If it were 50/50 fusion or taking the best of both cultures, that would be good. But when a chef trained in Chinese cuisine tries to make an American dish, their Chinese cuisine skills are going to somehow, even unintentionally, influence the dish. Their Mushroom Noodle Casserole might taste great to someone unfamiliar with American cuisine, but its not going to taste like mamas home cookin. Also, some narrow range between 50/50 and 100% culturepure just aint good because the culturepure is a superior dish. Thats my impression of their pork chop, tomato dishes, gravy dishes, and spaghetti. My other problems with their food is that their more American styled sauces taste artificial, probably because theyre factory made.  Users/CraigBrozinsky

20070722 12:51:18 nbsp i was recently reading the panda experiment, and i noticed the picture for DNCs orange chicken.  It looks fried?  ive ordered it there, and it has never been fried.  While panda is way worse than noodle, i actually like their orange chicken a little better.  I use to love noodle city but lately i think shanhai town, noodle express, and home town do a better job.  I still come back to noodle city for their great noodle salad though. Users/MattHh
  Funny you should mention that Im now more of a Shanghai Town regular, and I havent eaten DNCs orange chicken since the experiment, but got it a few days back. It was really good, but it seemed different from what Id remembered. Too bad I ate it too quickly to notice if it were fried or not :)  Users/CraigBrozinsky

20070809 08:15:18 nbsp
Why hasnt craig eaten at Delta of Venus?  There IS good food there. I like their Samosa. Users/JoshuaLee
  I only leave reviews if I can remember something specific about my dining experience. Its been several years since Ive eaten there, but I should definitely return so I can write something. At present all I can say is that I had a good chickenny thing there once, and that Im allergic to their curried tofu sandwich. Users/CraigBrozinsky

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module rabs_mcp
!
!***** February 2012, wide *****
!-----------------------------------------------------------------------------------------------------------------------
! This module provides an 'interface' to the routines of the MCP/MCBP components and the NJGRAH library as they are
! implemented in GRASP92. These parts of the program are considered to be 'obsolet' for the further development of the
! RABS package and will be replaced in the future by a more efficient treatment of 'angular integration'. At the present,
! however, these routines for calculating angular coefficients are required to generate continuum orbital in the mean-
! field of some given symmetry-adapted atomic state or for calculating Auger matrix elements, for instance.
! In this module, a CSF basis for a pair of CSF (two_csf) is defined into which all required quantum numbers are defined.
! Only the angular coefficients for one pair of CSF can be calculated by a call of the MCP routines. The quantum numbers
! of these CSF are 'used' by the (independent, non-module) procedures in rabs_mcp_adaptation to set the proper COMMON
! arrays of the MCP package. This somehow sophisticated looking procedure ensures however, that COMMON arrays and other
! features (which are obsolet in Fortran 90) need only to appear in a small number of procedures in rabs_mcp_adaptation.
!-----------------------------------------------------------------------------------------------------------------------
   !
   use rabs_angular
   use rabs_constant
   use rabs_csl
   implicit none
   !
   public  :: mcbp_generate_coefficients
                 ! Generates all (non-zero) MCBP coefficients for a given 'configuration scheme'.
   public  :: mcp_generate_coefficients
                 ! Generates all (non-zero) MCP coefficients for a given 'configuration scheme'.
   public  :: mcp_initialize_two_csf
                 ! Initializes an intermediate CSF scheme for two CSF in order to tranfer proper quantum numbers to the
                 ! MCP and NJGRAF components.
   public  :: mct_generate_coefficients
                 ! Generates all (non-zero) MCT coefficients for a given 'configuration scheme'.
   !
   ! Define data structure to allow for an adaptation of the GRASP92 procedures of the MCP and NJGRAF components.
   integer, save         :: csf_a_No, csf_b_No
   logical, save         :: mcp_common_need_initialization = .true.
   type(csf_basis), save :: two_csf
   !
   ! Define a structure to handle with MCP or MCBP coefficients
   type :: mcp_coefficient
      integer           :: r, s
      integer(kind=i2b) :: a, b, c, d, mu, nu
      real(kind=dp)     :: V
   end type mcp_coefficient
   !
   ! Define a structure to handle with MCT coefficients
   type, bind(c) :: mct_coefficient
      integer           :: r, s
      integer(kind=i2b) :: a, b, nu
      real(kind=dp)     :: T
   end type mct_coefficient
   !
   integer :: number_of_mcp_coefficients,  number_of_mcp_coefficients_max
   integer :: number_of_mcbp_coefficients, number_of_mcbp_coefficients_max
   !
   integer, bind(c) :: number_of_mct_coefficients
   !
   integer, parameter :: number_of_mct_coefficients_max = 1000
   !
   type (mcp_coefficient), dimension(:), allocatable :: mcp_list, mcp_list_save
   type (mcp_coefficient), dimension(:), allocatable :: mcbp_list
   !
   type (mct_coefficient), dimension(number_of_mct_coefficients_max), bind(c) :: mct_list
   !
   ! Define debug variable to support a few GRASP92 debug features
   integer, save :: debug_ibug1, debug_ibug2, debug_ibug3, debug_ibug4, debug_ibug5, debug_ibug6
   !
   logical :: debug_mcp_coefficients = .false.
   logical :: debug_mct_coefficients = .false.
   !
contains
   !
   subroutine cfp_1
   !--------------------------------------------------------------------------------------------------------------------
   ! Calculates
   !
   ! Calls:
   !--------------------------------------------------------------------------------------------------------------------
      !
   end subroutine cfp_1
   !
   !
   subroutine mcbp_generate_coefficients(rl,ru,sl,su,csf_set)
   !--------------------------------------------------------------------------------------------------------------------
   ! Generates the (non-zero) MCBP coefficients for all pairs of CSF in the range   r = rl,...,ru and  s = sl,...,su
   ! of the given 'configuration scheme'. The procedure applies an adaption of the MCBP component from  GRASP92 to the
   ! Fortran 90 standard. All coefficients for the full configuration scheme or for the selected part are stored in the
   ! array mcbp_list() and can be 'used' later from any other module of the RABS package..
   !
   ! Calls: mcp_initialize_two_csf(), breid(), breit(), rkco() [from GRASP92/MCP].
   !--------------------------------------------------------------------------------------------------------------------
      !
      integer, intent(in)         :: rl, ru, sl, su
      type(csf_basis), intent(in) :: csf_set
      !
      integer  :: incor, r, s
      external    breid, breit, rkco
      !
      ! Allocate fresh storage for the MCP coefficients
      if (allocated(mcbp_list))  deallocate( mcbp_list )
      number_of_mcbp_coefficients = 0;   number_of_mcbp_coefficients_max = 1000
      allocate( mcbp_list(1:number_of_mcbp_coefficients_max) )
      !
      if (rl < 1   .or.  ru > csf_set%nocsf   .or.  sl < 1   .or.  su > csf_set%nocsf)  then
         stop "mcbp_generate_coefficients(): program stop A."
      end if
      !
      debug_ibug1 = 0
      !
      ! Accumulate contributions from the two-electron Coulomb operator
      if (debug_mcp_coefficients) then
         write (99,*)
         write (99,*) "                                     kt"
         write (99,*) "                                    v  (abcd)"
         write (99,*) "  r  s   a    b    c    d   k  t     rs"
         debug_ibug1 = 1
      end if
      !
      do  r = rl,ru
         do  s = sl,su
            call mcp_initialize_two_csf(r,s,csf_set)
            !
            incor = 1
            call rkco (r,s,breit,breid,incor)
         end do
      end do
      !
   end subroutine mcbp_generate_coefficients
   !
   !
   subroutine mcp_generate_coefficients(rl,ru,sl,su,csf_set)
   !--------------------------------------------------------------------------------------------------------------------
   ! Generates the (non-zero) MCP coefficients for all pairs of CSF in the range   r = rl,...,ru and  s = sl,...,su  of
   ! the given 'configuration scheme'. The procedure applies an adaption of the MCP component from GRASP92 to the
   ! Fortran 90 standard. All coefficients for the full configuration scheme or for the selected part are stored in the
   ! array mcp_list() and can be 'used' later from any other module of the RABS package..
   !
   ! Calls: mcp_initialize_two_csf(), rkco() [from GRASP92/MCP].
   !--------------------------------------------------------------------------------------------------------------------
      !
      integer, intent(in)         :: rl, ru, sl, su
      type(csf_basis), intent(in) :: csf_set
      !
      integer  :: incor, r, s
      external    cor, cord, rkco
      !
      ! Allocate fresh storage for the MCP coefficients
      if (allocated(mcp_list))  deallocate( mcp_list )
      number_of_mcp_coefficients = 0;   number_of_mcp_coefficients_max = 1000
      allocate( mcp_list(1:number_of_mcp_coefficients_max) )
      !
      if (rl < 1   .or.  ru > csf_set%nocsf   .or.  sl < 1   .or.  su > csf_set%nocsf)  then
         stop "mcp_generate_coefficients(): program stop A."
      end if
      !
      debug_ibug1 = 0
      !
      ! Accumulate contributions from the two-electron Coulomb operator
      if (debug_mcp_coefficients) then
         write (99,*)
         write (99,*) "                                    k"
         write (99,*) "                                   V  (abcd)"
         write (99,*) "  r  s   a    c    b    d   k       rs"
         debug_ibug1 = 1
      end if
      !
      do  r = rl,ru
         do  s = sl,su
            call mcp_initialize_two_csf(r,s,csf_set)
            !
            incor = 1
            call rkco(r,s,cor,cord,incor)
         end do
      end do
      !
   end subroutine mcp_generate_coefficients
   !
   !
   subroutine mcp_initialize_two_csf(r,s,csf_set)
   !--------------------------------------------------------------------------------------------------------------------
   ! Initializes some storage of quantum numbers for the calculation MCP coefficients for the pair r, s of CSF.
   !
   ! Calls: mcp_initialize_grasp92().
   !--------------------------------------------------------------------------------------------------------------------
      !
      integer, intent(in)         :: r, s
      type(csf_basis), intent(in) :: csf_set
      !
      integer :: nw
      !
      call deallocate_csf_basis(two_csf)
      !
      csf_a_No = r;   csf_b_No = s
      nw = csf_set%nwshells
      two_csf%nwshells = csf_set%nwshells
      two_csf%nocsf    = 2
      allocate( two_csf%csf(1:2), two_csf%subshell(1:csf_set%nwshells) )
      two_csf%subshell(1:nw)%n     = csf_set%subshell(1:nw)%n
      two_csf%subshell(1:nw)%kappa = csf_set%subshell(1:nw)%kappa
      allocate( two_csf%csf(1)%occupation(1:csf_set%nwshells),   two_csf%csf(1)%seniority(1:csf_set%nwshells),   &
                two_csf%csf(1)%subshellJ(1:csf_set%nwshells),    two_csf%csf(1)%subshellX(1:csf_set%nwshells),   &
                two_csf%csf(2)%occupation(1:csf_set%nwshells),   two_csf%csf(2)%seniority(1:csf_set%nwshells),   &
                two_csf%csf(2)%subshellJ(1:csf_set%nwshells),    two_csf%csf(2)%subshellX(1:csf_set%nwshells) )
      two_csf%csf(1)%totalJ           = csf_set%csf(r)%totalJ
      two_csf%csf(1)%parity           = csf_set%csf(r)%parity
      two_csf%csf(1)%occupation(1:nw) = csf_set%csf(r)%occupation(1:nw)
      two_csf%csf(1)%seniority(1:nw)  = csf_set%csf(r)%seniority(1:nw)
      two_csf%csf(1)%subshellJ(1:nw)  = csf_set%csf(r)%subshellJ(1:nw)
      two_csf%csf(1)%subshellX(1:nw)  = csf_set%csf(r)%subshellX(1:nw)
      !
      two_csf%csf(2)%totalJ           = csf_set%csf(s)%totalJ
      two_csf%csf(2)%parity           = csf_set%csf(s)%parity
      two_csf%csf(2)%occupation(1:nw) = csf_set%csf(s)%occupation(1:nw)
      two_csf%csf(2)%seniority(1:nw)  = csf_set%csf(s)%seniority(1:nw)
      two_csf%csf(2)%subshellJ(1:nw)  = csf_set%csf(s)%subshellJ(1:nw)
      two_csf%csf(2)%subshellX(1:nw)  = csf_set%csf(s)%subshellX(1:nw)
      !
      call mcp_initialize_grasp92()
      !
   end subroutine mcp_initialize_two_csf
   !
   !
   subroutine mct_generate_coefficients(rl,ru,sl,su,csf_set,iopar,rank)
   !--------------------------------------------------------------------------------------------------------------------
   ! Generates the (non-zero) MCT coefficients for all pairs of CSF in the range   r = rl,...,ru and  s = sl,...,su  of
   ! the given 'configuration scheme'. The procedure applies an adaption of the MCT component from GRASP92 to the
   ! Fortran 90 standard. All coefficients for the full configuration scheme or for the selected part are stored in the
   ! array mct_list() and can be 'used' later from any other module of the RABS package..
   !
   ! Calls: mcp_initialize_two_csf(), rkco() [from GRASP92/MCP].
   !--------------------------------------------------------------------------------------------------------------------
      !
      integer, intent(in)           :: rl, ru, sl, su, iopar, rank
      type(csf_basis), intent(in)   :: csf_set
      !
      integer                       :: ia, ib, r, s, ncf, nw_grasp92, pntriq
      real(kind=dp), dimension(149) :: tshell
      !
      common/orb2/ ncf, nw_grasp92, pntriq
      !
      number_of_mct_coefficients = 0
      !
      if (rl < 1   .or.  ru > csf_set%nocsf   .or.  sl < 1   .or.  su > csf_set%nocsf)  then
         stop "mct_generate_coefficients(): program stop A."
      end if
      !
      debug_ibug1 = 0
      !
      ! Accumulate contributions from the two-electron Coulomb operator
      if (debug_mct_coefficients) then
         write (99,*)
         write (99,*) "                          k"
         write (99,*) "                         T  (ab)"
         write (99,*) "  r  s   a    b   k       rs"
         debug_ibug1 = 1
      end if
      !
      do  r = rl,ru
         do  s = sl,su
            call mcp_initialize_two_csf(r,s,csf_set)
            !
            call tnsrjj (rank,iopar,r,s,ia,ib,tshell)
            if (ia .ne. 0) then
               if (ia .eq. ib) then
                  do 2  ia = 1,nw_grasp92
                     if (abs (tshell(ia)) .gt.  eps10) then
                        number_of_mct_coefficients = number_of_mct_coefficients + 1
                        if(number_of_mct_coefficients > number_of_mct_coefficients_max) then
                           stop "mct_generate_coefficients(): program stop B."
                        end if
                        !
                        mct_list(number_of_mct_coefficients)%r  = r
                        mct_list(number_of_mct_coefficients)%s  = s
                        mct_list(number_of_mct_coefficients)%a  = ia
                        mct_list(number_of_mct_coefficients)%b  = ia
                        mct_list(number_of_mct_coefficients)%nu = rank
                        mct_list(number_of_mct_coefficients)%T  = tshell(ia)
                     endif
                2 continue
               else
                  if (abs (tshell(1)) .gt. eps10) then
                     number_of_mct_coefficients = number_of_mct_coefficients + 1
                     if(number_of_mct_coefficients > number_of_mct_coefficients_max) then
                        stop "mct_generate_coefficients(): program stop C."
                     end if
                     !
                     mct_list(number_of_mct_coefficients)%r  = r
                     mct_list(number_of_mct_coefficients)%s  = s
                     mct_list(number_of_mct_coefficients)%a  = ia
                     mct_list(number_of_mct_coefficients)%b  = ib
                     mct_list(number_of_mct_coefficients)%nu = rank
                     mct_list(number_of_mct_coefficients)%T  = tshell(1)
                  endif
               endif
            endif
         end do
      end do
      !
   end subroutine mct_generate_coefficients
   !
end module rabs_mcp


text: Hi Mr. SOS, Welcome to the Wiki!  Im Users/BrentLaabs  a real person typing this, and this is your userpage.  I was also MiniLuv in the Students for an Orwellian Society back when I went to UCD.  Id like to note that Wiki Community/Organizational accounts organizational accounts are frowned upon in this here wiki, and as an Orwellian, as I assume you are, you certainly would not want to unfollow any goodthink rules that Big Brother has given us.  Unless youre with the terrorists, that is.  If you have any questions, feel free to ask me on my page.   Users/BrentLaabs