/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/c20a/plot_c20a.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/eddy_flux/c20a/plot_c20a.m

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


Revision 1.7 - (hide annotations) (download)
Fri Mar 4 03:21:58 2005 UTC (20 years, 4 months ago) by edhill
Branch: MAIN
Changes since 1.6: +8 -9 lines
 o fix zonal average

1 edhill 1.1 %=======================================================
2     %
3 edhill 1.7 % $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/eddy_flux/c20a/plot_c20a.m,v 1.6 2005/03/03 21:05:24 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_20a" or "c20a" 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     % ssh eddy
24     % cd /r/r0/edhill/eddy_stats/eddy_flux_c20a
25    
26     % matlab -nojvm
27     % matlab -nojvm -nodisplay
28    
29     clear all
30     close all
31    
32    
33 edhill 1.2 %==================================================================
34     % Read the tile00?.mitgrid files
35 edhill 1.1 gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
36     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
37    
38 edhill 1.2 ne = 510;
39     nep1 = ne + 1;
40 edhill 1.1 iface = 1;
41     for iface = 1:6
42     fname = sprintf('grid/tile%03d.mitgrid', iface);
43     gid = fopen(fname, 'r', 'ieee-be');
44 edhill 1.2 tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
45 edhill 1.1 fclose(gid);
46     % surf(tmp(:,:,1)), view(2), shading interp
47     % for jj = 1:length(gvars)
48     for jj = 1:7
49     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
50     [gvars{jj}], iface, jj);
51     eval(comm);
52     end
53     end
54     % surf(XC(:,:,1)), view(2), shading interp
55     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
56 edhill 1.2 % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
57 edhill 1.1 % surf(YC(:,:,1)), view(2), shading interp
58     % surf(XG(:,:,1)), view(2), shading interp
59     % surf(YG(:,:,1)), view(2), shading interp
60 edhill 1.2 is = [1:ne];
61     vs = { 'XC','YC','DXF','DYF','RA' };
62     for i = 1:length(vs)
63     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
64 edhill 1.1 end
65    
66 edhill 1.2 delR = [ ...
67     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
68     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
69     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
70     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
71     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
72     341.50,364.50,387.50,410.50,433.50,456.50 ];
73     R = cumsum(delR) - 0.5*delR;
74    
75     %==================================================================
76     % Project fields to lower-res 1-degree Lat-Lon and write
77     % as NetCDF for viewing with Ingrid
78     %
79 edhill 1.1 % !echo `ls -1 ave__92_99 | sed -e 's|\.ave||g' \
80     % | sed -e "s|^|'|g" | sed -e "s|$|',|g"`
81 edhill 1.2 fields_3d = { ...
82     'DRHODR', 'RHOANOSQ', 'RHOAnoma', 'SALT', ...
83     'SALTSQ', ...
84     'THETA', 'THETASQ', 'URHOMASS', ...
85 edhill 1.1 'USLTMASS', 'UTHMASS', 'UVEL', 'UVELMASS', 'UVELSQ', ...
86     'UV_VEL_Z', 'VRHOMASS', 'VSLTMASS', 'VTHMASS', ...
87     'VVEL', 'VVELMASS', 'VVELSQ', 'WRHOMASS', 'WSLTMASS', ...
88     'WTHMASS', 'WU_VEL', 'WVELMASS', 'WVELSQ', 'WV_VEL' };
89 edhill 1.2 fields_2d = { ...
90     'ETANSQ', 'SFLUX', 'SRELAX', 'TAUX', 'TAUY', 'TFLUX', ...
91     'TICE', 'TRELAX' };
92     % fid = fopen('ave__92_99/DRHODR.ave','r','ieee-be');
93     % tmp = fread(gid,inf,'real*4',0,'ieee-be');
94     % fclose(fid);
95    
96     ne = 510;
97     nf = 6;
98     nz = 50;
99     nslab = ne*ne*nf;
100     adir = 'ave__92_99';
101     lat = [-90:90];
102     lon = [0:360];
103    
104     % ! rm -f cube_20a_*.nc
105    
106     ifld = 5;
107     fields = union(fields_2d, fields_3d);
108     for ifld = 1:length(fields)
109    
110 edhill 1.4 ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
111 edhill 1.2 if ismember(fields{ifld},fields_2d)
112     ir = [ 1 ];
113     end
114    
115     disp([ ' ' fields{ifld} ' :' ]);
116    
117     nc = netcdf(['cube_20a_' fields{ifld} '.nc'], 'clobber');
118     nc.reference = [ 'Results from Dimitris Menemenlis''' ...
119     '"cube 20a" integrations' ];
120     nc.author = 'Ed Hill <eh3@mit.edu>';
121     nc.date = 'Feb 24, 2005';
122     nc('X') = length(lon);
123     nc('Y') = length(lat);
124     nc('Z') = length(ir);
125     nc{'X'} = 'X';
126     nc{'Y'} = 'Y';
127     nc{'Z'} = 'Z';
128     nc{'X'}.uniquename = 'X';
129     nc{'X'}.long_name = 'longitude';
130     nc{'X'}.gridtype = ncint(1);
131     nc{'X'}.units = 'degree_east';
132     nc{'Y'}.uniquename = 'Y';
133     nc{'Y'}.long_name = 'latitude';
134     nc{'Y'}.gridtype = ncint(0);
135     nc{'Y'}.units = 'degree_north';
136     nc{'Z'}.uniquename = 'Z';
137     nc{'Z'}.long_name = 'depth';
138     nc{'Z'}.gridtype = ncint(0);
139     nc{'Z'}.units = 'm';
140     nc{'X'}(:) = lon;
141     nc{'Y'}(:) = lat;
142     nc{'Z'}(:) = R(ir);
143    
144     fname = sprintf('%s/%s.ave',adir,fields{ifld});
145     fid = fopen(fname,'r','ieee-be');
146    
147     id = fields{ifld};
148     nc{ id } = { 'Z' 'Y' 'X' };
149     nc{ id }.missing_value = ncdouble(NaN);
150     nc{ id }.FillValue_ = ncdouble(0.0);
151    
152     ii = 1;
153     for ii = 1:length(ir)
154    
155     iz = ir(ii);
156     disp(sprintf(' iz = %g',iz));
157    
158     fseek(fid,nslab*4*(iz-1),'bof');
159     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
160     tr = permute(reshape(tmp,[ 510 6 510 ]),[1 3 2]);
161     % surf(tr(:,:,1)), view(2), shading interp
162     xc360 = XC + 180;
163     trn = tr;
164     trn(find(tr == 0.0)) = NaN;
165     clear tmp tr
166     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
167     v = ll_regrid(xc360,YC,trn,lon,lat);
168     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
169    
170     nc{ id }(ii,:,:) = permute(v,[2 1]);
171    
172     end
173    
174     fclose(fid);
175     nc = close(nc);
176 edhill 1.4
177 edhill 1.2 end
178 edhill 1.1
179 edhill 1.3 % ! ncdump cube_20a.nc | more
180     % ! rm -rf output_netCDF ; mkdir output_netCDF
181     % ! scp output_netCDF/*.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_20a/
182    
183    
184    
185    
186     %=======================================================
187     % Compute [uvw]'[tsb]'
188 edhill 1.1
189 edhill 1.5 clear all
190     close all
191    
192 edhill 1.6 gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
193     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
194 edhill 1.5 ne = 510;
195     nep1 = ne + 1;
196     iface = 1;
197     for iface = 1:6
198     fname = sprintf('grid/tile%03d.mitgrid', iface);
199     gid = fopen(fname, 'r', 'ieee-be');
200     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
201     fclose(gid);
202     % surf(tmp(:,:,1)), view(2), shading interp
203     % for jj = 1:length(gvars)
204     for jj = 1:7
205     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
206     [gvars{jj}], iface, jj);
207     eval(comm);
208     end
209     end
210     % surf(XC(:,:,1)), view(2), shading interp
211     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
212     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
213     % surf(YC(:,:,1)), view(2), shading interp
214     % surf(XG(:,:,1)), view(2), shading interp
215     % surf(YG(:,:,1)), view(2), shading interp
216     is = [1:ne];
217     vs = { 'XC','YC','DXF','DYF','RA' };
218     for i = 1:length(vs)
219     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
220     end
221    
222     delR = [ ...
223     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
224     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
225     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
226     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
227     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
228     341.50,364.50,387.50,410.50,433.50,456.50 ];
229     R = cumsum(delR) - 0.5*delR;
230    
231     n1 = ne - 1;
232     dux = zeros(size(XC));
233     duy = zeros(size(XC));
234     dvx = zeros(size(XC));
235     dvy = zeros(size(XC));
236     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
237     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
238     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
239     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
240     dux = dux + 360*double(dux < 180);
241     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
242     duy = duy + 360*double(duy < 180);
243     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
244     dvx = dvx + 360*double(dvx < 180);
245     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
246     dvy = dvy + 360*double(dvy < 180);
247     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
248     llux = dux ./ sqrt(dux.^2 + duy.^2);
249     lluy = duy ./ sqrt(dux.^2 + duy.^2);
250     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
251     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
252 edhill 1.1
253 edhill 1.3 lpath = 'ave__92_04/';
254     u__id = fopen( [lpath 'UVEL.ave'], 'r', 'ieee-be'); % 1
255     v__id = fopen( [lpath 'VVEL.ave'], 'r', 'ieee-be'); % 2
256     %w__id = fopen( [lpath 'WVEL.ave'], 'r', 'ieee-be');
257     u2_id = fopen( [lpath 'UVELSQ.ave'], 'r', 'ieee-be'); % 3
258     v2_id = fopen( [lpath 'VVELSQ.ave'], 'r', 'ieee-be'); % 4
259     w2_id = fopen( [lpath 'WVELSQ.ave'], 'r', 'ieee-be'); % 5
260     um_id = fopen( [lpath 'UVELMASS.ave'], 'r', 'ieee-be'); % 6
261     vm_id = fopen( [lpath 'VVELMASS.ave'], 'r', 'ieee-be'); % 7
262     wm_id = fopen( [lpath 'WVELMASS.ave'], 'r', 'ieee-be'); % 8
263     t__id = fopen( [lpath 'THETA.ave'], 'r', 'ieee-be'); % 9
264     t2_id = fopen( [lpath 'THETASQ.ave'], 'r', 'ieee-be'); % 10
265     s__id = fopen( [lpath 'SALT.ave'], 'r', 'ieee-be'); % 11
266     s2_id = fopen( [lpath 'SALTSQ.ave'], 'r', 'ieee-be'); % 12
267     b__id = fopen( [lpath 'RHOAnoma.ave'], 'r', 'ieee-be'); % 13
268     b2_id = fopen( [lpath 'RHOANOSQ.ave'], 'r', 'ieee-be'); % 14
269     ut_id = fopen( [lpath 'UTHMASS.ave'], 'r', 'ieee-be'); % 15
270     vt_id = fopen( [lpath 'VTHMASS.ave'], 'r', 'ieee-be'); % 16
271     wt_id = fopen( [lpath 'WTHMASS.ave'], 'r', 'ieee-be'); % 17
272     us_id = fopen( [lpath 'USLTMASS.ave'], 'r', 'ieee-be'); % 18
273     vs_id = fopen( [lpath 'VSLTMASS.ave'], 'r', 'ieee-be'); % 19
274     ws_id = fopen( [lpath 'WSLTMASS.ave'], 'r', 'ieee-be'); % 20
275     ub_id = fopen( [lpath 'URHOMASS.ave'], 'r', 'ieee-be'); % 21
276     vb_id = fopen( [lpath 'VRHOMASS.ave'], 'r', 'ieee-be'); % 22
277     wb_id = fopen( [lpath 'WRHOMASS.ave'], 'r', 'ieee-be'); % 23
278     dr_id = fopen( [lpath 'DRHODR.ave'], 'r', 'ieee-be'); % 24
279     iids = [ u__id v__id u2_id v2_id w2_id um_id vm_id wm_id ...
280     t__id t2_id s__id s2_id b__id b2_id ...
281     ut_id vt_id wt_id us_id vs_id ws_id ...
282     ub_id vb_id wb_id dr_id ];
283 edhill 1.6
284     % ! rm -rf primes_92_04 ; mkdir primes_92_04
285 edhill 1.3 opath = 'primes_92_04/';
286     up2___id = fopen([opath 'up2'], 'wb', 'ieee-be');
287     vp2___id = fopen([opath 'vp2'], 'wb', 'ieee-be');
288     wp2___id = fopen([opath 'wp2'], 'wb', 'ieee-be');
289     tp2___id = fopen([opath 'tp2'], 'wb', 'ieee-be');
290     sp2___id = fopen([opath 'sp2'], 'wb', 'ieee-be');
291     bp2___id = fopen([opath 'bp2'], 'wb', 'ieee-be');
292     uptp__id = fopen([opath 'uptp'], 'wb', 'ieee-be');
293     vptp__id = fopen([opath 'vptp'], 'wb', 'ieee-be');
294     wptp__id = fopen([opath 'wptp'], 'wb', 'ieee-be');
295     upsp__id = fopen([opath 'upsp'], 'wb', 'ieee-be');
296     vpsp__id = fopen([opath 'vpsp'], 'wb', 'ieee-be');
297     wpsp__id = fopen([opath 'wpsp'], 'wb', 'ieee-be');
298     upbp__id = fopen([opath 'upbp'], 'wb', 'ieee-be');
299     vpbp__id = fopen([opath 'vpbp'], 'wb', 'ieee-be');
300     wpbp__id = fopen([opath 'wpbp'], 'wb', 'ieee-be');
301 edhill 1.5 vpbpdzid = fopen([opath 'vpbp_dbdz'], 'wb', 'ieee-be');
302 edhill 1.6 str___id = fopen([opath 'stress'], 'wb', 'ieee-be');
303 edhill 1.5
304 edhill 1.3 comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
305     '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
306     readslab = inline(comm,'id','nslab','ne');
307    
308     ne = 510;
309     nslab = 6 * ne * ne;
310     nztot = 50;
311     iz = 1;
312     for iz = 1:nztot
313 edhill 1.1
314 edhill 1.6 disp(sprintf(' iz = %d',iz));
315 edhill 1.3 offset = (iz - 1)*nslab*4;
316     for iid = 1:length(iids)
317     fseek(iids(iid), offset, 'bof');
318     end
319     u = readslab(u__id,nslab,ne); u2 = readslab(u2_id,nslab,ne);
320     v = readslab(v__id,nslab,ne); v2 = readslab(v2_id,nslab,ne);
321     um = readslab(um_id,nslab,ne);
322     vm = readslab(vm_id,nslab,ne);
323 edhill 1.4 wm = readslab(wm_id,nslab,ne);
324     % surf(v2(:,:,1)), view(2), shading interp
325     if (iz < nztot)
326     wmp1 = readslab(wm_id,nslab,ne);
327     else
328     wmp1 = zeros(size(wm));
329     end
330     wmp05 = (wm + wmp1)/2.0;
331     w2 = readslab(w2_id,nslab,ne);
332 edhill 1.3 t = readslab(t__id,nslab,ne); t2 = readslab(t2_id,nslab,ne);
333     s = readslab(s__id,nslab,ne); s2 = readslab(s2_id,nslab,ne);
334     b = readslab(b__id,nslab,ne); b2 = readslab(b2_id,nslab,ne);
335    
336     % "simple squared" quantities
337     up2 = u2 - u.^2;
338     vp2 = v2 - v.^2;
339     wp2 = w2 - wm.^2;
340     fwrite(up2___id, up2, 'real*4');
341     fwrite(vp2___id, vp2, 'real*4');
342     fwrite(wp2___id, wp2, 'real*4');
343     clear up2 vp2 wp2 u2 v2 w2
344     tp2 = t2 - t.^2;
345     sp2 = s2 - s.^2;
346     bp2 = b2 - b.^2;
347     fwrite(tp2___id, tp2, 'real*4');
348     fwrite(sp2___id, sp2, 'real*4');
349     fwrite(bp2___id, bp2, 'real*4');
350     clear tp2 sp2 bp2 t2 s2 b2
351     ut = readslab(ut_id,nslab,ne); vt = readslab(vt_id,nslab,ne);
352     us = readslab(us_id,nslab,ne); vs = readslab(vs_id,nslab,ne);
353     ub = readslab(ub_id,nslab,ne); vb = readslab(vb_id,nslab,ne);
354 edhill 1.5 drhodr = readslab(dr_id,nslab,ne);
355 edhill 1.3 [ tonu, tonv ] = mass_on_u_v(t);
356     [ sonu, sonv ] = mass_on_u_v(s);
357     [ bonu, bonv ] = mass_on_u_v(b);
358     uptp = ut - u .* tonu; vptp = vt - v .* tonv;
359     upsp = us - u .* sonu; vpsp = vs - v .* sonv;
360     upbp = ub - u .* bonu; vpbp = vb - v .* bonv;
361 edhill 1.5
362 edhill 1.6 % llupbp = upbp .* llux + vpbp .* llvx;
363 edhill 1.5 llvpbp = upbp .* lluy + vpbp .* llvy;
364 edhill 1.6 if iz > 1
365     % ave_llupbp = (llupbp + old_llupbp)/2.0;
366     ave_llvpbp = (llvpbp + old_llvpbp)/2.0;
367     tmp = drhodr;
368     tmp(find(tmp == 0.0)) = 1.0;
369     vpbpdbdz = ave_llvpbp ./ tmp;
370     fwrite(vpbpdzid, vpbpdbdz, 'real*4');
371     fac = 1000 * 2*pi/(24*3600);
372     stress = fac * sin(pi*YC/180) .* vpbpdbdz;
373     fwrite(str___id, stress, 'real*4');
374     end
375     % old_llupbp = llupbp;
376     old_llvpbp = llvpbp;
377 edhill 1.5
378 edhill 1.3 fwrite(uptp__id, uptp, 'real*4');
379     fwrite(vptp__id, vptp, 'real*4');
380     fwrite(upsp__id, upsp, 'real*4');
381     fwrite(vpsp__id, vpsp, 'real*4');
382     fwrite(upbp__id, upbp, 'real*4');
383     fwrite(vpbp__id, vpbp, 'real*4');
384 edhill 1.5 clear uptp vptp upsp vpsp upbp vpbp vpbpdbdz
385 edhill 1.4 wt = readslab(wt_id,nslab,ne); wptp = wt - wmp05 .* t;
386     ws = readslab(ws_id,nslab,ne); wpsp = ws - wmp05 .* s;
387     wb = readslab(wb_id,nslab,ne); wpbp = wb - wmp05 .* b;
388     fwrite(wptp__id, wptp, 'real*4');
389     fwrite(wpsp__id, wpsp, 'real*4');
390     fwrite(wpbp__id, wpbp, 'real*4');
391 edhill 1.3 end
392 edhill 1.6
393     ne = 510;
394     nz = 1;
395     nr = 50;
396     nrm1 = nr - 1;
397     nlat = 181; nlatm1 = nlat - 1;
398     % save indicies for zonal averages
399     hvals = linspace(-90,90,nlat);
400     i = 2;
401     for i = 2:nlat
402     inds = find(hvals(i-1)<YC & YC<hvals(i));
403     comm = sprintf('inds%04d = uint32(inds);',i-1);
404     eval(comm);
405     end
406    
407     comm = [ 'reshape(fread( id ,nslab,''real*4'',0,' ...
408     '''ieee-be''),[ne ne 6]);' ];
409     readcubelev = inline(comm,'id','nslab','ne');
410    
411     % Zonally average vpbp
412     acc = zeros(nlatm1, nrm1);
413     num = zeros(size(acc));
414     zid = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
415     for iz = 1:50,
416     disp(sprintf('iz = %d',iz));
417     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
418     vpbp = readcubelev(zid,nslab,ne);
419     % surf(log(vpbp(:,:,1))), view(2), shading interp
420     for jj = 1:nlatm1
421     eval( sprintf('clear inds; inds = inds%04d;',jj) );
422     tmp = vpbp(inds);
423     nzinds = find(tmp ~= 0.0);
424     num(jj,iz) = length(nzinds);
425     acc(jj,iz) = sum(tmp(nzinds));
426     end
427     end
428     fclose(zid);
429     zvpbp = acc ./ num;
430     % surf(log(zvpbp')), view(2), shading interp
431     % save primes_92_04/za_vpbp.mat zvpbp
432    
433     % zonally average vpbp_dbdz
434     acc = zeros(nlatm1, nrm1);
435     num = zeros(size(acc));
436     zid = fopen('primes_92_04/vpbp_dbdz', 'r', 'ieee-be');
437     iz = 1;
438 edhill 1.7 for iz = 1:49,
439 edhill 1.6 disp(sprintf('iz = %d',iz));
440     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
441     vpbpdbdz = readcubelev(zid,nslab,ne);
442     for jj = 1:nlatm1
443     eval( sprintf('clear inds; inds = inds%04d;',jj) );
444     tmp = vpbpdbdz(inds);
445     nzinds = find(tmp ~= 0.0);
446     num(jj,iz) = length(nzinds);
447     acc(jj,iz) = sum(tmp(nzinds));
448     end
449     end
450     fclose(zid);
451     za_vpbp_dbdz = acc ./ num;
452     % save primes_92_04/za_vpbp_dbdz.mat za_vpbp_dbdz
453 edhill 1.3
454 edhill 1.6 % zonally average stress
455     acc = zeros(nlatm1, nrm1);
456     num = zeros(size(acc));
457 edhill 1.7 zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
458 edhill 1.6 iz = 1;
459     for iz = 1:49,
460     disp(sprintf('iz = %d',iz));
461 edhill 1.7 fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
462     stress = readcubelev(zid,nslab,ne);
463 edhill 1.6 for jj = 1:nlatm1
464     eval( sprintf('clear inds; inds = inds%04d;',jj) );
465     tmp = stress(inds);
466     nzinds = find(tmp ~= 0.0);
467     num(jj,iz) = length(nzinds);
468     acc(jj,iz) = sum(tmp(nzinds));
469     end
470     end
471     fclose(zid);
472     stress = acc ./ num;
473 edhill 1.7 % save primes_92_04/stress.mat stress
474    
475    
476 edhill 1.4
477     clear all
478     close all
479 edhill 1.3
480 edhill 1.4 %==================================================================
481     % Read the tile00?.mitgrid files
482     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
483     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
484    
485     ne = 510;
486     nep1 = ne + 1;
487     iface = 1;
488     for iface = 1:6
489     fname = sprintf('grid/tile%03d.mitgrid', iface);
490     gid = fopen(fname, 'r', 'ieee-be');
491     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
492     fclose(gid);
493     % surf(tmp(:,:,1)), view(2), shading interp
494     % for jj = 1:length(gvars)
495     for jj = 1:7
496     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
497     [gvars{jj}], iface, jj);
498     eval(comm);
499     end
500     end
501     % surf(XC(:,:,1)), view(2), shading interp
502     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
503     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
504     % surf(YC(:,:,1)), view(2), shading interp
505     % surf(XG(:,:,1)), view(2), shading interp
506     % surf(YG(:,:,1)), view(2), shading interp
507     is = [1:ne];
508     vs = { 'XC','YC','DXF','DYF','RA' };
509     for i = 1:length(vs)
510     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
511     end
512    
513     delR = [ ...
514     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
515     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
516     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
517     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
518     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
519     341.50,364.50,387.50,410.50,433.50,456.50 ];
520     R = cumsum(delR) - 0.5*delR;
521    
522     n1 = ne - 1;
523     dux = zeros(size(XC));
524     duy = zeros(size(XC));
525     dvx = zeros(size(XC));
526     dvy = zeros(size(XC));
527     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
528     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
529     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
530     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
531     dux = dux + 360*double(dux < 180);
532     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
533     duy = duy + 360*double(duy < 180);
534     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
535     dvx = dvx + 360*double(dvx < 180);
536     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
537     dvy = dvy + 360*double(dvy < 180);
538     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
539     llux = dux ./ sqrt(dux.^2 + duy.^2);
540     lluy = duy ./ sqrt(dux.^2 + duy.^2);
541     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
542     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
543    
544     %==================================================================
545     % Project fields to lower-res 1-degree Lat-Lon and write
546     % as NetCDF for viewing with Ingrid
547     %
548     ne = 510;
549     nf = 6;
550     nz = 50;
551     nslab = ne*ne*nf;
552     adir = 'primes_92_04';
553     lat = [-90:90];
554     lon = [0:360];
555     ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
556     xc360 = XC + 180;
557    
558     % ! rm -f cube_20a_primes.nc
559     nc = netcdf(['cube_20a_primes.nc'], 'clobber');
560     nc.reference = [ 'Results from Dimitris Menemenlis''' ...
561     '"cube 20a" integrations' ];
562     nc.author = 'Ed Hill <eh3@mit.edu>';
563     nc.date = 'Feb 28, 2005';
564     nc('X') = length(lon);
565     nc('Y') = length(lat);
566     nc('Z') = length(ir);
567     nc{'X'} = 'X';
568     nc{'Y'} = 'Y';
569     nc{'Z'} = 'Z';
570     nc{'X'}.uniquename = 'X';
571     nc{'X'}.long_name = 'longitude';
572     nc{'X'}.gridtype = ncint(1);
573     nc{'X'}.units = 'degree_east';
574     nc{'Y'}.uniquename = 'Y';
575     nc{'Y'}.long_name = 'latitude';
576     nc{'Y'}.gridtype = ncint(0);
577     nc{'Y'}.units = 'degree_north';
578     nc{'Z'}.uniquename = 'Z';
579     nc{'Z'}.long_name = 'depth';
580     nc{'Z'}.gridtype = ncint(0);
581     nc{'Z'}.units = 'm';
582     nc{'X'}(:) = lon;
583     nc{'Y'}(:) = lat;
584     nc{'Z'}(:) = R(ir);
585    
586 edhill 1.6 f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'} };
587 edhill 1.4
588     ifg = 1;
589     for ifg = 1:length(f_s_3d)
590    
591     acell = f_s_3d{ifg};
592     tname = acell{1};
593     disp([ ' ' tname ' :' ]);
594    
595     fname = sprintf('%s/%s',adir,tname);
596     fid = fopen(fname,'r','ieee-be');
597    
598     id = tname;
599     nc{ id } = { 'Z' 'Y' 'X' };
600     nc{ id }.missing_value = ncdouble(NaN);
601     nc{ id }.FillValue_ = ncdouble(0.0);
602    
603     ii = 1;
604     for ii = 1:length(ir)
605    
606     iz = ir(ii);
607     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
608    
609     fseek(fid,nslab*4*(iz-1),'bof');
610     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
611     tr = reshape(tmp,[ 510 510 6 ]);
612     % surf(tr(:,:,1)), view(2), shading interp
613     trn = tr;
614     trn(find(tr == 0.0)) = NaN;
615     clear tmp tr
616     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
617     v = ll_regrid(xc360,YC,trn,lon,lat);
618     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
619    
620     nc{ id }(ii,:,:) = permute(v,[2 1]);
621    
622     end
623    
624     fclose(fid);
625    
626     end
627    
628     id = 'sum_up2_vp2';
629     nc{ id } = { 'Z' 'Y' 'X' };
630     nc{ id }.missing_value = ncdouble(NaN);
631     nc{ id }.FillValue_ = ncdouble(0.0);
632     fidu = fopen(sprintf('%s/%s',adir,'up2'),'r','ieee-be');
633     fidv = fopen(sprintf('%s/%s',adir,'vp2'),'r','ieee-be');
634     for ii = 1:length(ir)
635     iz = ir(ii);
636     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
637     fseek(fidu,nslab*4*(iz-1),'bof');
638     fseek(fidv,nslab*4*(iz-1),'bof');
639     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
640     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
641     trnu = tru; trnu(find(tru == 0.0)) = NaN;
642     trnv = trv; trnv(find(trv == 0.0)) = NaN;
643     clear tmp tru trv
644     lluv = ll_regrid(xc360,YC,trnu+trnv,lon,lat);
645     nc{ id }(ii,:,:) = permute(lluv,[2 1]);
646     end
647     fclose(fidu);
648     fclose(fidv);
649    
650 edhill 1.5 f_v_3d = { {'up2','vp2'}, ...
651     {'uptp','vptp'}, {'upsp','vpsp'}, {'upbp','vpbp'} };
652 edhill 1.4 for ip = 1:length(f_v_3d)
653     cell = f_v_3d{ip};
654     idu = cell{1};
655     idv = cell{2};
656 edhill 1.5 disp([' ' idu ' ' idv]);
657 edhill 1.4 nc{ idu } = { 'Z' 'Y' 'X' };
658     nc{ idu }.missing_value = ncdouble(NaN);
659     nc{ idu }.FillValue_ = ncdouble(0.0);
660     nc{ idv } = { 'Z' 'Y' 'X' };
661     nc{ idv }.missing_value = ncdouble(NaN);
662     nc{ idv }.FillValue_ = ncdouble(0.0);
663     fidu = fopen(sprintf('%s/%s',adir,idu),'r','ieee-be');
664     fidv = fopen(sprintf('%s/%s',adir,idv),'r','ieee-be');
665     for ii = 1:length(ir)
666     iz = ir(ii);
667     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
668     fseek(fidu,nslab*4*(iz-1),'bof');
669     fseek(fidv,nslab*4*(iz-1),'bof');
670     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
671     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
672     trnu = tru; trnu(find(tru == 0.0)) = NaN;
673     trnv = trv; trnv(find(trv == 0.0)) = NaN;
674     clear tmp tru trv
675     llru = trnu .* llux + trnv .* llvx;
676     llrv = trnu .* lluy + trnv .* llvy;
677     llu = ll_regrid(xc360,YC,llru,lon,lat);
678     llv = ll_regrid(xc360,YC,llrv,lon,lat);
679     nc{ idu }(ii,:,:) = permute(llu,[2 1]);
680     nc{ idv }(ii,:,:) = permute(llv,[2 1]);
681     end
682     end
683     fclose(fidu);
684     fclose(fidv);
685 edhill 1.5
686    
687 edhill 1.4
688     nc = close(nc);
689    
690     % ! scp cube_20a_primes.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_20a/primes/
691     % ! mv cube_20a_primes.nc primes_92_04
692    
693 edhill 1.1
694    
695    
696    
697    
698 edhill 1.5 %------- UNUSED -----------------------
699 edhill 1.1
700    
701     %=======================================================
702     % Calculate : (v'B')/(dB/dz)
703    
704     bid = fopen( 'b-ave-1992-2003', 'r', 'ieee-be');
705     bm1id = fopen('bm1-ave-1992-2003', 'r', 'ieee-be');
706    
707     iz = 1;
708     offset = (iz - 1)*nps*4;
709     fseek(upbpid, offset, 'bof');
710     fseek(vpbpid, offset, 'bof');
711     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
712     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
713     llvpbp = upbp .* lluy + vpbp .* llvy;
714     llvpbpm1 = llvpbp;
715     ne = 510; tx = 85; ty = 85; nt = 216;
716    
717     iz = 2;
718     iz = 15;
719     iz = 16;
720     for iz = 2:50,
721     izm1 = iz - 1;
722     disp(sprintf('iz = %d',iz));
723     % iz = 21;
724     offset = (iz - 1)*nps*4;
725     fseek(upbpid, offset, 'bof');
726     fseek(vpbpid, offset, 'bof');
727     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
728     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
729    
730     % v'B' on ll coords
731     llvpbp = upbp .* lluy + vpbp .* llvy;
732    
733     % B and Bm1
734     fseek(bid, offset, 'bof');
735     b = reshape(fread(bid, nps, 'real*4'),ne,ne,6);
736     offm1 = (iz - 2)*nps*4;
737     fseek(bm1id, offm1, 'bof');
738     bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);
739    
740     % ( v'B' )/( dB/dz )
741     dbdz = (b - bm1)/(R(iz) - R(izm1));
742     ind0 = find(dbdz==0.0);
743     dbdz(ind0) = 1.0;
744     rdbdz = 1.0 / (dbdz);
745     rdbdz(ind0) = 0.0;
746     vpbpdbdz = 0.5*(llvpbp + llvpbpm1) .* rdbdz;
747    
748     % Write the results
749     mid = fopen('vpbpdbdz-1992-2003', 'a', 'ieee-be');
750     fwrite(mid, vpbpdbdz, 'real*4');
751     fclose(mid);
752    
753     % Plot results
754     clow = [ -10 ]; % -20;
755     chigh = [ 10 ]; % 20;
756     ll = zeros(6*510,510);
757     for i = 1:6
758     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
759     yb = 1; ye = 510;
760     ll(xb:xe,yb:ye) = vpbpdbdz(:,:,i);
761     end
762     shift=-1;
763     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
764     title([ '(v''B'')/(dB/dz) at ' ...
765     sprintf('%g',Rmid(iz)) ...
766     'm depth on the 510x510x6 cubesphere for 1992--2003 ["cube5"]']);
767     % print('-painters', '-dpng', '-r650', ...
768     % ['vpTpdTdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
769     print('-painters', '-dpng', '-r150', ...
770     ['vpBpdBdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
771    
772    
773     % Stress: calc, write, and plot
774     stress = 1000 * (2*pi/(24*3600)*sin(pi*yc/180));
775     stress = stress .* vpbpdbdz;
776     sid = fopen('stress-b-1992-2003', 'a', 'ieee-be');
777     fwrite(sid, stress, 'real*4');
778     fclose(sid);
779     clow = []; % [ -1 ];
780     chigh = []; % [ 1 ];
781     ll = zeros(6*510,510);
782     for i = 1:6
783     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
784     yb = 1; ye = 510;
785     ll(xb:xe,yb:ye) = stress(:,:,i);
786     end
787     shift=-1;
788     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
789     title([ 'Stress at ' ...
790     sprintf('%g',Rmid(iz)) ...
791     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
792     % print('-painters', '-dpng', '-r650', ...
793     % ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
794     print('-painters', '-dpng', '-r150', ...
795     ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
796    
797    
798     % Next level
799     llvpbpm1 = llvpbp;
800    
801     end
802    
803    
804     %=======================================================
805     % Zonally average vpbpdbdz and stress
806    
807     clear all ; close all
808     tx = 85;
809     ty = 85;
810     nt = 216;
811     cx = 510;
812     cy = 510;
813     nz = 1;
814     ne = 510;
815     nps = ne * ne * 6;
816    
817     % XCS YCS XGS YGS
818     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
819     for in = 1:4
820     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
821     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
822     tx*nt,ty), ne, tx );
823     fclose(uid);
824     eval([lower(fnam(in,:)) ' = phi;']);
825     a = zeros(6*510,510);
826     for i = 1:6
827     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
828     yb = 1; ye = 510;
829     a(xb:xe,yb:ye) = phi(:,:,i);
830     end
831     eval([lower(fnam(in,:)) 's = a;'])
832     end
833     xcs = xcs + -360.0*(xcs > 180.0);
834     xgs = xgs + -360.0*(xgs > 180.0);
835     clear phi a
836    
837     nz = 1;
838     nr = 50;
839     nrm1 = nr - 1;
840     nlat = 181; nlatm1 = nlat - 1;
841    
842     hvals = linspace(-90,90,nlat);
843     % save indicies for zonal averages
844     i = 2;
845     for i = 2:nlat
846     inds = find(hvals(i-1)<yg & yg<hvals(i));
847     comm = sprintf('inds%04d = uint32(inds);',i-1);
848     eval(comm);
849     end
850    
851     % Zonally average vpbp
852     acc = zeros(nlatm1, nrm1);
853     num = zeros(size(acc));
854     ne = 510;
855     nps = ne * ne * 6;
856     zid = fopen('vpbpll-1992-2003', 'r', 'ieee-be');
857     iz = 1;
858     for iz = 1:50,
859     disp(sprintf('iz = %d',iz));
860     offset = (iz - 1)*nps*4;
861     fseek(zid, offset, 'bof');
862     vpbp = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
863    
864     for jj = 1:nlatm1
865     eval( sprintf('clear inds; inds = inds%04d;',jj) );
866     tmp = vpbp(inds);
867     nzinds = find(tmp ~= 0.0);
868     num(jj,iz) = length(nzinds);
869     acc(jj,iz) = sum(tmp(nzinds));
870     end
871     end
872     fclose(zid);
873     zvpbp = acc ./ num;
874     % save zvpbp zvpbp
875    
876     % zonally average vpbpdbdz
877     acc = zeros(nlatm1, nrm1);
878     num = zeros(size(acc));
879     ne = 510;
880     nps = ne * ne * 6;
881     zid = fopen('vpbpdbdz-1992-2003', 'r', 'ieee-be');
882     iz = 1;
883     for iz = 1:49,
884     disp(sprintf('iz = %d',iz));
885     offset = (iz - 1)*nps*4;
886     fseek(zid, offset, 'bof');
887     vpbpdbdz = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
888    
889     for jj = 1:nlatm1
890     eval( sprintf('clear inds; inds = inds%04d;',jj) );
891     tmp = vpbpdbdz(inds);
892     nzinds = find(tmp ~= 0.0);
893     num(jj,iz) = length(nzinds);
894     acc(jj,iz) = sum(tmp(nzinds));
895     end
896     end
897     fclose(zid);
898     zvpbpdbdz = acc ./ num;
899     % save zvpbpdbdz zvpbpdbdz
900    
901     % zonally average stress
902     acc = zeros(nlatm1, nrm1);
903     num = zeros(size(acc));
904     ne = 510;
905     nps = ne * ne * 6;
906     zid = fopen('stress-b-1992-2003', 'r', 'ieee-be');
907     iz = 1;
908     for iz = 1:49,
909     disp(sprintf('iz = %d',iz));
910     offset = (iz - 1)*nps*4;
911     fseek(zid, offset, 'bof');
912     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
913     jj = 22;
914     for jj = 1:nlatm1
915     eval( sprintf('clear inds; inds = inds%04d;',jj) );
916     tmp = stress(inds);
917     nzinds = find(tmp ~= 0.0);
918     num(jj,iz) = length(nzinds);
919     acc(jj,iz) = sum(tmp(nzinds));
920     end
921     end
922     fclose(zid);
923     stress = acc ./ num;
924     % save stress stress
925    
926     %------- UNUSED -----------------------
927     % zonally average u
928     acc = zeros(nlatm1, 50);
929     num = zeros(size(acc));
930     ne = 510;
931     nps = ne * ne * 6;
932     zid = fopen('Ull_ave_1994--2003', 'r', 'ieee-be');
933     iz = 1;
934     for iz = 1:50,
935     disp(sprintf('iz = %d',iz));
936     offset = (iz - 1)*nps*4;
937     fseek(zid, offset, 'bof');
938     U = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
939     jj = 22;
940     for jj = 1:nlatm1
941     eval( sprintf('clear inds; inds = inds%04d;',jj) );
942     tmp = U(inds);
943     nzinds = find(tmp ~= 0.0);
944     num(jj,iz) = length(nzinds);
945     acc(jj,iz) = sum(tmp(nzinds));
946     end
947     end
948     fclose(zid);
949     zonalu = acc ./ num;
950     % save zonalu zonalu
951     %------- UNUSED -----------------------
952    
953     % zonally average B
954     acc = zeros(nlatm1, 50);
955     num = zeros(size(acc));
956     ne = 510;
957     nps = ne * ne * 6;
958     zid = fopen('b-ave-1992-2003', 'r', 'ieee-be');
959     iz = 1;
960     for iz = 1:50,
961     disp(sprintf('iz = %d',iz));
962     offset = (iz - 1)*nps*4;
963     fseek(zid, offset, 'bof');
964     B = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
965     jj = 22;
966     for jj = 1:nlatm1
967     eval( sprintf('clear inds; inds = inds%04d;',jj) );
968     tmp = B(inds);
969     nzinds = find(tmp ~= 0.0);
970     num(jj,iz) = length(nzinds);
971     acc(jj,iz) = sum(tmp(nzinds));
972     end
973     end
974     fclose(zid);
975     zonalB = acc ./ num;
976     % save zonalB zonalB
977    
978     %------- UNUSED -----------------------
979     % Vertical gradient of zonally averaged B (dB/dz)
980     load zonalT
981     load allR
982     nz = size(zonalT,1);
983     dzaTdz = zeros(nz,size(zonalT,2)-1);
984     for iz = 2:50
985     dzaTdz(:,iz-1) = (zonalT(:,iz) - zonalT(:,iz-1)) ./ (R(iz) - R(iz-1));
986     end
987     % save dzaTdz dzaTdz
988     %------- UNUSED -----------------------
989    
990    
991     delR = [
992     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
993     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
994     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
995     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
996     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
997     341.50,364.50,387.50,410.50,433.50,456.50 ];
998     R = cumsum(delR) - 0.5*delR;
999     Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
1000     hmid = [ -89.5:1:89.5 ];
1001     % save allR delR R Rmid hmid
1002    
1003     %------- UNUSED -----------------------
1004     surf(hmid,-Rmid,num'), view(2), shading interp, colorbar
1005     caxis([ -1 1 ])
1006     caxis('manual')
1007     surf(hmid,-Rmid,mu'), view(2), shading interp, colorbar
1008     title('Zonally Averaged Stress for 1994--2003 [cube5]')
1009     xlabel('Latitude [deg]')
1010     zlabel('Depth [m]')
1011     axis([ -90 90 -5500 200 ])
1012    
1013     % print -dps -painters t_001.ps
1014     print('-painters', '-dpng', '-r150', 'za_str_94--03_r150.png')
1015     print('-painters', '-dpng', '-r650', 'za_str_94--03_r650.png')
1016    
1017     % plot and print bar(U)
1018     load zonalu
1019     surf(hmid,-R,zonalu'), view(2), shading interp, colorbar
1020     contour(hmid,-R,mu'), colorbar
1021     ac = 25;
1022     ac = [-0.3:.005:.3];
1023     [c,h] = contourf(hmid,-R,mu',ac); colorbar
1024     hold on, contour(hmid,-R,mu',ac), hold off
1025     title('Zonally Averaged U for 1994--2003 [cube5]')
1026     xlabel('Latitude [deg]')
1027     zlabel('Depth [m]')
1028     axis([ -90 90 -6000 200 ])
1029     print('-painters', '-dpng', '-r150', 'Ull_94--03_r150.png')
1030     print('-painters', '-dpng', '-r650', 'Ull_94--03_r650.png')
1031     %------- UNUSED -----------------------
1032    
1033    

  ViewVC Help
Powered by ViewVC 1.1.22