| 1 |
function [RHOP,RHOIS,RHOR] = density(TN,SN,GDEPT,ZREF) |
| 2 |
% [RHOP,RHOIS,RHOR] = density(Tpot,S,P,PREF) |
| 3 |
% |
| 4 |
% Adapted from OPA, subroutine EOS.F (bferron@ifremer.fr) |
| 5 |
% |
| 6 |
% ROUTINE EOS |
| 7 |
% ******************* |
| 8 |
% |
| 9 |
% PURPOSE : |
| 10 |
% ------- |
| 11 |
% COMPUTE THE POTENTIAL DENSITY RHOP, IN SITU DENSITY RHOIS, AND, IF A |
| 12 |
% REFERENCE LEVEL PREF IS SPECIFIED, THE DENSITY RHOR REFERENCED TO PREF |
| 13 |
% FROM POTENTIAL TEMPERATURE AND SALINITY FIELDS. |
| 14 |
% |
| 15 |
% IMPORTANT NOTE: This version minimize the deviations of in-situ |
| 16 |
% density compared to the UNESCO equation within |
| 17 |
% [0-4500]dbars!!!!!! It is the version used in OPA |
| 18 |
% (differences < 2e-4 kg/m3) |
| 19 |
% Below 6000 dbar , differences are > 1e-3 kg/m-3 |
| 20 |
% |
| 21 |
% METHOD : |
| 22 |
% ------ |
| 23 |
% JACKETT AND McDOUGALL (1994) EQUATION OF STATE |
| 24 |
% ********************* |
| 25 |
% THE IN SITU DENSITY IS COMPUTED DIRECTLY AS A FUNCTION OF |
| 26 |
% POTENTIAL TEMPERATURE RELATIVE TO THE SURFACE, SALT AND PRESSURE |
| 27 |
% |
| 28 |
% |
| 29 |
% RHO(T,S,P) |
| 30 |
% |
| 31 |
% WITH PRESSURE P DECIBARS |
| 32 |
% POTENTIAL TEMPERATURE T DEG CELSIUS |
| 33 |
% SALINITY S PSU |
| 34 |
% DENSITY RHO KG/M**3 |
| 35 |
% |
| 36 |
% CHECK VALUE: RHO = 1.04183326 KG/M**3 FOR P=10000 DBAR, |
| 37 |
% T = 40 DEG CELCIUS, S=40 PSU |
| 38 |
% |
| 39 |
|
| 40 |
% |
| 41 |
% JACKETT AND McDOUGALL (1994) FORMULATION |
| 42 |
%-------------------------------------------- |
| 43 |
% |
| 44 |
%... potential temperature, salinity and depth |
| 45 |
ZT = TN; |
| 46 |
ZS = SN; |
| 47 |
if size(TN,2)>1 & size(TN,1)>1 & size(TN,1)==length(GDEPT) |
| 48 |
ZH = GDEPT*ones(1,size(ZT,2)); % comented 13/04/99 for WHP Thorpe calc. |
| 49 |
elseif size(TN,2)>1 & size(TN,1)>1 & size(TN,2)==length(GDEPT) |
| 50 |
ZH = ones(size(ZT,1),1) *GDEPT'; |
| 51 |
elseif sum(size(TN))==2 |
| 52 |
ZH = ones(size(TN))*GDEPT; |
| 53 |
else |
| 54 |
ZH = GDEPT; % Added 13/04/99 for WHP Thorpe calc. |
| 55 |
end |
| 56 |
|
| 57 |
%... square root salinity |
| 58 |
ZSR= sqrt(ZS); |
| 59 |
%... compute density pure water at atm pressure |
| 60 |
ZR1= ((((6.536332E-9*ZT-1.120083E-6).*ZT+1.001685E-4).*ZT ... |
| 61 |
-9.095290E-3).*ZT+6.793952E-2).*ZT+999.842594; |
| 62 |
%... seawater density atm pressure |
| 63 |
ZR2= (((5.3875E-9*ZT-8.2467E-7).*ZT+7.6438E-5).*ZT ... |
| 64 |
-4.0899E-3).*ZT+0.824493; |
| 65 |
ZR3= (-1.6546E-6*ZT+1.0227E-4).*ZT-5.72466E-3; |
| 66 |
ZR4= 4.8314E-4; |
| 67 |
% |
| 68 |
%... potential density (referenced to the surface) |
| 69 |
RHOP= (ZR4*ZS + ZR3.*ZSR + ZR2).*ZS + ZR1; |
| 70 |
% |
| 71 |
%... add the compression terms |
| 72 |
ZE = (-3.508914E-8*ZT-1.248266E-8).*ZT-2.595994E-6; |
| 73 |
ZBW= ( 1.296821E-6*ZT-5.782165E-9).*ZT+1.045941E-4; |
| 74 |
ZB = ZBW + ZE .* ZS; |
| 75 |
% |
| 76 |
ZD = -2.042967E-2; |
| 77 |
ZC = (-7.267926E-5*ZT+2.598241E-3).*ZT+0.1571896; |
| 78 |
ZAW= ((5.939910E-6*ZT+2.512549E-3).*ZT-0.1028859).*ZT-4.721788; |
| 79 |
ZA = ( ZD*ZSR + ZC).*ZS + ZAW; |
| 80 |
% |
| 81 |
ZB1= (-0.1909078*ZT+7.390729).*ZT-55.87545; |
| 82 |
ZA1= ((2.326469E-3*ZT+1.553190).*ZT-65.00517).*ZT+1044.077; |
| 83 |
ZKW= (((-1.361629E-4*ZT-1.852732E-2).*ZT-30.41638).*ZT ... |
| 84 |
+2098.925).*ZT+190925.6; |
| 85 |
ZK0= (ZB1.*ZSR + ZA1).*ZS + ZKW; |
| 86 |
% |
| 87 |
%... in situ density |
| 88 |
RHOIS = RHOP ./ (1.0-ZH./(ZK0-ZH.*(ZA-ZH.*ZB))); |
| 89 |
|
| 90 |
%... density referenced to level ZREF |
| 91 |
if nargin==4 |
| 92 |
RHOR = RHOP ./ (1.0-ZREF./(ZK0-ZREF.*(ZA-ZREF.*ZB))); |
| 93 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph