| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph