/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/plot_all_10.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/eddy_flux/plot_all_10.m

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


Revision 1.1 - (hide annotations) (download)
Thu Aug 19 05:32:33 2004 UTC (20 years, 11 months ago) by edhill
Branch: MAIN
 o latest analysis commands

1 edhill 1.1 %=======================================================
2     %
3     % $Id: $
4     %
5     % Ed Hill
6     %
7    
8     % The following are the MatLAB commands used to create the various
9     % plots related to eddy fluxes using average velocities and densities
10     % (called bouyancy or "b" in many of the variables) from Dimitris'
11     % "cube_5" integration.
12    
13     % Groups of commands contained within the following "UNUSED" comments
14     % we not used in this analyssis. They are "left over" from the
15     % previous temperature-based calculations and have been kept for
16     % reference purposes only.
17    
18     %------- UNUSED -----------------------
19     % ...Commands...
20     %------- UNUSED -----------------------
21    
22    
23     matlab -nojvm
24     matlab -nojvm -nodisplay
25    
26     %=======================================================
27     % Compute 10-yr averages from the ten 1-yr averages
28    
29     clear all ; close all
30     ne = 510;
31     ntot = ne*ne*6 * 50;
32     by = 1994;
33     ey = 2003;
34     vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
35     'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
36     nm = 1;
37     for nm = 1:size(vnam,1)
38     sum = zeros(ntot,1);
39     eval(['pat = ''' sprintf('%s',vnam(nm,:)) '_tave_%4.4d'';']);
40     disp(sprintf('pat = "%s"',pat));
41     iy = by;
42     for iy = by:ey
43     disp([' opening: ' sprintf(pat,iy)]);
44     fid = fopen(sprintf(pat,iy), 'r', 'ieee-be');
45     a = fread(fid, ntot, 'real*4');
46     fclose(fid);
47     sum = sum + a;
48     end
49     a = sum ./ (ey-by+1);
50     ido = fopen(sprintf([pat '-%4.4d'],by,ey), 'w', 'ieee-be');
51     fwrite(ido, a, 'real*4');
52     fclose(ido);
53     end
54    
55     %=======================================================
56     % Compute the "primes"
57    
58     clear all; close all
59    
60     vnam = [ 'U_'; 'V_'; 'T_'; 'S_'; 'B_'; 'B1'; ...
61     'UT'; 'VT'; 'US'; 'VS'; 'UB'; 'VB' ];
62     nztot = 50;
63     ne = 510;
64     nps = 6 * ne * ne;
65     iz = 2;
66     iz = 1;
67     for iz = 1:nztot
68    
69     disp(sprintf('iz = %d', iz));
70    
71     offset = (iz - 1)*nps*4;
72     nm = 1;
73     for nm = 1:size(vnam,1)
74     eval(['pat = ''' sprintf('%s',vnam(nm,:)) ...
75     '_tave_1994-2003'';']);
76     fid = fopen(pat, 'r', 'ieee-be');
77     fseek(fid, offset, 'bof');
78     a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
79     eval([lower(vnam(nm,:)) ' = a;' ]);
80     fclose(fid);
81     end
82     clear a
83     % surf(u_(:,:,1)), view(2), shading interp, axis equal
84     % surf(v_(:,:,1)), view(2), shading interp, axis equal
85     % surf(t_(:,:,1)), view(2), shading interp, axis equal
86     % s_(find(s_==0.0)) = NaN;
87     % surf(s_(:,:,1)), view(2), shading interp, axis equal
88     % b_(find(b_==0.0)) = NaN;
89     % surf(b_(:,:,1)), view(2), shading interp, axis equal
90     % b1(find(b1==0.0)) = NaN;
91     % surf(b1(:,:,1)), view(2), shading interp, axis equal
92     % surf(ut(:,:,1)), view(2), shading interp, axis equal
93     % surf(vt(:,:,1)), view(2), shading interp, axis equal
94     % surf(us(:,:,1)), view(2), shading interp, axis equal
95     % surf(vs(:,:,1)), view(2), shading interp, axis equal
96     % surf(ub(:,:,1)), view(2), shading interp, axis equal
97     % surf(vb(:,:,1)), view(2), shading interp, axis equal
98    
99     % tmask = double(t_ ~= 0.0);
100    
101     % Fill the (ne+1)*(ne+1) grid with T,S,B values
102     ni = ne; nip1 = ni + 1;
103     nj = ne; njp1 = nj + 1;
104     tcubep1 = zeros( [ nip1 njp1 6 ] );
105     scubep1 = zeros( [ nip1 njp1 6 ] );
106     bcubep1 = zeros( [ nip1 njp1 6 ] );
107     for icf = 1:6
108     tcubep1(2:nip1,2:njp1,icf) = t_(:,:,icf);
109     scubep1(2:nip1,2:njp1,icf) = s_(:,:,icf);
110     bcubep1(2:nip1,2:njp1,icf) = b_(:,:,icf);
111     end
112    
113     % Do the upwind-edge T exchanges
114     tcubep1(1,2:njp1,1) = t_(ni:-1:1,nj,5); % -
115     tcubep1(2:nip1,1,1) = t_(1:ni,nj,6); % -
116     tcubep1(1,2:njp1,2) = t_(ni,1:nj,1); % -
117     tcubep1(2:nip1,1,2) = t_(ni,nj:-1:1,6); % -
118     tcubep1(1,2:njp1,3) = t_(ni:-1:1,nj,1); % -
119     tcubep1(2:nip1,1,3) = t_(1:ni,nj,2); % -
120     tcubep1(1,2:njp1,4) = t_(ni,1:nj,3); % -
121     tcubep1(2:nip1,1,4) = t_(ni,nj:-1:1,2); % -
122     tcubep1(1,2:njp1,5) = t_(ni:-1:1,nj,3); % -
123     tcubep1(2:nip1,1,5) = t_(1:ni,nj,4); % -
124     tcubep1(1,2:njp1,6) = t_(ni,1:nj,5); % -
125     tcubep1(2:nip1,1,6) = t_(ni,nj:-1:1,4); % -
126    
127     % Do the upwind-edge S exchanges
128     scubep1(1,2:njp1,1) = s_(ni:-1:1,nj,5); % -
129     scubep1(2:nip1,1,1) = s_(1:ni,nj,6); % -
130     scubep1(1,2:njp1,2) = s_(ni,1:nj,1); % -
131     scubep1(2:nip1,1,2) = s_(ni,nj:-1:1,6); % -
132     scubep1(1,2:njp1,3) = s_(ni:-1:1,nj,1); % -
133     scubep1(2:nip1,1,3) = s_(1:ni,nj,2); % -
134     scubep1(1,2:njp1,4) = s_(ni,1:nj,3); % -
135     scubep1(2:nip1,1,4) = s_(ni,nj:-1:1,2); % -
136     scubep1(1,2:njp1,5) = s_(ni:-1:1,nj,3); % -
137     scubep1(2:nip1,1,5) = s_(1:ni,nj,4); % -
138     scubep1(1,2:njp1,6) = s_(ni,1:nj,5); % -
139     scubep1(2:nip1,1,6) = s_(ni,nj:-1:1,4); % -
140    
141     % Do the upwind-edge B exchanges
142     bcubep1(1,2:njp1,1) = b_(ni:-1:1,nj,5); % -
143     bcubep1(2:nip1,1,1) = b_(1:ni,nj,6); % -
144     bcubep1(1,2:njp1,2) = b_(ni,1:nj,1); % -
145     bcubep1(2:nip1,1,2) = b_(ni,nj:-1:1,6); % -
146     bcubep1(1,2:njp1,3) = b_(ni:-1:1,nj,1); % -
147     bcubep1(2:nip1,1,3) = b_(1:ni,nj,2); % -
148     bcubep1(1,2:njp1,4) = b_(ni,1:nj,3); % -
149     bcubep1(2:nip1,1,4) = b_(ni,nj:-1:1,2); % -
150     bcubep1(1,2:njp1,5) = b_(ni:-1:1,nj,3); % -
151     bcubep1(2:nip1,1,5) = b_(1:ni,nj,4); % -
152     bcubep1(1,2:njp1,6) = b_(ni,1:nj,5); % -
153     bcubep1(2:nip1,1,6) = b_(ni,nj:-1:1,4); % -
154    
155     % Get T values on the U,V grid points
156     masku = 1.0 - abs(diff(double(tcubep1(2:nip1,:,:) == 0.0),1,2));
157     maskv = 1.0 - abs(diff(double(tcubep1(:,2:nip1,:) == 0.0),1,1));
158     diffu = 0.5*diff(tcubep1(2:nip1,:,:),1,2);
159     diffv = 0.5*diff(tcubep1(:,2:nip1,:),1,1);
160     tonu = t_(:,:,:) + masku.*diffu;
161     tonv = t_(:,:,:) + maskv.*diffv;
162    
163     % Get S values on the U,V grid points
164     diffu = 0.5*diff(scubep1(2:nip1,:,:),1,2);
165     diffv = 0.5*diff(scubep1(:,2:nip1,:),1,1);
166     sonu = s_(:,:,:) + masku.*diffu;
167     sonv = s_(:,:,:) + maskv.*diffv;
168    
169     % Get B values on the U grid points
170     diffu = 0.5*diff(bcubep1(2:nip1,:,:),1,2);
171     diffv = 0.5*diff(bcubep1(:,2:nip1,:),1,1);
172     bonu = b_(:,:,:) + masku.*diffu;
173     bonv = b_(:,:,:) + maskv.*diffv;
174    
175     % Compute the primes
176     uptp = ut - u_ .* tonu;
177     vptp = vt - v_ .* tonv;
178     upsp = us - u_ .* sonu;
179     vpsp = vs - v_ .* sonv;
180     upbp = ub - u_ .* bonu;
181     vpbp = vb - v_ .* bonv;
182    
183     % Write the results
184     for f = [ 't' 's' 'b' ]
185     for uv = [ 'u' 'v' ]
186     % disp([' writing: ' uv f]);
187     ido = fopen(['out/' uv 'p' f 'p_1994-2003'], 'a', 'ieee-be');
188     eval(['fwrite(ido, ' uv 'p' f 'p, ''real*4'');']);
189     fclose(ido);
190     end
191     end
192     % f = [ 't_'; 's_'; 'b_'; 'b1'; 'u_'; 'v_'; ];
193     % for fi = 1:size(f,1)
194     % disp([' writing: ' f(fi,:) ]);
195     % ido = fopen(['out/' f(fi,:) '_1994-2003'], 'a', 'ieee-be');
196     % eval(['fwrite(ido, ' f(fi,:) ', ''real*4'');']);
197     % fclose(ido);
198     % end
199    
200     end
201    
202     %=======================================================
203     % Plot u'B' and v'B' on a Lat-Lon grid
204    
205     clear all, close all
206     tx = 85;
207     ty = 85;
208     nt = 216;
209     cx = 510;
210     cy = 510;
211     nz = 1;
212     ne = 510;
213     nps = ne * ne * 6;
214    
215     % XCS YCS XGS YGS
216     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
217     for in = 1:4
218     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
219     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
220     tx*nt,ty), ne, tx );
221     fclose(uid);
222     eval([lower(fnam(in,:)) ' = phi;']);
223     a = zeros(6*510,510);
224     for i = 1:6
225     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
226     yb = 1; ye = 510;
227     a(xb:xe,yb:ye) = phi(:,:,i);
228     end
229     eval([lower(fnam(in,:)) 's = a;'])
230     end
231     xcs = xcs + -360.0*(xcs > 180.0);
232     xgs = xgs + -360.0*(xgs > 180.0);
233     clear phi a
234    
235     %=======================================================
236     % Calculation rotations for Lat-Lon vector alignment
237    
238     ne = 510;
239     n1 = ne - 1;
240     dux = zeros(size(xg));
241     duy = zeros(size(xg));
242     dvx = zeros(size(xg));
243     dvy = zeros(size(xg));
244     dux(1:n1,:,:) = diff(xg,1,1); dux(ne,:,:) = dux(n1,:,:);
245     dvx(:,1:n1,:) = diff(xg,1,2); dvx(:,ne,:) = dvx(:,n1,:);
246     duy(1:n1,:,:) = diff(yg,1,1); duy(ne,:,:) = duy(n1,:,:);
247     dvy(:,1:n1,:) = diff(yg,1,2); dvy(:,ne,:) = dvy(:,n1,:);
248     dux = dux + 360*double(dux < 180);
249     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
250     duy = duy + 360*double(duy < 180);
251     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
252     dvx = dvx + 360*double(dvx < 180);
253     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
254     dvy = dvy + 360*double(dvy < 180);
255     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
256     llux = dux ./ sqrt(dux.^2 + duy.^2);
257     lluy = duy ./ sqrt(dux.^2 + duy.^2);
258     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
259     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
260    
261     % surf(xg(:,:,2)), view(2), axis equal, shading interp, colorbar
262     % surf(dux(:,:,2)), view(2), axis equal, shading interp, colorbar
263     % surf(duy(:,:,2)), view(2), axis equal, shading interp, colorbar
264    
265     delR = [
266     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
267     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
268     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
269     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
270     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
271     341.50,364.50,387.50,410.50,433.50,456.50 ];
272     R = cumsum(delR) - 0.5*delR;
273     Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
274    
275     ne = 510;
276     nps = ne * ne * 6;
277    
278     clow = [ -1.0 ]; % -0.2;
279     chigh = [ 1.0 ]; % 0.2;
280    
281     iz = 1;
282     iz = 15;
283     iz = 25;
284     for iz = [ 1:50 ]
285    
286     disp(sprintf('iz = %d',iz));
287     offset = (iz - 1)*nps*4;
288     % Load the "primes"
289     for f = [ 't' 's' 'b' ]
290     for uv = [ 'u' 'v' ]
291     % disp([' reading: ' uv f]);
292     fid = fopen(['out/' uv 'p' f 'p_1994-2003'],'r','ieee-be');
293     fseek(fid, offset, 'bof');
294     a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
295     fclose(fid);
296     eval([uv 'p' f 'p = a;']);
297     end
298     end
299    
300     % Get LatLon-oriented vectors at the CubeSphere points
301     lluptp = uptp .* llux + vptp .* llvx;
302     llvptp = uptp .* lluy + vptp .* llvy;
303     llupsp = upsp .* llux + vpsp .* llvx;
304     llvpsp = upsp .* lluy + vpsp .* llvy;
305     llupbp = upbp .* llux + vpbp .* llvx;
306     llvpbp = upbp .* lluy + vpbp .* llvy;
307     % surf(xg(:,:,2),yg(:,:,2),lluptp(:,:,2)), view(2),shading interp,colorbar
308     % surf(xg(:,:,2),yg(:,:,2),llupbp(:,:,2)), view(2),shading interp,colorbar
309    
310     % save the LatLon-rotated "primes"
311     for f = [ 't' 's' 'b' ]
312     for uv = [ 'u' 'v' ]
313     % disp([' writing: ' uv f]);
314     fid = fopen(['out/ll' uv 'p' f 'p_1994-2003'], ...
315     'a','ieee-be');
316     eval(['fwrite(fid,ll' uv 'p' f 'p,''real*4'');']);
317     fclose(fid);
318     end
319     end
320    
321     % Look at Lat-Lon projections of the "primes"
322     clim = [ 0.6 0.15 0.0015 ];
323     fi = 0;
324     for f = [ 't' 's' 'b' ]
325    
326     fi = fi + 1;
327     eval(['llu = up' f 'p;']);
328     eval(['llv = vp' f 'p;']);
329    
330     pllu = zeros(6*510,510);
331     pllv = zeros(6*510,510);
332     for i = 1:6
333     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
334     yb = 1; ye = 510;
335     pllu(xb:xe,yb:ye) = llu(:,:,i);
336     pllv(xb:xe,yb:ye) = llv(:,:,i);
337     end
338     pllu(find(pllu == 0.0)) = NaN;
339     pllv(find(pllv == 0.0)) = NaN;
340     shift=-1;
341     for uv = [ 'u' 'v' ]
342     eval(['a = pll' uv ';']);
343     clo = -clim(fi);
344     chi = clim(fi);
345     grph_CS(a,xcs,ycs,xgs,ygs,[clo],[chi],shift);
346     title([ uv '''' upper(f) ''' at ' ...
347     sprintf('%g',R(iz)) ...
348     'm depth on the 510x510x6 cubesphere ' ...
349     'for 1994--2003 ["cube5"]']);
350     print('-painters', '-dpng', '-r150', ...
351     [ uv 'p' upper(f) 'p_' ...
352     sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])
353     end
354    
355     end
356     end
357    
358    
359     %=======================================================
360     % Calculate : (v'B')/(dB/dz)
361    
362     bid = fopen('B__tave_1994-2003', 'r', 'ieee-be');
363     bm1id = fopen('B1_tave_1994-2003', 'r', 'ieee-be');
364     upbpid = fopen('out/upbp_1994-2003', 'r', 'ieee-be');
365     vpbpid = fopen('out/vpbp_1994-2003', 'r', 'ieee-be');
366    
367     iz = 1;
368     offset = (iz - 1)*nps*4;
369     fseek(upbpid, offset, 'bof');
370     fseek(vpbpid, offset, 'bof');
371     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
372     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
373     llvpbp = upbp .* lluy + vpbp .* llvy;
374     llvpbpm1 = llvpbp;
375     ne = 510; tx = 85; ty = 85; nt = 216;
376    
377     iz = 2;
378     iz = 15;
379     iz = 16;
380     for iz = 2:50,
381     izm1 = iz - 1;
382     disp(sprintf('iz = %d',iz));
383    
384     offset = (iz - 1)*nps*4;
385     fseek(upbpid, offset, 'bof');
386     fseek(vpbpid, offset, 'bof');
387     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
388     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
389    
390     % v'B' on ll coords
391     llvpbp = upbp .* lluy + vpbp .* llvy;
392    
393     % B and Bm1
394     fseek(bid, offset, 'bof');
395     b = reshape(fread(bid, nps, 'real*4'),ne,ne,6);
396     offm1 = (iz - 2)*nps*4;
397     fseek(bm1id, offm1, 'bof');
398     bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);
399    
400     % ( v'B' )/( dB/dz )
401     dbdz = (b - bm1)/(R(iz) - R(izm1));
402     ind0 = find(dbdz==0.0);
403     dbdz(ind0) = 1.0;
404     rdbdz = 1.0 / (dbdz);
405     rdbdz(ind0) = 0.0;
406     vpbpdbdz = 0.5*(llvpbp + llvpbpm1) .* rdbdz;
407    
408     % Write the results
409     mid = fopen('out/vpbpdbdz_1994-2003', 'a', 'ieee-be');
410     fwrite(mid, vpbpdbdz, 'real*4');
411     fclose(mid);
412    
413     % Plot results
414     clo = [ -0.0015 ]; % -20;
415     chi = [ 0.0015 ]; % 20;
416     ll = zeros(6*510,510);
417     for i = 1:6
418     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
419     yb = 1; ye = 510;
420     ll(xb:xe,yb:ye) = vpbpdbdz(:,:,i);
421     end
422     shift=-1;
423     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
424     title([ '(v''B'')/(dB/dz) at ' ...
425     sprintf('%g',Rmid(iz)) ...
426     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
427     % print('-painters', '-dpng', '-r650', ...
428     % ['vpTpdTdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
429     print('-painters', '-dpng', '-r150', ...
430     ['vpBpdBdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
431    
432    
433     % Stress: calc, write, and plot
434     stress = 1000 * (2*pi/(24*3600)*sin(pi*yc/180));
435     stress = stress .* vpbpdbdz;
436     sid = fopen('out/stress_1994-2003', 'a', 'ieee-be');
437     fwrite(sid, stress, 'real*4');
438     fclose(sid);
439     clim = 0.15;
440     clo = [ -clim ]; % [ -1 ];
441     chi = [ clim ]; % [ 1 ];
442     ll = zeros(6*510,510);
443     for i = 1:6
444     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
445     yb = 1; ye = 510;
446     ll(xb:xe,yb:ye) = stress(:,:,i);
447     end
448     shift=-1;
449     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
450     title([ 'Stress at ' ...
451     sprintf('%g',Rmid(iz)) ...
452     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
453     % print('-painters', '-dpng', '-r650', ...
454     % ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
455     print('-painters', '-dpng', '-r150', ...
456     ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
457    
458    
459     % Next level
460     llvpbpm1 = llvpbp;
461    
462     end
463    
464    
465     %=======================================================
466     % Zonally average vpbpdbdz and stress
467    
468     clear all ; close all
469     tx = 85;
470     ty = 85;
471     nt = 216;
472     cx = 510;
473     cy = 510;
474     nz = 1;
475     ne = 510;
476     nps = ne * ne * 6;
477    
478     % XCS YCS XGS YGS
479     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
480     for in = 1:4
481     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
482     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
483     tx*nt,ty), ne, tx );
484     fclose(uid);
485     eval([lower(fnam(in,:)) ' = phi;']);
486     a = zeros(6*510,510);
487     for i = 1:6
488     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
489     yb = 1; ye = 510;
490     a(xb:xe,yb:ye) = phi(:,:,i);
491     end
492     eval([lower(fnam(in,:)) 's = a;'])
493     end
494     xcs = xcs + -360.0*(xcs > 180.0);
495     xgs = xgs + -360.0*(xgs > 180.0);
496     clear phi a
497    
498     nz = 1;
499     nr = 50;
500     nrm1 = nr - 1;
501     nlat = 181; nlatm1 = nlat - 1;
502    
503     hvals = linspace(-90,90,nlat);
504     % save indicies for zonal averages
505     i = 2;
506     for i = 2:nlat
507     inds = find(hvals(i-1)<yg & yg<hvals(i));
508     comm = sprintf('inds%04d = uint32(inds);',i-1);
509     eval(comm);
510     end
511    
512     % Zonally average Lat-Lon [uv]p[tsb]p
513     for f = [ 't' 's' 'b' ]
514     for uv = [ 'u' 'v' ]
515     acc = zeros(nlatm1, nrm1);
516     num = zeros(size(acc));
517     ne = 510;
518     nps = ne * ne * 6;
519     fname = [ 'out/ll' uv 'p' f 'p_1994-2003' ];
520     disp(['fname = "' fname '"']);
521     zid = fopen(fname, 'r', 'ieee-be');
522     iz = 1;
523     for iz = 1:50,
524     disp(sprintf(' iz = %d',iz));
525     offset = (iz - 1)*nps*4;
526     fseek(zid, offset, 'bof');
527     var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
528    
529     for jj = 1:nlatm1
530     eval( sprintf('clear inds; inds = inds%04d;',jj) );
531     tmp = var(inds);
532     nzinds = find(tmp ~= 0.0);
533     num(jj,iz) = length(nzinds);
534     acc(jj,iz) = sum(tmp(nzinds));
535     end
536     end
537     fclose(zid);
538     eval(['z' uv 'p' f 'p = acc ./ num;']);
539     eval(['save z' uv 'p' f 'p z' uv 'p' f 'p']);
540     end
541     end
542    
543     % zonally average T,S,B
544     for f = [ 't' 's' 'b' 'u' 'v' ]
545     acc = zeros(nlatm1, 50);
546     num = zeros(size(acc));
547     ne = 510;
548     nps = ne * ne * 6;
549     fname = [upper(f) '__tave_1994-2003'];
550     disp(fname);
551     zid = fopen(fname, 'r', 'ieee-be' );
552     iz = 1;
553     for iz = 1:50,
554     disp(sprintf(' iz = %d',iz));
555     offset = (iz - 1)*nps*4;
556     fseek(zid, offset, 'bof');
557     var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
558     jj = 22;
559     for jj = 1:nlatm1
560     eval( sprintf('clear inds; inds = inds%04d;',jj) );
561     tmp = var(inds);
562     nzinds = find(tmp ~= 0.0);
563     num(jj,iz) = length(nzinds);
564     acc(jj,iz) = sum(tmp(nzinds));
565     end
566     end
567     fclose(zid);
568     eval(['zonal' upper(f) ' = acc ./ num;']);
569     eval(['save zonal' upper(f) ' zonal' upper(f)]);
570     end
571    
572     % zonally average vpbpdbdz
573     acc = zeros(nlatm1, nrm1);
574     num = zeros(size(acc));
575     ne = 510;
576     nps = ne * ne * 6;
577     zid = fopen('out/vpbpdbdz_1994-2003', 'r', 'ieee-be');
578     iz = 1;
579     for iz = 1:49,
580     disp(sprintf('iz = %d',iz));
581     offset = (iz - 1)*nps*4;
582     fseek(zid, offset, 'bof');
583     vpbpdbdz = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
584    
585     for jj = 1:nlatm1
586     eval( sprintf('clear inds; inds = inds%04d;',jj) );
587     tmp = vpbpdbdz(inds);
588     nzinds = find(tmp ~= 0.0);
589     num(jj,iz) = length(nzinds);
590     acc(jj,iz) = sum(tmp(nzinds));
591     end
592     end
593     fclose(zid);
594     zvpbpdbdz = acc ./ num;
595     save zvpbpdbdz zvpbpdbdz
596    
597     % zonally average stress
598     acc = zeros(nlatm1, nrm1);
599     num = zeros(size(acc));
600     ne = 510;
601     nps = ne * ne * 6;
602     zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
603     iz = 1;
604     for iz = 1:49,
605     disp(sprintf('iz = %d',iz));
606     offset = (iz - 1)*nps*4;
607     fseek(zid, offset, 'bof');
608     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
609     jj = 22;
610     for jj = 1:nlatm1
611     eval( sprintf('clear inds; inds = inds%04d;',jj) );
612     tmp = stress(inds);
613     nzinds = find(tmp ~= 0.0);
614     num(jj,iz) = length(nzinds);
615     acc(jj,iz) = sum(tmp(nzinds));
616     end
617     end
618     fclose(zid);
619     stress = acc ./ num;
620     save stress stress
621    
622     %------- UNUSED -----------------------
623     % zonally average u
624     acc = zeros(nlatm1, 50);
625     num = zeros(size(acc));
626     ne = 510;
627     nps = ne * ne * 6;
628     zid = fopen('Ull_ave_1994--2003', 'r', 'ieee-be');
629     iz = 1;
630     for iz = 1:50,
631     disp(sprintf('iz = %d',iz));
632     offset = (iz - 1)*nps*4;
633     fseek(zid, offset, 'bof');
634     U = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
635     jj = 22;
636     for jj = 1:nlatm1
637     eval( sprintf('clear inds; inds = inds%04d;',jj) );
638     tmp = U(inds);
639     nzinds = find(tmp ~= 0.0);
640     num(jj,iz) = length(nzinds);
641     acc(jj,iz) = sum(tmp(nzinds));
642     end
643     end
644     fclose(zid);
645     zonalu = acc ./ num;
646     % save zonalu zonalu
647    
648     % Vertical gradient of zonally averaged B (dB/dz)
649     load zonalT
650     load allR
651     nz = size(zonalT,1);
652     dzaTdz = zeros(nz,size(zonalT,2)-1);
653     for iz = 2:50
654     dzaTdz(:,iz-1) = (zonalT(:,iz) - zonalT(:,iz-1)) ./ (R(iz) - R(iz-1));
655     end
656     % save dzaTdz dzaTdz
657     %------- UNUSED -----------------------
658    
659    
660     delR = [
661     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
662     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
663     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
664     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
665     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
666     341.50,364.50,387.50,410.50,433.50,456.50 ];
667     R = cumsum(delR) - 0.5*delR;
668     Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
669     hmid = [ -89.5:1:89.5 ];
670     save allR delR R Rmid hmid
671    
672     %------- UNUSED -----------------------
673     surf(hmid,-Rmid,num'), view(2), shading interp, colorbar
674     caxis([ -1 1 ])
675     caxis('manual')
676     surf(hmid,-Rmid,mu'), view(2), shading interp, colorbar
677     title('Zonally Averaged Stress for 1994--2003 [cube5]')
678     xlabel('Latitude [deg]')
679     zlabel('Depth [m]')
680     axis([ -90 90 -5500 200 ])
681    
682     % print -dps -painters t_001.ps
683     print('-painters', '-dpng', '-r150', 'za_str_94--03_r150.png')
684     print('-painters', '-dpng', '-r650', 'za_str_94--03_r650.png')
685    
686     % plot and print bar(U)
687     load zonalu
688     surf(hmid,-R,zonalu'), view(2), shading interp, colorbar
689     contour(hmid,-R,mu'), colorbar
690     ac = 25;
691     ac = [-0.3:.005:.3];
692     [c,h] = contourf(hmid,-R,mu',ac); colorbar
693     hold on, contour(hmid,-R,mu',ac), hold off
694     title('Zonally Averaged U for 1994--2003 [cube5]')
695     xlabel('Latitude [deg]')
696     zlabel('Depth [m]')
697     axis([ -90 90 -6000 200 ])
698     print('-painters', '-dpng', '-r150', 'Ull_94--03_r150.png')
699     print('-painters', '-dpng', '-r650', 'Ull_94--03_r650.png')
700     %------- UNUSED -----------------------
701    
702    
703     %=======================================================
704     % Put it all in NetCDF
705    
706     clear all ; close all
707    
708     load allR
709    
710     load zonalU ; i1 = 'U'; u1 = 'm/s';
711     load zonalV ; i2 = 'V'; u2 = 'm/s';
712     load zonalT ; i3 = 'T'; u3 = 'deg C';
713     load zonalS ; i4 = 'S'; u4 = 'PSU';
714     load zonalB ; i5 = 'B'; u5 = 'bouyancy';
715     load zuptp ; i10 = 'uT'; u10 = 'm*(deg C)/s';
716     load zvptp ; i11 = 'vT'; u11 = 'm*(deg C)/s';
717     load zupsp ; i12 = 'uS'; u12 = 'm*(PSU)/s';
718     load zvpsp ; i13 = 'vS'; u13 = 'm*(PSU)/s';
719     load zupbp ; i14 = 'uB'; u14 = 'm/s';
720     load zvpbp ; i15 = 'vB'; u15 = 'm/s';
721     load zvpbpdbdz ; i20 = 'zvpbpdbdz'; u20 = '1/s';
722     load stress ; i21 = 'stress'; u21 = 'N/(m^2)';
723    
724     !rm -f Zave_94-03.nc
725     nc = netcdf('Zave_94-03.nc', 'clobber');
726     nc.reference = 'Zonal averages from Dimitris "cube 5" integrations';
727     nc.author = 'Ed Hill <eh3@mit.edu>';
728     nc.date = 'Aug 19, 2004';
729     nc('X') = length(hmid);
730     nc('Y') = length(R);
731     nc('Ym1') = length(R) - 1;
732     nc{'X'} = 'X';
733     nc{'Y'} = 'Y';
734     nc{'Ym1'} = 'Ym1';
735     for ii = [ 1:5 10:15 ]
736     eval(sprintf('nc{ i%d } = { ''Y'', ''X'' };',ii));
737     end
738     nc{ i20 } = { 'Ym1', 'X' };
739     nc{ i21 } = { 'Ym1', 'X' };
740     nc{'X'}.uniquename = 'X';
741     nc{'X'}.long_name = 'latitude';
742     nc{'X'}.gridtype = ncint(0);
743     nc{'X'}.units = 'degree_north';
744     nc{'Y'}.uniquename = 'Y';
745     nc{'Y'}.long_name = 'depth';
746     nc{'Y'}.gridtype = ncint(1);
747     nc{'Y'}.units = 'm';
748     nc{'Ym1'}.uniquename = 'Ym1';
749     nc{'Ym1'}.long_name = 'depth';
750     nc{'Ym1'}.gridtype = ncint(1);
751     nc{'Ym1'}.units = 'm';
752     for ii = [ 1:5 10:15 20:21 ]
753     eval(sprintf('nc{ i%d }.units = u%d;',ii,ii));
754     eval(sprintf('nc{ i%d }.long_name = i%d;',ii,ii));
755     eval(sprintf('nc{ i%d }.missing_value = ncdouble(NaN);',ii));
756     eval(sprintf('nc{ i%d }.FillValue_ = ncdouble(0.);',ii));
757     end
758     nc{'X'}(:) = hmid;
759     nc{'Y'}(:) = -R;
760     nc{'Ym1'}(:) = -Rmid;
761     nc{ i1 }(:) = permute(zonalU,[ 2 1 ]);
762     nc{ i2 }(:) = permute(zonalV,[ 2 1 ]);
763     nc{ i3 }(:) = permute(zonalT,[ 2 1 ]);
764     nc{ i4 }(:) = permute(zonalS,[ 2 1 ]);
765     nc{ i5 }(:) = permute(zonalB,[ 2 1 ]);
766     nc{ i10 }(:) = permute(zuptp,[ 2 1 ]);
767     nc{ i11 }(:) = permute(zvptp,[ 2 1 ]);
768     nc{ i12 }(:) = permute(zupsp,[ 2 1 ]);
769     nc{ i13 }(:) = permute(zvpsp,[ 2 1 ]);
770     nc{ i14 }(:) = permute(zupbp,[ 2 1 ]);
771     nc{ i15 }(:) = permute(zvpbp,[ 2 1 ]);
772     nc{ i20 }(:) = permute(zvpbpdbdz,[ 2 1 ]);
773     nc{ i21 }(:) = permute(stress,[ 2 1 ]);
774     nc = close(nc);
775    
776     % !scp Zave_94-03.nc ingrid.mit.edu:INGRID_PEOPLE/EH3/eddy_flux/cube_5/

  ViewVC Help
Powered by ViewVC 1.1.22