0
|
1 |
#======================================================================
|
|
2 |
# L I B E O S 8 3 . P L
|
|
3 |
# doc: Mon Mar 8 08:22:05 1999
|
|
4 |
# dlm: Thu Oct 6 11:44:00 2011
|
|
5 |
# (c) 1999 A.M. Thurnherr
|
|
6 |
# uE-Info: 1 0 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# Perl Implementation of UNESCO Eqn of State 1983
|
|
10 |
|
|
11 |
# Notes:
|
|
12 |
# - copied from eos83calc.c which was in turn...
|
|
13 |
# - copied from pexec_v4 (incl comments)
|
|
14 |
# - only subset of functions implemented
|
|
15 |
# - some attempt at cleaning up the code has been made
|
|
16 |
# - pressures in dbar throughout
|
|
17 |
# - no temperature scale assumed; set PARAM ITS=90|68
|
|
18 |
# - T90*1.00024=T68
|
|
19 |
# - check values calculated with T68
|
|
20 |
# - no conductivity unit assumed; set PARAM cond.unit=S/p|mS/cm
|
|
21 |
|
|
22 |
# HISTORY:
|
|
23 |
# Mar 08, 1999: - translated by hand from krc
|
|
24 |
# - added &alpha(), &beta(), &Rrho(), &TurnerAngle()
|
|
25 |
# Mar 13, 1999: - cosmetic changes
|
|
26 |
# Mar 14, 1999: - BUG NAN instead of NaN
|
|
27 |
# Mar 21, 1999: - make $sigmaR optional in &sVolAnom()
|
|
28 |
# Mar 31, 1999: - alias &potemp() = &theta(); &podens() = &sigma()
|
|
29 |
# - added &temp() to calc in-situ from potemp
|
|
30 |
# Sep 18, 1999: - parameter typechecking
|
|
31 |
# Aug 28, 2000: - added PARAM T68 check
|
|
32 |
# Aug 29, 2000: - forced temp_scale param check --- set to T68 or T90
|
|
33 |
# Sep 25, 2000: - changed temp_scale to 68 or 90 (easier check)
|
|
34 |
# - check for temp_scale during loading
|
|
35 |
# Nov 07, 2000: - added &dynHt()
|
|
36 |
# - strictified
|
|
37 |
# Nov 13, 2000: - removed temp_scale check during loading (STUPID,
|
|
38 |
# because PARAMs are only available after header
|
|
39 |
# is read)
|
|
40 |
# Feb 28, 2001: - added &rho(S,T,P)
|
|
41 |
# Mar 28, 2001: - changed Rrho() to use podens
|
|
42 |
# - optimized &theta() for P == Pref
|
|
43 |
# Apr 2, 2001: - &TurnerAngle() disabled pending adaptation to new Rrho()
|
|
44 |
# Apr 3, 2001: - updated &TurnerAngle()
|
|
45 |
# Jul 17, 2001: - cosmetics
|
|
46 |
# - added &grav(), &potErgAnom()
|
|
47 |
# Jul 18, 2001: - cosmetics
|
|
48 |
# Jul 20, 2001: - changed temp_scale to ITS
|
|
49 |
# Nov 26, 2001: - added &g(), &f()
|
|
50 |
# Dec 26, 2005: - update notes
|
|
51 |
# Jul 1, 2006: - Version 3.3 [HISTORY]
|
|
52 |
# - nan must be quoted for use strict
|
|
53 |
# Nov 5, 2006: - added K15toSalin
|
|
54 |
# Oct 18, 2007: - adapted to new &antsFunUsage()
|
|
55 |
# - changed order of &salin() params
|
|
56 |
# Oct 19, 2007: - continued
|
|
57 |
# - removed &grav(), podens(), potemp()
|
|
58 |
# Dec 1, 2007: - made theta(), sigma(), rho() return nan on nan input
|
|
59 |
# Dec 21, 2007: - BUG: grav() was still used
|
|
60 |
# Jan 20, 2008: - BUG: theta(), sigma(), rho() still generated error on nan in
|
|
61 |
# - made depth(), dynht(), potErgAnom() return nan on nan in
|
|
62 |
# Jan 4, 2011: - maded salin(), sVel() return nan on nan in
|
|
63 |
# Oct 6, 2011: - added %cond.unit (analogous to %ITS)
|
|
64 |
|
|
65 |
require "$ANTS/libvec.pl";
|
|
66 |
use strict;
|
|
67 |
|
|
68 |
#======================================================================
|
|
69 |
# PART I: stuff taken from PEXEC
|
|
70 |
#======================================================================
|
|
71 |
|
|
72 |
{ # BEGIN STATIC SCOPE
|
|
73 |
|
|
74 |
my($TCONV);
|
|
75 |
|
|
76 |
sub TCONV()
|
|
77 |
{
|
|
78 |
unless (defined($TCONV)) {
|
|
79 |
my($ITS) = &antsRequireParam('ITS');
|
|
80 |
if ($ITS == 68) {
|
|
81 |
$TCONV = 1;
|
|
82 |
} elsif ($ITS == 90) {
|
|
83 |
$TCONV = 1.00024;
|
|
84 |
} else {
|
|
85 |
croak("$0: illegal PARAM-value ITS=$ITS\n");
|
|
86 |
}
|
|
87 |
}
|
|
88 |
return $TCONV;
|
|
89 |
}
|
|
90 |
|
|
91 |
} # END STATIC SCOPE
|
|
92 |
|
|
93 |
#----------------------------------------------------------------------
|
|
94 |
# ADIABATIC TEMPERATURE GRADIENT DEG C/BAR
|
|
95 |
# REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
|
|
96 |
# CHECK VALUE: ATGR80=3.255976E-3 FOR S=40 NSU,T=40 DEG C,
|
|
97 |
# PIN=10000 DECIBARS
|
|
98 |
# NB: real check value appears to be 3.255976E-4; this makes sense
|
|
99 |
# as check value of potential temperature below is correct. Also
|
|
100 |
# note the 0.1 factor on return.
|
|
101 |
#----------------------------------------------------------------------
|
|
102 |
|
|
103 |
{ my(@fc);
|
|
104 |
sub adiaTempGrad(@)
|
|
105 |
{
|
|
106 |
my($S,$T,$P) = &antsFunUsage(3,'fff','[salin, temp, press(db)]',
|
|
107 |
\@fc,'salin','temp','press',@_);
|
|
108 |
my($DS);
|
|
109 |
my($TCONV) = &TCONV();
|
|
110 |
|
|
111 |
$T *= $TCONV; # use T68
|
|
112 |
$P *= 0.1;
|
|
113 |
$DS = $S - 35.0;
|
|
114 |
return 0.1 * ((((-2.1687E-13*$T+1.8676E-11)*$T-4.6206E-10)*$P
|
|
115 |
+((2.7759E-10*$T-1.1351E-8)*$DS+((-5.4481E-12*$T
|
|
116 |
+8.733E-10)*$T-6.7795E-8)*$T+1.8741E-6))*$P
|
|
117 |
+(-4.2393E-7*$T+1.8932E-5)*$DS
|
|
118 |
+((6.6228E-9*$T-6.836E-7)*$T+8.5258E-5)*$T+3.5803E-4) / $TCONV;
|
|
119 |
}
|
|
120 |
}
|
|
121 |
|
|
122 |
#----------------------------------------------------------------------
|
|
123 |
# TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR
|
|
124 |
# USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE
|
|
125 |
# AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.
|
|
126 |
# REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408
|
|
127 |
# FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491
|
|
128 |
# CHECK VALUE: PTMP83 =36.89072 FOR S=40 NSU,T=40 DEG C,
|
|
129 |
# P(MEASURED PRESSURE)=10000 DECIBARS,Pref=0 DECIBARS
|
|
130 |
#----------------------------------------------------------------------
|
|
131 |
|
|
132 |
{ my(@fc);
|
|
133 |
sub theta(@)
|
|
134 |
{
|
|
135 |
my($S,$T,$P,$Pref) =
|
|
136 |
&antsFunUsage(4,'....','[salin, temp, press(db),] refpress(db)',
|
|
137 |
\@fc,'salin','temp','press',undef,@_);
|
|
138 |
return 'nan' unless numberp($S) && numberp($T) && numberp($P);
|
|
139 |
my($H,$XK,$Q);
|
|
140 |
my($TCONV) = &TCONV();
|
|
141 |
|
|
142 |
return $T if ($P == $Pref);
|
|
143 |
$T *= $TCONV; # use T68
|
|
144 |
$H = $Pref - $P;
|
|
145 |
$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;
|
|
146 |
$T += 0.5 * $XK;
|
|
147 |
$Q = $XK;
|
|
148 |
$P += 0.5 * $H;
|
|
149 |
$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;
|
|
150 |
$T += 0.29289322 * ($XK-$Q);
|
|
151 |
$Q = 0.58578644*$XK + 0.121320344*$Q;
|
|
152 |
$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;
|
|
153 |
$T += 1.707106781 * ($XK-$Q);
|
|
154 |
$Q = 3.414213562*$XK - 4.121320344*$Q;
|
|
155 |
$P += 0.5*$H;
|
|
156 |
$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;
|
|
157 |
return ($T + ($XK-2.0*$Q)/6.0)/$TCONV;
|
|
158 |
}
|
|
159 |
}
|
|
160 |
|
|
161 |
# can't easily do default fields because theta field name is not unique
|
|
162 |
sub temp(@)
|
|
163 |
{
|
|
164 |
my($S,$T,$P,$Pref) =
|
|
165 |
&antsFunUsage(4,"ffff","salin, potemp, press(db), refpress(db)",@_);
|
|
166 |
return &theta($S,$T,$Pref,$P);
|
|
167 |
}
|
|
168 |
|
|
169 |
#----------------------------------------------------------------------
|
|
170 |
# Calculation of specific volume anomaly
|
|
171 |
# Gill (1982), p.215: - specific volume = inverse of density
|
|
172 |
# - SVA referenced to same press, T=0, S=35
|
|
173 |
#
|
|
174 |
# *************************************************
|
|
175 |
# *** USES EQUATION OF STATE FOR SEA WATER 1980 ***
|
|
176 |
# *************************************************
|
|
177 |
#
|
|
178 |
# Check value. SVAN83=981.30210E-8, SIGMA=59.82037
|
|
179 |
# for S=40, T=40, P=10000
|
|
180 |
#
|
|
181 |
# NB: - this one was a bitch to translate, used f2c and it shows :-)
|
|
182 |
# - the 1e-8 factor in the check value is wrong; yup: see gill p.215
|
|
183 |
# (typical units are 1e-8m^3kg^-1)
|
|
184 |
# - check val appears to be slightly off, I get 981.30187319561
|
|
185 |
#
|
|
186 |
#----------------------------------------------------------------------
|
|
187 |
|
|
188 |
sub sVolAnom(@)
|
|
189 |
{
|
|
190 |
my($S,$T,$P,$sigmaR) =
|
|
191 |
&antsFunUsage(-3,'fff','salin, temp, press(db)[, ref to sigma var]',@_);
|
|
192 |
my($r3500) = 1028.1063;
|
|
193 |
my($r4) = 4.8314e-4;
|
|
194 |
my($dr350) = 28.106331;
|
|
195 |
|
|
196 |
my($temp0,$temp1,$temp2,$temp3,$temp4);
|
|
197 |
my($dvan,$dr35p,$p,$t,$dk,$pk,$sr,$gam,$sig,$rk35,$sva,$v350p);
|
|
198 |
|
|
199 |
my($TCONV) = &TCONV();
|
|
200 |
|
|
201 |
$t = $T*$TCONV; # use T68
|
|
202 |
$p = $P / 10.;
|
|
203 |
# $sr = sqrt(abs($S));
|
|
204 |
$sr = sqrt($S);
|
|
205 |
$temp3 = (((($t * 6.536332e-9 - 1.120083e-6) * $t +
|
|
206 |
1.001685e-4) * $t - .00909529) * $t + .06793952) * $t
|
|
207 |
- 28.263737;
|
|
208 |
$temp2 = ((($t * 5.3875e-9 - 8.2467e-7) * $t + 7.6438e-5)
|
|
209 |
* $t - .0040899) * $t + .824493;
|
|
210 |
$temp1 = ($t * -1.6546e-6 + 1.0227e-4) * $t - .00572466;
|
|
211 |
$sig = ($r4 * $S + $temp1 * $sr + $temp2) * $S + $temp3;
|
|
212 |
$v350p = 1. / $r3500;
|
|
213 |
$sva = -$sig * $v350p / ($r3500 + $sig);
|
|
214 |
$$sigmaR = $sig + $dr350 if (defined($sigmaR));
|
|
215 |
return $sva * 1e8 if ($p == 0.0);
|
|
216 |
|
|
217 |
$temp0 = ($t * 9.1697e-10 + 2.0816e-8) * $t - 9.9348e-7;
|
|
218 |
$temp1 = ($t * 5.2787e-8 - 6.12293e-6) * $t + 3.47718e-5;
|
|
219 |
$temp1 = $temp1 + $temp0 * $S;
|
|
220 |
$temp0 = 1.91075e-4;
|
|
221 |
$temp2 = ($t * -1.6078e-6 - 1.0981e-5) * $t + .0022838;
|
|
222 |
$temp3 = (($t * -5.77905e-7 + 1.16092e-4) * $t + .00143713) * $t - .1194975;
|
|
223 |
$temp3 = ($temp0 * $sr + $temp2) * $S + $temp3;
|
|
224 |
$temp0 = ($t * -5.3009e-4 + .016483) * $t + .07944;
|
|
225 |
$temp2 = (($t * -6.167e-5 + .0109987) * $t - .603459) * $t + 54.6746;
|
|
226 |
$temp4 = ((($t * -5.155288e-5 + .01360477) * $t - 2.327105) * $t + 148.4206) * $t - 1930.06;
|
|
227 |
$temp4 = ($temp0 * $sr + $temp2) * $S + $temp4;
|
|
228 |
$dk = ($temp1 * $p + $temp3) * $p + $temp4;
|
|
229 |
$rk35 = ($p * 5.03217e-5 + 3.359406) * $p + 21582.27;
|
|
230 |
$gam = $p / $rk35;
|
|
231 |
$pk = 1. - $gam;
|
|
232 |
$sva = $sva * $pk + ($v350p + $sva) * $p * $dk / ($rk35 * ($rk35 + $dk));
|
|
233 |
$v350p *= $pk;
|
|
234 |
$dr35p = $gam / $v350p;
|
|
235 |
$dvan = $sva / ($v350p * ($v350p + $sva));
|
|
236 |
$$sigmaR = $dr350 + $dr35p - $dvan if (defined($sigmaR));
|
|
237 |
return $sva * 1e8;
|
|
238 |
}
|
|
239 |
|
|
240 |
#----------------------------------------------------------------------
|
|
241 |
# Dynamic Height (from dyht83.F)
|
|
242 |
# Usage:
|
|
243 |
# - dynHt(salin,temp,press[,idx])
|
|
244 |
# Check Values from PEXEC:
|
|
245 |
# - dynHt(40,18,20) -> -.01879
|
|
246 |
# - dynHt(??,16,50) -> -.03194
|
|
247 |
# - dynHt(??,14,100) -> -.00255
|
|
248 |
# NB: - check values 2 & 3 do not appear to check out
|
|
249 |
# - checked against WHOI-supplied BBTRE data indicates max
|
|
250 |
# diff of 1e-3 dyn.m at 5141m
|
|
251 |
# - use different idx to use multiple times in single program
|
|
252 |
# (c.f. [gshear])
|
|
253 |
# - behaves the same as the pexec version with PFIRST=0
|
|
254 |
# - dynamic height seems to be defined (1e-5 factor) to allow
|
|
255 |
# station distance to be in km and velocities in cm/s
|
|
256 |
#----------------------------------------------------------------------
|
|
257 |
|
|
258 |
{ # BEGIN static scope
|
|
259 |
|
|
260 |
my(@lastP,@lastAnom,@ht);
|
|
261 |
|
|
262 |
sub dynHt(@)
|
|
263 |
{
|
|
264 |
my($S,$T,$P,$idx) =
|
|
265 |
&antsFunUsage(-3,"...","salin, temp, press(db)[, <idx>]",@_);
|
|
266 |
return 'nan' unless numbersp($S,$T,$P);
|
|
267 |
my($anom) = sVolAnom($S,$T,$P) * 1e-5;
|
|
268 |
|
|
269 |
if (!defined($ht[$idx])) { # first call
|
|
270 |
$ht[$idx] = $P * $anom;
|
|
271 |
} else { # successive calls
|
|
272 |
croak("$0: pressure not increasing monotonically ($lastP[$idx] -> $P)\n")
|
|
273 |
if ($P < $lastP[$idx]);
|
|
274 |
$ht[$idx] += ($P-$lastP[$idx]) * ($anom+$lastAnom[$idx])/2;
|
|
275 |
}
|
|
276 |
$lastP[$idx] = $P; $lastAnom[$idx] = $anom;
|
|
277 |
|
|
278 |
return $ht[$idx];
|
|
279 |
}
|
|
280 |
|
|
281 |
} # END static scope
|
|
282 |
|
|
283 |
#----------------------------------------------------------------------
|
|
284 |
# Local Gravity (from grav83.F)
|
|
285 |
# Usage:
|
|
286 |
# - g(press,lat)
|
|
287 |
# Check Values:
|
|
288 |
# - g(10000,30) = 9.804160
|
|
289 |
#----------------------------------------------------------------------
|
|
290 |
|
|
291 |
{ my(@fc);
|
|
292 |
sub g(@)
|
|
293 |
{
|
|
294 |
my($P,$lat) = &antsFunUsage(2,'ff','[press(db), lat]',\@fc,'press','%lat',@_);
|
|
295 |
|
|
296 |
my($x) = sin($lat/57.29578) ** 2;
|
|
297 |
my($g) = 9.780318*(1.0+(5.2788e-3+2.36e-5*$x)*$x); # global variation
|
|
298 |
return $g + 1.092e-6 * $P; # pressure correction
|
|
299 |
}
|
|
300 |
}
|
|
301 |
|
|
302 |
#----------------------------------------------------------------------
|
|
303 |
# Coriolis frequency
|
|
304 |
# Usage:
|
|
305 |
# - f(%lat)
|
|
306 |
#----------------------------------------------------------------------
|
|
307 |
|
|
308 |
{ my(@fc);
|
|
309 |
sub f(@)
|
|
310 |
{
|
|
311 |
my($lat) = &antsFunUsage(1,'f','[lat]',\@fc,'%lat',@_);
|
|
312 |
my($Omega) = 7.292e-5; # Gill (1982)
|
|
313 |
return 2 * $Omega * sin(rad($lat));
|
|
314 |
}
|
|
315 |
}
|
|
316 |
|
|
317 |
#----------------------------------------------------------------------
|
|
318 |
# Potential Energy Anomaly (from pean83.F)
|
|
319 |
# Usage:
|
|
320 |
# - potErgAnom(salin,temp,press,refPress,lat)
|
|
321 |
# Check Values from PEXEC: (NB: sequence of calls)
|
|
322 |
# - potErgAnom(40,18, 20,0,49.2235) = -1.91462
|
|
323 |
# - potErgAnom(38,16, 50,0,49.2235) = -4.31092
|
|
324 |
# - potErgAnom(36,14,100,0,49.2235) = 24.82922
|
|
325 |
# Units:
|
|
326 |
# - Gill (1982) p.45: m^2/s^2 = J/kg
|
|
327 |
# NB:
|
|
328 |
# - geopotential Phi is PEA wrt zero pressure
|
|
329 |
# - PE of volume is volume integral of rho*PEA (c.f. Gill, p.80)
|
|
330 |
# - check value lat calculated so that g(0,lat) ~ 9.81
|
|
331 |
# - IF Pref<P ON 1ST CALL, HEIGHT WILL BE CALCULATED FROM LEVEL
|
|
332 |
# Pref. HENCE IF Pref = 0.0, FROM SEA SURFACE DOWN.
|
|
333 |
# - different to pexec, this version sets values for P<Pref to nan
|
|
334 |
#----------------------------------------------------------------------
|
|
335 |
|
|
336 |
{ # BEGIN STATIC SCOPE
|
|
337 |
|
|
338 |
my($gzero,$H,$lastZ,$lastP);
|
|
339 |
|
|
340 |
sub potErgAnom(@)
|
|
341 |
{
|
|
342 |
my($S,$T,$P,$Pref,$lat) =
|
|
343 |
&antsFunUsage(5,".....","salin, temp, press(db), refpress(db), lat",@_);
|
|
344 |
return 'nan' unless numbersp($S,$T,$P,$Pref,$lat);
|
|
345 |
|
|
346 |
return 'nan' if ($P < $Pref);
|
|
347 |
|
|
348 |
$gzero = g(0,$lat) unless (defined($gzero)); # 1st time
|
|
349 |
my($g) = $gzero + 1.113e-4*$P;
|
|
350 |
|
|
351 |
my($anom) = sVolAnom($S,$T,$P) * 1e-3;
|
|
352 |
my($Z) = $anom * $P/$g;
|
|
353 |
|
|
354 |
if (defined($H)) { # not 1st time
|
|
355 |
croak("$0: pressure not increasing monotonically\n")
|
|
356 |
if ($P < $lastP);
|
|
357 |
$H += ($Z+$lastZ) * ($P-$lastP)*0.5;
|
|
358 |
} else {
|
|
359 |
$H = (($anom*$Pref/$g)+$Z)*($P-$Pref)*0.5;
|
|
360 |
}
|
|
361 |
|
|
362 |
$lastP = $P; $lastZ = $Z;
|
|
363 |
return $H;
|
|
364 |
}
|
|
365 |
|
|
366 |
} # END static scope
|
|
367 |
|
|
368 |
#----------------------------------------------------------------------
|
|
369 |
# Density at pressure P (P - Measured pressure)
|
|
370 |
# (PREF - Reference pressure)
|
|
371 |
# Equation of state for seawater proposed by JPOTS 1980
|
|
372 |
# References
|
|
373 |
# Millero et al 1980, Deep Sea Res.,27A,255-264
|
|
374 |
# Jpots Ninth Report 1978, Tenth Report 1980
|
|
375 |
# Units:
|
|
376 |
# Pressure P Decibars
|
|
377 |
# Temperature T Deg Celcius (IPTS-68)
|
|
378 |
# Salinity S NSU (IPSS-78)
|
|
379 |
# Density RHO KG/M**3
|
|
380 |
# Spec. Vol. EOS80 M**3/KG
|
|
381 |
# check value. 43.331642 for P=10000,PREF=5000,T=40,S=40
|
|
382 |
# NB: check value appears to be 42.33164!!!
|
|
383 |
#----------------------------------------------------------------------
|
|
384 |
|
|
385 |
{ my(@fc);
|
|
386 |
sub sigma(@)
|
|
387 |
{
|
|
388 |
my($S,$T,$P,$Pref) =
|
|
389 |
&antsFunUsage(4,'....','[salin, temp, press(db),] refpress(db)',
|
|
390 |
\@fc,'salin','temp','press',undef,@_);
|
|
391 |
return 'nan' unless numberp($S) && numberp($T) && numberp($P);
|
|
392 |
my($sig);
|
|
393 |
&sVolAnom($S,&theta($S,$T,$P,$Pref),$Pref,\$sig);
|
|
394 |
return $sig;
|
|
395 |
}
|
|
396 |
}
|
|
397 |
|
|
398 |
{ my(@fc);
|
|
399 |
sub rho(@)
|
|
400 |
{
|
|
401 |
my($S,$T,$P) = &antsFunUsage(3,'...','[salin, temp, press(db)]',
|
|
402 |
\@fc,'salin','temp','press',@_);
|
|
403 |
return 'nan' unless numberp($S) && numberp($T) && numberp($P);
|
|
404 |
return 1000 + &sigma($S,$T,$P,$P);
|
|
405 |
}
|
|
406 |
}
|
|
407 |
|
|
408 |
#----------------------------------------------------------------------
|
|
409 |
# FUNCTION TO CONVERT CONDUCTIVITY TO SALINITY ACCORDING TO THE
|
|
410 |
# ALGORITHMS RECOMMMENDED BY JPOTS USING THE 1978 PRACTICAL
|
|
411 |
# SALINITY SCALE (IPSS-78) AND IPTS-68 FOR TEMPERATURE.
|
|
412 |
# C1535=CONDUCTIVITY (FROM CULKIN&SMITH,1980,J.OCEAN.ENG.VOL.5
|
|
413 |
# PP22-23) = 42.9140 AT PRES=0.0,
|
|
414 |
# SALINITY 35 NSU AND TEMPERATURE 15 DEG CELSIUS (IPTS-68)
|
|
415 |
# PRESSURE=DECIBARS
|
|
416 |
# RETURNS ZERO FOR CND<.0005
|
|
417 |
# CHECKVALUES....
|
|
418 |
# SAL83 =40.000000 FOR CND=81.025545,T=40 DEG C,PIN=10000
|
|
419 |
# DECIBARS
|
|
420 |
# WRITTEN BY N FOFONOFF; REVISED OCT 6 1980
|
|
421 |
# FCN SAL83, XR=SQRT(RT)
|
|
422 |
# DERIVATIVE WRT XR ; DSAL/DXR
|
|
423 |
# RT35
|
|
424 |
# C,B,A, POLYNOMIALS
|
|
425 |
# NB: - done using f2c (couldn't be bothered)
|
|
426 |
# - removed SAL->COND
|
|
427 |
#----------------------------------------------------------------------
|
|
428 |
|
|
429 |
{ my(@fc);
|
|
430 |
my($cond_scale); # 1 for mS/cm, 10 for S/m
|
|
431 |
sub salin(@)
|
|
432 |
{
|
|
433 |
my($C,$T,$P) = &antsFunUsage(3,'...','[cond, temp, press(db)]',
|
|
434 |
\@fc,'cond','temp','press',@_);
|
|
435 |
return 'nan' unless numberp($C) && numberp($T) && numberp($P);
|
|
436 |
my($r__,$dt,$rt,$c1535);
|
|
437 |
my($TCONV) = &TCONV();
|
|
438 |
|
|
439 |
unless (defined($cond_scale)) { # deal with different conductivity units
|
|
440 |
my($cu) = &antsRequireParam('cond.unit');
|
|
441 |
if ($cu eq 'S/m') { $cond_scale = 10; }
|
|
442 |
elsif ($cu eq 'mS/cm') { $cond_scale = 1; }
|
|
443 |
else { croak("$0: illegal PARAM-value cond.unit=$cu\n"); }
|
|
444 |
}
|
|
445 |
$C *= $cond_scale;
|
|
446 |
|
|
447 |
return 0.0 if ($C <= 5e-4); # zero salinity trap
|
|
448 |
$T *= $TCONV; # use T68 scale
|
|
449 |
$c1535 = 42.914;
|
|
450 |
$dt = $T - 15.0;
|
|
451 |
$P *= .1; # convert pressure to bars
|
|
452 |
$r__ = $C / $c1535; # convert cond to salin
|
|
453 |
$rt = $r__ / ((((($T * 1.0031e-9 - 6.9698e-7) * $T +
|
|
454 |
1.104259e-4) * $T + .0200564) * $T + .6766097) * ((($P
|
|
455 |
* 3.989e-12 - 6.37e-8) * $P + 2.07e-4) * $P / (
|
|
456 |
($T * 4.464e-4 + .03426) * $T + 1. + ($T *
|
|
457 |
-.003107 + .4215) * $r__) + 1.));
|
|
458 |
$rt = sqrt(abs($rt));
|
|
459 |
return (((($rt * 2.7081 - 7.0261) * $rt + 14.0941) *
|
|
460 |
$rt + 25.3851) * $rt - .1692) * $rt + .008 +
|
|
461 |
$dt / ($dt * .0162 + 1.) * ((((($rt * -.0144 +
|
|
462 |
.0636) * $rt - .0375) * $rt - .0066) * $rt -
|
|
463 |
.0056) * $rt + 5e-4);
|
|
464 |
}
|
|
465 |
}
|
|
466 |
|
|
467 |
#----------------------------------------------------------------------
|
|
468 |
# K15toSalin; see http://ioc.unesco.org/oceanteacher/OceanTeacher2/01_GlobOcToday/02_CollDta/02_OcDtaFunda/03_T&SScales/TemperatureAndSalinityScales.htm
|
|
469 |
#----------------------------------------------------------------------
|
|
470 |
|
|
471 |
sub K15toSalin(@)
|
|
472 |
{
|
|
473 |
my($K15) = &antsFunUsage(1,'f','K15',@_);
|
|
474 |
return 0.0080 - 0.1692 * $K15**0.5
|
|
475 |
+ 25.3851 * $K15**1.0
|
|
476 |
+ 14.0941 * $K15**1.5
|
|
477 |
- 7.0261 * $K15**2.0
|
|
478 |
+ 2.7081 * $K15**2.5;
|
|
479 |
}
|
|
480 |
|
|
481 |
#----------------------------------------------------------------------
|
|
482 |
# DEPTH IN METERS FROM PRESSURE IN DECIBARS USING
|
|
483 |
# SAUNDERS AND FOFONOFF'S METHOD.
|
|
484 |
# DEEP SEA RES., 1976,23,109-111.
|
|
485 |
# FORMULA REFITTED FOR EOS80
|
|
486 |
# CHECK VALUE. 9712.654 M FOR PIN=10000 DECIBARS,LATITUDE=30 DEG
|
|
487 |
# ..CONVERT PRESSURE TO BARS
|
|
488 |
# NB: results are closer to Saunders, JPO 11, 1981 than ref given
|
|
489 |
# above. They are not exact, either, however.
|
|
490 |
#----------------------------------------------------------------------
|
|
491 |
|
|
492 |
{ my(@fc);
|
|
493 |
sub depth(@)
|
|
494 |
{
|
|
495 |
my($P,$lat) = &antsFunUsage(2,'..','[press(db), lat]',\@fc,'press','%lat',@_);
|
|
496 |
return 'nan' unless numbersp($P,$lat);
|
|
497 |
my($x) = sin($lat/57.29578);
|
|
498 |
$P *= 0.1;
|
|
499 |
$x = $x*$x;
|
|
500 |
return ((((-1.82E-11*$P+2.279E-7)*$P-2.2512E-3)*$P+97.2659)*$P) /
|
|
501 |
(9.780318*(1.0+(5.2788E-3+2.36E-5*$x)*$x) + 1.092E-5*$P);
|
|
502 |
}
|
|
503 |
}
|
|
504 |
|
|
505 |
#----------------------------------------------------------------------
|
|
506 |
# Convert depth to pressure using 2nd order polynomial
|
|
507 |
# From file pdepth.F (pexec)
|
|
508 |
# This is based on Saunders, JPO 11, 1981
|
|
509 |
#----------------------------------------------------------------------
|
|
510 |
|
|
511 |
{ my(@fc);
|
|
512 |
sub press(@)
|
|
513 |
{
|
|
514 |
my($d,$lat) = &antsFunUsage(2,'ff','[depth, lat]',\@fc,'depth','%lat',@_);
|
|
515 |
my($c1) = 5.92E-3 + 5.25E-3 * sin($lat/57.29578) ** 2;
|
|
516 |
my($c2) = 2.21E-6;
|
|
517 |
my($press) = (1 - $c1 - sqrt((1 - $c1)**2 - 4*$d*$c2)) / (2*$c2);
|
|
518 |
my($derr) = abs(&depth($press,$lat) - $d);
|
|
519 |
|
|
520 |
&antsInfo("WARNING (libEOS83.pl): %.1gm depth error due to pressure " .
|
|
521 |
"approximation", $derr) if ($derr >= 1);
|
|
522 |
return $press;
|
|
523 |
}
|
|
524 |
}
|
|
525 |
|
|
526 |
#----------------------------------------------------------------------
|
|
527 |
# SOUND SPEED SEAWATER (CHEN & MILLERO 1977, JASA,62,1129-1135)
|
|
528 |
# SPEED IN M/S, P IN DECIBARS, T DEG C (IPTS-68), S NSU(IPSS-78)
|
|
529 |
# CHECK VALUE : 1731.9954 M/S FOR PIN=10000, T=40, C,S=40
|
|
530 |
# NB: - f2c used (because of fortran `equivalence')
|
|
531 |
#----------------------------------------------------------------------
|
|
532 |
|
|
533 |
{ my(@fc);
|
|
534 |
sub sVel(@)
|
|
535 |
{
|
|
536 |
my($S,$T,$P) = &antsFunUsage(3,'...','[salin, temp, press(db)]',
|
|
537 |
\@fc,'salin','temp','press',@_);
|
|
538 |
return 'nan' unless numberp($S) && numberp($T) && numberp($P);
|
|
539 |
my($temp0,$temp1,$temp2,$temp3);
|
|
540 |
my($a,$b,$c__,$d__,$sr);
|
|
541 |
my($TCONV) = &TCONV();
|
|
542 |
|
|
543 |
$T *= $TCONV; # T68
|
|
544 |
$P *= .1; # CONVERT PRESSURE TO BARS
|
|
545 |
$sr = sqrt(abs($S)); # S**2 TERM
|
|
546 |
$d__ = .001727 - $P * 7.9836e-6;
|
|
547 |
$temp1 = $T * 1.7945e-7 + 7.3637e-5; # S**3/2 TERM
|
|
548 |
$temp0 = -.01922 - $T * 4.42e-5;
|
|
549 |
$b = $temp0 + $temp1 * $P;
|
|
550 |
$temp3 = ($T * -3.389e-13 + 6.649e-12) * $T + 1.1e-10; # S**1 TERM
|
|
551 |
$temp2 = (($T * 7.988e-12 - 1.6002e-10) * $T + 9.1041e-9) * $T - 3.9064e-7;
|
|
552 |
$temp1 = ((($T * -2.0122e-10 + 1.0507e-8) * $T - 6.4885e-8) * $T
|
|
553 |
- 1.258e-5) * $T + 9.4742e-5;
|
|
554 |
$temp0 = ((($T * -3.21e-8 + 2.006e-6) * $T + 7.164e-5) * $T - .01262)
|
|
555 |
* $T + 1.389;
|
|
556 |
$a = (($temp3 * $P + $temp2) * $P + $temp1) * $P + $temp0;
|
|
557 |
$temp3 = ($T * -2.3643e-12 + 3.8504e-10) * $T - 9.7729e-9; # S**0 TERM
|
|
558 |
$temp2 = ((($T * 1.0405e-12 - 2.5335e-10) * $T + 2.5974e-8) * $T
|
|
559 |
- 1.7107e-6) * $T + 3.126e-5;
|
|
560 |
$temp1 = ((($T * -6.1185e-10 + 1.3621e-7) * $T - 8.1788e-6) * $T
|
|
561 |
+ 6.8982e-4) * $T + .153563;
|
|
562 |
$temp0 = (((($T * 3.1464e-9 - 1.478e-6) * $T + 3.342e-4) * $T - .0580852)
|
|
563 |
* $T + 5.03711) * $T + 1402.388;
|
|
564 |
$c__ = (($temp3 * $P + $temp2) * $P + $temp1) * $P + $temp0;
|
|
565 |
return $c__ + ($a + $b * $sr + $d__ * $S) * $S; # SOUND SPEED RETURN
|
|
566 |
}
|
|
567 |
}
|
|
568 |
|
|
569 |
#======================================================================
|
|
570 |
# PART II: Homegrown Stuff
|
|
571 |
#======================================================================
|
|
572 |
|
|
573 |
#----------------------------------------------------------------------
|
|
574 |
# &alpha(S,T,P) linear thermal expansion coefficient
|
|
575 |
# Notes: - use temperature interval of 0.2 degrees
|
|
576 |
# - depth instead of pressure ok
|
|
577 |
# Check Value: - alpha(35,2,6000) = 0.000209447279306853
|
|
578 |
#----------------------------------------------------------------------
|
|
579 |
|
|
580 |
sub alpha(@)
|
|
581 |
{
|
|
582 |
my($S,$T,$P) = &antsFunUsage(3,"fff","salin, temp, press | depth",@_);
|
|
583 |
return (&sigma($S,$T-.1,$P,$P) - &sigma($S,$T+.1,$P,$P))
|
|
584 |
/ (.2 * (1000+&sigma($S,$T,$P,$P)));
|
|
585 |
}
|
|
586 |
|
|
587 |
#----------------------------------------------------------------------
|
|
588 |
# &beta(S,T,P) linear haline contraction coefficient
|
|
589 |
# Notes: - use salinity interval of 0.02 psu
|
|
590 |
# - depth instead of pressure ok
|
|
591 |
# Check Value: - beta(35,2,6000) = 0.000718513652714652 (should check Gill)
|
|
592 |
#----------------------------------------------------------------------
|
|
593 |
|
|
594 |
sub beta(@)
|
|
595 |
{
|
|
596 |
my($S,$T,$P) = &antsFunUsage(3,"fff","salin, temp, press | depth",@_);
|
|
597 |
return (&sigma($S+0.01,$T,$P,$P) - &sigma($S-0.01,$T,$P,$P))
|
|
598 |
/ (.02 * (1000+&sigma($S,$T,$P,$P)));
|
|
599 |
}
|
|
600 |
|
|
601 |
#----------------------------------------------------------------------
|
|
602 |
# &Rrho(midS,S0,S1,midT,T0,T1,midP,P0,P1) stability ratio
|
|
603 |
# Notes: - depth instead of pressure ok
|
|
604 |
# Check Value: - Rrho(35.05,35.1,35,4.27,4.82,3.72,2100,2100,2100) = 2.222
|
|
605 |
#----------------------------------------------------------------------
|
|
606 |
|
|
607 |
sub alphaDT(@)
|
|
608 |
{
|
|
609 |
my($S,$T0,$T1,$P,$P0,$P1) = @_;
|
|
610 |
# &antsFunUsage(6,"ffffff","S, T0, T1, P, P0, P1",@_);
|
|
611 |
my($sgn) = ($P1 > $P0) ? 1 : -1;
|
|
612 |
return $sgn * (&sigma($S,$T1,$P1,$P) - &sigma($S,$T0,$P0,$P));
|
|
613 |
}
|
|
614 |
|
|
615 |
sub betaDS(@)
|
|
616 |
{
|
|
617 |
my($S0,$S1,$T,$P,$P0,$P1) = @_;
|
|
618 |
# &antsFunUsage(6,"ffffff","S0, S1, T, P, P0, P1",@_);
|
|
619 |
my($sgn) = ($P1 > $P0) ? 1 : -1;
|
|
620 |
return $sgn * (&sigma($S0,$T,$P0,$P) - &sigma($S1,$T,$P1,$P));
|
|
621 |
}
|
|
622 |
|
|
623 |
sub Rrho(@)
|
|
624 |
{
|
|
625 |
my($S,$S0,$S1,$T,$T0,$T1,$P,$P0,$P1) =
|
|
626 |
&antsFunUsage(9,"fffffffff","S, S0, S1, T, T0, T1, P, P0, P1",@_);
|
|
627 |
my($aDT) = &alphaDT($S,$T0,$T1,$P,$P0,$P1);
|
|
628 |
my($bDS) = &betaDS($S0,$S1,$T,$P,$P0,$P1);
|
|
629 |
return $bDS == 0 ? 'nan' : $aDT / $bDS;
|
|
630 |
}
|
|
631 |
|
|
632 |
#----------------------------------------------------------------------
|
|
633 |
# &TurnerAngle(midS,S0,S1,midT,T0,T1,midP,P0,P1) Turner Angle
|
|
634 |
# Notes: - c.f. &Rrho()
|
|
635 |
# - -45<Tu<45 doubly stable
|
|
636 |
# - Tu<-90, Tu>90 statically unstable
|
|
637 |
# - 45<Tu<90 fingering regime
|
|
638 |
# - -90<Tu<-45 diffusive regime
|
|
639 |
# Check Value: TurnerAngle(35.05,35.1,35,4.27,4.82,3.72,2100,2100,2100.01) = 69.2296333539803
|
|
640 |
#----------------------------------------------------------------------
|
|
641 |
|
|
642 |
sub TurnerAngle(@)
|
|
643 |
{
|
|
644 |
my($S,$S0,$S1,$T,$T0,$T1,$P,$P0,$P1) =
|
|
645 |
&antsFunUsage(9,"fffffffff","S, S0, S1, T, T0, T1, P, P0, P1",@_);
|
|
646 |
my($aDT) = &alphaDT($S,$T0,$T1,$P,$P0,$P1);
|
|
647 |
my($bDS) = &betaDS($S0,$S1,$T,$P,$P0,$P1);
|
|
648 |
return atan2($aDT+$bDS,$aDT-$bDS)*57.29577951;
|
|
649 |
}
|
|
650 |
|
|
651 |
1;
|
|
652 |
|