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.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