/* ******************************************************************* * * * Copyright (c) L-DGO/MIT/JGOFS * * * * * * File : svan.c * * * * Purpose : * * * * Version Number : 1.1 * * * * Revision History : * * * * Date Developer * * ---- --------- * * * * Sat Oct 17 1992 10.00 Glenn Flierl * * * * * ******************************************************************* */ /* #include */ #include /* float atg(); float theta(); float svan(); main() { float delta,sigma; delta=svan(40.0,40.0,10000.0,&sigma); printf("delta %e sigma %e\n",delta,sigma); printf("atg %e\n",atg(40.0,40.0,10000.0)); printf("theta %e\n",theta(40.0,40.0,10000.0,0.0)); } */ float svan(s,t,p0,sigma) float s,t,p0; float *sigma; { /* C CHECK VALUE SVAN=981.30210 AND SIGMA=59.82037,S=40.0 ,T=40 P=10000 DB C ****************************************************** C SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION C OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE. C REFERENCES C MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264 C MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629. C BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981) C UNITS: C PRESSURE P0 DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C SALINITY S (IPSS-78) C SPEC. VOL. ANA. SVAN M**3/KG *1.0E-8 C DENSITY ANA. SIGMA KG/M**3 C ******************************************************* */ float p,sig,r1,r2,r3,b1,k,k35; float v350p,sva; double sr; /* C ******************** C DATA */ #define r3500 1028.1063 #define r4 4.8314e-4 #define dr350 28.106331 /* C R4 IS REFERED TO AS C IN MILLERO AND POISSON 1981 C CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY. */ p=p0/10.0; sr = sqrt(fabs((double)s)); /* C ********************************************************* C PURE WATER DENSITY AT ATMOSPHERIC PRESSURE C BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537. C */ r1 = ((((6.536332e-9*t-1.120083e-6)*t+1.001685e-4)*t -9.095290e-3)*t+6.793952e-2)*t-28.263737; /* C SEAWATER DENSITY ATM PRESS. C COEFFICIENTS INVOLVING SALINITY C R2 = A IN NOTATION OF MILLERO AND POISSON 1981 */ r2 = (((5.3875e-9*t-8.2467e-7)*t+7.6438e-5)*t-4.0899e-3)*t +8.24493e-1 ; /* C R3 = B IN NOTATION OF MILLERO AND POISSON 1981 */ r3 = (-1.6546e-6*t+1.0227e-4)*t-5.72466e-3; /* C INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER */ sig = (r4*s + r3*sr + r2)*s + r1; /* C SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE */ v350p = 1.0/r3500; sva = -sig*v350p/(r3500+sig); *sigma=sig+dr350; /* C SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS */ if(p==0.0) return(sva*1.0e+8); /* C ****************************************************************** C ****** NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ******** C ****************************************************************** C MILLERO, ET AL , 1980 DSR 27A, PP 255-264 C CONSTANT NOTATION FOLLOWS ARTICLE C******************************************************** C COMPUTE COMPRESSION TERMS */ b1 = (9.1697e-10*t+2.0816e-8)*t-9.9348e-7; r3 = (5.2787e-8*t-6.12293e-6)*t+3.47718e-5; r3 = r3 + b1*s; /* C */ b1 = 1.91075e-4; r2 = (-1.6078e-6*t-1.0981e-5)*t+2.2838e-3; r1 = ((-5.77905e-7*t+1.16092e-4)*t+1.43713e-3)*t-0.1194975; r1 = (b1*sr + r2)*s + r1; /* C C COEFFICIENT K0 */ b1 = (-5.3009e-4*t+1.6483e-2)*t+7.944e-2; r2 = ((-6.1670e-5*t+1.09987e-2)*t-0.603459)*t+54.6746; k = (((-5.155288e-5*t+1.360477e-2)*t-2.327105)*t+148.4206)*t-1930.06; k = (b1*sr + r2)*s + k; /* C EVALUATE PRESSURE POLYNOMIAL C *********************************************** C K EQUALS THE SECANT BULK MODULUS OF SEAWATER C DK=K(S,T,P)-K(35,0,P) C K35=K(35,0,P) C *********************************************** */ b1 = (r3*p + r1)*p + k; k35 = (5.03217e-5*p+3.359406)*p+21582.27; r2=p/k35; r3 = 1.0 - r2; sva = sva*r3 + (v350p+sva)*p*b1/(k35*(k35+b1)); /* c SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS */ v350p = v350p*r3; /* C **************************************************** C COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3 C 1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS C 2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C , PRES. VARIATION C 3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY C ******************************************************************** C CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78), C T = 40 DEG C, P0= 10000 DECIBARS. C ******************************************************* */ r2=r2/v350p; r1=sva/(v350p*(v350p+sva)); *sigma=dr350+r2-r1; return(sva*1.0e+8); } float atg(s,t,p) float s,t,p; { float ds,atg1; /* c units adiabatic temp grad atg deg. c/decibar c check atg=3.255976e-4 for s=40 t=40 p0=10000 */ ds=s-35.0; atg1=(((-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; return (atg1); } float theta(s,t0,p0,pr) float s,t0,p0,pr; { float p,t,h,xk,q; /* c units ref press pr decibars c pot temp theta deg. C c check theta=36.89073, s=40,t=40,p=10000,pr=0 */ p=p0; t=t0; h=pr-p; xk=h*atg(s,t,p); t=t+0.5*xk; q=xk; p=p+0.5*h; xk=h*atg(s,t,p); t=t+0.29289322*(xk-q); q=0.58578644*xk+0.121320344*q; xk=h*atg(s,t,p); t=t+1.707106781*(xk-q); q=3.414213562*xk-4.121320344*q; p=p+0.5*h; xk=h*atg(s,t,p); return(t+(xk-2.0*q)/6.0); }