1 |
% [F,A,D,CROP] = diagWALIN(FLAG,C1,C2,Qnet,Snet,Classes,dA) |
2 |
% |
3 |
% DESCRIPTION: |
4 |
% Compute the transformation rate of a surface outcrop class (potential |
5 |
% density or SST) from surface net heat flux Qnet and salt flux Snet |
6 |
% according to the Walin theory. |
7 |
% |
8 |
% INPUTS: |
9 |
% FLAG : Can either be: 0, 1 or 2 |
10 |
% 0: Outcrop field is surface potential density computed |
11 |
% from C1=SST and C2=SSS |
12 |
% 1: Outcrop field is surface potential density given by C1 |
13 |
% 2: Outcrop field is SST and potential density is computed |
14 |
% from C1=SST and C2=SSS |
15 |
% C1,C2 : Depends on option FLAG: |
16 |
% - FLAG = 0 : |
17 |
% C1 : Sea surface temperature (degC) |
18 |
% C2 : Sea surface salinity (PSU) |
19 |
% - FLAG = 1 : |
20 |
% C1 : Surface potential density (kg/m3) |
21 |
% C2 : Not used |
22 |
% - FLAG = 2 : |
23 |
% C1 : Sea surface temperature (degC) |
24 |
% C2 : Sea surface salinity (PSU) |
25 |
% Qnet : Downward net surface heat flux (W/m2) |
26 |
% Snet : Downward net surface salt flux (kg/m2/s) -> |
27 |
% ie, Snet = rho*beta*SSS*(E-P) |
28 |
% Classes : Range of outcrops to explore (eg: [20:.1:30] for potential density) |
29 |
% lon,lat : axis |
30 |
% dA : Matrix of grid surface elements (m2) centered in (lon,lat) of Ci |
31 |
% |
32 |
% |
33 |
% OUTPUTS: |
34 |
% F(3,:) : Transformation rate (m3/s) (from 1:Qnet, 2:Snet and 3:Total) |
35 |
% A : Surface of each outcrops |
36 |
% D(3,:,:) : Maps of density flux (kg/m2/s) from 1:Qnet, 2:Snet and 3:Total |
37 |
% CROP(:,:) : Map of the surface field used to compute outcrop's contours |
38 |
% |
39 |
% |
40 |
% NOTES: |
41 |
% - Fields are of the format: C(LAT,LON) |
42 |
% - The potential density is computed with the equation of state routine from |
43 |
% the MITgcm called densjmd95.m |
44 |
% (see: http://mitgcm.org/cgi-bin/viewcvs.cgi/MITgcm_contrib/gmaze_pv/subfct/densjmd95.m) |
45 |
% - Snet may be filled of NaN if not available, its F component won't computed |
46 |
% |
47 |
% |
48 |
% AUTHOR: |
49 |
% Guillaume Maze / MIT 2006 |
50 |
% |
51 |
% HISTORY: |
52 |
% - Revised: 06/28/2007 |
53 |
% * Add option do directly give the pot. density as input |
54 |
% * Add options do take SST as outcrop |
55 |
% - Created: 06/22/2007 |
56 |
% |
57 |
% REFERENCES: |
58 |
% Walin G. 1982: On the relation between sea-surface |
59 |
% heat flow and thermal circulation in the ocean. Tellus N24 |
60 |
% |
61 |
|
62 |
% The routine is not optimized for speed but for clarity, that's why we |
63 |
% compute buoyancy fluxes, etc... |
64 |
% |
65 |
% TO DO: |
66 |
% - Fix signs in density fluxes to be correct albeit consistent with F right now |
67 |
% - Create options for non regular CLASS |
68 |
% - Create options to also compute the formation rate M |
69 |
% - Create options to compute an error bar |
70 |
% - Create check of inputs section |
71 |
|
72 |
function varargout = diagWALIN(FLAG,C1,C2,QNET,SNET,CLASS,dA) |
73 |
|
74 |
|
75 |
% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC |
76 |
% Variables: |
77 |
nlat = size(C1,1); |
78 |
nlon = size(C1,2); |
79 |
CLASS = CLASS(:); |
80 |
|
81 |
% Determine surface fields from which we'll take outcrops contours: |
82 |
switch FLAG |
83 |
|
84 |
case {0,2} % Need to compute SIGMA THETA |
85 |
SST = C1; |
86 |
SSS = C2; |
87 |
if FLAG == 0 % Outcrop is SIGMA THETA: |
88 |
OUTCROP = ST; |
89 |
ST = densjmd95(SSS,SST,zeros(nlat,nlon)) - 1000; % Real surface (depth = 0) |
90 |
%dpt = -5; ST = densjmd95(SSS,SST,(0.09998*9.81*dpt)*ones(nlat,nlon)) - 1000; % Model surface |
91 |
elseif FLAG == 2 % Outcrop is SST: |
92 |
OUTCROP = SST; |
93 |
if length(find(isnan(SSS)==1)) == nlat*nlon |
94 |
ST = ones(nlat,nlon).*1035; |
95 |
else |
96 |
ST = densjmd95(SSS,SST,zeros(nlat,nlon)) - 1000; |
97 |
end |
98 |
end |
99 |
|
100 |
case 1 |
101 |
ST = C1; % Potential density |
102 |
OUTCROP = ST; |
103 |
end |
104 |
|
105 |
% Create a flag if we don't find salt flux: |
106 |
if length(find(isnan(SNET)==1)) == nlat*nlon |
107 |
do_ep = 0; |
108 |
else |
109 |
do_ep = 1; |
110 |
end |
111 |
|
112 |
% Physical constants: |
113 |
g = 9.81; % Gravity (m/s2) |
114 |
Cp = 3994; % Specific heat of sea water (J/K/kg) |
115 |
rho0 = 1035; % Density of reference (kg/m3) |
116 |
rho = ST+1000; % Density (kg/m3) |
117 |
% Thermal expansion coefficient (1/K) |
118 |
if exist('SST') & exist('SSS') & length(find(isnan(SSS)==1)) ~= nlat*nlon |
119 |
alpha = sw_alpha(SSS,SST,zeros(nlat,nlon)); |
120 |
else |
121 |
alpha = 2.*1e-4; |
122 |
end |
123 |
|
124 |
%ix=200;iy=100;[SST(iy,ix),SSS(iy,ix),QNET(iy,ix),SNET(iy,ix)] |
125 |
|
126 |
% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BUOYANCY FLUX: b |
127 |
% The buoyancy flux (m/s2*m/s=m2/s3) is computed as: |
128 |
% b = g/rho*( alpha/Cp*QNET - SNET ) |
129 |
% b = g/rho*alpha/Cp*QNET - g/rho*SNET |
130 |
% b = b_hf + b_ep |
131 |
% QNET the net heat flux (W/m2) and SNET the net salt flux (kg/m2/s) |
132 |
b_hf = g.*alpha./Cp.*QNET./rho; |
133 |
if do_ep==1, b_ep = -g*SNET./rho; else b_ep = zeros(nlat,nlon); end |
134 |
b = b_hf + do_ep*b_ep; |
135 |
|
136 |
|
137 |
% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DENSITY FLUX: bd |
138 |
% Buoyancy flux is transformed into density flux (kg/m3*m/s = kg/m2/s): |
139 |
% bd = - rho/g * b |
140 |
% with b the buoyancy flux |
141 |
bd_hf = - rho/g.*b_hf; |
142 |
bd_ep = - rho/g.*b_ep; |
143 |
bd = - rho/g.*b; |
144 |
|
145 |
%[bd_hf(iy,ix),bd_ep(iy,ix),bd(iy,ix)] |
146 |
|
147 |
% 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NET MASS FLUX INTEGRATED OVER OUTCROPS: Bd |
148 |
% The amount of mass water flux over an outcrop is computed as: |
149 |
% Bd = SUM_ij bd(i,j)*dA(i,j)*MASK(i,j,OUTCROP) |
150 |
% with MASK(i,j,OUTCROP) = 1 where OUTCROP(i,j)-dC/2 <= OUTCROP(i,j) < OUTCROP(i,j)+dC/2 |
151 |
% = 0 otherwise |
152 |
% Outcrops are defined with an increment of: |
153 |
dCROP = diff(CLASS(1:2)); |
154 |
|
155 |
switch FLAG |
156 |
case {0,1}, coef = 1; % Potential density as outcrops |
157 |
case 2, coef = 1./(alpha.*rho0); % SST as outcrops |
158 |
end %switch |
159 |
|
160 |
|
161 |
% Surface integral: |
162 |
for iC = 1 : length(CLASS) |
163 |
CROPc = CLASS(iC); |
164 |
mask = zeros(nlat,nlon); |
165 |
mask(find( (CROPc-dCROP/2 <= OUTCROP) & (OUTCROP < CROPc+dCROP/2) )) = 1; |
166 |
%if CROPc == 18,[CROPc-dCROP/2 CROPc+dCROP/2],global mask18,mask18=mask;end; |
167 |
Bd_hf(iC) = nansum(nansum(dA.*mask.*bd_hf.*coef,1),2); |
168 |
Bd_ep(iC) = nansum(nansum(dA.*mask.*bd_ep.*coef,1),2); |
169 |
Bd(iC) = nansum(nansum(dA.*mask.*bd.*coef,1),2); |
170 |
AA(iC) = nansum(nansum(dA.*mask,1),2); |
171 |
end %for iC |
172 |
|
173 |
|
174 |
% 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSFORMATION RATE: F |
175 |
% F is defined as the convergence/divergence of the integrated mass flux Bd. |
176 |
% F = Bd(CROP) / dCROP |
177 |
% where Bd is the mass flux over an outcrop. |
178 |
F_hf = Bd_hf./dCROP; |
179 |
F_ep = Bd_ep./dCROP; |
180 |
F = Bd./dCROP; |
181 |
|
182 |
|
183 |
|
184 |
% 5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS |
185 |
% Transformation rate: |
186 |
TRANSFORM_RATE(1,:) = F_hf; |
187 |
TRANSFORM_RATE(2,:) = F_ep; |
188 |
TRANSFORM_RATE(3,:) = F; |
189 |
|
190 |
% Density flux: |
191 |
DENSITY_FLUX(1,:,:) = bd_hf; |
192 |
DENSITY_FLUX(2,:,:) = bd_ep; |
193 |
DENSITY_FLUX(3,:,:) = bd; |
194 |
|
195 |
switch nargout |
196 |
case 1 |
197 |
varargout(1) = {TRANSFORM_RATE}; |
198 |
case 2 |
199 |
varargout(1) = {TRANSFORM_RATE}; |
200 |
varargout(2) = {AA}; |
201 |
case 3 |
202 |
varargout(1) = {TRANSFORM_RATE}; |
203 |
varargout(2) = {AA}; |
204 |
varargout(3) = {DENSITY_FLUX}; |
205 |
case 4 |
206 |
varargout(1) = {TRANSFORM_RATE}; |
207 |
varargout(2) = {AA}; |
208 |
varargout(3) = {DENSITY_FLUX}; |
209 |
varargout(4) = {OUTCROP}; |
210 |
end %switch |