/[MITgcm]/MITgcm_contrib/enderton/Diagnostics/DiagUtility/calc_ZonAv_CS.m
ViewVC logotype

Annotation of /MITgcm_contrib/enderton/Diagnostics/DiagUtility/calc_ZonAv_CS.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.2 - (hide annotations) (download)
Thu Aug 18 15:02:05 2005 UTC (19 years, 11 months ago) by enderton
Branch: MAIN
Changes since 1.1: +2 -2 lines
Handles 2D fields now.

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

  ViewVC Help
Powered by ViewVC 1.1.22