R/sw_gibbs.R

In marelac: Tools for Aquatic Sciences

## -----------------------------------------------------------------------------
## Seawater Specific Gibbs Free Energy and Derivatives up to Order 2
## -----------------------------------------------------------------------------

sw_gibbs <- function (S = 35, t = 25, p = P-1.013253, P = 1.013253,
                      dS = 0, dt = 0, dp = 0) {
  isnaS <- which (is.na(S))
  if (length(isnaS))
    S[isnaS] <- 0
  if (any (S<0))
    stop ("Salinity should be >= 0")
     
  Su <- 40.188617 # g/dkg
  tu <- 40        # dgC
  pu <- 1e8       # pa
  gu <- 1         # J/kg

  x2 <- S/Su
  x  <- sqrt(x2)

  y <- t/tu
  z <- p*1e5/pu

  mfac <- 1
  if (dt ==1)
    mfac <- mfac/tu

  if (dt==2)
    mfac <- mfac/tu/tu

  if (dp ==1)
    mfac<-mfac*1e-8

  if (dp ==2)
    mfac<-mfac*1e-16


  if (dS ==0 ) {  # the coefficients for pure water
    gjk <- matrix(data = c(
      0.101342743139674e3, 0.100015695367145e6,
     -0.254457654203630e4, 0.284517778446287e3,
     -0.333146754253611e2, 0.420263108803084e1,
     -0.546428511471039  , 0.590578347909402e1,
     -0.270983805184062e3, 0.776153611613101e3,
     -0.196512550881220e3, 0.289796526294175e2,
     -0.213290083518327e1, 0                  ,
     -0.123577859330390e5, 0.145503645404680e4,
     -0.756558385769359e3, 0.273479662323528e3,
     -0.555604063817218e2, 0.434420671917197e1,
      0                  , 0.736741204151612e3,
      -0.672507783145070e3, 0.499360390819152e3,
     -0.239545330654412e3, 0.488012518593872e2,
     -0.166307106208905e1, 0                  ,
     -0.148185936433658e3, 0.397968445406972e3,
     -0.301815380621876e3, 0.152196371733841e3,
     -0.263748377232802e2, 0                  ,
      0                  , 0.580259125842571e2,
     -0.194618310617595e3, 0.120520654902025e3,
     -0.552723052340152e2, 0.648190668077221e1,
      0                  , 0                  ,
     -0.189843846514172e2, 0.635113936641785e2,
     -0.222897317140459e2, 0.817060541818112e1,
      0                  , 0                  ,
      0                  , 0.305081646487967e1,
     -0.963108119393062e1, 0                  ,
      0                  , 0                  ,
      0                  , 0),
      nrow = 8, ncol = 7, byrow = TRUE)

    nr <- 8
    nc <- 7
    if (dt >= 1){
      nr  <- nr-1
      gjk <- gjk[-1,]
      for (j in 2:nr) gjk[j,] <- gjk[j,]*j
    }

    if (dt ==2) {
      nr  <- nr-1
      gjk <- gjk[-1,]
      for (j in 2:nr) gjk[j,] <- gjk[j,]*j
    }

    if (dp >= 1) {
      nc  <- nc-1
      gjk <- gjk[,-1]
      for (k in 2:nc) gjk[,k] <- gjk[,k]*k
    }
    if (dp ==2) {
      nc  <- nc-1
      gjk <- gjk[,-1]
      for (k in 2:nc) gjk[,k] <- gjk[,k]*k
    }

    Gpure <- 0
    for (j in 1:nr) {
      for ( k in 1:nc) Gpure <- Gpure + gjk[j,k]*y^(j-1)*z^(k-1)
    }
  } else Gpure <- 0   # dS != 0
  Gsea <- 0
  if ( any (S>0 ) ) { # the coefficients for seawater
    gijk <- array(dim=c(7,7,6),0)

    Carr <- matrix(data=c(
      1, 0, 0,  5812.81456626732,1, 1, 0,  851.226734946706,
      2, 0, 0,  1416.27648484197,3, 0, 0, -2432.14662381794,
      4, 0, 0,  2025.80115603697,5, 0, 0, -1091.66841042967,
      6, 0, 0,  374.601237877840,7, 0, 0, -48.5891069025409,
      2, 1, 0,  168.072408311545,3, 1, 0, -493.407510141682,
      4, 1, 0,  543.835333000098,5, 1, 0, -196.028306689776,
      6, 1, 0,  36.7571622995805,2, 2, 0,  880.031352997204,
      3, 2, 0, -43.0664675978042,4, 2, 0, -68.5572509204491,
      2, 3, 0, -225.267649263401,3, 3, 0, -10.0227370861875,
      4, 3, 0,  49.3667694856254,2, 4, 0,  91.4260447751259,
      3, 4, 0, 0.875600661808945,4, 4, 0, -17.1397577419788,
      2, 5, 0, -21.6603240875311,4, 5, 0,  2.49697009569508,
      2, 6, 0,  2.13016970847183,2, 0, 1, -3310.49154044839,
      3, 0, 1,  199.459603073901,4, 0, 1, -54.7919133532887,
      5, 0, 1,  36.0284195611086,2, 1, 1,  729.116529735046,
      3, 1, 1, -175.292041186547,4, 1, 1, -22.6683558512829,
      2, 2, 1, -860.764303783977,3, 2, 1,  383.058066002476,
      2, 3, 1,  694.244814133268,3, 3, 1, -460.319931801257,
      2, 4, 1, -297.728741987187,3, 4, 1,  234.565187611355,
      2, 0, 2,  384.794152978599,3, 0, 2, -52.2940909281335,
      4, 0, 2, -4.08193978912261,2, 1, 2, -343.956902961561,
      3, 1, 2,  83.1923927801819,2, 2, 2,  337.409530269367,
      3, 2, 2, -54.1917262517112,2, 3, 2, -204.889641964903,
      2, 4, 2,  74.7261411387560,2, 0, 3, -96.5324320107458,
      3, 0, 3,  68.0444942726459,4, 0, 3, -30.1755111971161,
      2, 1, 3,  124.687671116248,3, 1, 3, -29.4830643494290,
      2, 2, 3, -178.314556207638,3, 2, 3,  25.6398487389914,
      2, 3, 3,  113.561697840594,2, 4, 3, -36.4872919001588,
      2, 0, 4,  15.8408172766824,3, 0, 4, -3.41251932441282,
      2, 1, 4, -31.6569643860730,2, 2, 4,  44.2040358308000,
      2, 3, 4, -11.1282734326413,2, 0, 5, -2.62480156590992,
      2, 1, 5,  7.04658803315449,2, 2, 5, -7.92001547211682 ),
      ncol=4, byrow=TRUE)

    Carr[,c(2,3)] <- Carr[,c(2,3)]+1
    gijk[cbind(Carr[,1],Carr[,2],Carr[,3])] <- Carr[,4]

    nj <- 7
    nk <- 6
    ni <- 7

    if (dt >= 1) {
      nj   <- nj-1
      gijk <- gijk[,-1,]
      for (j in 2:nj) gijk[,j,] <- gijk[,j,]*j
    }
    if (dt == 2) {
      nj   <- nj-1
      gijk <- gijk[,-1,]
      for (j in 2:nj) gijk[,j,] <- gijk[,j,]*j
    }
    if (dp >= 1) {
      nk   <- nk-1
      gijk <- gijk[,,-1]
      for (k in 2:nk) gijk[,,k] <- gijk[,,k]*k
    }
    if (dp ==2) {
      nk   <- nk-1
      gijk <- gijk[,,-1]
      for (k in 2:nk) gijk[,,k] <- gijk[,,k]*k
    }
    x2 <- S/Su
    x  <- sqrt(x2)

    Gsea <- 0
    if (dS==0 & dt==0 ) {
      for (j in 1:nj){
        for (k in 1:nk) {
          it <- 0
          for (i in 2:ni) it <- it + gijk[i,j,k]*x^(i)
          Gsea <- Gsea + (gijk[1,j,k]*x2*log(x) +it) *y^(j-1) *z^(k-1)
        }
      }
    } else if (dS==0) {
      for (j in 1:nj) {
        for (k in 1:nk) {
          it <- 0
          for (i in 2:ni) it <- it + gijk[i,j,k]*x^(i)
          Gsea <- Gsea  + it*y^(j-1)*z^(k-1)
        }
      }
      if (dt == 1) Gsea <- Gsea + gijk[1,1,1]*x2*log(x)
      if (dp == 1) Gsea <- Gsea + gijk[1,2,1]*x2*log(x)

    } else if (dS == 1) {
      for (j in 1:nj) {
        for (k in 1:nk) {
          it <- 0
          for (i in 2:ni) it <- it + 0.5*i*gijk[i,j,k]/Su*x2^(0.5*i-1)
          Gsea <- Gsea  + it*y^(j-1)*z^(k-1)
        }
      }
      Gsea <- Gsea + (log(x) +0.5)/Su*(gijk[1,1,1]+gijk[1,2,1]*y)
      if (dt == 1) Gsea <- Gsea + gijk[1,2,1]*0.5*log(x)
    }
    Gsea[is.nan(Gsea)]<-0
  }
  Gibbs <- Gsea+Gpure
  if (length(isnaS))
    Gibbs[isnaS] <- NA

  return(Gibbs*mfac)
}