| 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 |
%nr_hfac=size(hFacC,3); |
| 116 |
if nBas < 0, hFacC=ones(ncx,nc,nr); end%else hFacC=rdmds([rac,'hFacC']); end |
| 117 |
mskC=min(1,hFacC); mskC=ceil(mskC); |
| 118 |
%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 |
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 |
Parent Directory
| 
Revision Log
| 
Revision Graph