R/esfofo.R

esfofo = function(Sal, Temp, Dep, Alt){
  #References the dpress function from this package (aqeco)
  DepP = dpress(Dep)

  #Uses ideal gas law with scale height of 8.4km
  AltP = exp(-(Alt/1000)/8.4)

  #Scale pressure to baseline of 0 at sea level and 0 depth
  Press = AltP + DepP - 1

  #Define all variables from Table 3 in Fofonoff (1985)
  A = function(t){
    AV = 999.842594 + 6.793952*10**(-2)*t - 9.095290*10**(-3)*t**2 +
      1.001685*10**(-4)*t**3 - 1.120083*10**(-6)*t**4 + 6.536332*10**(-9)*t**5
    return(AV)
  }

  B = function(t){
    BV = 8.24493*10**(-1) - 4.0899*10**(-3)*t + 7.6438*10**(-5)*t**2 -
      8.2467*10**(-7)*t**3 + 5.3875*10**(-9)*t**4
    return(BV)
  }

  C = function(t){
    CV = -5.72466*10**(-3) + 1.0227*10**(-4)*t - 1.6546*10**(-6)*t**2
    return(CV)
  }

  D <- 4.8314*10**(-4)

  E = function(t){
    EV = 19652.21 + 148.4206*t - 2.327105*t**2 + 1.360477*10**(-2)*t**3 -
      5.155288*10**(-5)*t**4
    return(EV)
  }

  FF = function(t){
    FV = 54.6746 - 0.603459*t + 1.09987*10**(-2)*t**2 - 6.1670*10**(-5)*t**3
    return(FV)
  }

  G = function(t){
    GV = 7.944*10**(-2) + 1.6483*10**(-2)*t - 5.3009*10**(-4)*t**2
    return(GV)
  }

  H = function(t){
    HV = 3.239908 + 1.43713*10**(-3)*t + 1.16092*10**(-4)*t**2 - 5.77907*10**(-7)*t**3
    return(HV)
  }

  I = function(t){
    IV = 2.2838*10**(-3) - 1.0981*10**(-5)*t - 1.6078*10**(-6)*t**2
    return(IV)
  }

  J <- 1.91075*10**(-4)

  M = function(t){
    MV = 8.50935*10**(-5) - 6.12293*10**(-6)*t + 5.2787*10**(-8)*t**2
    return(MV)
  }

  N = function(t){
    NV = -9.9348*10**(-7) + 2.0816*10**(-8)*t + 9.1697*10**(-10)*t**2
    return(NV)
  }

  #Eqn from caption on Table 3 of Fofonott (1985)
  kstp <- E(Temp) + FF(Temp)*Sal + G(Temp)*Sal**(3/2) +
    (H(Temp) + I(Temp)*Sal + J*Sal**(3/2))*Press + (M(Temp) + N(Temp)*Sal)*Press**2

  #Eqn from caption on Table 3 of Fofonott (1985)
  pst0 <- A(Temp) + B(Temp)*Sal + C(Temp)*Sal**(3/2) + D*Sal**2

  #Eqn from caption on Table 3 of Fofonott (1985)
  vst0 <- 1/pst0

  #Eqn from caption on Table 3 of Fofonott (1985)
  vstp <- vst0 * (1 - Press/kstp)

  #Reconvert from specific volume to pressure
  pstp <- 1/vstp

  return(pstp)
}
kevans27/aqeco documentation built on May 16, 2019, 4:08 a.m.