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.2 |
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.2 |
hFacC=reshape(hFacC,nPg,nr_hfac); mskC=reshape(mskC,nPg,nr_hfac); |
119 |
enderton |
1.1 |
a6c=reshape(arc,nPg,1); |
120 |
|
|
if nBas > 0, |
121 |
|
|
%- standard 3. Sector definition: |
122 |
|
|
%mskB=rdda([Rac,'maskC_bas.bin'],[nPg nBas],1,'real*4','b'); |
123 |
|
|
%- with Mediter in Indian Ocean Sector : |
124 |
|
|
mskB=rdda([Rac,'maskC_MdI.bin'],[nPg nBas],1,'real*4','b'); |
125 |
|
|
elseif nBas < 0, |
126 |
|
|
|
127 |
|
|
%landFrc=rdda([dir,'landFrc.cpl_FM.bin'],[nPg 1],1,'real*8','b'); |
128 |
|
|
landFrc=rdda([dir,'landFrc_cs32.zero.bin'],[nPg 1],1,'real*8','b'); |
129 |
|
|
|
130 |
|
|
mskB=zeros(nPg,abs(nBas)); |
131 |
|
|
mskB(:,1)=landFrc; |
132 |
|
|
if nBas==-2, mskB(:,2)=1-landFrc ; end |
133 |
|
|
end |
134 |
|
|
end |
135 |
|
|
|
136 |
|
|
if kpr == 1 | kpr == 2, |
137 |
|
|
if nBas < 0, mskZb=zeros(ny,1+abs(nBas)); else mskZb=zeros(ny,nr,1+abs(nBas)); end |
138 |
|
|
for nb=1:1+abs(nBas), |
139 |
|
|
if nb == 1, vv1=a6c ; else vv1=a6c.*mskB(:,nb-1); end |
140 |
|
|
if nBas < 0, |
141 |
|
|
var=vv1; |
142 |
|
|
for j=1:ny, |
143 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
144 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
145 |
|
|
mskZb(j,nb)=sum(vvLoc.*var(ijLoc)); |
146 |
|
|
end |
147 |
|
|
else |
148 |
|
|
for k=1:nr, |
149 |
|
|
var=vv1.*hFacC(:,k); |
150 |
|
|
for j=1:ny, |
151 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
152 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
153 |
|
|
mskZb(j,k,nb)=sum(vvLoc.*var(ijLoc)); |
154 |
|
|
end |
155 |
|
|
end |
156 |
|
|
end |
157 |
|
|
end |
158 |
|
|
|
159 |
|
|
if kwr == 1, |
160 |
|
|
if nBas < 0, namH='landFrc'; else namH='hFacC'; end |
161 |
|
|
namfil=[rac,namH,sufx02]; |
162 |
|
|
fid=fopen(namfil,'w','b'); fwrite(fid,mskZb,'real*8'); fclose(fid); |
163 |
|
|
%fprintf([' -> write mskZb to file: ',namfil,'\n']); |
164 |
|
|
end |
165 |
|
|
|
166 |
|
|
elseif krd == 1, |
167 |
|
|
if nBas < 0, |
168 |
|
|
namH='landFrc'; |
169 |
|
|
namfil=[rac,namH,sufx02]; |
170 |
|
|
mskZb=rdda(namfil,[ny 1+abs(nBas)],1,'real*8','b'); |
171 |
|
|
else |
172 |
|
|
namH='hFacC'; |
173 |
|
|
namfil=[rac,namH,sufx02]; |
174 |
|
|
mskZb=rdda(namfil,[ny nr 1+abs(nBas)],1,'real*8','b'); |
175 |
|
|
end |
176 |
|
|
%fprintf([' <- read mskZb from file: ',namfil,'\n']); |
177 |
|
|
end |
178 |
|
|
|
179 |
|
|
if kpr > 0 , |
180 |
|
|
%- area Zon.Av: |
181 |
|
|
arZ=zeros(ny,1); |
182 |
|
|
for j=1:ny, |
183 |
|
|
ijLoc=ijZav(j,1:npZav(j)); vvLoc=alZav(j,1:npZav(j))'; |
184 |
|
|
arZ(j)=sum(vvLoc.*a6c(ijLoc)); |
185 |
|
|
end |
186 |
|
|
fldZb=zeros(ny,nr,1+abs(nBas),nti); |
187 |
|
|
for nt=1:nti, |
188 |
|
|
fld1t=reshape(fld(:,:,:,nt),nPg,nr); |
189 |
|
|
for nb=1:1+abs(nBas), |
190 |
|
|
if nb == 1, vv1=a6c ; else vv1=a6c.*mskB(:,nb-1); end |
191 |
|
|
for k=1:nr, |
192 |
|
|
if nBas < 0, mskLoc=mskZb(:,nb); else mskLoc=mskZb(:,k,nb); end |
193 |
|
|
var=vv1.*hFacC(:,k); |
194 |
|
|
var=var.*fld1t(:,k); |
195 |
|
|
for j=1:ny, |
196 |
|
|
if mskLoc(j) > frcZmn*arZ(j), |
197 |
|
|
ijLoc=ijZav(j,1:npZav(j)); |
198 |
|
|
vvLoc=alZav(j,1:npZav(j))'; |
199 |
|
|
fldZb(j,k,nb,nt)=sum(vvLoc.*var(ijLoc))/mskLoc(j); |
200 |
|
|
else |
201 |
|
|
fldZb(j,k,nb,nt)=NaN; |
202 |
|
|
end |
203 |
|
|
end |
204 |
|
|
end |
205 |
|
|
end |
206 |
|
|
end |
207 |
|
|
end |
208 |
|
|
|
209 |
|
|
return |