| 1 |
function HT = calcHeatTransDirect(varargin) |
| 2 |
|
| 3 |
% HT = calcHeatTransDirect(d,g,time,flu,blkFile,[optional]); |
| 4 |
% HT = calcHeatTransDirect(d,g,time,flu,blkFile,'grav',9.81); |
| 5 |
% |
| 6 |
% Input arguements: |
| 7 |
% The incoming field data (d) and grid data (g) must be in a structured |
| 8 |
% array format (which is the format that comes from rdmnc): |
| 9 |
% d [Field data] hUtave,hVtave,uVeltave,vVeltave,Ttave,UTtave, |
| 10 |
% VTtave,(Stave,UStave,VStave for atm) |
| 11 |
% g [Grid data ] drF,dxG,dyG,dxC,dyC,HFacW,HFacS,rA |
| 12 |
% Other input parameters: |
| 13 |
% time (vec) Time levels to analyze ([] for all) |
| 14 |
% flu (str) 'O' or 'A' for ocean or atmosphere |
| 15 |
% blkFile (str) Broken line file (eg 'isoLat_cs32_59.mat') |
| 16 |
% Optional parameters: |
| 17 |
% 'grav' (num, default 9.81) Acceleration due to gravity |
| 18 |
% 'LhVap' (num, default 2501) Latent heat of vaporization |
| 19 |
% 'CpO' (num, default 3994) Specific heat capacity of water |
| 20 |
% 'RhoO' (num, default 1030) Density of sea water |
| 21 |
% 'CpA' (num, default 1004) Specific heat capacity of water |
| 22 |
% %'DiffKh' (num, default 800) Horizontal diffusivity |
| 23 |
% |
| 24 |
% Output: |
| 25 |
% HT_Out is a structured array with the following fields: |
| 26 |
% time ([nt,1]) Time axis |
| 27 |
% ylatHT ([ny,1]) Heat transport latitude axis |
| 28 |
% ylatSF ([ny-1,1]) Implied heating latitude axis |
| 29 |
% AreaZon ([ny,1]) Area in zonal bands |
| 30 |
% SenHT ([ny,nBas,nt,nFld]) Sensible heat transport |
| 31 |
% SenSF ([ny-1,nBas,nt,nFld]) Implied heating above |
| 32 |
% LatHT ([ny,nBas,nt,nFld]) Latent heat transport (atm only) |
| 33 |
% LatSF ([ny-1,nBas,nt,nFld]) Implied heating from aboce (atm only) |
| 34 |
% Currently, the routine is only configured to handle the global basin, |
| 35 |
% so nBas = 1 for the output. ny is defined by the broken line file used |
| 36 |
% for the cube calculation. nFld is the heat transport component: |
| 37 |
% nFld = 1 = Eulerian circulation HT |
| 38 |
% nFld = 2 = HT by time mean circulation |
| 39 |
% nFld = 3 = Residual [3=1-2] |
| 40 |
% nFld = 4 = HT by zonal mean circulation |
| 41 |
% nFld = 5 = Residual [5=2-4] |
| 42 |
% Ocn only: |
| 43 |
% If GM advective form: |
| 44 |
% nFld = 6 = HT by horizontal diffusion |
| 45 |
% If Skew flux form: |
| 46 |
% nFld = 6 = total (GM+Redi) eddy transport |
| 47 |
% nFld = 7 = GM (advective) eddy transport |
| 48 |
% |
| 49 |
% Description: |
| 50 |
% Calculation heat transport, and to degree possible, decomposition. |
| 51 |
% Heat transport is given in PW and the implied surface heating/flux |
| 52 |
% in W/m^2. The incoming data arrays are all assumed to be of the |
| 53 |
% dimensions [6*nc,nc,nr,nt]. |
| 54 |
% |
| 55 |
% Original Author: Jean-Michel Campin |
| 56 |
% Modifications: Daniel Enderton |
| 57 |
|
| 58 |
% Default constants |
| 59 |
LhVap = 2501; |
| 60 |
grav = 9.81; |
| 61 |
CpO = 3994; |
| 62 |
RhoO = 1030; |
| 63 |
CpA = 1004; |
| 64 |
kappa = 2/7; |
| 65 |
masking=0; |
| 66 |
|
| 67 |
PHIref = grav * [ 431.199 2105.998 5536.327 10440.915 18190.889]; |
| 68 |
%PHIref = grav * [ 173.888 529.598 896.724 1276.210 1669.052 ... |
| 69 |
% 2076.380 2499.554 2940.137 3399.937 3881.061 ... |
| 70 |
% 4385.993 4917.689 5480.084 6078.962 6722.700 ... |
| 71 |
% 7420.744 8183.151 9023.181 9959.131 11019.028... |
| 72 |
% 12253.003 13723.604 15599.238 18412.089 24539.570]; |
| 73 |
|
| 74 |
% Read input parameters. |
| 75 |
d = varargin{1}; |
| 76 |
g = varargin{2}; |
| 77 |
dE= varargin{3}; |
| 78 |
GM= varargin{4}; |
| 79 |
nt= varargin{5}; |
| 80 |
flu = varargin{6}; |
| 81 |
blkFile = varargin{7}; |
| 82 |
%for ii = 8:2:length(varargin) |
| 83 |
% temp = varargin{ii+1}; eval([varargin{ii},' = temp;']); |
| 84 |
%end |
| 85 |
if length(varargin) >= 8 & ~isempty(varargin{8}) |
| 86 |
mask = varargin{8}; |
| 87 |
masking = 1; |
| 88 |
disp('Masking is on...') |
| 89 |
end |
| 90 |
if length(varargin) == 9 |
| 91 |
PHIref= varargin{9}; |
| 92 |
end |
| 93 |
|
| 94 |
nBas = 0; |
| 95 |
if isequal(flu,'A'), nout = 5; end |
| 96 |
if isequal(flu,'O'), nout = 7; end |
| 97 |
flag_tot=0; |
| 98 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 99 |
% Prepare / reform incoming data % |
| 100 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 101 |
|
| 102 |
% Determine time indecies. |
| 103 |
%if isempty(time), time = d.iters_read_from_file; i_time = 1:length(time); |
| 104 |
%else [dump,i_time] = ismember(time,d.iters_read_from_file); end |
| 105 |
|
| 106 |
nc = size(g.HFacC,2); |
| 107 |
nr = size(g.HFacC,3); |
| 108 |
i_time=1:nt; |
| 109 |
|
| 110 |
if nr~=5 |
| 111 |
LhVap = 2500e3; |
| 112 |
grav = 9.79764; |
| 113 |
RDGAS = 287.04; |
| 114 |
CpA = RDGAS/kappa; |
| 115 |
end |
| 116 |
|
| 117 |
ac = reshape(g.rA ,[6*nc*nc, 1]); |
| 118 |
hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 119 |
hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]); |
| 120 |
dxc = reshape(g.dxC(1:6*nc,1:nc) ,[6*nc*nc, 1]); |
| 121 |
dyc = reshape(g.dyC(1:6*nc,1:nc) ,[6*nc*nc, 1]); |
| 122 |
dxg = reshape(g.dxG(1:6*nc,1:nc) ,[6*nc*nc, 1]); |
| 123 |
dyg = reshape(g.dyG(1:6*nc,1:nc) ,[6*nc*nc, 1]); |
| 124 |
drf = reshape(g.drF,[1,length(g.drF)]); |
| 125 |
Z = reshape(g.Z,[1,length(g.Z)]); |
| 126 |
|
| 127 |
%u = reshape(d.UVEL(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 128 |
%v = reshape(d.VVEL(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 129 |
ut = reshape(dE.UTHMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 130 |
vt = reshape(dE.VTHMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 131 |
hu = reshape(dE.UVELMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 132 |
hv = reshape(dE.VVELMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 133 |
t = d.THETA(1:6*nc,1:nc,:,i_time); |
| 134 |
if isequal(flu,'O') |
| 135 |
Kux = reshape(GM.GM_Kux(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 136 |
Kvy = reshape(GM.GM_Kvy(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 137 |
|
| 138 |
if isfield(GM,'GM_ubT') %%% if advective form of GM |
| 139 |
UbTgm = reshape(GM.GM_ubT(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 140 |
VbTgm = reshape(GM.GM_vbT(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 141 |
elseif isfield(GM,'GM_KuzTz') %%% if skew-flux form of GM |
| 142 |
UbTgm = reshape(GM.GM_KuzTz(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 143 |
VbTgm = reshape(GM.GM_KvzTz(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 144 |
elseif isfield(dE,'ADVx_TH') %%% if total advective field |
| 145 |
flag_tot = 1; |
| 146 |
UbTgm = reshape(dE.ADVx_TH(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 147 |
VbTgm = reshape(dE.ADVy_TH(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]); |
| 148 |
else |
| 149 |
UbTgm = zeros(6*nc*nc,nr,nt); |
| 150 |
VbTgm = zeros(6*nc*nc,nr,nt); |
| 151 |
end |
| 152 |
end |
| 153 |
|
| 154 |
if isequal(flu,'A') |
| 155 |
up = reshape(dE.UVELPHI(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 156 |
vp = reshape(dE.VVELPHI(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 157 |
uq = reshape(dE.USLTMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 158 |
vq = reshape(dE.VSLTMASS(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]); |
| 159 |
q = d.SALT(1:6*nc,1:nc,:,i_time); |
| 160 |
p = d.PHIHYD(1:6*nc,1:nc,:,i_time); |
| 161 |
ut = ut .* ( repmat(Z,[6*nc*nc 1 nt]) / 1e5 ).^kappa; |
| 162 |
vt = vt .* ( repmat(Z,[6*nc*nc 1 nt]) / 1e5 ).^kappa; |
| 163 |
t = t .* ( repmat(reshape(Z,[1 1 nr]),[6*nc nc 1 nt]) / 1e5 ).^kappa; |
| 164 |
end |
| 165 |
|
| 166 |
if masking == 1 |
| 167 |
hu = repmat(reshape(mask.maskW(:,:,1),6*nc*nc,1),[1 nr nt]) .* hu; |
| 168 |
hv = repmat(reshape(mask.maskS(:,:,1),6*nc*nc,1),[1 nr nt]) .* hv; |
| 169 |
ut = repmat(reshape(mask.maskW(:,:,1),6*nc*nc,1),[1 nr nt]) .* ut; |
| 170 |
vt = repmat(reshape(mask.maskS(:,:,1),6*nc*nc,1),[1 nr nt]) .* vt; |
| 171 |
if isequal(flu,'O') |
| 172 |
Kux = repmat(reshape(mask.maskW(:,:,1),6*nc*nc,1),[1 nr nt]) .* Kux; |
| 173 |
Kvy = repmat(reshape(mask.maskS(:,:,1),6*nc*nc,1),[1 nr nt]) .* Kvy; |
| 174 |
UbTgm = repmat(reshape(mask.maskW(:,:,1),6*nc*nc,1),[1 nr nt]) .* UbTgm; |
| 175 |
VbTgm = repmat(reshape(mask.maskS(:,:,1),6*nc*nc,1),[1 nr nt]) .* VbTgm; |
| 176 |
end |
| 177 |
end |
| 178 |
|
| 179 |
|
| 180 |
% Load broken line information. Compute (tracer point) cell area between |
| 181 |
% broken lines for each basin. There are nbkl broken lines and nbkl+1 |
| 182 |
% bands between broken lines. The variable "bkl_Zon" gives the zone number |
| 183 |
% (nbkl+1 total) for a given index between 0 and nbkl, that is, nbkl+1 |
| 184 |
% total zones. Comments block is for eventual addition of multiple basin |
| 185 |
% calculations. |
| 186 |
load(blkFile); |
| 187 |
ydim=length(bkl_Ylat); |
| 188 |
AreaZon=zeros(ydim+1,1+nBas); |
| 189 |
for j = 1:ydim+1 |
| 190 |
izon = find(bkl_Zon == j-1); |
| 191 |
AreaZon(j,1) = sum(ac(izon)); |
| 192 |
% for b = 1:nBas, |
| 193 |
% tmp = ac.*mskBc(:,b); |
| 194 |
% AreaZon(j,1+b) = sum(tmp(izon)); |
| 195 |
% end |
| 196 |
end |
| 197 |
|
| 198 |
% Latitute plotting information. Average latitude of a broken line |
| 199 |
% (ylatHF) is calculated from a mean value of the y value of all the edges. |
| 200 |
% The latitude at the surface flux points is a mean of the broken line mean |
| 201 |
% values. |
| 202 |
YlatAv=sum(bkl_Ysg,1)./(1+bkl_Npts'); %' |
| 203 |
ylatHT = [-90,YlatAv,90]; |
| 204 |
ylatSF = ( ylatHT(2:ydim+2) + ylatHT(1:ydim+1) )./2; |
| 205 |
|
| 206 |
% The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point |
| 207 |
% and -2/2 if it has a v point. Positive/negative values contribute |
| 208 |
% positively/negatively to northward heat transport (this depends on the |
| 209 |
% orientation of the cell). A zero value indicates an end of edges that |
| 210 |
% contribute to a broken line. The u and v information is parced into two |
| 211 |
% seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac). |
| 212 |
ufac = zeros([size(bkl_Flg),1+nBas]); |
| 213 |
vfac = zeros([size(bkl_Flg),1+nBas]); |
| 214 |
ufac(:,:,1) = rem(bkl_Flg,2); |
| 215 |
vfac(:,:,1) = fix(bkl_Flg/2); |
| 216 |
% for jl=1:ydim, |
| 217 |
% ie=bkl_Npts(jl); |
| 218 |
% for b=1:nBas, |
| 219 |
% ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1); |
| 220 |
% vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1); |
| 221 |
% end |
| 222 |
% end |
| 223 |
ufacabs = abs(ufac); |
| 224 |
vfacabs = abs(vfac); |
| 225 |
|
| 226 |
% Prepare mask(s). |
| 227 |
% ??? I temporarily took out the code to configure the masks beyond this |
| 228 |
% global one. Does this need to account for a ridge if present? |
| 229 |
mskG=ones(ydim+1,1+nBas); |
| 230 |
|
| 231 |
% Area factors. "ArW_Dif" and "ArS_Dif" the areas for the western and |
| 232 |
% southern edge of cells, respectively The "_Dif" suffix indicates the |
| 233 |
% areas used for the diffusivity because there is an extra "dxc" or "dyc" |
| 234 |
% from the gradient (here for computational efficiency reasons). The |
| 235 |
% division by "dxc" and "dyc" is associates with gradient of temperature. |
| 236 |
% |
| 237 |
% ??? Why is the land mask (hw and hs) only in the diffusive area term? |
| 238 |
ArW = dyg*reshape(drf,[1,length(drf)]); |
| 239 |
ArS = dxg*reshape(drf,[1,length(drf)]); |
| 240 |
ArW_Dif=hw.*((dyg./dxc)*reshape(drf,[1,length(drf)])); |
| 241 |
ArS_Dif=hs.*((dxg./dyc)*reshape(drf,[1,length(drf)])); |
| 242 |
|
| 243 |
% Compute the temperature and its gradient and the velocity points: |
| 244 |
% tbi/tdi = temperature between/difference i points (at u points) |
| 245 |
% tbj/tdj = temperature between/difference j points (at v points) |
| 246 |
% The cube is first split and the extra points are added to the files. |
| 247 |
% Then the means and differences are taken. Note that the division of dxc |
| 248 |
% and dyc for the gradient is not applied until later for computational |
| 249 |
% efficiency purposes. The arrays are then croped and reshaped to the |
| 250 |
% format for the other field variables, [6*nc*nc,nr]. Note that the |
| 251 |
% split_C_cub function adds a row/column of tracer points in front of the |
| 252 |
% first row/column of the tile matries. This, when the differences and |
| 253 |
% gradients are computed and cropped, the off indecies are selected from |
| 254 |
% [2:nc+1] rather than [1:nc]. (This was a bit mystifying to me). |
| 255 |
t6bi=zeros(nc,nc+1,nr,nt,6); t6di=t6bi; q6bi=t6bi; |
| 256 |
t6bj=zeros(nc+1,nc,nr,nt,6); t6dj=t6bj; q6bj=t6bj; |
| 257 |
t6t=split_C_cub(t); |
| 258 |
t6bi([1:nc],:,:,:,:) = ( t6t([1:nc],:,:,:,:) + t6t([2:nc+1],:,:,:,:) )./2; |
| 259 |
t6bj(:,[1:nc],:,:,:) = ( t6t(:,[1:nc],:,:,:) + t6t(:,[2:nc+1],:,:,:) )./2; |
| 260 |
t6di([1:nc],:,:,:,:) = ( t6t([1:nc],:,:,:,:) - t6t([2:nc+1],:,:,:,:) ); |
| 261 |
t6dj(:,[1:nc],:,:,:) = ( t6t(:,[1:nc],:,:,:) - t6t(:,[2:nc+1],:,:,:) ); |
| 262 |
tbi = t6bi([1:nc],[2:nc+1],:,:,:); |
| 263 |
tbj = t6bj([2:nc+1],[1:nc],:,:,:); |
| 264 |
tdi = t6di([1:nc],[2:nc+1],:,:,:); |
| 265 |
tdj = t6dj([2:nc+1],[1:nc],:,:,:); |
| 266 |
tbi=reshape(permute(tbi,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 267 |
tbj=reshape(permute(tbj,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 268 |
tdi=reshape(permute(tdi,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 269 |
tdj=reshape(permute(tdj,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 270 |
if isequal(flu,'A') |
| 271 |
q6t=split_C_cub(q); |
| 272 |
q6bi([1:nc],:,:,:,:) = (q6t([1:nc],:,:,:,:)+q6t([2:nc+1],:,:,:,:))./2; |
| 273 |
q6bj(:,[1:nc],:,:,:) = (q6t(:,[1:nc],:,:,:)+q6t(:,[2:nc+1],:,:,:))./2; |
| 274 |
qbi = q6bi([1:nc],[2:nc+1],:,:,:); |
| 275 |
qbj = q6bj([2:nc+1],[1:nc],:,:,:); |
| 276 |
qbi=reshape(permute(qbi,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 277 |
qbj=reshape(permute(qbj,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 278 |
|
| 279 |
p6t=split_C_cub(p); |
| 280 |
p6bi([1:nc],:,:,:,:) = (p6t([1:nc],:,:,:,:)+p6t([2:nc+1],:,:,:,:))./2; |
| 281 |
p6bj(:,[1:nc],:,:,:) = (p6t(:,[1:nc],:,:,:)+p6t(:,[2:nc+1],:,:,:))./2; |
| 282 |
pbi = p6bi([1:nc],[2:nc+1],:,:,:); |
| 283 |
pbj = p6bj([2:nc+1],[1:nc],:,:,:); |
| 284 |
pbi=reshape(permute(pbi,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 285 |
pbj=reshape(permute(pbj,[1,5,2,3,4]),[6*nc*nc,nr,nt]); |
| 286 |
end |
| 287 |
|
| 288 |
% Prepare output arrays. "nout" is the number of transport output fields. |
| 289 |
% It is currently hard-coded, but could eventually be an input parameters |
| 290 |
% to set which output fields are desired if some of then become |
| 291 |
% computationally expensive. |
| 292 |
% SenHT = Sensible heat transport |
| 293 |
% SenSF = Sensible implied surface flux |
| 294 |
% IntV = Integrated volume transport |
| 295 |
% IntT = Integrated temperature |
| 296 |
SenHT = zeros(ydim+2,1+nBas,nt,nout); |
| 297 |
SenSF = zeros(ydim+1,1+nBas,nt,nout); |
| 298 |
IntV = zeros(ydim,nr,1+nBas,nt); |
| 299 |
IntT = zeros(ydim,nr,1+nBas,nt); |
| 300 |
IntM = zeros(ydim,nr,1+nBas,nt); |
| 301 |
if isequal(flu,'A') |
| 302 |
LatHT = zeros(ydim+2,1+nBas,nt,nout); |
| 303 |
PotHT = zeros(ydim+2,1+nBas,nt,nout); |
| 304 |
LatSF = zeros(ydim+1,1+nBas,nt,nout); |
| 305 |
IntQ = zeros(ydim,nr,1+nBas,nt); |
| 306 |
IntP = zeros(ydim,nr,1+nBas,nt); |
| 307 |
end |
| 308 |
Psi = zeros(ydim+2,1+nBas,nt); |
| 309 |
|
| 310 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 311 |
% Make heat transport calculations % |
| 312 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 313 |
|
| 314 |
% Preparation for calculation of zonal average temperature. The |
| 315 |
% tempereature multiplied by the appropriate length scale ("tbi_temp", |
| 316 |
% "tbj_temp") is summed up ("IntT" in the next section) and |
| 317 |
% divided by the total length ("IntM", composed from summing "hw_temp", |
| 318 |
% "hs_temp"). |
| 319 |
hw_temp = zeros(size(hw)); |
| 320 |
hs_temp = zeros(size(hs)); |
| 321 |
for k=1:nr, |
| 322 |
hw_temp(:,k) = dyg.*hw(:,k); |
| 323 |
hs_temp(:,k) = dxg.*hs(:,k); |
| 324 |
end |
| 325 |
|
| 326 |
for it = 1:length(i_time) |
| 327 |
|
| 328 |
% uz / vz = Volume transport though cell faces (velocity times area). |
| 329 |
% Used for zonal mean volume transport (4). |
| 330 |
% utz1/vtz1 = Eulerian sensible heat transport through cell faces (1). |
| 331 |
% utz2/vtz2 = Sensible heat transport through cell faces by Eulerian |
| 332 |
% mean circulations (2). |
| 333 |
% dtx1/dty1 = Temperatude gradient at cell face times the area (when |
| 334 |
% multiplied by the diffusion will be the horizontal |
| 335 |
% diffusion heat transport) (6). |
| 336 |
uz = ArW.*hu(:,:,it); |
| 337 |
vz = ArS.*hv(:,:,it); |
| 338 |
utz1 = sum(ArW.*ut(:,:,it),2); |
| 339 |
vtz1 = sum(ArS.*vt(:,:,it),2); |
| 340 |
utz2 = sum(ArW.*hu(:,:,it).*tbi(:,:,it),2); |
| 341 |
vtz2 = sum(ArS.*hv(:,:,it).*tbj(:,:,it),2); |
| 342 |
if isequal(flu,'A') |
| 343 |
uqz1 = sum(ArW.*uq(:,:,it),2); |
| 344 |
vqz1 = sum(ArS.*vq(:,:,it),2); |
| 345 |
uqz2 = sum(ArW.*hu(:,:,it).*qbi(:,:,it),2); |
| 346 |
vqz2 = sum(ArS.*hv(:,:,it).*qbj(:,:,it),2); |
| 347 |
upz1 = sum(ArW.*(up(:,:,it)+hu(:,:,it).*repmat(PHIref,[6*nc*nc 1])),2); |
| 348 |
vpz1 = sum(ArS.*(vp(:,:,it)+hv(:,:,it).*repmat(PHIref,[6*nc*nc 1])),2); |
| 349 |
upz2 = sum(ArW.*hu(:,:,it).*(pbi(:,:,it)+repmat(PHIref,[6*nc*nc 1])),2); |
| 350 |
vpz2 = sum(ArS.*hv(:,:,it).*(pbj(:,:,it)+repmat(PHIref,[6*nc*nc 1])),2); |
| 351 |
end |
| 352 |
if isequal(flu,'O') |
| 353 |
dtx1 = sum(ArW_Dif.*Kux(:,:,it).*tdi(:,:,it),2); |
| 354 |
dty1 = sum(ArS_Dif.*Kvy(:,:,it).*tdj(:,:,it),2); |
| 355 |
Kuztz1 = sum(UbTgm(:,:,it),2); |
| 356 |
Kvztz1 = sum(VbTgm(:,:,it),2); |
| 357 |
end |
| 358 |
|
| 359 |
% Preparation for calculation of zonal average temperature. The |
| 360 |
% temperature multiplied by the appropriate length scale ("tbi_temp", |
| 361 |
% "tbj_temp") is summed up ("IntT" in the next section) and |
| 362 |
% divided by the total length ("IntM", composed from summing "hw_temp", |
| 363 |
% "hs_temp"). |
| 364 |
tbi_temp = hw_temp.*tbi(:,:,it); |
| 365 |
tbj_temp = hs_temp.*tbj(:,:,it); |
| 366 |
if isequal(flu,'A') |
| 367 |
qbi_temp = hw_temp.*qbi(:,:,it); |
| 368 |
qbj_temp = hs_temp.*qbj(:,:,it); |
| 369 |
pbi_temp = hw_temp.*(pbi(:,:,it)+repmat(PHIref,[6*nc*nc 1])); |
| 370 |
pbj_temp = hs_temp.*(pbj(:,:,it)+repmat(PHIref,[6*nc*nc 1])); |
| 371 |
end |
| 372 |
|
| 373 |
% Block 1: |
| 374 |
% With the vertical integral of heat transport calculated across cell |
| 375 |
% edges, the zonal integral (along the broken line) is computed to |
| 376 |
% determine the total northward heat transport. The first for loop is |
| 377 |
% over the individual broken lines (determining the northward HT at the |
| 378 |
% representative latitude). The second loop is over the basins. The |
| 379 |
% third loop is over the individual edges along a specific broken line. |
| 380 |
% Note that an individual cell can have two edges (u and v) that have a |
| 381 |
% HT contributions. Hence the variable containing indecies of cells |
| 382 |
% with edges along the broken lines, "bkl_IJuv", has some repeats. |
| 383 |
% Note that the variable "bkl_Npts" is the number of edges along a |
| 384 |
% given broken line. Note also that the latitude axis starts at 2 |
| 385 |
% because heat transport at extremes (latitute = -90/90) is zero by |
| 386 |
% definition. Recall that index (1) is total Eulerian transport, (2) |
| 387 |
% is from the mean circulations, and (3) is from the horizontal |
| 388 |
% diffusion. |
| 389 |
% |
| 390 |
% Block 2: |
| 391 |
% Here zonal average circulation and temperature / moisture is |
| 392 |
% calculated. The zonal average volume transport v (IntV) and t |
| 393 |
% (IntT/IntM) are computed first and are multiplied together at the |
| 394 |
% end. |
| 395 |
for jl=1:ydim |
| 396 |
ie=bkl_Npts(jl); |
| 397 |
for b=1:1+nBas |
| 398 |
ij=bkl_IJuv(1:ie,jl); |
| 399 |
% Block 1: |
| 400 |
SenHT(1+jl,b,it,1) = SenHT(1+jl,b,it,1) + ... |
| 401 |
sum(ufac(1:ie,jl,b).*utz1(ij) + ... |
| 402 |
vfac(1:ie,jl,b).*vtz1(ij)); |
| 403 |
SenHT(1+jl,b,it,2) = SenHT(1+jl,b,it,2) + ... |
| 404 |
sum(ufac(1:ie,jl,b).*utz2(ij) + ... |
| 405 |
vfac(1:ie,jl,b).*vtz2(ij)); |
| 406 |
if isequal(flu,'A') |
| 407 |
LatHT(1+jl,b,it,1) = LatHT(1+jl,b,it,1) + ... |
| 408 |
sum(ufac(1:ie,jl,b).*uqz1(ij) + ... |
| 409 |
vfac(1:ie,jl,b).*vqz1(ij)); |
| 410 |
LatHT(1+jl,b,it,2) = LatHT(1+jl,b,it,2) + ... |
| 411 |
sum(ufac(1:ie,jl,b).*uqz2(ij) + ... |
| 412 |
vfac(1:ie,jl,b).*vqz2(ij)); |
| 413 |
PotHT(1+jl,b,it,1) = PotHT(1+jl,b,it,1) + ... |
| 414 |
sum(ufac(1:ie,jl,b).*upz1(ij) + ... |
| 415 |
vfac(1:ie,jl,b).*vpz1(ij)); |
| 416 |
PotHT(1+jl,b,it,2) = PotHT(1+jl,b,it,2) + ... |
| 417 |
sum(ufac(1:ie,jl,b).*upz2(ij) + ... |
| 418 |
vfac(1:ie,jl,b).*vpz2(ij)); |
| 419 |
end |
| 420 |
if isequal(flu,'O') |
| 421 |
SenHT(1+jl,b,it,6) = SenHT(1+jl,b,it,6) + ... |
| 422 |
sum(ufac(1:ie,jl,b).*dtx1(ij) + ... |
| 423 |
vfac(1:ie,jl,b).*dty1(ij)); |
| 424 |
SenHT(1+jl,b,it,7) = SenHT(1+jl,b,it,7) + ... |
| 425 |
sum(ufac(1:ie,jl,b).*Kuztz1(ij) + ... |
| 426 |
vfac(1:ie,jl,b).*Kvztz1(ij)); |
| 427 |
end |
| 428 |
% Block 2: |
| 429 |
IntV(jl,:,b,it) = IntV(jl,:,b,it) ... |
| 430 |
+ ufac(1:ie,jl,b)'*uz(ij,:) ... |
| 431 |
+ vfac(1:ie,jl,b)'*vz(ij,:) ; |
| 432 |
IntT(jl,:,b,it) = IntT(jl,:,b,it) ... |
| 433 |
+ ufacabs(1:ie,jl,b)'*tbi_temp(ij,:) ... |
| 434 |
+ vfacabs(1:ie,jl,b)'*tbj_temp(ij,:); |
| 435 |
IntM(jl,:,b,it) = IntM(jl,:,b,it) ... |
| 436 |
+ ufacabs(1:ie,jl,b)'*hw_temp(ij,:) ... |
| 437 |
+ vfacabs(1:ie,jl,b)'*hs_temp(ij,:); |
| 438 |
if isequal(flu,'A') |
| 439 |
IntQ(jl,:,b,it) = IntQ(jl,:,b,it) ... |
| 440 |
+ ufacabs(1:ie,jl,b)'*qbi_temp(ij,:) ... |
| 441 |
+ vfacabs(1:ie,jl,b)'*qbj_temp(ij,:); |
| 442 |
IntP(jl,:,b,it) = IntP(jl,:,b,it) ... |
| 443 |
+ ufacabs(1:ie,jl,b)'*pbi_temp(ij,:) ... |
| 444 |
+ vfacabs(1:ie,jl,b)'*pbj_temp(ij,:); |
| 445 |
end |
| 446 |
end |
| 447 |
end |
| 448 |
|
| 449 |
% Prepare HT output: Calculate HT by zonal flows and tabulate |
| 450 |
% residuals. Also, the multiplicative constants (Cp,rho,grav,LhVap) are |
| 451 |
% applied here to put moisture and potential temperature fluxes in |
| 452 |
% terms of heat transports. |
| 453 |
tmp = IntV(:,:,:,it) .* IntT(:,:,:,it) ... |
| 454 |
./ change(IntM(:,:,:,it),'==',0,NaN); |
| 455 |
SenHT(2:ydim+1,:,it,4) = sum(change(tmp,'==',NaN,0),2); |
| 456 |
% SenHT(2:ydim+1,:,it,4) = sum( IntV(:,:,:,it) ... |
| 457 |
% .* IntT(:,:,:,it) ... |
| 458 |
% ./ IntM(:,:,:,it),2); |
| 459 |
SenHT(:,:,it,3) = SenHT(:,:,it,1) - SenHT(:,:,it,2); |
| 460 |
SenHT(:,:,it,5) = SenHT(:,:,it,2) - SenHT(:,:,it,4); |
| 461 |
if isequal(flu,'O') |
| 462 |
% SenHT(:,:,it,6) = DiffKh.*SenHT(:,:,it,6); |
| 463 |
SenHT(:,:,it,:) = (CpO*RhoO)*SenHT(:,:,it,:); |
| 464 |
if flag_tot==1 |
| 465 |
SenHT(:,:,it,7)=SenHT(:,:,it,7)-SenHT(:,:,it,1); |
| 466 |
end |
| 467 |
elseif isequal(flu,'A') |
| 468 |
SenHT(:,:,it,:) = (CpA./grav).*SenHT(:,:,it,:); |
| 469 |
LatHT(2:ydim+1,:,it,4) = sum( IntV(:,:,:,it) ... |
| 470 |
.* IntQ(:,:,:,it) ... |
| 471 |
./ IntM(:,:,:,it),2); |
| 472 |
LatHT(:,:,it,3) = LatHT(:,:,it,1) - LatHT(:,:,it,2); |
| 473 |
LatHT(:,:,it,5) = LatHT(:,:,it,2) - LatHT(:,:,it,4); |
| 474 |
LatHT(:,:,it,:) = (LhVap./grav).*LatHT(:,:,it,:); |
| 475 |
|
| 476 |
PotHT(2:ydim+1,:,it,4) = sum( IntV(:,:,:,it) ... |
| 477 |
.* IntP(:,:,:,it) ... |
| 478 |
./ IntM(:,:,:,it),2); |
| 479 |
PotHT(:,:,it,3) = PotHT(:,:,it,1) - PotHT(:,:,it,2); |
| 480 |
PotHT(:,:,it,5) = PotHT(:,:,it,2) - PotHT(:,:,it,4); |
| 481 |
PotHT(:,:,it,:) = 1/grav * PotHT(:,:,it,:); |
| 482 |
end |
| 483 |
Psi(2:ydim+1,:,it) = - sum( IntV(:,:,:,it) ,2 ); |
| 484 |
|
| 485 |
% Implied surface heat flux from heat transports (implied heating). |
| 486 |
% Tabulated as the difference in heat transports between two broken |
| 487 |
% lines divided by the zonal band area. |
| 488 |
mskG = reshape(mskG,(ydim+1)*(1+nBas),1); |
| 489 |
I = find(mskG==0); |
| 490 |
mskG = reshape(mskG,ydim+1,1+nBas); |
| 491 |
var = zeros(ydim+1,1+nBas); |
| 492 |
for n=1:min(nout,6), |
| 493 |
varT = SenHT([2:ydim+2],:,it,n) ... |
| 494 |
- SenHT([1:ydim+1],:,it,n); |
| 495 |
varT = reshape(varT,(ydim+1)*(1+nBas),1); varT(I)=NaN; |
| 496 |
varT = reshape(varT,ydim+1,1+nBas); |
| 497 |
SenSF([1:ydim+1],:,it,n) = varT./AreaZon; |
| 498 |
if isequal(flu,'A') |
| 499 |
for n=1:min(nout,6) |
| 500 |
varQ = LatHT([2:ydim+2],:,it,n) ... |
| 501 |
- LatHT([1:ydim+1],:,it,n); |
| 502 |
varQ = reshape(varQ,(ydim+1)*(1+nBas),1); varQ(I)=NaN; |
| 503 |
varQ = reshape(varQ,ydim+1,1+nBas); |
| 504 |
LatSF([1:ydim+1],:,it,n) = varQ./AreaZon; |
| 505 |
end |
| 506 |
end |
| 507 |
end |
| 508 |
end |
| 509 |
|
| 510 |
|
| 511 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 512 |
% Assign outputs, put in units of PW % |
| 513 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 514 |
|
| 515 |
SenHT = SenHT*1e-15; |
| 516 |
|
| 517 |
HT.time = i_time; |
| 518 |
HT.SenHT = reshape(SenHT,ydim+2,nt,nout); |
| 519 |
HT.SenSF = reshape(SenSF,ydim+1,nt,nout); |
| 520 |
HT.ylatHT = ylatHT; |
| 521 |
HT.ylatSF = ylatSF; |
| 522 |
HT.AreaZon = AreaZon; |
| 523 |
HT.PsiSurf = 1e-6*reshape(Psi,ydim+2,nt); |
| 524 |
|
| 525 |
if isequal(flu,'A') |
| 526 |
LatHT = LatHT*1e-15; |
| 527 |
PotHT = PotHT*1e-15; |
| 528 |
HT.LatHT = reshape(LatHT,ydim+2,nt,nout); |
| 529 |
HT.LatSF = reshape(LatSF,ydim+1,nt,nout); |
| 530 |
HT.PotHT = reshape(PotHT,ydim+2,nt,nout); |
| 531 |
end |
| 532 |
|
| 533 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 534 |
% Stuff that might need to be added back in later (for basins) % |
| 535 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 536 |
|
| 537 |
% Block for constructing mskG for different basins: |
| 538 |
% if nBas > 0, |
| 539 |
% mskBc=rdda([Rac,'maskC_bas.bin'],[6*nc*nc 3],1,'real*4','b'); |
| 540 |
% mskBw=rdda([Rac,'maskW_bas.bin'],[6*nc*nc 3],1,'real*4','b'); |
| 541 |
% mskBs=rdda([Rac,'maskS_bas.bin'],[6*nc*nc 3],1,'real*4','b'); |
| 542 |
% if nBas==2, |
| 543 |
% mskBc(:,2)=mskBc(:,2)+mskBc(:,3); |
| 544 |
% mskBw(:,2)=mskBw(:,2)+mskBw(:,3); |
| 545 |
% mskBs(:,2)=mskBs(:,2)+mskBs(:,3); |
| 546 |
% mskBc=min(1,mskBc); mskBw=min(1,mskBw); mskBs=min(1,mskBs); |
| 547 |
% end |
| 548 |
% %- load: np_Sep, ij_Sep, tp_Sep: |
| 549 |
% sep_lineF=[Rac,'sepBas_cs32_60']; |
| 550 |
% load(sep_lineF); |
| 551 |
% fprintf([' + bassin mask & Sep.line:',sep_lineF,' \n']); |
| 552 |
% %- compute mask for each bassin: ----------- |
| 553 |
% kMsep=1; |
| 554 |
% if kMsep, |
| 555 |
% mskW=1+min(1,ceil(hw(:,1))); |
| 556 |
% mskS=1+min(1,ceil(hs(:,1))); |
| 557 |
% for b=1:nBas, |
| 558 |
% bs=b; b1=1+bs; b2=2+rem(bs,nBas); |
| 559 |
% if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end |
| 560 |
% for j=1:ydim+1, |
| 561 |
% for i=1:np_Sep(bs,j), |
| 562 |
% ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i)); |
| 563 |
% if typ == 1, |
| 564 |
% mskG(j,b1)=mskG(j,b1)*mskW(ij); |
| 565 |
% mskG(j,b2)=mskG(j,b2)*mskW(ij); |
| 566 |
% elseif typ == 2, |
| 567 |
% mskG(j,b1)=mskG(j,b1)*mskS(ij); |
| 568 |
% mskG(j,b2)=mskG(j,b2)*mskS(ij); |
| 569 |
% end |
| 570 |
% end |
| 571 |
% end |
| 572 |
% end |
| 573 |
% mskG=2-min(2,mskG); |
| 574 |
% end |
| 575 |
% end |
| 576 |
|
| 577 |
% % Masking for different basin. |
| 578 |
% mskZ = ones(ydim+2,1+nBas); |
| 579 |
% mskZ([2:ydim+1],:) = mskG([1:ydim],:) + mskG([2:ydim+1],:); |
| 580 |
% mskZ = reshape(min(mskZ,1),(ydim+2)*(1+nBas),1); |
| 581 |
% I = find(mskZ == 0); |
| 582 |
% mskZ = reshape(mskZ,(ydim+2),1+nBas); |
| 583 |
% for n=1:min(6,nout) |
| 584 |
% var=SenHT(:,:,n); |
| 585 |
% var=reshape(var,(ydim+2)*(1+nBas),1); var(I) = NaN; |
| 586 |
% var=reshape(var,(ydim+2),1+nBas); SenHT(:,:,n) = var; |
| 587 |
% end |
| 588 |
|
| 589 |
% Makes sure that there are no zeros in area AreaZon: |
| 590 |
% - set AreaZon to 1 if = zero |
| 591 |
% AreaZon=reshape(AreaZon,(ydim+1)*(1+nBas),1); |
| 592 |
% [I]=find(AreaZon==0); AreaZon(I)=1; |
| 593 |
% AreaZon=reshape(AreaZon,ydim+1,1+nBas); |
Parent Directory
| 
Revision Log
| 
Revision Graph