1 | function [rho,bvf]=rho_eos(Tt,Ts,z_r,z_w,g,rho0) |
---|
2 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
3 | % |
---|
4 | % function [rho,bvf]=rho_eos(Tt,Ts,z_r,z_w,g,rho0) |
---|
5 | % |
---|
6 | % Computes density via Equation Of State (EOS) for seawater. |
---|
7 | % If so prescribed, non-linear EOS of Jackett and McDougall (1995) |
---|
8 | % is used. |
---|
9 | % |
---|
10 | % Tt potential temperature [deg Celsius]. |
---|
11 | % Ts salinity [PSU]. |
---|
12 | % Tz pressure/depth, [depth in meters and negative]. |
---|
13 | % |
---|
14 | % K0, K1 and K2 are the pressure polynomial coefficients for secant |
---|
15 | % bulk modulus, so that |
---|
16 | % |
---|
17 | % bulk = K0 - K1 * z + K2 * z**2 ; |
---|
18 | % |
---|
19 | % while rho1 is sea-water density [kg/m^3] at standard pressure |
---|
20 | % of 1 Atm, so that the density anomaly at in-sity pressure is |
---|
21 | % |
---|
22 | % rho = rho1 / (1 + z / bulk) - 1000 |
---|
23 | % |
---|
24 | % If so prescribed, it also computes the Brunt-Vaisala frequency |
---|
25 | % [1/s^2] at horizontal RHO-points and vertical W-points, |
---|
26 | % |
---|
27 | % bvf = - g/rho0 d(rho)/d(z). |
---|
28 | % |
---|
29 | % In computation of bvf the density anomaly difference is computed |
---|
30 | % by adiabatically lowering/rising the water parcel from RHO point |
---|
31 | % above/below to the W-point depth at "z_w". |
---|
32 | % |
---|
33 | % Reference: |
---|
34 | % |
---|
35 | % Jackett, D. R. and T. J. McDougall, 1995, Minimal Adjustment of |
---|
36 | % Hydrostatic Profiles to Achieve Static Stability, Journ of Atmos. |
---|
37 | % and Oceanic Techn., vol. 12, pp. 381-389. |
---|
38 | % |
---|
39 | % << This equation of state formulation has been derived by Jackett |
---|
40 | % and McDougall (1992), unpublished manuscript, CSIRO, Australia. |
---|
41 | % It computes in-situ density anomaly as a function of potential |
---|
42 | % temperature (Celsius) relative to the surface, salinity (PSU), |
---|
43 | % and depth (meters). It assumes no pressure variation along |
---|
44 | % geopotential surfaces, that is, depth and pressure are |
---|
45 | % interchangeable. >> |
---|
46 | % John Wilkin, 29 July 92 |
---|
47 | % |
---|
48 | % Check Values: T=3 C S=35.5 PSU Z=-5000 m rho=1050.3639165364 |
---|
49 | % |
---|
50 | % |
---|
51 | % Further Information: |
---|
52 | % http://www.brest.ird.fr/Roms_tools/ |
---|
53 | % |
---|
54 | % This file is part of ROMSTOOLS |
---|
55 | % |
---|
56 | % ROMSTOOLS is free software; you can redistribute it and/or modify |
---|
57 | % it under the terms of the GNU General Public License as published |
---|
58 | % by the Free Software Foundation; either version 2 of the License, |
---|
59 | % or (at your option) any later version. |
---|
60 | % |
---|
61 | % ROMSTOOLS is distributed in the hope that it will be useful, but |
---|
62 | % WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
63 | % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
64 | % GNU General Public License for more details. |
---|
65 | % |
---|
66 | % You should have received a copy of the GNU General Public License |
---|
67 | % along with this program; if not, write to the Free Software |
---|
68 | % Foundation, Inc., 59 Temple Place, Suite 330, Boston, |
---|
69 | % MA 02111-1307 USA |
---|
70 | % |
---|
71 | % Copyright (c) 2001-2006 by Pierrick Penven |
---|
72 | % e-mail:Pierrick.Penven@ird.fr |
---|
73 | % |
---|
74 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
---|
75 | [N,M,L]=size(Tt); |
---|
76 | |
---|
77 | A00=+19092.56;A01=+209.8925; |
---|
78 | A02=-3.041638;A03=-1.852732e-3;A04=-1.361629e-5;A10=104.4077; |
---|
79 | A11=-6.500517;A12=+0.1553190;A13=2.326469e-4;AS0=-5.587545; |
---|
80 | AS1=+0.7390729;AS2=-1.909078e-2;B00=+4.721788e-1;B01=+1.028859e-2; |
---|
81 | B02=-2.512549e-4;B03=-5.939910e-7;B10=-1.571896e-2;B11=-2.598241e-4; |
---|
82 | B12=+7.267926e-6;BS1=+2.042967e-3;E00=+1.045941e-5;E01=-5.782165e-10; |
---|
83 | E02=+1.296821e-7;E10=-2.595994e-7;E11=-1.248266e-9;E12=-3.508914e-9; |
---|
84 | |
---|
85 | QR=+999.842594;Q01=+6.793952e-2;Q02=-9.095290e-3; |
---|
86 | Q03=+1.001685e-4;Q04=-1.120083e-6;Q05=+6.536332e-9;Q10=+0.824493; |
---|
87 | Q11=-4.08990e-3;Q12=+7.64380e-5;Q13=-8.24670e-7;Q14=+5.38750e-9; |
---|
88 | QS0=-5.72466e-3;QS1=+1.02270e-4;QS2=-1.65460e-6;Q20=+4.8314e-4; |
---|
89 | |
---|
90 | |
---|
91 | sqrtTs=sqrt(Ts); |
---|
92 | |
---|
93 | K0=A00+Tt.*(A01+Tt.*(A02+Tt.*(A03+Tt.*A04)))... |
---|
94 | +Ts.*(A10+Tt.*(A11+Tt.*(A12+Tt.*A13))... |
---|
95 | +sqrtTs.*(AS0+Tt.*(AS1+Tt.*AS2))); |
---|
96 | K1=B00+Tt.*(B01+Tt.*(B02+Tt.*B03))... |
---|
97 | +Ts.*(B10+Tt.*(B11+Tt.*B12)+sqrtTs.*BS1); |
---|
98 | K2=E00+Tt.*(E01+Tt.*E02)... |
---|
99 | +Ts.*(E10+Tt.*(E11+Tt.*E12)); |
---|
100 | rho1=QR+Tt.*(Q01+Tt.*(Q02+Tt.*(Q03+Tt.*(Q04+Tt.*Q05))))... |
---|
101 | +Ts.*(Q10+Tt.*(Q11+Tt.*(Q12+Tt.*(Q13+Tt.*Q14)))... |
---|
102 | +sqrtTs.*(QS0+Tt.*(QS1+Tt.*QS2))+Ts.*Q20); |
---|
103 | rho=rho1./(1+0.1.*z_r./(K0-z_r.*(K1-z_r.*K2))); |
---|
104 | |
---|
105 | if nargin > 3 |
---|
106 | bvf=0.*z_w; |
---|
107 | cff=g./rho0; |
---|
108 | bvf(2:N,:,:)=-cff*(rho1(2:N,:,:)./... |
---|
109 | (1.+0.1.*z_w(2:N,:,:)./... |
---|
110 | ( K0(2:N,:,:)-z_w(2:N,:,:).*(K1(2:N,:,:)-z_w(2:N,:,:).*K2(2:N,:,:))))... |
---|
111 | -rho1(1:N-1,:,:)./( 1.+0.1.*z_w(2:N,:,:)./... |
---|
112 | ( K0(1:N-1,:,:)-z_w(2:N,:,:).*(K1(1:N-1,:,:)-z_w(2:N,:,:).*K2(1:N-1,:,:)))))... |
---|
113 | ./(z_r(2:N,:,:)-z_r(1:N-1,:,:)); |
---|
114 | else |
---|
115 | bvf=0; |
---|
116 | end |
---|
117 | return |
---|
118 | |
---|