Warning, /utils/matlab/cs_grid/calcHeatTransDirect.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 4cf7d6e4 on 2006-06-29 18:52:07 UTC
d43e8328b9 Dani*0001 function HT = calcHeatTransDirect(varargin)
0002
0003 % HT = calcHeatTransDirect(d,g,time,flu,blkFile,[optional]);
0004 % HT = calcHeatTransDirect(d,g,time,flu,blkFile,'grav',9.81);
0005 %
0006 % Input arguements:
0007 % The incoming field data (d) and grid data (g) must be in a structured
0008 % array format (which is the format that comes from rdmnc):
0009 % d [Field data] hUtave,hVtave,uVeltave,vVeltave,Ttave,UTtave,
0010 % VTtave,(Stave,UStave,VStave for atm)
0011 % g [Grid data ] drF,dxG,dyG,dxC,dyC,HFacW,HFacS,rA
0012 % Other input parameters:
0013 % time (vec) Time levels to analyze ([] for all)
0014 % flu (str) 'O' or 'A' for ocean or atmosphere
0015 % blkFile (str) Broken line file (eg 'isoLat_cs32_59.mat')
0016 % Optional parameters:
0017 % 'grav' (num, default 9.81) Acceleration due to gravity
0018 % 'LhVap' (num, default 2501) Latent heat of vaporization
0019 % 'CpO' (num, default 3994) Specific heat capacity of water
0020 % 'RhoO' (num, default 1035) Density of sea water
0021 % 'CpA' (num, default 1004) Specific heat capacity of water
0022 % 'DiffKh' (num, default 800) Horizontal diffusivity
0023 %
32ec7e18d0 Dani*0024 % Output:
d43e8328b9 Dani*0025 % HT_Out is a structured array with the following fields:
0026 % time ([nt,1]) Time axis
0027 % ylatHT ([ny,1]) Heat transport latitude axis
0028 % ylatSF ([ny-1,1]) Implied heating latitude axis
0029 % AreaZon ([ny,1]) Area in zonal bands
0030 % SenHT ([ny,nBas,nt,nFld]) Sensible heat transport
0031 % SenSF ([ny-1,nBas,nt,nFld]) Implied heating above
0032 % LatHT ([ny,nBas,nt,nFld]) Latent heat transport (atm only)
0033 % LatSF ([ny-1,nBas,nt,nFld]) Implied heating from aboce (atm only)
0034 % Currently, the routine is only configured to handle the global basin,
0035 % so nBas = 1 for the output. ny is defined by the broken line file used
0036 % for the cube calculation. nFld is the heat transport component:
0037 % nFld = 1 = Eulerian circulation HT
0038 % nFld = 2 = HT by mean circulation
0039 % nFld = 3 = Residual [3=1-2]
0040 % nFld = 4 = HT by zonal mean circulation
0041 % nFld = 5 = Residual [5=2-4]
0042 % nFld = 6 = HT by horizontal diffusion (ocn only)
32ec7e18d0 Dani*0043 %
0044 % Description:
0045 % Calculation heat transport, and to degree possible, decomposition.
d43e8328b9 Dani*0046 % Heat transport is given in PW and the implied surface heating/flux
32ec7e18d0 Dani*0047 % in W/m^2. The incoming data arrays are all assumed to be of the
d43e8328b9 Dani*0048 % dimensions [6*nc,nc,nr,nt].
0049 %
0050 % Original Author: Jean-Michel Campin
0051 % Modifications: Daniel Enderton
0052
0053 % Default constants (can be overriden).
0054 LhVap = 2501;
0055 grav = 9.81;
0056 CpO = 3994;
0057 RhoO = 1035;
0058 CpA = 1004;
0059 DiffKh = 800;
0060
0061 % Read input parameters.
0062 d = varargin{1};
0063 g = varargin{2};
0064 time = varargin{3};
0065 flu = varargin{4};
0066 blkFile = varargin{5};
0067 for ii = 6:2:length(varargin)
0068 temp = varargin{ii+1}; eval([varargin{ii},' = temp;']);
0069 end
32ec7e18d0 Dani*0070
d43e8328b9 Dani*0071 nBas = 0;
0072 nlfs = 1;
0073 if isequal(flu,'A'), nout = 5; end
0074 if isequal(flu,'O'), nout = 6; end
32ec7e18d0 Dani*0075 if nlfs ~= 1, error('Not ready to handle non-nonlinear free surface!'); end
0076 if nBas ~= 0, error('Not ready to handle multiple basins'); end
0077
0078
0079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0080 % Prepare / reform incoming data %
0081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0082
0083 % Determine time indecies.
0084 if isempty(time), time = d.T; i_time = 1:length(time);
0085 else [dump,i_time] = ismember(time,d.T); end
0086
0087 nc = size(g.HFacC,2);
0088 nr = size(g.HFacC,3);
0089 nt = length(time);
0090
0091 ac = reshape(g.rA ,[6*nc*nc, 1]);
0092 hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0093 hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0094 dxc = reshape(g.dxC(1:6*nc,1:nc) ,[6*nc*nc, 1]);
0095 dyc = reshape(g.dyC(1:6*nc,1:nc) ,[6*nc*nc, 1]);
0096 dxg = reshape(g.dxG(1:6*nc,1:nc) ,[6*nc*nc, 1]);
0097 dyg = reshape(g.dyG(1:6*nc,1:nc) ,[6*nc*nc, 1]);
0098 drf = reshape(g.drF,[1,length(g.drF)]);
0099
0100 u = reshape(d.uVeltave(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]);
0101 v = reshape(d.vVeltave(1:6*nc,1:nc,:,i_time),[6*nc*nc,nr,nt]);
0102 ut = reshape(d.UTtave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
0103 vt = reshape(d.VTtave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
4cf7d6e47a Dani*0104 if isequal(flu,'A'),
0105 uq = reshape(d.UStave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
0106 vq = reshape(d.VStave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
0107 end
32ec7e18d0 Dani*0108 hu = reshape(d.hUtave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
0109 hv = reshape(d.hVtave(1:6*nc,1:nc,:,i_time) ,[6*nc*nc,nr,nt]);
0110 t = d.Ttave(1:6*nc,1:nc,:,i_time);
4cf7d6e47a Dani*0111 if isequal(flu,'A'),
0112 q = d.Stave(1:6*nc,1:nc,:,i_time);
0113 end
32ec7e18d0 Dani*0114
0115 % Load broken line information. Compute (tracer point) cell area between
0116 % broken lines for each basin. There are nbkl broken lines and nbkl+1
0117 % bands between broken lines. The variable "bkl_Zon" gives the zone number
0118 % (nbkl+1 total) for a given index between 0 and nbkl, that is, nbkl+1
0119 % total zones. Comments block is for eventual addition of multiple basin
0120 % calculations.
0121 load(blkFile);
0122 ydim=length(bkl_Ylat);
0123 AreaZon=zeros(ydim+1,1+nBas);
0124 for j = 1:ydim+1
0125 izon = find(bkl_Zon == j-1);
0126 AreaZon(j,1) = sum(ac(izon));
0127 % for b = 1:nBas,
0128 % tmp = ac.*mskBc(:,b);
0129 % AreaZon(j,1+b) = sum(tmp(izon));
0130 % end
0131 end
0132
0133 % Latitute plotting information. Average latitude of a broken line
0134 % (ylatHF) is calculated from a mean value of the y value of all the edges.
0135 % The latitude at the surface flux points is a mean of the broken line mean
0136 % values.
92699e6cf5 Dani*0137 YlatAv=sum(bkl_Ysg,1)./(1+bkl_Npts'); %'
32ec7e18d0 Dani*0138 ylatHT = [-90,YlatAv,90];
0139 ylatSF = ( ylatHT(2:ydim+2) + ylatHT(1:ydim+1) )./2;
0140
0141 % The variable "bkl_Flg" is -1/1 if edge (on a given broken) has a u point
0142 % and -2/2 if it has a v point. Positive/negative values contribute
0143 % positively/negatively to northward heat transport (this depends on the
0144 % oreientation of the cell). A zero value indicates an end of edges that
0145 % contribute to a broken line. The u and v information is parced into two
0146 % seperate fields, ufac and vfac (-2/2 are reduced to -1/1 for vfac).
0147 ufac = zeros([size(bkl_Flg),1+nBas]);
0148 vfac = zeros([size(bkl_Flg),1+nBas]);
0149 ufac(:,:,1) = rem(bkl_Flg,2);
0150 vfac(:,:,1) = fix(bkl_Flg/2);
0151 % for jl=1:ydim,
0152 % ie=bkl_Npts(jl);
0153 % for b=1:nBas,
0154 % ufac(1:ie,jl,1+b)=mskBw(bkl_IJuv(1:ie,jl),b).*ufac(1:ie,jl,1);
0155 % vfac(1:ie,jl,1+b)=mskBs(bkl_IJuv(1:ie,jl),b).*vfac(1:ie,jl,1);
0156 % end
0157 % end
0158 ufacabs = abs(ufac);
0159 vfacabs = abs(vfac);
0160
0161 % Prepare mask(s).
0162 % ??? I temporarily took out the code to configure the masks beyond this
0163 % global one. Does this need to account for a ridge if present?
0164 mskG=ones(ydim+1,1+nBas);
0165
0166 % Area factors. "ArW_Dif" and "ArS_Dif" the areas for the western and
0167 % southern edge of cells, respectively The "_Dif" suffix indicates the
0168 % areas used for the diffusivity because there is an extra "dxc" or "dyc"
0169 % from the gradient (here for computational efficiency reasons). The
0170 % division by "dxc" and "dyc" is associates with gradient of temperature.
0171 %
0172 % ??? Why is the land mask (hw and hs) only in the diffusive area term?
0173 ArW = dyg*reshape(drf,[1,length(drf)]);
0174 ArS = dxg*reshape(drf,[1,length(drf)]);
0175 ArW_Dif=hw.*((dyg./dxc)*reshape(drf,[1,length(drf)]));
0176 ArS_Dif=hs.*((dxg./dyc)*reshape(drf,[1,length(drf)]));
0177
0178 % Compute the temperature and its gradient and the velocity points:
0179 % tbi/tdi = temperature between/difference i points (at u points)
0180 % tbj/tdj = temperature between/difference j points (at v points)
0181 % The cube is first split and the extra points are added to the files.
0182 % Then the means and differences are taken. Note that the division of dxc
0183 % and dyc for the gradient is not applied until later for computational
0184 % efficiency purposes. The arrays are then croped and reshaped to the
0185 % format for the other field variables, [6*nc*nc,nr]. Note that the
0186 % split_C_cub function adds a row/column of tracer points in front of the
0187 % first row/column of the tile matries. This, when the differences and
0188 % gradients are computed and cropped, the off indecies are selected from
0189 % [2:nc+1] rather than [1:nc]. (This was a bit mystifying to me).
64c5c5fbe0 Dani*0190 t6bi=zeros(nc,nc+1,nr,nt,6); t6di=t6bi; q6bi=t6bi;
0191 t6bj=zeros(nc+1,nc,nr,nt,6); t6dj=t6bj; q6bj=t6bj;
32ec7e18d0 Dani*0192 t6t=split_C_cub(t);
0193 t6bi([1:nc],:,:,:,:) = ( t6t([1:nc],:,:,:,:) + t6t([2:nc+1],:,:,:,:) )./2;
0194 t6bj(:,[1:nc],:,:,:) = ( t6t(:,[1:nc],:,:,:) + t6t(:,[2:nc+1],:,:,:) )./2;
0195 t6di([1:nc],:,:,:,:) = ( t6t([1:nc],:,:,:,:) - t6t([2:nc+1],:,:,:,:) );
0196 t6dj(:,[1:nc],:,:,:) = ( t6t(:,[1:nc],:,:,:) - t6t(:,[2:nc+1],:,:,:) );
0197 tbi = t6bi([1:nc],[2:nc+1],:,:,:);
0198 tbj = t6bj([2:nc+1],[1:nc],:,:,:);
0199 tdi = t6di([1:nc],[2:nc+1],:,:,:);
0200 tdj = t6dj([2:nc+1],[1:nc],:,:,:);
0201 tbi=reshape(permute(tbi,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
0202 tbj=reshape(permute(tbj,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
0203 tdi=reshape(permute(tdi,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
0204 tdj=reshape(permute(tdj,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
64c5c5fbe0 Dani*0205 if isequal(flu,'A')
0206 q6t=split_C_cub(q);
0207 q6bi([1:nc],:,:,:,:) = (q6t([1:nc],:,:,:,:)+q6t([2:nc+1],:,:,:,:))./2;
0208 q6bj(:,[1:nc],:,:,:) = (q6t(:,[1:nc],:,:,:)+q6t(:,[2:nc+1],:,:,:))./2;
0209 qbi = q6bi([1:nc],[2:nc+1],:,:,:);
0210 qbj = q6bj([2:nc+1],[1:nc],:,:,:);
0211 qbi=reshape(permute(qbi,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
0212 qbj=reshape(permute(qbj,[1,5,2,3,4]),[6*nc*nc,nr,nt]);
0213 end
32ec7e18d0 Dani*0214
0215 % Prepare output arrays. "nout" is the number of transport output fields.
0216 % It is currently hard-coded, but could eventually be an input parameters
0217 % to set which output fields are desired if some of then become
0218 % computationally expensive.
d43e8328b9 Dani*0219 % SenHT = Sensible heat transport
0220 % SenSF = Sensible implied surface flux
32ec7e18d0 Dani*0221 % IntV = Integrated volume transport
0222 % IntT = Integrated temperature
64c5c5fbe0 Dani*0223 SenHT = zeros(ydim+2,1+nBas,nt,nout);
0224 SenSF = zeros(ydim+1,1+nBas,nt,nout);
32ec7e18d0 Dani*0225 IntV = zeros(ydim,nr,1+nBas,nt);
0226 IntT = zeros(ydim,nr,1+nBas,nt);
0227 IntM = zeros(ydim,nr,1+nBas,nt);
64c5c5fbe0 Dani*0228 if isequal(flu,'A')
d43e8328b9 Dani*0229 LatHT = zeros(ydim+2,1+nBas,nt,nout);
92699e6cf5 Dani*0230 Lat = LatHT;
d43e8328b9 Dani*0231 LatSF = zeros(ydim+1,1+nBas,nt,nout);
0232 IntQ = zeros(ydim,nr,1+nBas,nt);
64c5c5fbe0 Dani*0233 end
32ec7e18d0 Dani*0234
0235
0236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0237 % Make heat transport calculations %
0238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0239
0240 % Preparation for calculation of zonal average temperature. The
0241 % tempereature multiplied by the appropriate length scale ("tbi_temp",
0242 % "tbj_temp") is summed up ("IntT" in the next section) and
0243 % divided by the total length ("IntM", composed from summing "hw_temp",
0244 % "hs_temp").
0245 hw_temp = zeros(size(hw));
0246 hs_temp = zeros(size(hs));
0247 for k=1:nr,
0248 hw_temp(:,k) = dyg.*hw(:,k);
0249 hs_temp(:,k) = dxg.*hs(:,k);
0250 end
0251
0252 for it = 1:length(time)
0253
0254 % uz / vz = Volume transport though cell faces (velocity times area).
0255 % Used for zonal mean volume transport (4).
0256 % utz1/vtz1 = Eulerian sensible heat transport through cell faces (1).
0257 % utz2/vtz2 = Sensible heat transport through cell faces by Eulerian
0258 % mean circulations (2).
0259 % dtx1/dty1 = Temperatude gradient at cell face times the area (when
0260 % multiplied by the diffusion will be the horizontal
0261 % diffusion heat transport) (6).
0262 uz = ArW.*hu(:,:,it);
0263 vz = ArS.*hv(:,:,it);
0264 utz1 = sum(ArW.*ut(:,:,it),2);
0265 vtz1 = sum(ArS.*vt(:,:,it),2);
0266 utz2 = sum(ArW.*hu(:,:,it).*tbi(:,:,it),2);
0267 vtz2 = sum(ArS.*hv(:,:,it).*tbj(:,:,it),2);
64c5c5fbe0 Dani*0268 if isequal(flu,'A')
0269 uqz1 = sum(ArW.*uq(:,:,it),2);
0270 vqz1 = sum(ArS.*vq(:,:,it),2);
0271 uqz2 = sum(ArW.*hu(:,:,it).*qbi(:,:,it),2);
0272 vqz2 = sum(ArS.*hv(:,:,it).*qbj(:,:,it),2);
0273 end
d43e8328b9 Dani*0274 dtx1 = sum(ArW_Dif.*tdi(:,:,it),2);
0275 dty1 = sum(ArS_Dif.*tdj(:,:,it),2);
32ec7e18d0 Dani*0276
0277 % Preparation for calculation of zonal average temperature. The
0278 % tempereature multiplied by the appropriate length scale ("tbi_temp",
0279 % "tbj_temp") is summed up ("IntT" in the next section) and
0280 % divided by the total length ("IntM", composed from summing "hw_temp",
0281 % "hs_temp").
0282 tbi_temp = hw_temp.*tbi(:,:,it);
0283 tbj_temp = hs_temp.*tbj(:,:,it);
64c5c5fbe0 Dani*0284 if isequal(flu,'A')
0285 qbi_temp = hw_temp.*qbi(:,:,it);
0286 qbj_temp = hs_temp.*qbj(:,:,it);
0287 end
32ec7e18d0 Dani*0288
0289 % Block 1:
0290 % With the vertical integral of heat transport calculated across cell
0291 % edges, the zonal integral (along the broken line) is computed to
d43e8328b9 Dani*0292 % determine the total northward heat transport. The first for loop is
0293 % over the individual broken lines (determining the northward HT at the
0294 % representative latitude). The second loop is over the basins. The
0295 % third loop is over the individual edges along a specific broken line.
0296 % Note that an individual cell can have two edges (u and v) that have a
0297 % HT contributions. Hence the variable containing indecies of cells
0298 % with edges along the broken lines, "bkl_IJuv", has some repeats.
0299 % Note that the variable "bkl_Npts" is the number of edges along a
0300 % given broken line. Note also that the latitude axis starts at 2
0301 % because heat transport at extremes (latitute = -90/90) is zero by
0302 % definition. Recall that index (1) is total Eulerian transport, (2)
0303 % is from the mean circulations, and (3) is from the horizontal
0304 % diffusion.
32ec7e18d0 Dani*0305 %
0306 % Block 2:
d43e8328b9 Dani*0307 % Here zonal average circulation and temperature / moisture is
32ec7e18d0 Dani*0308 % calculated. The zonal average volume transport v (IntV) and t
0309 % (IntT/IntM) are computed first and are multiplied together at the
0310 % end.
0311 for jl=1:ydim
0312 ie=bkl_Npts(jl);
0313 for b=1:1+nBas
92699e6cf5 Dani*0314 ij=bkl_IJuv(1:ie,jl);
0315 % Block 1:
0316 SenHT(1+jl,b,it,1) = SenHT(1+jl,b,it,1) + ...
0317 sum(ufac(1:ie,jl,b).*utz1(ij) + ...
0318 vfac(1:ie,jl,b).*vtz1(ij));
0319 SenHT(1+jl,b,it,2) = SenHT(1+jl,b,it,2) + ...
0320 sum(ufac(1:ie,jl,b).*utz2(ij) + ...
0321 vfac(1:ie,jl,b).*vtz2(ij));
0322 if isequal(flu,'A')
0323 LatHT(1+jl,b,it,1) = LatHT(1+jl,b,it,1) + ...
0324 sum(ufac(1:ie,jl,b).*uqz1(ij) + ...
0325 vfac(1:ie,jl,b).*vqz1(ij));
0326 LatHT(1+jl,b,it,2) = LatHT(1+jl,b,it,2) + ...
0327 sum(ufac(1:ie,jl,b).*uqz2(ij) + ...
0328 vfac(1:ie,jl,b).*vqz2(ij));
0329 end
0330 if isequal(flu,'O')
0331 SenHT(1+jl,b,it,6) = SenHT(1+jl,b,it,6) + ...
0332 sum(ufac(1:ie,jl,b).*dtx1(ij) + ...
0333 vfac(1:ie,jl,b).*dty1(ij));
0334 end
0335 % Block 2:
0336 IntV(jl,:,b,it) = IntV(jl,:,b,it) ...
0337 + ufac(1:ie,jl,b)'*uz(ij,:) ...
0338 + vfac(1:ie,jl,b)'*vz(ij,:) ;
0339 IntT(jl,:,b,it) = IntT(jl,:,b,it) ...
0340 + ufacabs(1:ie,jl,b)'*tbi_temp(ij,:) ...
0341 + vfacabs(1:ie,jl,b)'*tbj_temp(ij,:);
0342 IntM(jl,:,b,it) = IntM(jl,:,b,it) ...
0343 + ufacabs(1:ie,jl,b)'*hw_temp(ij,:) ...
0344 + vfacabs(1:ie,jl,b)'*hs_temp(ij,:);
0345 if isequal(flu,'A')
0346 IntQ(jl,:,b,it) = IntQ(jl,:,b,it) ...
0347 + ufacabs(1:ie,jl,b)'*qbi_temp(ij,:) ...
0348 + vfacabs(1:ie,jl,b)'*qbj_temp(ij,:);
32ec7e18d0 Dani*0349 end
0350 end
0351 end
0352
d43e8328b9 Dani*0353 % Prepare HT output: Calculate HT by zonal flows and tabulate
0354 % residuals. Also, the multipliative constants (Cp,rho,grav,LhVap) are
0355 % applied here to put moisture and potential temperature fluxes in
0356 % terms of heat transports.
64c5c5fbe0 Dani*0357 SenHT(2:ydim+1,:,it,4) = sum( IntV(:,:,:,it) ...
d43e8328b9 Dani*0358 .* IntT(:,:,:,it) ...
0359 ./ IntM(:,:,:,it),2);
64c5c5fbe0 Dani*0360 SenHT(:,:,it,3) = SenHT(:,:,it,1) - SenHT(:,:,it,2);
0361 SenHT(:,:,it,5) = SenHT(:,:,it,2) - SenHT(:,:,it,4);
d43e8328b9 Dani*0362 if isequal(flu,'O')
0363 SenHT(:,:,it,6) = DiffKh.*SenHT(:,:,it,6);
0364 SenHT(:,:,it,:) = (CpO*RhoO)*SenHT(:,:,it,:);
0365 elseif isequal(flu,'A')
0366 SenHT(:,:,it,:) = (CpA./grav).*SenHT(:,:,it,:);
64c5c5fbe0 Dani*0367 LatHT(2:ydim+1,:,it,4) = sum( IntV(:,:,:,it) ...
d43e8328b9 Dani*0368 .* IntQ(:,:,:,it) ...
64c5c5fbe0 Dani*0369 ./ IntM(:,:,:,it),2);
0370 LatHT(:,:,it,3) = LatHT(:,:,it,1) - LatHT(:,:,it,2);
0371 LatHT(:,:,it,5) = LatHT(:,:,it,2) - LatHT(:,:,it,4);
d43e8328b9 Dani*0372 LatHT(:,:,it,:) = (LhVap./grav).*LatHT(:,:,it,:);
64c5c5fbe0 Dani*0373 end
32ec7e18d0 Dani*0374
0375 % Implied surface heat flux from heat transports (implied heating).
d43e8328b9 Dani*0376 % Tabulated as the difference in heat transports between two broken
0377 % lines divided by the zonal band area.
32ec7e18d0 Dani*0378 mskG = reshape(mskG,(ydim+1)*(1+nBas),1);
0379 I = find(mskG==0);
0380 mskG = reshape(mskG,ydim+1,1+nBas);
0381 var = zeros(ydim+1,1+nBas);
0382 for n=1:min(nout,6),
64c5c5fbe0 Dani*0383 varT = SenHT([2:ydim+2],:,it,n) ...
0384 - SenHT([1:ydim+1],:,it,n);
0385 varT = reshape(varT,(ydim+1)*(1+nBas),1); varT(I)=NaN;
0386 varT = reshape(varT,ydim+1,1+nBas);
0387 SenSF([1:ydim+1],:,it,n) = varT./AreaZon;
d43e8328b9 Dani*0388 if isequal(flu,'A')
0389 for n=1:min(nout,6)
0390 varQ = LatHT([2:ydim+2],:,it,n) ...
0391 - LatHT([1:ydim+1],:,it,n);
0392 varQ = reshape(varQ,(ydim+1)*(1+nBas),1); varQ(I)=NaN;
0393 varQ = reshape(varQ,ydim+1,1+nBas);
0394 LatSF([1:ydim+1],:,it,n) = varQ./AreaZon;
0395 end
64c5c5fbe0 Dani*0396 end
32ec7e18d0 Dani*0397 end
0398 end
0399
0400
0401 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
d43e8328b9 Dani*0402 % Assign outputs, put in units of PW %
32ec7e18d0 Dani*0403 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0404
64c5c5fbe0 Dani*0405 SenHT = SenHT*1e-15;
32ec7e18d0 Dani*0406
d43e8328b9 Dani*0407 HT.time = time;
0408 HT.SenHT = SenHT;
0409 HT.SenSF = SenSF;
0410 HT.ylatHT = ylatHT;
0411 HT.ylatSF = ylatSF;
0412 HT.AreaZon = AreaZon;
64c5c5fbe0 Dani*0413
0414 if isequal(flu,'A')
0415 LatHT = LatHT*1e-15;
d43e8328b9 Dani*0416 HT.LatHT = LatHT;
0417 HT.LatSF = LatSF;
64c5c5fbe0 Dani*0418 end
32ec7e18d0 Dani*0419
0420
0421 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0422 % Stuff that might need to be added back in later (for basins) %
0423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0424
0425 % Block for constructing mskG for different basins:
0426 % if nBas > 0,
0427 % mskBc=rdda([Rac,'maskC_bas.bin'],[6*nc*nc 3],1,'real*4','b');
0428 % mskBw=rdda([Rac,'maskW_bas.bin'],[6*nc*nc 3],1,'real*4','b');
0429 % mskBs=rdda([Rac,'maskS_bas.bin'],[6*nc*nc 3],1,'real*4','b');
0430 % if nBas==2,
0431 % mskBc(:,2)=mskBc(:,2)+mskBc(:,3);
0432 % mskBw(:,2)=mskBw(:,2)+mskBw(:,3);
0433 % mskBs(:,2)=mskBs(:,2)+mskBs(:,3);
0434 % mskBc=min(1,mskBc); mskBw=min(1,mskBw); mskBs=min(1,mskBs);
0435 % end
0436 % %- load: np_Sep, ij_Sep, tp_Sep:
0437 % sep_lineF=[Rac,'sepBas_cs32_60'];
0438 % load(sep_lineF);
0439 % fprintf([' + bassin mask & Sep.line:',sep_lineF,' \n']);
0440 % %- compute mask for each bassin: -----------
0441 % kMsep=1;
0442 % if kMsep,
0443 % mskW=1+min(1,ceil(hw(:,1)));
0444 % mskS=1+min(1,ceil(hs(:,1)));
0445 % for b=1:nBas,
0446 % bs=b; b1=1+bs; b2=2+rem(bs,nBas);
0447 % if nBas == 2, bs=b+b-1; b1=2; b2=3 ; end
0448 % for j=1:ydim+1,
0449 % for i=1:np_Sep(bs,j),
0450 % ij=ij_Sep(bs,j,i); typ=abs(tp_Sep(bs,j,i));
0451 % if typ == 1,
0452 % mskG(j,b1)=mskG(j,b1)*mskW(ij);
0453 % mskG(j,b2)=mskG(j,b2)*mskW(ij);
0454 % elseif typ == 2,
0455 % mskG(j,b1)=mskG(j,b1)*mskS(ij);
0456 % mskG(j,b2)=mskG(j,b2)*mskS(ij);
0457 % end
0458 % end
0459 % end
0460 % end
0461 % mskG=2-min(2,mskG);
0462 % end
0463 % end
0464
d43e8328b9 Dani*0465 % % Masking for different basin.
0466 % mskZ = ones(ydim+2,1+nBas);
0467 % mskZ([2:ydim+1],:) = mskG([1:ydim],:) + mskG([2:ydim+1],:);
0468 % mskZ = reshape(min(mskZ,1),(ydim+2)*(1+nBas),1);
0469 % I = find(mskZ == 0);
0470 % mskZ = reshape(mskZ,(ydim+2),1+nBas);
0471 % for n=1:min(6,nout)
0472 % var=SenHT(:,:,n);
0473 % var=reshape(var,(ydim+2)*(1+nBas),1); var(I) = NaN;
0474 % var=reshape(var,(ydim+2),1+nBas); SenHT(:,:,n) = var;
0475 % end
0476
32ec7e18d0 Dani*0477 % Makes sure that there are no zeros in area AreaZon:
0478 % - set AreaZon to 1 if = zero
0479 % AreaZon=reshape(AreaZon,(ydim+1)*(1+nBas),1);
0480 % [I]=find(AreaZon==0); AreaZon(I)=1;
92699e6cf5 Dani*0481 % AreaZon=reshape(AreaZon,ydim+1,1+nBas);
