MITgcm_contrib/gmaze_pv/diagWALIN.m

% [F,A,D,CROP] = diagWALIN(FLAG,C1,C2,Qnet,Snet,Classes,dA)
% 
% DESCRIPTION:
% Compute the transformation rate of a surface outcrop class (potential
% density or SST) from surface net heat flux Qnet and salt flux Snet
% according to the Walin theory.
%
% INPUTS: 
% FLAG    : Can either be: 0, 1 or 2
%           0: Outcrop field is surface potential density computed 
%              from C1=SST and C2=SSS
%           1: Outcrop field is surface potential density given by C1
%           2: Outcrop field is SST and potential density is computed 
%              from C1=SST and C2=SSS
% C1,C2   : Depends on option FLAG:
%           - FLAG = 0 : 
%                        C1 : Sea surface temperature (degC) 
%                        C2 : Sea surface salinity (PSU)
%           - FLAG = 1 : 
%                        C1 : Surface potential density (kg/m3) 
%                        C2 : Not used
%           - FLAG = 2 : 
%                        C1 : Sea surface temperature (degC) 
%                        C2 : Sea surface salinity (PSU)
% Qnet    : Downward net surface heat flux (W/m2)
% Snet    : Downward net surface salt flux (kg/m2/s) -> 
%           ie, Snet = rho*beta*SSS*(E-P)
% Classes : Range of outcrops to explore (eg: [20:.1:30] for potential density)
% lon,lat : axis
% dA      : Matrix of grid surface elements (m2) centered in (lon,lat) of Ci
%
%
% OUTPUTS:
% F(3,:)    : Transformation rate (m3/s) (from 1:Qnet, 2:Snet and 3:Total)
% A         : Surface of each outcrops
% D(3,:,:)  : Maps of density flux (kg/m2/s) from 1:Qnet, 2:Snet and 3:Total
% CROP(:,:) : Map of the surface field used to compute outcrop's contours
%
%
% NOTES:
% - Fields are of the format: C(LAT,LON)
% - The potential density is computed with the equation of state routine from
%   the MITgcm called densjmd95.m 
%   (see: http://mitgcm.org/cgi-bin/viewcvs.cgi/MITgcm_contrib/gmaze_pv/subfct/densjmd95.m)
% - Snet may be filled of NaN if not available, its F component won't computed
%
%
% AUTHOR: 
% Guillaume Maze / MIT 2006
% 
% HISTORY:
% - Revised: 06/28/2007
%            * Add option do directly give the pot. density as input
%            * Add options do take SST as outcrop 
% - Created: 06/22/2007
%
% REFERENCES: 
% Walin G. 1982: On the relation between sea-surface 
% heat flow and thermal circulation in the ocean. Tellus N24
%

% The routine is not optimized for speed but for clarity, that's why we
% compute buoyancy fluxes, etc...
%
% TO DO: 
% - Fix signs in density fluxes to be correct albeit consistent with F right now
% - Create options for non regular CLASS
% - Create options to also compute the formation rate M
% - Create options to compute an error bar
% - Create check of inputs section

function varargout = diagWALIN(FLAG,C1,C2,QNET,SNET,CLASS,dA)


% 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPROC
% Variables:
nlat = size(C1,1);
nlon = size(C1,2);
CLASS = CLASS(:);
  
% Determine surface fields from which we'll take outcrops contours:
switch FLAG
  
 case {0,2} % Need to compute SIGMA THETA
  SST = C1;
  SSS = C2;
  if FLAG == 0     % Outcrop is SIGMA THETA:
     OUTCROP = ST;
     ST = densjmd95(SSS,SST,zeros(nlat,nlon)) - 1000;  % Real surface (depth = 0)
     %dpt = -5; ST = densjmd95(SSS,SST,(0.09998*9.81*dpt)*ones(nlat,nlon)) - 1000; % Model surface
  elseif FLAG == 2 % Outcrop is SST:
     OUTCROP = SST;
     if length(find(isnan(SSS)==1)) == nlat*nlon 
       ST = ones(nlat,nlon).*1035;
     else
       ST = densjmd95(SSS,SST,zeros(nlat,nlon)) - 1000;
     end
  end
  
 case 1
  ST = C1; % Potential density
  OUTCROP = ST;
end
  
% Create a flag if we don't find salt flux:
if length(find(isnan(SNET)==1)) == nlat*nlon
  do_ep = 0;
else
  do_ep = 1;
end

% Physical constants:
g = 9.81;        % Gravity (m/s2)
Cp = 3994;       % Specific heat of sea water (J/K/kg)
rho0 = 1035;     % Density of reference (kg/m3)
rho  = ST+1000;  % Density (kg/m3)
                 % Thermal expansion coefficient (1/K)
if exist('SST') & exist('SSS') & length(find(isnan(SSS)==1)) ~= nlat*nlon 
  alpha = sw_alpha(SSS,SST,zeros(nlat,nlon));
else
  alpha = 2.*1e-4; 
end

%ix=200;iy=100;[SST(iy,ix),SSS(iy,ix),QNET(iy,ix),SNET(iy,ix)]

% 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BUOYANCY FLUX: b
% The buoyancy flux (m/s2*m/s=m2/s3) is computed as:
% b = g/rho*( alpha/Cp*QNET - SNET )
% b = g/rho*alpha/Cp*QNET - g/rho*SNET
% b = b_hf + b_ep
% QNET the net heat flux (W/m2) and SNET the net salt flux (kg/m2/s) 
              b_hf =  g.*alpha./Cp.*QNET./rho;
if do_ep==1,  b_ep = -g*SNET./rho; else b_ep = zeros(nlat,nlon); end
                 b = b_hf + do_ep*b_ep;


% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DENSITY FLUX: bd
% Buoyancy flux is transformed into density flux (kg/m3*m/s = kg/m2/s):
% bd = - rho/g * b
% with b the buoyancy flux
             bd_hf = - rho/g.*b_hf; 
             bd_ep = - rho/g.*b_ep;
             bd    = - rho/g.*b;

%[bd_hf(iy,ix),bd_ep(iy,ix),bd(iy,ix)]
             
% 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% NET MASS FLUX INTEGRATED OVER OUTCROPS: Bd
% The amount of mass water flux over an outcrop is computed as:
% Bd = SUM_ij bd(i,j)*dA(i,j)*MASK(i,j,OUTCROP)
% with MASK(i,j,OUTCROP) = 1 where  OUTCROP(i,j)-dC/2 <=  OUTCROP(i,j) < OUTCROP(i,j)+dC/2
%                        = 0 otherwise
% Outcrops are defined with an increment of:
dCROP = diff(CLASS(1:2));

switch FLAG
 case {0,1}, coef = 1;                 % Potential density as outcrops
 case 2,     coef = 1./(alpha.*rho0);  % SST as outcrops
end %switch


% Surface integral:
for iC = 1 : length(CLASS)
  CROPc  = CLASS(iC);
  mask   = zeros(nlat,nlon);
  mask(find( (CROPc-dCROP/2 <= OUTCROP) & (OUTCROP < CROPc+dCROP/2) )) = 1;
  %if CROPc == 18,[CROPc-dCROP/2 CROPc+dCROP/2],global mask18,mask18=mask;end;
               Bd_hf(iC) = nansum(nansum(dA.*mask.*bd_hf.*coef,1),2);
               Bd_ep(iC) = nansum(nansum(dA.*mask.*bd_ep.*coef,1),2);
                  Bd(iC) = nansum(nansum(dA.*mask.*bd.*coef,1),2);
                  AA(iC) = nansum(nansum(dA.*mask,1),2);
end %for iC


% 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSFORMATION RATE: F
% F is defined as the convergence/divergence of the integrated mass flux Bd.
% F = Bd(CROP) / dCROP
% where Bd is the mass flux over an outcrop.
             F_hf = Bd_hf./dCROP;
             F_ep = Bd_ep./dCROP; 
             F    = Bd./dCROP;


% 5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% OUTPUTS
% Transformation rate:
TRANSFORM_RATE(1,:) = F_hf;
TRANSFORM_RATE(2,:) = F_ep;
TRANSFORM_RATE(3,:) = F;             

% Density flux:
DENSITY_FLUX(1,:,:) = bd_hf;
DENSITY_FLUX(2,:,:) = bd_ep;
DENSITY_FLUX(3,:,:) = bd;

switch nargout
 case 1
  varargout(1) = {TRANSFORM_RATE};
 case 2
  varargout(1) = {TRANSFORM_RATE};
  varargout(2) = {AA};
 case 3
  varargout(1) = {TRANSFORM_RATE};
  varargout(2) = {AA};
  varargout(3) = {DENSITY_FLUX};
 case 4
  varargout(1) = {TRANSFORM_RATE};
  varargout(2) = {AA};
  varargout(3) = {DENSITY_FLUX};
  varargout(4) = {OUTCROP};
end %switch
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 = rhobetaSSS*(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.099989.81dpt)*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/CpQNET - SNET )
129	% b = g/rhoalpha/CpQNET - 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