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	gmaze	1.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