Warning, /utils/matlab/cs_grid/calcZonalAvgCube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 6b4f212f on 2024-06-05 16:55:17 UTC
b3d495b63f Dani*0001 function [fldzon,mskzon,ylat,areazon] = calcZonalAvgCube(fld,ny,yc,ar,hc);
0002
0003 % [fldzon,mskzon,ylat,areazon] = calcZonalAvgCube(fld,ny,yc,ar,hc);
0004 %
6b4f212f67 Jean*0005 % Input arguments:
b3d495b63f Dani*0006 % fld Data to zonally average ([6*n,n,nr,nt], [6*n,n,nr], [6*n,n])
0007 % ny Number of evenly spaces latitude lines to zonally average onto
0008 % yc Cell centered latitude
0009 % ar Cell area
0010 % hc HFacC
0011 %
6b4f212f67 Jean*0012 % Output arguments:
b3d495b63f Dani*0013 % fldzon Zonally averaged field ([ny,nr,nt], [ny,nr], [ny])
0014 % mskzon Mask for zonal slice
0015 % ylat Latitude coordinates for zonal average data
0016 % areazon Area associated with zonal band
0017 %
0018 % Description:
0019 % Calculates zonal average of "fld".
0020 %
0021 % Original author: JMC
0022 % Modifications: Daniel Enderton
0023
0024 % Dimension information.
0025 nBas = 0;
0026 dims = size(fld);
f7a0cd83c0 Jean*0027 %if length(dims) > 4,
0028 % error(['To many dimensions (Max 4): ',mat2str(dims)]);
0029 %end
6b4f212f67 Jean*0030 n1h=dims(1); n2h=dims(2);
0031 if n1h == n2h*6, nc=n2h; elseif n1h*6 == n2h, nc=n1h;
0032 else
0033 error(['Error in CS-dimensions: ',...
0034 'n1h,n2h = ',int2str(n1h),', ',int2str(n2h)])
0035 end
0036 nr=1; nt=1; nDim=4;
f7a0cd83c0 Jean*0037 if length(dims) > 2,
0038 if length(size(hc)) == 3, nr=size(hc,3); end
0039 if dims(3) ~= nr, nr=1; nDim=3; end
0040 nt=prod(dims(3:end))/nr;
b3d495b63f Dani*0041 end
0042
0043 % Reshape data.
0044 if nr == 1, hc = hc(:,:,1); end % Account for nr = 1
f7a0cd83c0 Jean*0045 hc = reshape(hc, [6*nc*nc,nr]);
0046 yc = reshape(yc, [6*nc*nc,1] );
0047 ar = reshape(ar, [6*nc*nc,1] );
0048 fld= reshape(fld,[6*nc*nc,nr,nt]);
b3d495b63f Dani*0049
0050 % Define latitude axis (with regular spacing) to average to.
0051 dy=180/ny; ylat=-90+([1:ny]-0.5)*dy;
0052
0053 % Define minimum area fraction (relative to full latitude band):
0054 frcZmn=0.05; % ??? Why
0055
0056 % Compute zonal line weights.
0057 nMax = 6*nc*3; % Max number of cells near zonal line
0058 ijzav = zeros(ny,nMax); % Index of cells near band line (within dy)
0059 alzav = zeros(ny,nMax); % Area weights for cells near band line
0060 npzav = zeros(ny,1); % Number of cells close to band line
0061 for ij = 1:6*nc*nc % Loop over all cells
0062 del=(yc(ij)-ylat(1))/dy; % Latitude difference from bottom
0063 jz=1+floor(del); % Number of zonal lines away (round up)
0064 del=rem(del,1); % Fracional distance from upper line
0065 if jz < 1
0066 n = npzav(jz+1) + 1; % Contribution from cell NORTH of line
0067 ijzav(jz+1,n) = ij;
0068 alzav(jz+1,n) = 1; % Weight = 1 because all area accounted for
6b4f212f67 Jean*0069 npzav(jz+1) = n;
b3d495b63f Dani*0070 elseif jz >= ny
0071 n = npzav(jz) + 1; % Contribution from cell NORTH of line
0072 ijzav(jz,n) = ij;
0073 alzav(jz,n) = 1; % Weight = 1 because all area accounted for
0074 npzav(jz) = n;
0075 else
0076 n = npzav(jz) + 1; % Contribution from cell SOUTH of line
0077 ijzav(jz,n) = ij;
0078 alzav(jz,n) = 1 - del;
0079 npzav(jz) = n;
0080 n = npzav(jz+1) + 1; % Contribution from cell NORTH of line
0081 ijzav(jz+1,n) = ij;
0082 alzav(jz+1,n) = del;
0083 npzav(jz+1) = n;
0084 end
0085 end
0086 [NbMx,jM] = max(npzav); % Reduce size
0087 ijzav = ijzav(:,1:NbMx);
0088 alzav = alzav(:,1:NbMx);
0089
0090 mskC = ceil(hc); % Cell centered mask: 1 = any land, 0 = no land
0091
0092 % Basin stuff:
0093 % if nBas > 0
0094 % % standard 3. Sector definition:
0095 % % mskB=rdda([rac,'maskC_bas.bin'],[6*nc*nc nBas],1,'real*4','b');
0096 % % with Mediter in Indian Ocean Sector :
0097 % mskB = rdda([rac,'maskC_MdI.bin'],[6*nc*nc nBas],1,'real*4','b');
0098 % elseif nBas < 0
0099 % landFrc = rdda([rac,'landFrc_cs32.zero.bin'],[6*nc*nc 1],1,'real*8','b');
0100 % mskB = zeros(6*nc*nc,abs(nBas));
0101 % mskB(:,1) = landFrc;
0102 % if nBas == -2, mskB(:,2) = 1 - landFrc; end
0103 % end
0104
0105 % Basin stuff:
0106 if nBas < 0, mskzon = zeros(ny,1+abs(nBas));
0107 else, mskzon = zeros(ny,nr,1+abs(nBas)); end
0108 for nb = 1:1+abs(nBas)
0109 if nb == 1, vv1=ar; else vv1=ar.*mskB(:,nb-1); end
0110 if nBas < 0,
0111 var=vv1;
0112 for j=1:ny,
0113 ijLoc=ijzav(j,1:npzav(j));
0114 vvLoc=alzav(j,1:npzav(j))';
0115 mskzon(j,nb)=sum(vvLoc.*var(ijLoc));
0116 end
0117 else
0118 for k=1:nr,
6b4f212f67 Jean*0119 var=vv1.*hc(:,k);
b3d495b63f Dani*0120 for j=1:ny,
0121 ijLoc=ijzav(j,1:npzav(j));
0122 vvLoc=alzav(j,1:npzav(j))';
0123 mskzon(j,k,nb)=sum(vvLoc.*var(ijLoc));
0124 end
0125 end
0126 end
0127 end
0128
0129 % Area associated with zonal lines:
0130 areazon=zeros(ny,1);
0131 for j=1:ny,
0132 ijLoc = ijzav(j,1:npzav(j));
0133 vvLoc = alzav(j,1:npzav(j))';
0134 areazon(j) = sum(vvLoc.*ar(ijLoc));
0135 end
0136
0137 % Compute zonal average:
0138 fldzon = zeros(ny,nr,1+abs(nBas),nt);
0139 for it = 1:nt,
f7a0cd83c0 Jean*0140 fld1t = fld(:,:,it);
b3d495b63f Dani*0141 for nb = 1:1+abs(nBas),
0142 if nb == 1, vv1 = ar;
0143 else vv1 = ar.*mskB(:,nb-1); end
0144 for k = 1:nr,
0145 if nBas < 0, mskLoc=mskzon(:,nb);
0146 else mskLoc=mskzon(:,k,nb); end
0147 var = vv1.*hc(:,k);
0148 var = var.*fld1t(:,k);
0149 for j = 1:ny
0150 if mskLoc(j) > frcZmn*areazon(j),
0151 ijLoc=ijzav(j,1:npzav(j));
0152 vvLoc=alzav(j,1:npzav(j))';
0153 fldzon(j,k,nb,it)=sum(vvLoc.*var(ijLoc))/mskLoc(j);
0154 else
0155 fldzon(j,k,nb,it)=NaN;
0156 end
0157 end
0158 end
0159 end
0160 end
f7a0cd83c0 Jean*0161 if length(dims) > nDim,
0162 fldzon=reshape(fldzon,[ny nr 1+abs(nBas) dims(nDim:end)]);
0163 end
0164 fldzon = squeeze(fldzon);
