Home > omviz > swstate.m

swstate

PURPOSE ^

SWSTATE State equation for seawater

SYNOPSIS ^

function [SVAN,SIGMA]=swstate(S,T,P0,ftn);

DESCRIPTION ^

 SWSTATE State equation for seawater

        [SVAN,SIGMA]=SWSTATE(S,T,P) returns the specific volume
        anomaly SVAN (m^3/kg*1e-8) and the density anomaly SIGMA (kg/m^3)
        given the salinity S (ppt), temperature T (deg C) and pressure
        P (dbars).

        [dVdT,dRdT]=SWSTATE(S,T,P,'dT') returns derivatives w.r.t.
        temperature of the volume and density.

        [dVdS,dRdS]=SWSTATE(S,T,P,'dS') returns derivatives w.r.t.
        salinity.

        [dVdP,dRdP]=SWSTATE(S,T,P,'dP') returns derivatives w.r.t.
        pressure.

        All elements can be scalars, vectors, or matrices but should be
        the same size.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [SVAN,SIGMA]=swstate(S,T,P0,ftn);
0002 % SWSTATE State equation for seawater
0003 %
0004 %        [SVAN,SIGMA]=SWSTATE(S,T,P) returns the specific volume
0005 %        anomaly SVAN (m^3/kg*1e-8) and the density anomaly SIGMA (kg/m^3)
0006 %        given the salinity S (ppt), temperature T (deg C) and pressure
0007 %        P (dbars).
0008 %
0009 %        [dVdT,dRdT]=SWSTATE(S,T,P,'dT') returns derivatives w.r.t.
0010 %        temperature of the volume and density.
0011 %
0012 %        [dVdS,dRdS]=SWSTATE(S,T,P,'dS') returns derivatives w.r.t.
0013 %        salinity.
0014 %
0015 %        [dVdP,dRdP]=SWSTATE(S,T,P,'dP') returns derivatives w.r.t.
0016 %        pressure.
0017 %
0018 %        All elements can be scalars, vectors, or matrices but should be
0019 %        the same size.
0020 
0021 %Notes: RP (WHOI 2/dec/91)
0022 %
0023 %  This stuff is directly copied from the UNESCO algorithms, with some
0024 %  minor changes to make it Matlab compatible (like adding ";" and changing
0025 %  "*" to ".*" when necessary.
0026 %
0027 %      RP  (WHOI 3/dec/91)
0028 %
0029 %  Added first derivative calculations.
0030 
0031 derivT=0;
0032 derivS=0;
0033 derivP=0;
0034 
0035 % Rich Signell modification.  If no pressure is specified,
0036 % assume pressure is 0 (sea surface).
0037 %
0038 if (nargin==2),     
0039    P0=zeros(size(T));
0040 end
0041 
0042 if (nargin==4),
0043    if     (ftn=='dT'), derivT=1;
0044    elseif (ftn=='dS'), derivS=1;
0045    elseif (ftn=='dP'), derivP=1;
0046    else error('swstate: Unrecognized option!');
0047    end;
0048 end;
0049 
0050 % ******************************************************
0051 % SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
0052 % OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
0053 % REFERENCES
0054 % MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
0055 % MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
0056 % BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
0057 % UNITS:
0058 %       PRESSURE        P0       DECIBARS
0059 %       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
0060 %       SALINITY        S        (IPSS-78)
0061 %       SPEC. VOL. ANA. SVAN     M**3/KG *1.0E-8
0062 %       DENSITY ANA.    SIGMA    KG/M**3
0063 % ******************************************************************
0064 % CHECK VALUE: SVAN=981.3021 E-8 M**3/KG.  FOR S = 40 (IPSS-78) ,
0065 % T = 40 DEG C, P0= 10000 DECIBARS.
0066 % CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78) ,
0067 % T = 40 DEG C, P0= 10000 DECIBARS.
0068 %HECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS.
0069 %        DR/DP                  DR/DT                 DR/DS
0070 %       DRV(1,7)              DRV(2,3)             DRV(1,8)
0071 %
0072 % FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION
0073 %
0074 %       3.46969238E-3       -.43311722           .705110777
0075 %
0076 % EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80
0077 %
0078 %       3.4696929E-3        -.4331173            .7051107
0079 %
0080 % (RP...I think this ---------^^^^^^ should be -.4431173!);
0081 
0082 
0083 % *******************************************************
0084 % DATA
0085    R3500=1028.1063;
0086    R4=4.8314E-4;
0087    DR350=28.106331;
0088 
0089 % CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
0090       P=P0/10.;
0091       SAL=S;
0092       SR = sqrt(abs(S));
0093 % *********************************************************
0094 % PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
0095 %   BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
0096 %
0097       R1 = ((((6.536332E-9*T-1.120083E-6).*T+1.001685E-4).*T ...
0098             -9.095290E-3).*T+6.793952E-2).*T-28.263737;
0099 % SEAWATER DENSITY ATM PRESS.
0100 %  COEFFICIENTS INVOLVING SALINITY
0101       R2 = (((5.3875E-9*T-8.2467E-7).*T+7.6438E-5).*T-4.0899E-3).*T+8.24493E-1; 
0102       R3 = (-1.6546E-6*T+1.0227E-4).*T-5.72466E-3;
0103 %  INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
0104       SIG = (R4*S + R3.*SR + R2).*S + R1;
0105 % SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
0106       V350P = 1.0/R3500;
0107       SVA = -SIG*V350P./(R3500+SIG);
0108       SIGMA=SIG+DR350;
0109       V0 = 1.0./(1000.0 + SIGMA);
0110 %  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
0111       SVAN=SVA*1.0E+8;
0112 
0113 if (derivS),               % These are derivatives for (S,T,0).
0114       R4S=9.6628E-4;
0115       RHO1 = 1000.0 + SIGMA;
0116 
0117       RHOS=R4S*SAL+1.5.*R3.*SR+R2;  
0118       V0S=-RHOS./(RHO1.*RHO1);  
0119 elseif (derivT),
0120       R1 =(((3.268166E-8*T-4.480332E-6).*T+3.005055E-4).*T...
0121           -1.819058E-2).*T+6.793952E-2;
0122       R2 = ((2.155E-8*T-2.47401E-6).*T+1.52876E-4).*T-4.0899E-3;
0123       R3 = -3.3092E-6*T+1.0227E-4;
0124       RHO1 = 1000.0 + SIGMA;
0125 
0126       RHOT = (R3.*SR + R2).*SAL + R1;    
0127       V0T = -RHOT./(RHO1.*RHO1);
0128 end;
0129       
0130 % ******************************************************************
0131 % ******  NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
0132 % ******************************************************************
0133 %        MILLERO, ET AL , 1980 DSR 27A, PP 255-264
0134 %               CONSTANT NOTATION FOLLOWS ARTICLE
0135 %********************************************************
0136 % COMPUTE COMPRESSION TERMS
0137       E = (9.1697E-10*T+2.0816E-8).*T-9.9348E-7;
0138       BW = (5.2787E-8*T-6.12293E-6).*T+3.47718E-5;
0139       B = BW + E.*S;    % Bulk Modulus (almost)
0140 %  CORRECT B FOR ANAMOLY BIAS CHANGE
0141       Bout = B + 5.03217E-5;
0142 
0143 if (derivS),
0144       DBDS=E;
0145 elseif (derivT),
0146       BW = 1.05574E-7*T-6.12293E-6;
0147       E = 1.83394E-9*T +2.0816E-8;
0148       BT = BW + E.*SAL;
0149 end;
0150 %
0151       D = 1.91075E-4;
0152       C = (-1.6078E-6*T-1.0981E-5).*T+2.2838E-3;
0153       AW = ((-5.77905E-7*T+1.16092E-4).*T+1.43713E-3).*T-0.1194975;
0154       A = (D*SR + C).*S + AW;    
0155 %  CORRECT A FOR ANAMOLY BIAS CHANGE
0156       Aout = A + 3.3594055;
0157 
0158 if (derivS),
0159       DADS=2.866125E-4*SR+C;
0160 elseif (derivT),
0161       C = -3.2156E-6*T -1.0981E-5;
0162       AW = (-1.733715E-6*T+2.32184E-4).*T+1.43713E-3;
0163       AT = C.*SAL + AW;
0164 end;
0165            
0166       B1 = (-5.3009E-4*T+1.6483E-2).*T+7.944E-2;
0167       A1 = ((-6.1670E-5*T+1.09987E-2).*T-0.603459).*T+54.6746;
0168       KW = (((-5.155288E-5*T+1.360477E-2).*T-2.327105).*T+148.4206).*T-1930.06;
0169       K0 = (B1.*SR + A1).*S + KW;
0170 
0171 if (derivS),
0172       K0S=1.5*B1.*SR+A1;
0173       KS=(DBDS.*P+DADS).*P+K0S;
0174 elseif (derivT),
0175       B1 = -1.06018E-3*T+1.6483E-2;
0176       % APRIL 9 1984 CORRECT A1 BIAS FROM -.603457 !!!
0177       A1 = (-1.8501E-4*T+2.19974E-2).*T-0.603459;
0178       KW = ((-2.0621152E-4*T+4.081431E-2).*T-4.65421).*T+148.4206;
0179       K0T = (B1.*SR+A1).*SAL + KW;
0180       KT = (BT.*P + AT).*P + K0T;
0181 end;
0182 
0183 
0184 % EVALUATE PRESSURE POLYNOMIAL
0185 % ***********************************************
0186 %   K EQUALS THE SECANT BULK MODULUS OF SEAWATER
0187 %   DK=K(S,T,P)-K(35,0,P)
0188 %  K35=K(35,0,P)
0189 % ***********************************************
0190       DK = (B.*P + A).*P + K0;
0191       K35  = (5.03217E-5*P+3.359406).*P+21582.27;
0192       GAM=P./K35;
0193       PK = 1.0 - GAM;
0194       SVA = SVA.*PK + (V350P+SVA).*P.*DK./(K35.*(K35+DK));
0195 %  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
0196       SVAN=SVA*1.0E+8;      % Volume anomaly
0197       V350P = V350P.*PK;
0198 %  ****************************************************
0199 % COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
0200 %  1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS
0201 %  2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C ,  PRES. VARIATION
0202 %  3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
0203 % ********************************************************************
0204 % CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78),
0205 % T = 40 DEG C, P0= 10000 DECIBARS.
0206 % *******************************************************
0207       DR35P=GAM./V350P;
0208       DVAN=SVA./(V350P.*(V350P+SVA));
0209       SIGMA=DR350+DR35P-DVAN;  % Density anomaly
0210 
0211       K=K35+DK;
0212       VP=1.0-P./K;
0213       V = (1.) ./(SIGMA+1000.0);
0214 
0215 if (derivS),
0216       VS=V0S.*VP+V0.*P.*KS./(K.*K);
0217 
0218       SVAN=VS;              % dVdS
0219       SIGMA=-VS./(V.*V);    % dRdS
0220 elseif (derivT),
0221       VT = V0T.*VP + V0.*P.*KT./(K.*K);
0222 
0223       SVAN=VT;              % dVdT
0224       SIGMA=-VT./(V.*V);    % dRdT
0225 elseif (derivP),
0226       DKDP = 2.0*Bout.*P + Aout;   
0227 % CORRECT DVDP TO PER DECIBAR BY MULTIPLE *.1
0228       DVDP = -.1*V0.*(1.0 - P.*DKDP./K)./K;
0229 
0230       SVAN=DVDP;            % dVdP
0231       SIGMA=-DVDP./(V.*V);  % dRdP
0232 end;
0233 
0234 
0235 
0236

Generated on Thu 01-Dec-2005 08:50:35 by m2html © 2003