module biogas

  implicit none
  private
  
  real(8), parameter, private :: cvrt = 44.66_8  ! convert from ml/L to uM

  public :: oxygensat, pistonvel

contains

real(8) function oxygensat(temperature,salinity)
  !-----------------------------------------------------------------------------
  ! Function to calculate oxygen saturation in seawater
  ! Implemented by Henrik Andersson <han@dmi.dk>
  ! 2006-09-22
  !
  ! Code based on formula in:
  ! Garcia H. E., and L. I. Gordon. 1992._8
  ! Oxygen solubility in seawater: Better fitting equations
  !  Limnol. Oceanogr. 37:1307-1312._8
  !
  !-----------------------------------------------------------------------------
  ! Input: Temperature (C), Salinity (psu)
  ! Output: Oxygen solubility at atmospheric pressure in micromolar = mmol/m^3
  !-----------------------------------------------------------------------------
  
  !Dummy arguments
  real(8), intent(in) :: temperature, salinity

  !Local variables
  real(8) :: Ts, alpha, beta
  real(8), parameter :: A0 = 2.00856_8
  real(8), parameter :: A1 = 3.22400_8
  real(8), parameter :: A2 = 3.99063_8
  real(8), parameter :: A3 = 4.80299_8
  real(8), parameter :: A4 = 9.78188e-01_8
  real(8), parameter :: A5 = 1.71069_8
  real(8), parameter :: B0 = -6.24097e-03_8
  real(8), parameter :: B1 = -6.93498e-03_8
  real(8), parameter :: B2 = -6.90358e-03_8
  real(8), parameter :: B3 = -4.29155e-03_8
  real(8), parameter :: C0 = -3.11680e-07_8
  real(8), parameter :: C2K = 273.15_8        ! Celcius to Kelvin
  real(8), parameter :: C25toK = 25.0_8 + C2K ! 25 degrees C to K

  ! Ts is scaled temperature, a help variable
! FIXME: Is it sane to compare 298.15K to actual water temperature in C ?
  Ts = log( (C25toK - temperature) / (C2K + temperature) )

  alpha = salinity*((B0 + (B1 + (B2 + B3*Ts)*Ts)*Ts) + C0*salinity)

  beta  = A0 + (A1 + (A2 + (A3 + (A4 + A5*Ts)*Ts)*Ts)*Ts)*Ts

  oxygensat = exp(beta+alpha)*cvrt

end function oxygensat

real(8) function pistonvel(windspeed,ts)
  !-----------------------------------------------------------------------------
  ! Function to calculate the exchange velocity of oxygen between the atm. and 
  ! the ocean.
  !
  ! The formulation is taken from:
  ! Schimdt Number based on Wanninkhof, 1992, JGR 97, 7373-7382._8
  !-----------------------------------------------------------------------------

  real(8), intent(in) :: windspeed, ts

  !-Local variables
  real(8) :: alpha
  real(8), parameter :: A0 = 1953.4_8
  real(8), parameter :: A1 = 128.0_8
  real(8), parameter :: A2 = 3.9918_8
  real(8), parameter :: A3 = 0.050091_8
  real(8), parameter :: B0 = 660.0_8
  real(8), parameter :: C0 = 0.31_8
  real(8), parameter :: C1 = 0.24_8

  alpha = A0 - (A1 - (A2 - A3*ts)*ts)*ts
  alpha = sqrt(B0/alpha)
  pistonvel = C0*alpha*windspeed*windspeed*C1

end function pistonvel

end module biogas


