| 1 |
enderton |
1.1 |
function [fldZb,mskZb,ylat]=... |
| 2 |
|
|
calc_ZonAv_CS(fld,kpr,kwr,nBas,xcs,ycs,xcg,ycg,arc,datadir,hFacC); |
| 3 |
|
|
if (nargin < 4), nBas=0; end |
| 4 |
|
|
krd=1; |
| 5 |
|
|
|
| 6 |
|
|
%kpr=1 ; kwr=1; krd=1 ; %- comput weight + zonAv_mask + zonAv(fld) |
| 7 |
|
|
%kpr=2 ; kwr=1; krd=1 ; %- read weight but comput zonAv_mask + zonAv(fld); |
| 8 |
|
|
%kpr=3 ; kwr=1; krd=1 ; %- read weight & zonAv_mask but comput zonAv(fld); |
| 9 |
|
|
%kpr=3 ; kwr=1; krd=0 ; %- read nothing except field, comput zonAv(fld); |
| 10 |
|
|
%nBas=0: Global ZonAv only ; nBas=-1: atmos-> Glob+Land ; =-2: idem +Ocean ; |
| 11 |
|
|
%nBas>0 : ocean: Glob+(nbas)bassins |
| 12 |
|
|
|
| 13 |
|
|
%nBas=-2; |
| 14 |
|
|
%kpr=3 ; kwr=1; krd=1; |
| 15 |
|
|
% dim3=0; |
| 16 |
|
|
|
| 17 |
|
|
% Rac <- bassin mask , cs_grid & ZonAv weight |
| 18 |
|
|
% dir <- hFacC & zonAv_mask |
| 19 |
|
|
Rac = '/Users/enderton/Research/CoupledGFDExperiments/cs_grid/'; |
| 20 |
|
|
dir = datadir; |
| 21 |
|
|
|
| 22 |
|
|
%fprintf('calc_ZonAv_CS:'); |
| 23 |
|
|
dims=size(fld); |
| 24 |
|
|
if length(dims) > 4, |
| 25 |
|
|
fprintf('error: too many dimensions!\n'); |
| 26 |
|
|
dims; |
| 27 |
|
|
return |
| 28 |
|
|
end |
| 29 |
|
|
if length(dims) < 4, nti=1 ; else nti=dims(4); end |
| 30 |
|
|
if length(dims) < 3, nr=1 ; else nr=dims(3); end |
| 31 |
|
|
fld=reshape(fld,dims(1),dims(2),nr,nti); |
| 32 |
|
|
%fprintf('assume 3rd dim = nr = %i \n',nr); |
| 33 |
|
|
|
| 34 |
|
|
%- define latitude Axis (with regular spacing) to average to: |
| 35 |
|
|
ny=64 ; yyMx=90; |
| 36 |
|
|
yyMn=-yyMx ; ddy=(yyMx-yyMn)/ny; |
| 37 |
|
|
ylat=yyMn+([1:ny]-0.5)*ddy; |
| 38 |
|
|
|
| 39 |
|
|
%- define minimum area fraction (relative to full latitude band): |
| 40 |
|
|
frcZmn=0.05; |
| 41 |
|
|
sufx02=['_Zav_',int2str(ny)]; |
| 42 |
|
|
|
| 43 |
|
|
%------------------------------------ |
| 44 |
|
|
rac=Rac ; |
| 45 |
|
|
[ncx nc]=size(arc); nPg=ncx*nc; |
| 46 |
|
|
|
| 47 |
|
|
if kpr == 1, |
| 48 |
|
|
%- compute weight used for zonal average |
| 49 |
|
|
% (no bassin, no mask at this level) |
| 50 |
|
|
|
| 51 |
|
|
x6c=reshape(xcs,nPg,1); |
| 52 |
|
|
y6c=reshape(ycs,nPg,1); |
| 53 |
|
|
|
| 54 |
|
|
%--- |
| 55 |
|
|
nMax=ncx*3; |
| 56 |
|
|
ijzav=zeros(ny,nMax); |
| 57 |
|
|
alzav=zeros(ny,nMax); |
| 58 |
|
|
npzav=zeros(ny,1); |
| 59 |
|
|
|
| 60 |
|
|
|
| 61 |
|
|
for ij=1:nPg, |
| 62 |
|
|
yloc=y6c(ij); |
| 63 |
|
|
del=(yloc-ylat(1))/ddy; |
| 64 |
|
|
jz=1+floor(del); del=rem(del,1); |
| 65 |
|
|
if jz < 1, |
| 66 |
|
|
jj=jz+1; |
| 67 |
|
|
n=npzav(jj)+1; |
| 68 |
|
|
ijzav(jj,n)=ij; |
| 69 |
|
|
alzav(jj,n)=1; |
| 70 |
|
|
npzav(jj)=n; |
| 71 |
|
|
elseif jz >= ny, |
| 72 |
|
|
jj=jz; |
| 73 |
|
|
n=npzav(jj)+1; |
| 74 |
|
|
ijzav(jj,n)=ij; |
| 75 |
|
|
alzav(jj,n)=1; |
| 76 |
|
|
npzav(jj)=n; |
| 77 |
|
|
else |
| 78 |
|
|
jj=jz; |
| 79 |
|
|
n=npzav(jj)+1; |
| 80 |
|
|
ijzav(jj,n)=ij; |
| 81 |
|
|
alzav(jj,n)=1-del; |
| 82 |
|
|
npzav(jj)=n; |
| 83 |
|
|
jj=jz+1; |
| 84 |
|
|
n=npzav(jj)+1; |
| 85 |
|
|
ijzav(jj,n)=ij; |
| 86 |
|
|
alzav(jj,n)=del; |
| 87 |
|
|
npzav(jj)=n; |
| 88 |
|
|
end |
| 89 |
|
|
end |
| 90 |
|
|
|
| 91 |
|
|
%-- reduce size: |
| 92 |
|
|
[NbMx,jM]=max(npzav); |
| 93 |
|
|
%fprintf('Total Nb: %i \n',sum(npzav)); |
| 94 |
|
|
%fprintf('Max Nb for lat(j=%i)=%3.2f : %i \n',jM,ylat(jM),NbMx); |
| 95 |
|
|
|
| 96 |
|
|
ijZav=ijzav(:,1:NbMx); |
| 97 |
|
|
alZav=alzav(:,1:NbMx); |
| 98 |
|
|
npZav=npzav; |
| 99 |
|
|
if kwr == 1, |
| 100 |
|
|
namfil=['zonAvWeigth_cs32_',int2str(ny)]; |
| 101 |
|
|
save(namfil,'npZav','alZav','ijZav'); |
| 102 |
|
|
%fprintf(['save npZav alZav & ijZav to file: ',namfil,'.mat \n']); |
| 103 |
|
|
end |
| 104 |
|
|
else |
| 105 |
|
|
if krd > 0, |
| 106 |
|
|
namfil=['zonAvWeigth_cs32_',int2str(ny)]; |
| 107 |
|
|
fprintf(['read npZav alZav & ijZav from file: ',namfil,'\n']); |
| 108 |
|
|
load([Rac,namfil]); |
| 109 |
|
|
end |
| 110 |
|
|
end |
| 111 |
|
|
|
| 112 |
|
|
if kpr > 0 & krd > 0, |
| 113 |
|
|
rac=dir; |
| 114 |
|
|
%if nBas < 0, hFacC=ones(ncx,nc,dim3); else hFacC=rdmds([rac,'hFacC']); end |
| 115 |
enderton |
1.3 |
%nr_hfac=size(hFacC,3); |
| 116 |
enderton |
1.1 |
if nBas < 0, hFacC=ones(ncx,nc,nr); end%else hFacC=rdmds([rac,'hFacC']); end |
| 117 |
|
|
mskC=min(1,hFacC); mskC=ceil(mskC); |
| 118 |
enderton |
1.3 |
%hFacC=reshape(hFacC,nPg,nr_hfac); mskC=reshape(mskC,nPg,nr_hfac); |
| 119 |
|
|
if nr == 1, hFacC = hFacC(:,:,1); mskC = mskC(:,:,1); end % Cludge by Daniel. |
| 120 |
|
|
hFacC=reshape(hFacC,nPg,nr); mskC=reshape(mskC,nPg,nr); |
| 121 |
enderton |
1.1 |
a6c=reshape(arc,nPg,1); |
| 122 |
|
|
if nBas > 0, |
| 123 |
|
|
%- standard 3. Sector definition: |
| 124 |
|
|
%mskB=rdda([Rac,'maskC_bas.bin'],[nPg nBas],1,'real*4','b'); |
| 125 |
|
|
%- with Mediter in Indian Ocean Sector : |
| 126 |
|
|
mskB=rdda([Rac,'maskC_MdI.bin'],[nPg nBas],1,'real*4','b'); |
| 127 |
|
|
elseif nBas < 0, |
| 128 |
|
|
|
| 129 |
|
|
%landFrc=rdda([dir,'landFrc.cpl_FM.bin'],[nPg 1],1,'real*8','b'); |
| 130 |
|
|
landFrc=rdda([dir,'landFrc_cs32.zero.bin'],[nPg 1],1,'real*8','b'); |
| 131 |
|
|
|
| 132 |
|
|
mskB=zeros(nPg,abs(nBas)); |
| 133 |
|
|
mskB(:,1)=landFrc; |
| 134 |
|
|
if nBas==-2, mskB(:,2)=1-landFrc ; end |
| 135 |
|
|
end |
| 136 |
|
|
end |
| 137 |
|
|
|
| 138 |
|
|
if kpr == 1 | kpr == 2, |
| 139 |
|
|
if nBas < 0, mskZb=zeros(ny,1+abs(nBas)); else mskZb=zeros(ny,nr,1+abs(nBas)); end |
| 140 |
|
|
for nb=1:1+abs(nBas), |
| 141 |
|
|
if nb == 1, vv1=a6c ; else vv1=a6c.*mskB(:,nb-1); end |
| 142 |
|
|
if nBas < 0, |
| 143 |
|
|
var=vv1; |
| 144 |
|
|
for j=1:ny, |
| 145 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
| 146 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
| 147 |
|
|
mskZb(j,nb)=sum(vvLoc.*var(ijLoc)); |
| 148 |
|
|
end |
| 149 |
|
|
else |
| 150 |
|
|
for k=1:nr, |
| 151 |
|
|
var=vv1.*hFacC(:,k); |
| 152 |
|
|
for j=1:ny, |
| 153 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
| 154 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
| 155 |
|
|
mskZb(j,k,nb)=sum(vvLoc.*var(ijLoc)); |
| 156 |
|
|
end |
| 157 |
|
|
end |
| 158 |
|
|
end |
| 159 |
|
|
end |
| 160 |
|
|
|
| 161 |
|
|
if kwr == 1, |
| 162 |
|
|
if nBas < 0, namH='landFrc'; else namH='hFacC'; end |
| 163 |
|
|
namfil=[rac,namH,sufx02]; |
| 164 |
|
|
fid=fopen(namfil,'w','b'); fwrite(fid,mskZb,'real*8'); fclose(fid); |
| 165 |
|
|
%fprintf([' -> write mskZb to file: ',namfil,'\n']); |
| 166 |
|
|
end |
| 167 |
|
|
|
| 168 |
|
|
elseif krd == 1, |
| 169 |
|
|
if nBas < 0, |
| 170 |
|
|
namH='landFrc'; |
| 171 |
|
|
namfil=[rac,namH,sufx02]; |
| 172 |
|
|
mskZb=rdda(namfil,[ny 1+abs(nBas)],1,'real*8','b'); |
| 173 |
|
|
else |
| 174 |
|
|
namH='hFacC'; |
| 175 |
|
|
namfil=[rac,namH,sufx02]; |
| 176 |
|
|
mskZb=rdda(namfil,[ny nr 1+abs(nBas)],1,'real*8','b'); |
| 177 |
|
|
end |
| 178 |
|
|
%fprintf([' <- read mskZb from file: ',namfil,'\n']); |
| 179 |
|
|
end |
| 180 |
|
|
|
| 181 |
|
|
if kpr > 0 , |
| 182 |
|
|
%- area Zon.Av: |
| 183 |
|
|
arZ=zeros(ny,1); |
| 184 |
|
|
for j=1:ny, |
| 185 |
|
|
ijLoc=ijZav(j,1:npZav(j)); vvLoc=alZav(j,1:npZav(j))'; |
| 186 |
|
|
arZ(j)=sum(vvLoc.*a6c(ijLoc)); |
| 187 |
|
|
end |
| 188 |
|
|
fldZb=zeros(ny,nr,1+abs(nBas),nti); |
| 189 |
|
|
for nt=1:nti, |
| 190 |
|
|
fld1t=reshape(fld(:,:,:,nt),nPg,nr); |
| 191 |
|
|
for nb=1:1+abs(nBas), |
| 192 |
|
|
if nb == 1, vv1=a6c ; else vv1=a6c.*mskB(:,nb-1); end |
| 193 |
|
|
for k=1:nr, |
| 194 |
|
|
if nBas < 0, mskLoc=mskZb(:,nb); else mskLoc=mskZb(:,k,nb); end |
| 195 |
|
|
var=vv1.*hFacC(:,k); |
| 196 |
|
|
var=var.*fld1t(:,k); |
| 197 |
|
|
for j=1:ny, |
| 198 |
|
|
if mskLoc(j) > frcZmn*arZ(j), |
| 199 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
| 200 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
| 201 |
|
|
fldZb(j,k,nb,nt)=sum(vvLoc.*var(ijLoc))/mskLoc(j); |
| 202 |
|
|
else |
| 203 |
|
|
fldZb(j,k,nb,nt)=NaN; |
| 204 |
|
|
end |
| 205 |
|
|
end |
| 206 |
|
|
end |
| 207 |
|
|
end |
| 208 |
|
|
end |
| 209 |
|
|
end |
| 210 |
|
|
|
| 211 |
|
|
return |