function ATG=adiabat(S,T,P);
% ADIABAT Computes the adiabatic temperature gradient (called by POTENTMP).
%
%       ADIABAT(S,T,P) returns the gradient. Units are:
%   
%       PRESSURE        P        DECIBARS
%       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
%       SALINITY        S        (IPSS-78)
%       ADIABATIC       ATG      DEG. C/DECIBAR
%
% REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408

% Notes: RP 29/Nov/91
%
% I have modified the FORTRAN code to make it Matlab compatible, but
% no numbers have been changed. In certain places "*" has been replaced
% with ".*" to allow vectorization.
%
% This routine is called by potentmp.m, and is called ATG in the
% UNESCO algorithms.

%C CHECKVALUE: ATG=3.255976E-4 C/DBAR FOR S=40 (IPSS-78),
%C T=40 DEG C,P0=10000 DECIBARS
%       IMPLICIT REAL*8 (A-H,O-Z)

      DS = S - 35.0 ;
      ATG = (((-2.1687E-16*T+1.8676E-14).*T-4.6206E-13).*P+...
((2.7759E-12*T-1.1351E-10).*DS+((-5.4481E-14*T+8.733E-12).*T-6.7795E-10).*T+...
1.8741E-8)).*P+(-4.2393E-8*T+1.8932E-6).*DS+((6.6228E-10*T-6.836E-8).*T+...
8.5258E-6).*T+3.5803E-5;
% aou.m                                            by:  Edward T Peltzer, MBARI
%                                                  revised:  2007 Apr 26.
%
% CALCULATE OXYGEN CONCENTRATION AT SATURATION
%
% Source:  The solubility of nitrogen, oxygen and argon in water and
%         seawater - Weiss (1970) Deep Sea Research V17(4): 721-735.
%
% Molar volume of oxygen at STP obtained from NIST website on the
%          thermophysical properties of fluid systems:
%
%          http://webbook.nist.gov/chemistry/fluid/
%
%
% CALCULATE AOU BY DIFFERENCE:
%
%         AOU (umol/kg) = sat O2 (umol/kg) - obs o2 (umol/kg).
%
%
% Input:       S = Salinity (pss-78)
%              T = Potential Temp (deg C)
%              O2 = Meas'd Oxygen Conc (umol/kg)
%
% Output:      Apparant Oxygen Utilization (umol/kg).
%
%                        AOU = aou(S,T,O2).

function [AOU]=aou(S,T,O2)


% DEFINE CONSTANTS, ETC FOR SATURATION CALCULATION

%    The constants -177.7888, etc., are used for units of ml O2/kg.

  T1 = (T + 273.15) ./ 100;

  OSAT = -177.7888 + 255.5907 ./ T1 + 146.4813 .* log(T1) - 22.2040 .* T1;
  OSAT = OSAT + S .* (-0.037362 + T1 .* (0.016504 - 0.0020564 .* T1));
  OSAT = exp(OSAT);


% CONVERT FROM ML/KG TO UM/KG

  OSAT = OSAT * 1000 ./ 22.392;


% CALCULATE AOU

  AOU = OSAT - O2;
% densatp.m                                      by:  Edward T Peltzer, MBARI
%                                                  revised:  29 Jan 98.
%
% CALCULATE THE DENSITY OF SEAWATER AT A GIVEN S, T and P.
% Equation of State is from Millero & Poisson (1981) DSR V28: 625-629.
%
% INPUT:       	Salinity (S) in g/kg or pss.
%		Temperature (T) in degrees C.
%		Pressure (P) in decibar.
%
% OUTPUT:	Density [rhp] in g/cc:
%
%			rhp = densatp(S,T,P).

function [rhp]=densatp(S,T,P)

% DEFINE CONSTANTS FOR EQUATION OF STATE

  R0=+9.99842594E2;
  R1=+6.793952E-2;
  R2=-9.095290E-3;
  R3=+1.001685E-4;
  R4=-1.120083E-6;
  R5=+6.536332E-9;

  A0=+8.24493E-1;
  A1=-4.0899E-3;
  A2=+7.6438E-5;
  A3=-8.2467E-7;
  A4=+5.3875E-9;

  B0=-5.72466E-3;
  B1=+1.0227E-4;
  B2=-1.6546E-6;

  C=+4.8314E-4;

% CALCULATE RHO

    SR=sqrt(S);

    RHO0=R0+T.*(R1+T.*(R2+T.*(R3+T.*(R4+T.*R5))));
    A=A0+T.*(A1+T.*(A2+T.*(A3+T.*A4)));
    B=B0+T.*(B1+T.*B2);

    RHO=RHO0+S.*(A+B.*SR+C.*S);

% DEFINE CONSTANTS FOR SECANT BULK MODULUS

  D0=+1.965221E+4;
  D1=+1.484206E+2;
  D2=-2.327105E+0;
  D3=+1.360477E-2;
  D4=-5.155288E-5;

  E0=+5.46746E+1;
  E1=-6.03459E-1;
  E2=+1.09987E-2;
  E3=-6.1670E-5;

  F0=+7.944E-2;
  F1=+1.6483E-2;
  F2=-5.3009E-4;

  G0=+3.239908E+0;
  G1=+1.43713E-3;
  G2=+1.16092E-4;
  G3=-5.77905E-7;

  H0=+2.2838E-3;
  H1=-1.0981E-5;
  H2=-1.6078E-6;
  H3=+1.91075E-4;

  I0=+8.50935E-5;
  I1=-6.12293E-6;
  I2=+5.2787E-8;

  J0=-9.9348E-7;
  J1=+2.0816E-8;
  J2=+9.1697E-10;

% CORRECT P IN DECI-BARS TO BARS

    Pc=P./10.;

% CALCULATE K

    H=H0+T.*(H1+T.*H2)+H3.*SR;
    J=J0+T.*(J1+T.*J2);

    K1=D0+T.*(D1+T.*(D2+T.*(D3+T.*D4)));
    K2=E0+T.*(E1+T.*(E2+T.*E3));
    K3=F0+T.*(F1+T.*F2);
    K4=G0+T.*(G1+T.*(G2+T.*G3))+S.*H;
    K5=I0+T.*(I1+T.*I2)+S.*J;

    K=K1+S.*(K2+SR.*K3)+Pc.*(K4+Pc.*K5);

% CORRECT FOR PRESSURE

    RHP=RHO.*(1./(1-Pc./K));

% CONVERT KG/M3 TO g/cc

    rhp = RHP ./ 1000;
% density.m                                      by:  Edward T Peltzer, MBARI
%                                                  revised:  29 Jan 98.
%
% CALCULATE THE DENSITY OF SEAWATER AT A GIVEN S & T, P = 1 ATM
% Equation of State is from Millero & Poisson (1981) DSR V28: 625-629.
%
% INPUT:       	Salinity (S) in g/kg or pss.
%		Temperature (T) in degrees C.
%
% OUTPUT:	Density [rho] in g/cc:
%
%			rho = density(S,T).

function [rho]=density(S,T)

% DEFINE CONSTANTS FOR EQUATION OF STATE

  R0=+9.99842594E2;
  R1=+6.793952E-2;
  R2=-9.095290E-3;
  R3=+1.001685E-4;
  R4=-1.120083E-6;
  R5=+6.536332E-9;

  A0=+8.24493E-1;
  A1=-4.0899E-3;
  A2=+7.6438E-5;
  A3=-8.2467E-7;
  A4=+5.3875E-9;

  B0=-5.72466E-3;
  B1=+1.0227E-4;
  B2=-1.6546E-6;

  C=+4.8314E-4;

% CALCULATE RHO

    RHO0=R0+T.*(R1+T.*(R2+T.*(R3+T.*(R4+T.*R5))));

    A=A0+T.*(A1+T.*(A2+T.*(A3+T.*A4)));
    B=B0+T.*(B1+T.*B2);
    RHO=RHO0+S.*(A+B.*sqrt(S)+C.*S);

% CONVERT KG/M3 TO g/cc

    rho = RHO ./ 1000;
function DEPTH=depth(P,LAT);
% DEPTH   Computes depth given the pressure at some latitude
%         D=DEPTH(P,LAT) gives the depth D (m) for a pressure P
%             (dbars) at some latitude LAT (degrees).
%
%         Fofonoff and Millard (1982). UNESCO Tech Paper #44.
%
% Notes:  (ETP3, MBARI)
%         This algorithm was originally compiled by RP @ WHOI.
%         It was copied from the UNESCO technical report.
%         The algorithm was endorsed by SCOR Working Group 51.
%         The equations were originally developed by Saunders
%             and Fofonoff (1976).  DSR 23: 109-111.
%         The parameters were re-fit for the 1980 equation of
%             state for seawater (EOS80).
%
% CHECKVALUE: D = 9712.653 M FOR P=10000 DECIBARS, LAT=30 DEG
%
% CALCULATON ASSUMES STD OCEAN: T = 0 DEG C; S = 35 (IPSS-78)


  X = sin(LAT/57.29578);
  X = X.*X;

% GR= GRAVITY VARIATION WITH LAT: ANON (1970) BULLETIN GEODESIQUE

  GR = 9.780318*(1.0+(5.2788E-3+2.36E-5*X).*X) + 1.092E-6.*P;

  DEPTH = (((-1.82E-15*P+2.279E-10).*P-2.2512E-5).*P+9.72659).*P;
  DEPTH = DEPTH./GR;
% n2sat.m                                          by:  Edward T Peltzer, MBARI
%                                                  revised:  2007 Apr 25.
%
% CALCULATE NITROGEN CONCENTRATION AT SATURATION
%
% Source:  The solubility of nitrogen, oxygen and argon in water and
%          seawater - Weiss (1970) Deep Sea Research V17(4): 721-735.
%
% Molar volume of nitrogen at STP obtained from NIST website on the
%          thermophysical properties of fluid systems:
%
%          http://webbook.nist.gov/chemistry/fluid/
%
%
% Input:       S = Salinity (0/00)
%              T = Temp (deg C)
%
% Output:      Nitrogen saturation at one atmosphere (umol/kg).
%
%                        N2 = n2sat(S,T).

function [N2] = n2sat(S,T)


% DEFINE CONSTANTS, ETC FOR SATURATION CALCULATION

%    The constants -177.0212, etc., are used for units of ml N2/kg.

  T1 = (T + 273.15) ./ 100;

  NSAT = -177.0212 + 254.6078 ./ T1 + 146.3611 .* log(T1) - 22.0933 .* T1;
  NSAT = NSAT + S .* (-0.054052 + T1 .* (0.027266 - 0.0038430 .* T1));
  NSAT = exp(NSAT);


% CONVERT FROM ML/KG TO UM/KG

  N2 = NSAT * 1000 ./ 22.404;
% o2sat.m                                          by:  Edward T Peltzer, MBARI
%                                                  revised:  2007 Apr 26.
%
% CALCULATE OXYGEN CONCENTRATION AT SATURATION
%
% Source:  The solubility of nitrogen, oxygen and argon in water and
%         seawater - Weiss (1970) Deep Sea Research V17(4): 721-735.
%
% Molar volume of oxygen at STP obtained from NIST website on the
%          thermophysical properties of fluid systems:
%
%          http://webbook.nist.gov/chemistry/fluid/
%
%
% Input:       S = Salinity (pss-78)
%              T = Temp (deg C)
%
% Output:      Oxygen saturation at one atmosphere (umol/kg).
%
%                        O2 = o2sat(S,T).

function [O2] = o2sat(S,T)


% DEFINE CONSTANTS, ETC FOR SATURATION CALCULATION

%    The constants -177.7888, etc., are used for units of ml O2/kg.

  T1 = (T + 273.15) ./ 100;

  OSAT = -177.7888 + 255.5907 ./ T1 + 146.4813 .* log(T1) - 22.2040 .* T1;
  OSAT = OSAT + S .* (-0.037362 + T1 .* (0.016504 - 0.0020564 .* T1));
  OSAT = exp(OSAT);


% CONVERT FROM ML/KG TO UM/KG

  O2 = OSAT * 1000 ./ 22.392;
function THETA=potentmp(S,T0,P0,PR);
% POTENTMP Potential temperature from in-situ measurements.
%          THETA=POTENTMP(S,TO,PO,PR) computes the potential temp.
%          at a reference pressure PR (dbars) corresponding to the 
%          salinity S (ppt) and temperature TO (deg C) at pressure 
%          PO (dbars). The formula has been copied from the UNESCO
%          algorithms (See comments for details).
%
%          if PR is omitted, it is assumed to be zero.
%
%          S, T0, P0 can be vectors or matrices. If S (T) is a scalar,
%          and T (S) is a vector or matrix, we assume that the scalar
%          value corresponds to all values in the vector. 
%  
%          If P0 is a vector, and one of S,T0 is a matrix, we assume
%          that elements of P0 correspond with all row entries in the
%          matrices.
%
% REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
%      FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491


if (nargin==3), PR=0; end;

[NS,MS]=size(S);
[NT,MT]=size(T0);
if ((NS==1 & MS==1) & (NT ~=1 | MT ~=1) ),
   S=S*ones(NT,MT);
   [NS,MS]=size(S);
elseif ((NT==1 & MT==1) & (NS ~=1 | MS ~=1) ),
   T0=T0*ones(NS,MS);
   [NT,MT]=size(T0);
end;

[NP,MP]=size(P0);
if (MP==1), P0=P0*ones(1,MT); 
elseif (NP==1), P0=P0'*ones(1,MT); end;



% Notes: RP 29/Nov/91
%
% I have modified the FORTRAN code to make it Matlab compatible, but
% no numbers have been changed. In certain places "*" has been replaced
% with ".*" to allow vectorization.
%
% This routine calls adiabat.m (renamed from ATG).

%C ***********************************
%C TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR
%C USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE
%C AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.
%C UNITS:      
%C       PRESSURE        P0       DECIBARS
%C       TEMPERATURE     T0       DEG CELSIUS (IPTS-68)
%C       SALINITY        S        (IPSS-78)
%C       REFERENCE PRS   PR       DECIBARS
%C       POTENTIAL TMP.  THETA    DEG CELSIUS 
%C CHECKVALUE: THETA= 36.89073 C,S=40 (IPSS-78),T0=40 DEG C,
%C P0=10000 DECIBARS,PR=0 DECIBARS
%C             
%C      SET-UP INTERMEDIATE TEMPERATURE AND PRESSURE VARIABLES
%      IMPLICIT REAL*8 (A-H,O-Z)
      P=P0;
      T=T0;
%C**************
      H = PR - P;
      XK = H.*adiabat(S,T,P) ;
      T = T + 0.5*XK;
      Q = XK  ;
      P = P + 0.5*H ;
      XK = H.*adiabat(S,T,P) ;
      T = T + 0.29289322*(XK-Q) ;
      Q = 0.58578644*XK + 0.121320344*Q ;
      XK = H.*adiabat(S,T,P) ;
      T = T + 1.707106781*(XK-Q);
      Q = 3.414213562*XK - 4.121320344*Q;
      P = P + 0.5*H ;
      XK = H.*adiabat(S,T,P) ;
      THETA = T + (XK-2.0*Q)/6.0;
function P80=pressure(DPTH,XLAT);
% PRESSURE Computes pressure given the depth at some latitude
%          P=PRESSURE(D,LAT) gives the pressure P (dbars) at a depth D (m)
%          at some latitude LAT (degrees).
%
%          This probably works best in mid-latitude oceans, if anywhere!
%
%          Ref: Saunders, "Practical Conversion of Pressure to Depth",
%              J. Phys. Oceanog., April 1981.
% 

%Notes: RP (WHOI) 2/Dec/91
%         I copied this directly from the UNESCO algorithms.


% CHECK VALUE: P80=7500.004 DBARS;FOR LAT=30 DEG., DEPTH=7321.45 METERS

      PLAT=abs(XLAT*pi/180.);
      D=sin(PLAT);
      C1=5.92E-3+(D.*D)*5.25E-3;
      P80=((1-C1)-sqrt(((1-C1).^2)-(8.84E-6*DPTH)))/4.42E-6;
% salinity.m                                       by:  Edward T Peltzer, MBARI
%                                                  revised:  2007 Apr 28.
%
% SALINITY CALCULATION (pss) FROM CONDUCTIVITY.
%
% Reference: Fofonff, P. and Millard, R.C. Jr.  Unesco 1983.
%    Algorithms for computation of fundamental properties of seawater.
%    Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp.
%
% Input:       	C   - in situ conductivity  (S/m)
%               T   - temperature in degree C
%               P   - pressure in dbars  (not SI)
%               C15 - conductivity at 15 deg C
%                     MBARI default value = 4.2914
%
% Internal Variables:
%
%		R   - ratio of in situ conductivity to standard conductivity
%		rt  - ratio of conductivity of standard sea water S=35,T=t to S=35,T=15
%		Rp  - ratio of in situ conductivity to same at P=0
%		Rt  - ratio of sample conductivity to standard conductivity at T=t
%
%		S   - salinity at T=15
%		dS  - delta correction for T not equal 15
%
% Output:      
%              Salt = salinity(C,T,P).


function [Salt] = salinity(C,T,P)


% CHECK INPUTS ARE SAME DIMENSIONS

[mc,nc] = size(C);
[mt,nt] = size(T);
[mp,np] = size(P);

if ~(mc==mt | mc==mp | nc==nt | nc==np)
  error('sw_salt.m: cndr,T,P must all have the same dimensions')
end


% DEFINE CONSTANTS, ETC FOR CALCULATION

  C15= 4.2914;

  a0=  0.008;
  a1= -0.1692;
  a2= 25.3851;
  a3= 14.0941;
  a4= -7.0261;
  a5=  2.7081;

  b0=  0.0005;
  b1= -0.0056;
  b2= -0.0066;
  b3= -0.0375;
  b4=  0.0636;
  b5= -0.0144;

  c0=  0.6766097;
  c1=  2.00564e-2;
  c2=  1.104259e-4;
  c3= -6.9698e-7;
  c4=  1.0031e-9;

  d1=  3.426e-2;
  d2=  4.464e-4;
  d3=  4.215e-1;
  d4= -3.107e-3;

% The e# coefficients reflect the use of pressure in dbar
%   rather that in Pascals (SI).

  e1=  2.07e-5;
  e2= -6.37e-10;
  e3=  3.989e-15;

  k= 0.0162;


% Calculate internal variables

R = C ./ C15;
rt = c0+(c1+(c2+(c3+c4.*T).*T).*T).*T;
Rp = 1.0 + (e1+(e2+e3.*P).*P).*P ./ (1.0 + (d1+d2.*T).*T + (d3+d4.*T).*R);
Rt = R ./ Rp ./ rt;
sqrt_Rt = sqrt(Rt);


% Calculate salinity

Salt = a0 + (a1+(a3+a5.*Rt).*Rt).*sqrt_Rt + (a2+a4.*Rt).*Rt;

dS = b0 + (b1+(b3+b5.*Rt).*Rt).*sqrt_Rt + (b2+b4.*Rt).*Rt;
dS = dS .* (T-15) ./ (1+k.*(T-15));

Salt = Salt + dS;
% sbec2sal.m                                       by:  Edward T Peltzer, MBARI
%                                                  revised:  2007 Apr 28.
%
% SALINITY CALCULATION (pss) FROM CONDUCTIVITY.
%
% Reference: Fofonff, P. and Millard, R.C. Jr.  Unesco 1983.
%    Algorithms for computation of fundamental properties of seawater.
%    Unesco Tech. Pap. in Mar. Sci., No. 44, 53 pp.
%
% Input:       	C   - in situ conductivity  (mS/cm)
%               T   - temperature in degree C
%               P   - pressure in dbars  (not SI)
%               C15 - conductivity at 15 deg C
%                     SeaBird default value = 42.914
%
% Internal Variables:
%
%		R   - ratio of in situ conductivity to standard conductivity
%		rt  - ratio of conductivity of standard sea water S=35,T=t to S=35,T=15
%		Rp  - ratio of in situ conductivity to same at P=0
%		Rt  - ratio of sample conductivity to standard conductivity at T=t
%
%		S   - salinity at T=15
%		dS  - delta correction for T not equal 15
%
% Output:      
%              Salt = salinity(C,T,P).


function [Salt] = salinity(C,T,P)


% CHECK INPUTS ARE SAME DIMENSIONS

[mc,nc] = size(C);
[mt,nt] = size(T);
[mp,np] = size(P);

if ~(mc==mt | mc==mp | nc==nt | nc==np)
  error('sw_salt.m: cndr,T,P must all have the same dimensions')
end


% DEFINE CONSTANTS, ETC FOR CALCULATION

  C15= 42.914;

  a0=  0.008;
  a1= -0.1692;
  a2= 25.3851;
  a3= 14.0941;
  a4= -7.0261;
  a5=  2.7081;

  b0=  0.0005;
  b1= -0.0056;
  b2= -0.0066;
  b3= -0.0375;
  b4=  0.0636;
  b5= -0.0144;

  c0=  0.6766097;
  c1=  2.00564e-2;
  c2=  1.104259e-4;
  c3= -6.9698e-7;
  c4=  1.0031e-9;

  d1=  3.426e-2;
  d2=  4.464e-4;
  d3=  4.215e-1;
  d4= -3.107e-3;

% The e# coefficients reflect the use of pressure in dbar
%   rather that in Pascals (SI).

  e1=  2.07e-5;
  e2= -6.37e-10;
  e3=  3.989e-15;

  k= 0.0162;


% Calculate internal variables

R = C ./ C15;
rt = c0+(c1+(c2+(c3+c4.*T).*T).*T).*T;
Rp = 1.0 + (e1+(e2+e3.*P).*P).*P ./ (1.0 + (d1+d2.*T).*T + (d3+d4.*T).*R);
Rt = R ./ Rp ./ rt;
sqrt_Rt = sqrt(Rt);


% Calculate salinity

Salt = a0 + (a1+(a3+a5.*Rt).*Rt).*sqrt_Rt + (a2+a4.*Rt).*Rt;

dS = b0 + (b1+(b3+b5.*Rt).*Rt).*sqrt_Rt + (b2+b4.*Rt).*Rt;
dS = dS .* (T-15) ./ (1+k.*(T-15));

Salt = Salt + dS;
