CND=1.888091,T=40 DEG C,P=10000 ! C=28.48 ! t=5.86 ! p=0. ! s=sal(c,t,p) ! print *,'sal=',s ! end C C-------------------------------------------------------- C C THIS FUNCTION IS A 'FRONT END' TO THE SALINITY ROUTINE C GIVEN IN UNESCO TECHNICAL PAPES IN MARINE SCIENCE, NO 44. C C THE CODE FOR THESE ROUTINES WAS GENEROUSLY SUPPLIED C BY WOOD'S HOLE. C C THE QUANTITY WHICH IS USUALLY KNOWN AS DEPTH IN OUR C SOFTWARE IS, IN FACT, PRESSURE. C ******************************************************************* C UNITS: C PRESSURE P DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C CONDUCTIVITY CND RATIO mS/cm (M=0) C SALINITY SAL78 (PSS-78) (M=0) C CHECKVALUES: C SAL78=1.888091 :CND= 40.0000,T=40 DEG C,P= 10000 DECIBARS: M= 1 C SAL78=40.00000 :CND=1.888091,T=40 DEG C,P=10000 DECIBARS: M=0 C------------------------------------------------------- C FUNCTION SAL(C,T,P) C C C35150 IS THE CONDUCTIVITY OF STANDARD SEAWATER AT 35 PRACTICAL C SALINITY UNITS, 15 DEGREES, AND ATMOSPHERIC PRESSURE C SAVE DATA C35150 / 42.914 / C IF(ABS(C).GT. 0.01)THEN C C ALLOW FOR EFFECTS OF TEMPERATURE AND PRESSURE ON THE C CONDUCTIVITY CELL C G=(C*(1.-6.5E-06*(T-2.8)+1.5E-8*(P-3000.)))/C35150 SAL=SAL78(G,T,P,0) ELSE C ERROR - CONDUCTIVITY EFFECTIVELY ZERO C SAL=0.0 ENDIF RETURN END C C----------------------------------------------------- C C FUNCTION TO CALCULATE CONDUCTIVITY FROM S,T AND P C USING THE ABOVE-MENTIONED ROUTINES C C----------------------------------------------------- C FUNCTION CONDUC(S,T,P) C SAVE DATA C35150 / 42.914 / C G=SAL78(S,T,P,1) C CONDUC=G*C35150/(1.-6.5E-06*(T-2.8)+1.5E-8*(P-3000.)) C RETURN END C C-------------------------------------------------- C C FUNCTION TO CALCULATE SIGMA-T, GIVEN THE SALINITY AND C TEMPERATURE. C C--------------------------------------------------- C FUNCTION SIGMAT(S,T) C SV=SVAN(S,T,0.,SIG) SIGMAT = SIG C RETURN END C C============================================================================== C C START OF UNESCO ROUTINES C C============================================================================== C SEPT. 28 1983 C ADD TF(S,P) FREEZING PT. C WHOI CTD GROUP DISK FILE SPEC=BLUE::CTDA:PHYPROPSW.FOR C C TITLE: ALGORITHMS FOR OCEANOGRAPHIC COMPUTATIONS C N FOFONOFF & R MILLARD C SAL78 FCN ********** MAR 28 1983 ******* REAL FUNCTION SAL78(CND,T,P,M) Czaj CND=Conductivity / C(35,15,0) while unit of C is ms/cm, when m=0 czaj CND=salinity unit is PPT, when M=1 C ******************************************************************** C THE CONDUCTIVITY RATIO (CND) = 1.0000000 FOR SALINITY = 35 PSS-78 C TEMPERATURE = 15.0 DEG. CELSIUS , AND ATMOSPHERIC PRESSURE. C******************************************************************** C C FUNCTION TO CONVERT CONDUCTIVITY RATIO TO SALINITY (M = 0) C SALINITY TO CONDUCTIVITY RATIO (M = 1,CND BECOMES INPUT SALINITY) C***************************************************************** C REFERENCES: ALSO LOCATED IN UNESCO REPORT # 37 1981 C PRACTICAL SALINITY SCALE 1978: E.L. LEWIS IEEE OCEAN ENG. JAN. 1980 C ******************************************************************* C UNITS: C PRESSURE P DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C CONDUCTIVITY CND RATIO (M=0) C SALINITY SAL78 (PSS-78) (M=0) C CHECKVALUES: C SAL78=1.888091 :CND= 40.0000,T=40 DEG C,P= 10000 DECIBARS: M= 1 C SAL78=40.00000 :CND=1.888091,T=40 DEG C,P=10000 DECIBARS: M=0 C******************************************************************* C SAL78 RATIO: RETURNS ZERO FOR CONDUCTIVITY RATIO: < 0.0005 C SAL78: RETURNS ZERO FOR SALINITY: < 0.02 C **************************************************************** C INTERNAL FUNCTIONS C **************************************************************** C PRACTICAL SALINITY SCALE 1978 DEFINITION WITH TEMPERATURE CORRECTION C XT=T-15.0 : XR=SQRT(RT) SAL(XR,XT) =((((2.7081*XR-7.0261)*XR+14.0941)*XR+25.3851)*XR X-0.1692)* XR+0.0080 X +(XT/(1.0+0.0162*XT))*(((((-0.0144*XR+ X 0.0636)*XR-0.0375)*XR-0.0066)*XR-0.0056)*XR+0.0005) C DSAL(XR,XT) FUNCTION FOR DERIVATIVE OF SAL(XR,XT) WITH XR. DSAL(XR,XT) =((((13.5405*XR-28.1044)*XR+42.2823)*XR+50.7702)*XR X -0.1692)+(XT/(1.0+0.0162*XT))*((((-0.0720*XR+0.2544)*XR X -0.1125)*XR-0.0132)*XR-0.0056) C FUNCTION RT35 : C(35,T,0)/C(35,15,0) VARIATION WITH TEMPERATURE C WITH TEMPERATURE. RT35(XT) = (((1.0031E-9*XT-6.9698E-7)*XT+1.104259E-4)*XT X + 2.00564E-2)*XT + 0.6766097 C POLYNOMIALS OF RP: C(S,T,P)/C(S,T,0) VARIATION WITH PRESSURE C C(XP) POLYNOMIAL CORRESPONDS TO A1-A3 CONSTANTS: LEWIS 1980 C(XP) = ((3.989E-15*XP-6.370E-10)*XP+2.070E-5)*XP B(XT) = (4.464E-4*XT+3.426E-2)*XT + 1.0 C A(XT) POLYNOMIAL CORRESPONDS TO B3 AND B4 CONSTANTS: LEWIS 1980 A(XT) = -3.107E-3*XT + 0.4215 C******************************************************************* C ZERO SALINITY/CONDUCTIVITY TRAP SAL78=0.0 ! IF((M.EQ.0).AND.(CND.LE.5E-4)) RETURN ! IF((M.EQ.1).AND.(CND.LE.0.02)) RETURN C *********************************************** DT = T - 15.0 C SELECT BRANCH FOR SALINITY (M=0) OR CONDUCTIVITY (M=1) IF(M .NE. 1)THEN C ************************************************ C CONVERT CONDUCTIVITY TO SALINITY R = CND RT = R/(RT35(T)*(1.0 + C(P)/(B(T) + A(T)*R))) RT = SQRT(ABS(RT)) SAL78 = SAL(RT,DT) !zaj extended to low salinity in estuarine according to EPA's formulation 09/12/2006 a0=0.008 b0=0.0005 X=400.0*RT Y=100.0*RT F=DT/(1.+0.0162*DT) DELTS=a0/(1+1.5*X+X*X)+b0*F/(1+Y**0.5+Y**1.5) SAL78=SAL78-DELTS C ********* END OF CONDUCTIVITY TO SALINITY SECTION *********** C ******************************************************* C INVERT SALINITY TO CONDUCTIVITY BY THE C NEWTON-RAPHSON ITERATIVE METHOD. C ****************************************************** C FIRST APPROXIMATION ELSE RT = SQRT(CND/35.0) SI = SAL(RT,DT) N = 0 C C ITERATION LOOP BEGINS HERE WITH A MAXIMUM OF 10 CYCLES C 15 RT = RT + (CND - SI)/DSAL(RT,DT) SI = SAL(RT,DT) N = N + 1 DELS = ABS(SI - CND) IF((DELS.GT.1.0E-4).AND.(N.LT.10))GO TO 15 C C ****************************END OF ITERATION LOOP ******** C C COMPUTE CONDUCTIVITY RATIO RTT = RT35(T)*RT*RT AT = A(T) BT = B(T) CP = C(P) CP = RTT*(CP + BT) BT = BT - RTT*AT C C SOLVE QUADRATIC EQUATION FOR R: R=RT35*RT*(1+C/AR+B) C R = SQRT(ABS(BT*BT + 4.0*AT*CP)) - BT C CONDUCTIVITY RETURN SAL78 = 0.5*R/AT ENDIF RETURN END REAL FUNCTION SVAN(S,T,P0,SIGMA) C MODIFIED RCM 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 ****************************************************************** C CHECK VALUE: SVAN=981.3021 E-8 M**3/KG. FOR S = 40 (IPSS-78) , C T = 40 DEG C, P0= 10000 DECIBARS. C CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78) , C T = 40 DEG C, P0= 10000 DECIBARS. C ******************************************************* REAL P,T,S,SIG,SR,R1,R2,R3,R4 REAL A,B,C,D,E,A1,B1,AW,BW,K,K0,KW,K35 C EQUIV EQUIVALENCE (E,D,B1),(BW,B,R3),(C,A1,R2) EQUIVALENCE (AW,A,R1),(KW,K0,K) C ******************** C DATA DATA R3500,R4/1028.1063,4.8314E-4/ DATA DR350/28.106331/ SAVE 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. SR = SQRT(ABS(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 X-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 X+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 SVAN=SVA*1.0E+8 IF(P.EQ.0.0) RETURN 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 E = (9.1697E-10*T+2.0816E-8)*T-9.9348E-7 BW = (5.2787E-8*T-6.12293E-6)*T+3.47718E-5 B = BW + E*S C D = 1.91075E-4 C = (-1.6078E-6*T-1.0981E-5)*T+2.2838E-3 AW = ((-5.77905E-7*T+1.16092E-4)*T+1.43713E-3)*T X-0.1194975 A = (D*SR + C)*S + AW C B1 = (-5.3009E-4*T+1.6483E-2)*T+7.944E-2 A1 = ((-6.1670E-5*T+1.09987E-2)*T-0.603459)*T+54.6746 KW = (((-5.155288E-5*T+1.360477E-2)*T-2.327105)*T X+148.4206)*T-1930.06 K0 = (B1*SR + A1)*S + KW 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 *********************************************** DK = (B*P + A)*P + K0 K35 = (5.03217E-5*P+3.359406)*P+21582.27 GAM=P/K35 PK = 1.0 - GAM SVA = SVA*PK + (V350P+SVA)*P*DK/(K35*(K35+DK)) C SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS SVAN=SVA*1.0E+8 V350P = V350P*PK 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 ******************************************************* DR35P=GAM/V350P DVAN=SVA/(V350P*(V350P+SVA)) SIGMA=DR350+DR35P-DVAN RETURN END REAL FUNCTION DEPTH(P,LAT) C ******************************** C DEPTH IN METERS FROM PRESSURE IN DECIBARS USING C SAUNDERS AND FOFONOFF'S METHOD. C DEEP-SEA RES., 1976,23,109-111. C FORMULA REFITTED FOR 1980 EQUATION OF STATE C UNITS: C PRESSURE P DECIBARS C LATITUDE LAT DEGREES C DEPTH DEPTH METERS C CHECKVALUE: DEPTH = 9712.653 M FOR P=10000 DECIBARS, LATITUDE=30 DEG C ABOVE FOR STANDARD OCEAN: T=0 DEG. CELSUIS ; S=35 (IPSS-78) C REAL LAT C X = SIN(LAT/57.29578) C************************** X = X*X C GR= GRAVITY VARIATION WITH LATITUDE: 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 RETURN END REAL FUNCTION TF(S,P) C FUNCTION TO COMPUTE THE FREEZING POINT OF SEAWATER C C REFERENCE: UNESCO TECH. PAPERS IN THE MARINE SCIENCE NO. 28. 1978 C EIGHTH REPORT JPOTS C ANNEX 6 FREEZING POINT OF SEAWATER F.J. MILLERO PP.29-35. C C UNITS: C PRESSURE P DECIBARS C SALINITY S PSS-78 C TEMPERATURE TF DEGREES CELSIUS C FREEZING PT. C************************************************************ C CHECKVALUE: TF= -2.588567 DEG. C FOR S=40.0, P=500. DECIBARS TF=(-.0575+1.710523E-3*SQRT(ABS(S))-2.154996E-4*S)*S-7.53E-4*P RETURN END REAL FUNCTION CPSW(S,T,P0) C **************************** C UNITS: C PRESSURE P0 DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C SALINITY S (IPSS-78) C SPECIFIC HEAT CPSW J/(KG DEG C) C******************************************************** C REF: MILLERO ET AL,1973,JGR,78,4499-4507 C MILLERO ET AL, UNESCO REPORT NO. 38 1981 PP. 99-188. C PRESSURE VARIATION FROM LEAST SQUARES POLYNOMIAL C DEVELOPED BY FOFONOFF 1980. C CHECK VALUE: CPSW = 3849.500 J/(KG DEG. C) FOR S = 40 (IPSS-78), C T = 40 DEG C, P0= 10000 DECIBARS C SCALE PRESSURE TO BARS P=P0/10. C************************** C SQRT SALINITY FOR FRACTIONAL TERMS SR = SQRT(ABS(S)) C SPECIFIC HEAT CP0 FOR P=0 (MILLERO ET AL ,UNESCO 1981) A = (-1.38385E-3*T+0.1072763)*T-7.643575 B = (5.148E-5*T-4.07718E-3)*T+0.1770383 C = (((2.093236E-5*T-2.654387E-3)*T+0.1412855)*T X -3.720283)*T+4217.4 CP0 = (B*SR + A)*S + C C CP1 PRESSURE AND TEMPERATURE TERMS FOR S = 0 A = (((1.7168E-8*T+2.0357E-6)*T-3.13885E-4)*T+1.45747E-2)*T X -0.49592 B = (((2.2956E-11*T-4.0027E-9)*T+2.87533E-7)*T-1.08645E-5)*T X +2.4931E-4 C = ((6.136E-13*T-6.5637E-11)*T+2.6380E-9)*T-5.422E-8 CP1 = ((C*P+B)*P+A)*P C CP2 PRESSURE AND TEMPERATURE TERMS FOR S > 0 A = (((-2.9179E-10*T+2.5941E-8)*T+9.802E-7)*T-1.28315E-4)*T X +4.9247E-3 B = (3.122E-8*T-1.517E-6)*T-1.2331E-4 A = (A+B*SR)*S B = ((1.8448E-11*T-2.3905E-9)*T+1.17054E-7)*T-2.9558E-6 B = (B+9.971E-8*SR)*S C = (3.513E-13*T-1.7682E-11)*T+5.540E-10 C = (C-1.4300E-12*T*SR)*S CP2 = ((C*P+B)*P+A)*P C SPECIFIC HEAT RETURN CPSW = CP0 + CP1 + CP2 RETURN END REAL FUNCTION ATG(S,T,P) C **************************** C ADIABATIC TEMPERATURE GRADIENT DEG C PER DECIBAR C REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408 C UNITS: C PRESSURE P DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C SALINITY S (IPSS-78) C ADIABATIC ATG DEG. C/DECIBAR C CHECKVALUE: ATG=3.255976E-4 C/DBAR FOR S=40 (IPSS-78), C T=40 DEG C,P0=10000 DECIBARS DS = S - 35.0 ATG = (((-2.1687E-16*T+1.8676E-14)*T-4.6206E-13)*P X+((2.7759E-12*T-1.1351E-10)*DS+((-5.4481E-14*T X+8.733E-12)*T-6.7795E-10)*T+1.8741E-8))*P X+(-4.2393E-8*T+1.8932E-6)*DS X+((6.6228E-10*T-6.836E-8)*T+8.5258E-6)*T+3.5803E-5 RETURN END REAL FUNCTION THETA_F(S,T0,P0,PR) 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 REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408 C FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491 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 P=P0 T=T0 C************** 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) THETA_F = T + (XK-2.0*Q)/6.0 RETURN END REAL FUNCTION SVEL(S,T,P0) C ******************************* C SOUND SPEED SEAWATER CHEN AND MILLERO 1977,JASA,62,1129-1135 C UNITS: C PRESSURE P0 DECIBARS C TEMPERATURE T DEG CELSIUS (IPTS-68) C SALINITY S (IPSS-78) C SOUND SPEED SVEL METERS/SECOND C CHECKVALUE: SVEL=1731.995 M/S, S=40 (IPSS-78),T=40 DEG C,P=10000 DBAR C EQUIVALENCE (A0,B0,C0),(A1,B1,C1),(A2,C2),(A3,C3) C C SCALE PRESSURE TO BARS P=P0/10. C************************** SR = SQRT(ABS(S)) C S**2 TERM D = 1.727E-3 - 7.9836E-6*P C S**3/2 TERM B1 = 7.3637E-5 +1.7945E-7*T B0 = -1.922E-2 -4.42E-5*T B = B0 + B1*P C S**1 TERM A3 = (-3.389E-13*T+6.649E-12)*T+1.100E-10 A2 = ((7.988E-12*T-1.6002E-10)*T+9.1041E-9)*T-3.9064E-7 A1 = (((-2.0122E-10*T+1.0507E-8)*T-6.4885E-8)*T-1.2580E-5)*T X +9.4742E-5 A0 = (((-3.21E-8*T+2.006E-6)*T+7.164E-5)*T-1.262E-2)*T X +1.389 A = ((A3*P+A2)*P+A1)*P+A0 C S**0 TERM C3 = (-2.3643E-12*T+3.8504E-10)*T-9.7729E-9 C2 = (((1.0405E-12*T-2.5335E-10)*T+2.5974E-8)*T-1.7107E-6)*T X +3.1260E-5 C1 = (((-6.1185E-10*T+1.3621E-7)*T-8.1788E-6)*T+6.8982E-4)*T X +0.153563 C0 = ((((3.1464E-9*T-1.47800E-6)*T+3.3420E-4)*T-5.80852E-2)*T X +5.03711)*T+1402.388 C = ((C3*P+C2)*P+C1)*P+C0 C SOUND SPEED RETURN SVEL = C + (A+B*SR+D*S)*S RETURN END