/[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.3 - (hide annotations) (download)
Fri Aug 20 20:36:49 2004 UTC (20 years, 11 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +301 -38 lines
 o more fixes from David and John
   - also added a SST-binned averaging scheme for U, upbp, stress

1 edhill 1.1 %=======================================================
2     %
3 edhill 1.3 % $Id: plot_all_10.m,v 1.2 2004/08/20 03:06:58 edhill Exp $
4 edhill 1.1 %
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 edhill 1.2 % Save and print the Lat-Lon U,V values
359     % !rm -f out/ll[uv]_1994-2003
360     iz = 1;
361     for iz = [ 1:50 ]
362    
363     disp(sprintf('iz = %d',iz));
364     offset = (iz - 1)*nps*4;
365     % Load the U,V values
366     for uv = [ 'u' 'v' ]
367     fid = fopen([upper(uv) '__tave_1994-2003'],'r','ieee-be');
368     fseek(fid, offset, 'bof');
369     a = reshape(fread(fid, nps, 'real*4'),ne,ne,6);
370     fclose(fid);
371     eval([uv ' = a;']);
372     end
373    
374     % Get LatLon-oriented vectors at the CubeSphere points
375     llu = u .* llux + v .* llvx;
376     llv = u .* lluy + v .* llvy;
377    
378     % save the LatLon-rotated U,V
379     for uv = [ 'u' 'v' ]
380     % disp([' writing: ' uv f]);
381     fid = fopen(['out/ll' uv '_1994-2003'],'a','ieee-be');
382     eval(['fwrite(fid,ll' uv ',''real*4'');']);
383     fclose(fid);
384     end
385    
386     clim = 1.0;
387     pllu = zeros(6*510,510);
388     pllv = zeros(6*510,510);
389     for i = 1:6
390     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
391     yb = 1; ye = 510;
392     pllu(xb:xe,yb:ye) = llu(:,:,i);
393     pllv(xb:xe,yb:ye) = llv(:,:,i);
394     end
395     pllu(find(pllu == 0.0)) = NaN;
396     pllv(find(pllv == 0.0)) = NaN;
397     shift=-1;
398     for uv = [ 'u' 'v' ]
399     eval(['a = pll' uv ';']);
400     clo = -clim;
401     chi = clim;
402     grph_CS(a,xcs,ycs,xgs,ygs,[clo],[chi],shift);
403     title([ uv ' at ' ...
404     sprintf('%g',R(iz)) ...
405     'm depth on the 510x510x6 cubesphere ' ...
406     'for 1994--2003 ["cube5"]']);
407     print('-painters', '-dpng', '-r150', ...
408     [ uv '_' sprintf('%02d_%07.1f',iz,R(iz)) 'm_150.png'])
409     end
410    
411     end
412    
413 edhill 1.1 %=======================================================
414     % Calculate : (v'B')/(dB/dz)
415 edhill 1.3
416 edhill 1.1 bid = fopen('B__tave_1994-2003', 'r', 'ieee-be');
417     bm1id = fopen('B1_tave_1994-2003', 'r', 'ieee-be');
418     upbpid = fopen('out/upbp_1994-2003', 'r', 'ieee-be');
419     vpbpid = fopen('out/vpbp_1994-2003', 'r', 'ieee-be');
420    
421     iz = 1;
422     offset = (iz - 1)*nps*4;
423     fseek(upbpid, offset, 'bof');
424     fseek(vpbpid, offset, 'bof');
425     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
426     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
427     llvpbp = upbp .* lluy + vpbp .* llvy;
428     llvpbpm1 = llvpbp;
429     ne = 510; tx = 85; ty = 85; nt = 216;
430    
431     iz = 2;
432 edhill 1.3 iz = 9;
433     iz = 10;
434 edhill 1.2 % !rm out/stress_1994-2003 out/vpbpdbdz_1994-2003
435 edhill 1.1 for iz = 2:50,
436     izm1 = iz - 1;
437     disp(sprintf('iz = %d',iz));
438    
439     offset = (iz - 1)*nps*4;
440     fseek(upbpid, offset, 'bof');
441     fseek(vpbpid, offset, 'bof');
442     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
443     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
444    
445     % v'B' on ll coords
446     llvpbp = upbp .* lluy + vpbp .* llvy;
447    
448     % B and Bm1
449     fseek(bid, offset, 'bof');
450     b = reshape(fread(bid, nps, 'real*4'),ne,ne,6);
451 edhill 1.2 fseek(bm1id, offset, 'bof');
452 edhill 1.1 bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);
453    
454     % ( v'B' )/( dB/dz )
455 edhill 1.2 dbdz = (bm1 - b)/(R(iz) - R(izm1));
456 edhill 1.1 ind0 = find(dbdz==0.0);
457     dbdz(ind0) = 1.0;
458     rdbdz = 1.0 / (dbdz);
459     rdbdz(ind0) = 0.0;
460     vpbpdbdz = 0.5*(llvpbp + llvpbpm1) .* rdbdz;
461 edhill 1.3
462 edhill 1.1 % Write the results
463     mid = fopen('out/vpbpdbdz_1994-2003', 'a', 'ieee-be');
464     fwrite(mid, vpbpdbdz, 'real*4');
465     fclose(mid);
466    
467     % Plot results
468 edhill 1.2 clo = [ -10.0 ]; % -20;
469     chi = [ 10.0 ]; % 20;
470 edhill 1.1 ll = zeros(6*510,510);
471     for i = 1:6
472     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
473     yb = 1; ye = 510;
474     ll(xb:xe,yb:ye) = vpbpdbdz(:,:,i);
475     end
476     shift=-1;
477     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
478     title([ '(v''B'')/(dB/dz) at ' ...
479     sprintf('%g',Rmid(iz)) ...
480     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
481     % print('-painters', '-dpng', '-r650', ...
482     % ['vpTpdTdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
483     print('-painters', '-dpng', '-r150', ...
484     ['vpBpdBdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
485    
486     % Stress: calc, write, and plot
487 edhill 1.2 stress = 1000 * (2*2*pi/(24*3600)*sin(pi*yc/180));
488 edhill 1.1 stress = stress .* vpbpdbdz;
489     sid = fopen('out/stress_1994-2003', 'a', 'ieee-be');
490     fwrite(sid, stress, 'real*4');
491     fclose(sid);
492 edhill 1.3 clim = 1.0; % clim = 0.2
493 edhill 1.1 clo = [ -clim ]; % [ -1 ];
494     chi = [ clim ]; % [ 1 ];
495     ll = zeros(6*510,510);
496     for i = 1:6
497     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
498     yb = 1; ye = 510;
499     ll(xb:xe,yb:ye) = stress(:,:,i);
500     end
501     shift=-1;
502     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clo],[chi],shift);
503     title([ 'Stress at ' ...
504     sprintf('%g',Rmid(iz)) ...
505     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
506     % print('-painters', '-dpng', '-r650', ...
507     % ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
508     print('-painters', '-dpng', '-r150', ...
509 edhill 1.3 ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150_clim_0.2.png'])
510 edhill 1.1
511    
512     % Next level
513     llvpbpm1 = llvpbp;
514    
515     end
516    
517    
518     %=======================================================
519     % Zonally average vpbpdbdz and stress
520    
521     clear all ; close all
522     tx = 85;
523     ty = 85;
524     nt = 216;
525     cx = 510;
526     cy = 510;
527     nz = 1;
528     ne = 510;
529     nps = ne * ne * 6;
530    
531     % XCS YCS XGS YGS
532     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
533     for in = 1:4
534     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
535     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
536     tx*nt,ty), ne, tx );
537     fclose(uid);
538     eval([lower(fnam(in,:)) ' = phi;']);
539     a = zeros(6*510,510);
540     for i = 1:6
541     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
542     yb = 1; ye = 510;
543     a(xb:xe,yb:ye) = phi(:,:,i);
544     end
545     eval([lower(fnam(in,:)) 's = a;'])
546     end
547     xcs = xcs + -360.0*(xcs > 180.0);
548     xgs = xgs + -360.0*(xgs > 180.0);
549     clear phi a
550    
551     nz = 1;
552     nr = 50;
553     nrm1 = nr - 1;
554     nlat = 181; nlatm1 = nlat - 1;
555    
556     hvals = linspace(-90,90,nlat);
557     % save indicies for zonal averages
558     i = 2;
559     for i = 2:nlat
560     inds = find(hvals(i-1)<yg & yg<hvals(i));
561     comm = sprintf('inds%04d = uint32(inds);',i-1);
562     eval(comm);
563     end
564    
565     % Zonally average Lat-Lon [uv]p[tsb]p
566     for f = [ 't' 's' 'b' ]
567     for uv = [ 'u' 'v' ]
568     acc = zeros(nlatm1, nrm1);
569     num = zeros(size(acc));
570     ne = 510;
571     nps = ne * ne * 6;
572     fname = [ 'out/ll' uv 'p' f 'p_1994-2003' ];
573     disp(['fname = "' fname '"']);
574     zid = fopen(fname, 'r', 'ieee-be');
575     iz = 1;
576     for iz = 1:50,
577     disp(sprintf(' iz = %d',iz));
578     offset = (iz - 1)*nps*4;
579     fseek(zid, offset, 'bof');
580     var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
581    
582     for jj = 1:nlatm1
583     eval( sprintf('clear inds; inds = inds%04d;',jj) );
584     tmp = var(inds);
585     nzinds = find(tmp ~= 0.0);
586     num(jj,iz) = length(nzinds);
587     acc(jj,iz) = sum(tmp(nzinds));
588     end
589     end
590     fclose(zid);
591     eval(['z' uv 'p' f 'p = acc ./ num;']);
592     eval(['save z' uv 'p' f 'p z' uv 'p' f 'p']);
593     end
594     end
595    
596     % zonally average T,S,B
597 edhill 1.3 for f = [ 't' 's' 'b' ]
598 edhill 1.1 acc = zeros(nlatm1, 50);
599     num = zeros(size(acc));
600     ne = 510;
601     nps = ne * ne * 6;
602     fname = [upper(f) '__tave_1994-2003'];
603     disp(fname);
604     zid = fopen(fname, 'r', 'ieee-be' );
605     iz = 1;
606     for iz = 1:50,
607     disp(sprintf(' iz = %d',iz));
608     offset = (iz - 1)*nps*4;
609     fseek(zid, offset, 'bof');
610     var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
611     jj = 22;
612     for jj = 1:nlatm1
613     eval( sprintf('clear inds; inds = inds%04d;',jj) );
614     tmp = var(inds);
615     nzinds = find(tmp ~= 0.0);
616     num(jj,iz) = length(nzinds);
617     acc(jj,iz) = sum(tmp(nzinds));
618     end
619     end
620     fclose(zid);
621     eval(['zonal' upper(f) ' = acc ./ num;']);
622     eval(['save zonal' upper(f) ' zonal' upper(f)]);
623     end
624    
625 edhill 1.3 % zonally average Lat-Lon u,v
626     for f = [ 'u' 'v' ]
627     acc = zeros(nlatm1, 50);
628     num = zeros(size(acc));
629     ne = 510;
630     nps = ne * ne * 6;
631     fname = ['out/ll' f '_1994-2003'];
632     disp(fname);
633     zid = fopen(fname, 'r', 'ieee-be' );
634     iz = 1;
635     for iz = 1:50,
636     disp(sprintf(' iz = %d',iz));
637     offset = (iz - 1)*nps*4;
638     fseek(zid, offset, 'bof');
639     var = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
640     jj = 22;
641     for jj = 1:nlatm1
642     eval( sprintf('clear inds; inds = inds%04d;',jj) );
643     tmp = var(inds);
644     nzinds = find(tmp ~= 0.0);
645     num(jj,iz) = length(nzinds);
646     acc(jj,iz) = sum(tmp(nzinds));
647     end
648     end
649     fclose(zid);
650     eval(['zonal' upper(f) ' = acc ./ num;']);
651     eval(['save zonal' upper(f) ' zonal' upper(f)]);
652     end
653    
654 edhill 1.1 % zonally average vpbpdbdz
655     acc = zeros(nlatm1, nrm1);
656     num = zeros(size(acc));
657     ne = 510;
658     nps = ne * ne * 6;
659     zid = fopen('out/vpbpdbdz_1994-2003', 'r', 'ieee-be');
660     iz = 1;
661     for iz = 1:49,
662     disp(sprintf('iz = %d',iz));
663     offset = (iz - 1)*nps*4;
664     fseek(zid, offset, 'bof');
665     vpbpdbdz = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
666    
667     for jj = 1:nlatm1
668     eval( sprintf('clear inds; inds = inds%04d;',jj) );
669     tmp = vpbpdbdz(inds);
670     nzinds = find(tmp ~= 0.0);
671     num(jj,iz) = length(nzinds);
672     acc(jj,iz) = sum(tmp(nzinds));
673     end
674     end
675     fclose(zid);
676     zvpbpdbdz = acc ./ num;
677     save zvpbpdbdz zvpbpdbdz
678    
679     % zonally average stress
680     acc = zeros(nlatm1, nrm1);
681     num = zeros(size(acc));
682     ne = 510;
683     nps = ne * ne * 6;
684     zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
685     iz = 1;
686     for iz = 1:49,
687     disp(sprintf('iz = %d',iz));
688     offset = (iz - 1)*nps*4;
689     fseek(zid, offset, 'bof');
690     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
691     jj = 22;
692     for jj = 1:nlatm1
693     eval( sprintf('clear inds; inds = inds%04d;',jj) );
694     tmp = stress(inds);
695     nzinds = find(tmp ~= 0.0);
696     num(jj,iz) = length(nzinds);
697     acc(jj,iz) = sum(tmp(nzinds));
698     end
699     end
700     fclose(zid);
701     stress = acc ./ num;
702     save stress stress
703    
704 edhill 1.3
705     % zonally average stress IN A SPECIFIC LONGITUDE RANGE
706     hvals = linspace(-90,90,nlat);
707     % save indicies for zonal averages
708     i = 2;
709     for i = 2:nlat
710     inds = find(hvals(i-1)<yg & yg<hvals(i) & 210.0<xg & xg<230.0);
711     comm = sprintf('inds%04d = uint32(inds);',i-1);
712     eval(comm);
713     end
714     acc = zeros(nlatm1, nrm1);
715 edhill 1.1 num = zeros(size(acc));
716     ne = 510;
717     nps = ne * ne * 6;
718 edhill 1.3 zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
719 edhill 1.1 iz = 1;
720 edhill 1.3 for iz = 1:49,
721 edhill 1.1 disp(sprintf('iz = %d',iz));
722     offset = (iz - 1)*nps*4;
723     fseek(zid, offset, 'bof');
724 edhill 1.3 stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
725 edhill 1.1 jj = 22;
726     for jj = 1:nlatm1
727     eval( sprintf('clear inds; inds = inds%04d;',jj) );
728 edhill 1.3 tmp = stress(inds);
729 edhill 1.1 nzinds = find(tmp ~= 0.0);
730     num(jj,iz) = length(nzinds);
731     acc(jj,iz) = sum(tmp(nzinds));
732     end
733     end
734     fclose(zid);
735 edhill 1.3 s_210_230 = acc ./ num;
736     save s_210_230 s_210_230
737    
738     %------- UNUSED -----------------------
739 edhill 1.1 % Vertical gradient of zonally averaged B (dB/dz)
740     load zonalT
741     load allR
742     nz = size(zonalT,1);
743     dzaTdz = zeros(nz,size(zonalT,2)-1);
744     for iz = 2:50
745     dzaTdz(:,iz-1) = (zonalT(:,iz) - zonalT(:,iz-1)) ./ (R(iz) - R(iz-1));
746     end
747     % save dzaTdz dzaTdz
748     %------- UNUSED -----------------------
749    
750    
751     delR = [
752     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
753     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
754     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
755     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
756     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
757     341.50,364.50,387.50,410.50,433.50,456.50 ];
758     R = cumsum(delR) - 0.5*delR;
759     Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
760     hmid = [ -89.5:1:89.5 ];
761     save allR delR R Rmid hmid
762    
763     %------- UNUSED -----------------------
764     surf(hmid,-Rmid,num'), view(2), shading interp, colorbar
765     caxis([ -1 1 ])
766     caxis('manual')
767     surf(hmid,-Rmid,mu'), view(2), shading interp, colorbar
768     title('Zonally Averaged Stress for 1994--2003 [cube5]')
769     xlabel('Latitude [deg]')
770     zlabel('Depth [m]')
771     axis([ -90 90 -5500 200 ])
772    
773     % print -dps -painters t_001.ps
774     print('-painters', '-dpng', '-r150', 'za_str_94--03_r150.png')
775     print('-painters', '-dpng', '-r650', 'za_str_94--03_r650.png')
776     %------- UNUSED -----------------------
777    
778    
779     %=======================================================
780     % Put it all in NetCDF
781    
782     clear all ; close all
783    
784     load allR
785    
786     load zonalU ; i1 = 'U'; u1 = 'm/s';
787     load zonalV ; i2 = 'V'; u2 = 'm/s';
788     load zonalT ; i3 = 'T'; u3 = 'deg C';
789     load zonalS ; i4 = 'S'; u4 = 'PSU';
790 edhill 1.3 load zonalB ; i5 = 'B'; u5 = 'm/(s^2)';
791 edhill 1.1 load zuptp ; i10 = 'uT'; u10 = 'm*(deg C)/s';
792     load zvptp ; i11 = 'vT'; u11 = 'm*(deg C)/s';
793     load zupsp ; i12 = 'uS'; u12 = 'm*(PSU)/s';
794     load zvpsp ; i13 = 'vS'; u13 = 'm*(PSU)/s';
795 edhill 1.3 load zupbp ; i14 = 'uB'; u14 = 'm^2/(s^3)';
796     load zvpbp ; i15 = 'vB'; u15 = 'm^2/(s^3)';
797     load zvpbpdbdz ; i20 = 'zvpbpdbdz'; u20 = 'm^2/s';
798 edhill 1.1 load stress ; i21 = 'stress'; u21 = 'N/(m^2)';
799 edhill 1.3 load s_210_230 ; i30 = 'stress_210_230'; u30 = 'N/(m^2)';
800 edhill 1.1
801     !rm -f Zave_94-03.nc
802     nc = netcdf('Zave_94-03.nc', 'clobber');
803     nc.reference = 'Zonal averages from Dimitris "cube 5" integrations';
804     nc.author = 'Ed Hill <eh3@mit.edu>';
805     nc.date = 'Aug 19, 2004';
806     nc('X') = length(hmid);
807     nc('Y') = length(R);
808     nc('Ym1') = length(R) - 1;
809     nc{'X'} = 'X';
810     nc{'Y'} = 'Y';
811     nc{'Ym1'} = 'Ym1';
812     for ii = [ 1:5 10:15 ]
813     eval(sprintf('nc{ i%d } = { ''Y'', ''X'' };',ii));
814     end
815     nc{ i20 } = { 'Ym1', 'X' };
816     nc{ i21 } = { 'Ym1', 'X' };
817 edhill 1.3 nc{ i30 } = { 'Ym1', 'X' };
818 edhill 1.1 nc{'X'}.uniquename = 'X';
819     nc{'X'}.long_name = 'latitude';
820     nc{'X'}.gridtype = ncint(0);
821     nc{'X'}.units = 'degree_north';
822     nc{'Y'}.uniquename = 'Y';
823     nc{'Y'}.long_name = 'depth';
824     nc{'Y'}.gridtype = ncint(1);
825     nc{'Y'}.units = 'm';
826     nc{'Ym1'}.uniquename = 'Ym1';
827     nc{'Ym1'}.long_name = 'depth';
828     nc{'Ym1'}.gridtype = ncint(1);
829     nc{'Ym1'}.units = 'm';
830 edhill 1.3 for ii = [ 1:5 10:15 20:21 30 ]
831 edhill 1.1 eval(sprintf('nc{ i%d }.units = u%d;',ii,ii));
832     eval(sprintf('nc{ i%d }.long_name = i%d;',ii,ii));
833     eval(sprintf('nc{ i%d }.missing_value = ncdouble(NaN);',ii));
834     eval(sprintf('nc{ i%d }.FillValue_ = ncdouble(0.);',ii));
835     end
836     nc{'X'}(:) = hmid;
837     nc{'Y'}(:) = -R;
838     nc{'Ym1'}(:) = -Rmid;
839     nc{ i1 }(:) = permute(zonalU,[ 2 1 ]);
840     nc{ i2 }(:) = permute(zonalV,[ 2 1 ]);
841     nc{ i3 }(:) = permute(zonalT,[ 2 1 ]);
842     nc{ i4 }(:) = permute(zonalS,[ 2 1 ]);
843     nc{ i5 }(:) = permute(zonalB,[ 2 1 ]);
844     nc{ i10 }(:) = permute(zuptp,[ 2 1 ]);
845     nc{ i11 }(:) = permute(zvptp,[ 2 1 ]);
846     nc{ i12 }(:) = permute(zupsp,[ 2 1 ]);
847     nc{ i13 }(:) = permute(zvpsp,[ 2 1 ]);
848     nc{ i14 }(:) = permute(zupbp,[ 2 1 ]);
849     nc{ i15 }(:) = permute(zvpbp,[ 2 1 ]);
850     nc{ i20 }(:) = permute(zvpbpdbdz,[ 2 1 ]);
851     nc{ i21 }(:) = permute(stress,[ 2 1 ]);
852 edhill 1.3 nc{ i30 }(:) = permute(s_210_230,[ 2 1 ]);
853 edhill 1.1 nc = close(nc);
854    
855     % !scp Zave_94-03.nc ingrid.mit.edu:INGRID_PEOPLE/EH3/eddy_flux/cube_5/
856 edhill 1.3
857    
858     %=======================================================
859     % Create a new NetCDF file containing the
860     % temperature-averaged fields: stress, U, v'b'
861    
862    
863     clear all ; close all
864     tx = 85;
865     ty = 85;
866     nt = 216;
867     cx = 510;
868     cy = 510;
869     nz = 1;
870     ne = 510;
871     nps = ne * ne * 6;
872    
873     % XCS YCS XGS YGS
874     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
875     for in = 1:4
876     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
877     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
878     tx*nt,ty), ne, tx );
879     fclose(uid);
880     eval([lower(fnam(in,:)) ' = phi;']);
881     a = zeros(6*510,510);
882     for i = 1:6
883     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
884     yb = 1; ye = 510;
885     a(xb:xe,yb:ye) = phi(:,:,i);
886     end
887     eval([lower(fnam(in,:)) 's = a;'])
888     end
889     xcs = xcs + -360.0*(xcs > 180.0);
890     xgs = xgs + -360.0*(xgs > 180.0);
891     clear phi a
892    
893     nr = 50;
894     nrm1 = nr - 1;
895    
896     % Get the SST
897     f = [ 't' ];
898     ne = 510;
899     nps = ne * ne * 6;
900     fname = [upper(f) '__tave_1994-2003'];
901     zid = fopen(fname, 'r', 'ieee-be' );
902     iz = 1;
903     offset = (iz - 1)*nps*4;
904     fseek(zid, offset, 'bof');
905     sst = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
906     fclose(zid);
907     sst(find(sst==0.0)) = NaN;
908     % surf(sst(:,:,1)), shading interp, view(2)
909     % [ squeeze(min(min(sst)))' ; squeeze(max(max(sst)))' ]'
910     % [ squeeze(min(min(yg)))' ; squeeze(max(max(yg)))' ]'
911    
912     % Get the SST-based indicies
913     ntemp = 51;
914     tval = linspace(-3,33,ntemp);
915     for i = 2:length(tval)
916     % Northern Hemisphere (yg > -10)
917     inds = find(yg>-10.0 & tval(i-1)<sst & sst<tval(i));
918     eval(sprintf('ni%04d = uint32(inds);',i-1));
919     % Southern Hemisphere (yg < 10)
920     inds = find(yg<10.0 & tval(i-1)<sst & sst<tval(i));
921     eval(sprintf('si%04d = uint32(inds);',i-1));
922     end
923     tmid = tval(1:(length(tval)-1)) + diff(tval);
924     save tval tval tmid
925    
926     % "temperature-average" the stress field
927     nlatm1 = length(tval) - 1;
928     accn = zeros(nlatm1, nrm1);
929     accs = zeros(nlatm1, nrm1);
930     nnum = zeros(size(nacc));
931     snum = zeros(size(nacc));
932     ne = 510;
933     nps = ne * ne * 6;
934     zid = fopen('out/stress_1994-2003', 'r', 'ieee-be');
935     iz = 1;
936     for iz = 1:48,
937     disp(sprintf('iz = %d',iz));
938     offset = (iz - 1)*nps*4;
939     fseek(zid, offset, 'bof');
940     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
941     jj = 22;
942     for jj = 1:nlatm1
943     eval( sprintf('clear ni; ni = ni%04d;',jj) );
944     eval( sprintf('clear si; si = si%04d;',jj) );
945     ntmp = stress(ni);
946     stmp = stress(si);
947     nz_ni = find(ntmp ~= 0.0);
948     nz_si = find(stmp ~= 0.0);
949     nnum(jj,iz) = length(nz_ni);
950     snum(jj,iz) = length(nz_si);
951     nacc(jj,iz) = sum(ntmp(nz_ni));
952     sacc(jj,iz) = sum(stmp(nz_si));
953     end
954     end
955     fclose(zid);
956     ntstress = nacc ./ nnum;
957     ststress = sacc ./ snum;
958     save t_stress ntstress ststress
959    
960    
961     % "temperature-average" the U field
962     nlatm1 = length(tval) - 1;
963     accn = zeros(nlatm1, nrm1);
964     accs = zeros(nlatm1, nrm1);
965     nnum = zeros(size(nacc));
966     snum = zeros(size(nacc));
967     ne = 510;
968     nps = ne * ne * 6;
969     zid = fopen('out/llu_1994-2003', 'r', 'ieee-be');
970     iz = 1;
971     for iz = 1:49,
972     disp(sprintf('iz = %d',iz));
973     offset = (iz - 1)*nps*4;
974     fseek(zid, offset, 'bof');
975     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
976     jj = 22;
977     for jj = 1:nlatm1
978     eval( sprintf('clear ni; ni = ni%04d;',jj) );
979     eval( sprintf('clear si; si = si%04d;',jj) );
980     ntmp = stress(ni);
981     stmp = stress(si);
982     nz_ni = find(ntmp ~= 0.0);
983     nz_si = find(stmp ~= 0.0);
984     nnum(jj,iz) = length(nz_ni);
985     snum(jj,iz) = length(nz_si);
986     nacc(jj,iz) = sum(ntmp(nz_ni));
987     sacc(jj,iz) = sum(stmp(nz_si));
988     end
989     end
990     fclose(zid);
991     ntu = nacc ./ nnum;
992     stu = sacc ./ snum;
993     save t_u ntu stu
994    
995    
996     % "temperature-average" the llvpbp field
997     nlatm1 = length(tval) - 1;
998     accn = zeros(nlatm1, nrm1);
999     accs = zeros(nlatm1, nrm1);
1000     nnum = zeros(size(nacc));
1001     snum = zeros(size(nacc));
1002     ne = 510;
1003     nps = ne * ne * 6;
1004     zid = fopen('out/llvpbp_1994-2003', 'r', 'ieee-be');
1005     iz = 1;
1006     for iz = 1:49,
1007     disp(sprintf('iz = %d',iz));
1008     offset = (iz - 1)*nps*4;
1009     fseek(zid, offset, 'bof');
1010     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
1011     jj = 22;
1012     for jj = 1:nlatm1
1013     eval( sprintf('clear ni; ni = ni%04d;',jj) );
1014     eval( sprintf('clear si; si = si%04d;',jj) );
1015     ntmp = stress(ni);
1016     stmp = stress(si);
1017     nz_ni = find(ntmp ~= 0.0);
1018     nz_si = find(stmp ~= 0.0);
1019     nnum(jj,iz) = length(nz_ni);
1020     snum(jj,iz) = length(nz_si);
1021     nacc(jj,iz) = sum(ntmp(nz_ni));
1022     sacc(jj,iz) = sum(stmp(nz_si));
1023     end
1024     end
1025     fclose(zid);
1026     ntllvpbp = nacc ./ nnum;
1027     stllvpbp = sacc ./ snum;
1028     save t_llvpbp ntllvpbp stllvpbp
1029    
1030    
1031     clear all ; close all
1032    
1033     load allR
1034     load tval
1035     load t_stress;
1036     i1 = 'nt_stress'; u1 = 'N/(m^2)';
1037     i2 = 'st_stress'; u2 = 'N/(m^2)';
1038     load t_llvpbp
1039     i3 = 'nt_vpbp'; u3 = 'm^2/(s^3)';
1040     i4 = 'st_vpbp'; u4 = 'm^2/(s^3)';
1041     load t_u
1042     i5 = 'nt_u'; u5 = 'deg C';
1043     i6 = 'st_u'; u6 = 'deg C';
1044    
1045     !rm -f SSTave_94-03.nc
1046     nc = netcdf('SSTave_94-03.nc', 'clobber');
1047     nc.reference = 'SST-based averages from Dimitris "cube 5" integrations';
1048     nc.author = 'Ed Hill <eh3@mit.edu>';
1049     nc.date = 'Aug 20, 2004';
1050     nc('X') = length(tmid);
1051     nc('Y') = length(Rmid);
1052     nc('Ym1') = length(Rmid) - 1;
1053     Rmidm1 = Rmid([1:(length(Rmid) - 1)]);
1054     nc{'X'} = 'X';
1055     nc{'Y'} = 'Y';
1056     nc{'Ym1'} = 'Ym1';
1057     for ii = [ 1:2 ]
1058     eval(sprintf('nc{ i%d } = { ''Ym1'', ''X'' };',ii));
1059     end
1060     for ii = [ 3:6 ]
1061     eval(sprintf('nc{ i%d } = { ''Y'', ''X'' };',ii));
1062     end
1063     nc{'X'}.uniquename = 'X';
1064     nc{'X'}.long_name = 'SST';
1065     nc{'X'}.gridtype = ncint(0);
1066     nc{'X'}.units = 'deg C';
1067     nc{'Y'}.uniquename = 'Y';
1068     nc{'Y'}.long_name = 'depth';
1069     nc{'Y'}.gridtype = ncint(1);
1070     nc{'Y'}.units = 'm';
1071     nc{'Ym1'}.uniquename = 'Ym1';
1072     nc{'Ym1'}.long_name = 'depth';
1073     nc{'Ym1'}.gridtype = ncint(1);
1074     nc{'Ym1'}.units = 'm';
1075     for ii = [ 1:6 ]
1076     eval(sprintf('nc{ i%d }.units = u%d;',ii,ii));
1077     eval(sprintf('nc{ i%d }.long_name = i%d;',ii,ii));
1078     eval(sprintf('nc{ i%d }.missing_value = ncdouble(NaN);',ii));
1079     eval(sprintf('nc{ i%d }.FillValue_ = ncdouble(0.);',ii));
1080     end
1081     nc{'X'}(:) = tmid;
1082     nc{'Y'}(:) = -Rmid;
1083     nc{'Ym1'}(:) = -Rmidm1;
1084     nc{ i1 }(:) = permute(ntstress,[ 2 1 ]);
1085     nc{ i2 }(:) = permute(ststress,[ 2 1 ]);
1086     nc{ i3 }(:) = permute(ntllvpbp,[ 2 1 ]);
1087     nc{ i4 }(:) = permute(stllvpbp,[ 2 1 ]);
1088     nc{ i5 }(:) = permute(ntu,[ 2 1 ]);
1089     nc{ i6 }(:) = permute(stu,[ 2 1 ]);
1090     nc = close(nc);
1091    
1092     % !scp SSTave_94-03.nc ingrid.mit.edu:INGRID_PEOPLE/EH3/eddy_flux/cube_5/

  ViewVC Help
Powered by ViewVC 1.1.22