/[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.6 - (hide annotations) (download)
Thu Mar 3 21:05:24 2005 UTC (20 years, 5 months ago) by edhill
Branch: MAIN
Changes since 1.5: +105 -14 lines
 o add stress plus numerous small fixes

1 edhill 1.1 %=======================================================
2     %
3 edhill 1.5 % $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/eddy_flux/c20a/plot_c20a.m,v 1.4 2005/03/01 18:00:12 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     for iz = 1:50,
439     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     zid = fopen('', 'r', 'ieee-be');
458     iz = 1;
459     for iz = 1:49,
460     disp(sprintf('iz = %d',iz));
461     offset = (iz - 1)*nps*4;
462     fseek(zid, offset, 'bof');
463     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
464     jj = 22;
465     for jj = 1:nlatm1
466     eval( sprintf('clear inds; inds = inds%04d;',jj) );
467     tmp = stress(inds);
468     nzinds = find(tmp ~= 0.0);
469     num(jj,iz) = length(nzinds);
470     acc(jj,iz) = sum(tmp(nzinds));
471     end
472     end
473     fclose(zid);
474     stress = acc ./ num;
475     % save stress stress
476    
477 edhill 1.4
478     clear all
479     close all
480 edhill 1.3
481 edhill 1.4 %==================================================================
482     % Read the tile00?.mitgrid files
483     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
484     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
485    
486     ne = 510;
487     nep1 = ne + 1;
488     iface = 1;
489     for iface = 1:6
490     fname = sprintf('grid/tile%03d.mitgrid', iface);
491     gid = fopen(fname, 'r', 'ieee-be');
492     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
493     fclose(gid);
494     % surf(tmp(:,:,1)), view(2), shading interp
495     % for jj = 1:length(gvars)
496     for jj = 1:7
497     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
498     [gvars{jj}], iface, jj);
499     eval(comm);
500     end
501     end
502     % surf(XC(:,:,1)), view(2), shading interp
503     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
504     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
505     % surf(YC(:,:,1)), view(2), shading interp
506     % surf(XG(:,:,1)), view(2), shading interp
507     % surf(YG(:,:,1)), view(2), shading interp
508     is = [1:ne];
509     vs = { 'XC','YC','DXF','DYF','RA' };
510     for i = 1:length(vs)
511     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
512     end
513    
514     delR = [ ...
515     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
516     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
517     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
518     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
519     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
520     341.50,364.50,387.50,410.50,433.50,456.50 ];
521     R = cumsum(delR) - 0.5*delR;
522    
523     n1 = ne - 1;
524     dux = zeros(size(XC));
525     duy = zeros(size(XC));
526     dvx = zeros(size(XC));
527     dvy = zeros(size(XC));
528     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
529     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
530     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
531     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
532     dux = dux + 360*double(dux < 180);
533     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
534     duy = duy + 360*double(duy < 180);
535     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
536     dvx = dvx + 360*double(dvx < 180);
537     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
538     dvy = dvy + 360*double(dvy < 180);
539     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
540     llux = dux ./ sqrt(dux.^2 + duy.^2);
541     lluy = duy ./ sqrt(dux.^2 + duy.^2);
542     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
543     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
544    
545     %==================================================================
546     % Project fields to lower-res 1-degree Lat-Lon and write
547     % as NetCDF for viewing with Ingrid
548     %
549     ne = 510;
550     nf = 6;
551     nz = 50;
552     nslab = ne*ne*nf;
553     adir = 'primes_92_04';
554     lat = [-90:90];
555     lon = [0:360];
556     ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
557     xc360 = XC + 180;
558    
559     % ! rm -f cube_20a_primes.nc
560     nc = netcdf(['cube_20a_primes.nc'], 'clobber');
561     nc.reference = [ 'Results from Dimitris Menemenlis''' ...
562     '"cube 20a" integrations' ];
563     nc.author = 'Ed Hill <eh3@mit.edu>';
564     nc.date = 'Feb 28, 2005';
565     nc('X') = length(lon);
566     nc('Y') = length(lat);
567     nc('Z') = length(ir);
568     nc{'X'} = 'X';
569     nc{'Y'} = 'Y';
570     nc{'Z'} = 'Z';
571     nc{'X'}.uniquename = 'X';
572     nc{'X'}.long_name = 'longitude';
573     nc{'X'}.gridtype = ncint(1);
574     nc{'X'}.units = 'degree_east';
575     nc{'Y'}.uniquename = 'Y';
576     nc{'Y'}.long_name = 'latitude';
577     nc{'Y'}.gridtype = ncint(0);
578     nc{'Y'}.units = 'degree_north';
579     nc{'Z'}.uniquename = 'Z';
580     nc{'Z'}.long_name = 'depth';
581     nc{'Z'}.gridtype = ncint(0);
582     nc{'Z'}.units = 'm';
583     nc{'X'}(:) = lon;
584     nc{'Y'}(:) = lat;
585     nc{'Z'}(:) = R(ir);
586    
587 edhill 1.6 f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'} };
588 edhill 1.4
589     ifg = 1;
590     for ifg = 1:length(f_s_3d)
591    
592     acell = f_s_3d{ifg};
593     tname = acell{1};
594     disp([ ' ' tname ' :' ]);
595    
596     fname = sprintf('%s/%s',adir,tname);
597     fid = fopen(fname,'r','ieee-be');
598    
599     id = tname;
600     nc{ id } = { 'Z' 'Y' 'X' };
601     nc{ id }.missing_value = ncdouble(NaN);
602     nc{ id }.FillValue_ = ncdouble(0.0);
603    
604     ii = 1;
605     for ii = 1:length(ir)
606    
607     iz = ir(ii);
608     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
609    
610     fseek(fid,nslab*4*(iz-1),'bof');
611     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
612     tr = reshape(tmp,[ 510 510 6 ]);
613     % surf(tr(:,:,1)), view(2), shading interp
614     trn = tr;
615     trn(find(tr == 0.0)) = NaN;
616     clear tmp tr
617     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
618     v = ll_regrid(xc360,YC,trn,lon,lat);
619     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
620    
621     nc{ id }(ii,:,:) = permute(v,[2 1]);
622    
623     end
624    
625     fclose(fid);
626    
627     end
628    
629     id = 'sum_up2_vp2';
630     nc{ id } = { 'Z' 'Y' 'X' };
631     nc{ id }.missing_value = ncdouble(NaN);
632     nc{ id }.FillValue_ = ncdouble(0.0);
633     fidu = fopen(sprintf('%s/%s',adir,'up2'),'r','ieee-be');
634     fidv = fopen(sprintf('%s/%s',adir,'vp2'),'r','ieee-be');
635     for ii = 1:length(ir)
636     iz = ir(ii);
637     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
638     fseek(fidu,nslab*4*(iz-1),'bof');
639     fseek(fidv,nslab*4*(iz-1),'bof');
640     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
641     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
642     trnu = tru; trnu(find(tru == 0.0)) = NaN;
643     trnv = trv; trnv(find(trv == 0.0)) = NaN;
644     clear tmp tru trv
645     lluv = ll_regrid(xc360,YC,trnu+trnv,lon,lat);
646     nc{ id }(ii,:,:) = permute(lluv,[2 1]);
647     end
648     fclose(fidu);
649     fclose(fidv);
650    
651 edhill 1.5 f_v_3d = { {'up2','vp2'}, ...
652     {'uptp','vptp'}, {'upsp','vpsp'}, {'upbp','vpbp'} };
653 edhill 1.4 for ip = 1:length(f_v_3d)
654     cell = f_v_3d{ip};
655     idu = cell{1};
656     idv = cell{2};
657 edhill 1.5 disp([' ' idu ' ' idv]);
658 edhill 1.4 nc{ idu } = { 'Z' 'Y' 'X' };
659     nc{ idu }.missing_value = ncdouble(NaN);
660     nc{ idu }.FillValue_ = ncdouble(0.0);
661     nc{ idv } = { 'Z' 'Y' 'X' };
662     nc{ idv }.missing_value = ncdouble(NaN);
663     nc{ idv }.FillValue_ = ncdouble(0.0);
664     fidu = fopen(sprintf('%s/%s',adir,idu),'r','ieee-be');
665     fidv = fopen(sprintf('%s/%s',adir,idv),'r','ieee-be');
666     for ii = 1:length(ir)
667     iz = ir(ii);
668     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
669     fseek(fidu,nslab*4*(iz-1),'bof');
670     fseek(fidv,nslab*4*(iz-1),'bof');
671     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
672     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
673     trnu = tru; trnu(find(tru == 0.0)) = NaN;
674     trnv = trv; trnv(find(trv == 0.0)) = NaN;
675     clear tmp tru trv
676     llru = trnu .* llux + trnv .* llvx;
677     llrv = trnu .* lluy + trnv .* llvy;
678     llu = ll_regrid(xc360,YC,llru,lon,lat);
679     llv = ll_regrid(xc360,YC,llrv,lon,lat);
680     nc{ idu }(ii,:,:) = permute(llu,[2 1]);
681     nc{ idv }(ii,:,:) = permute(llv,[2 1]);
682     end
683     end
684     fclose(fidu);
685     fclose(fidv);
686 edhill 1.5
687    
688 edhill 1.4
689     nc = close(nc);
690    
691     % ! scp cube_20a_primes.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_20a/primes/
692     % ! mv cube_20a_primes.nc primes_92_04
693    
694 edhill 1.1
695    
696    
697    
698    
699 edhill 1.5 %------- UNUSED -----------------------
700 edhill 1.1
701    
702     %=======================================================
703     % Calculate : (v'B')/(dB/dz)
704    
705     bid = fopen( 'b-ave-1992-2003', 'r', 'ieee-be');
706     bm1id = fopen('bm1-ave-1992-2003', 'r', 'ieee-be');
707    
708     iz = 1;
709     offset = (iz - 1)*nps*4;
710     fseek(upbpid, offset, 'bof');
711     fseek(vpbpid, offset, 'bof');
712     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
713     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
714     llvpbp = upbp .* lluy + vpbp .* llvy;
715     llvpbpm1 = llvpbp;
716     ne = 510; tx = 85; ty = 85; nt = 216;
717    
718     iz = 2;
719     iz = 15;
720     iz = 16;
721     for iz = 2:50,
722     izm1 = iz - 1;
723     disp(sprintf('iz = %d',iz));
724     % iz = 21;
725     offset = (iz - 1)*nps*4;
726     fseek(upbpid, offset, 'bof');
727     fseek(vpbpid, offset, 'bof');
728     upbp = reshape(fread(upbpid, nps, 'real*4'),ne,ne,6);
729     vpbp = reshape(fread(vpbpid, nps, 'real*4'),ne,ne,6);
730    
731     % v'B' on ll coords
732     llvpbp = upbp .* lluy + vpbp .* llvy;
733    
734     % B and Bm1
735     fseek(bid, offset, 'bof');
736     b = reshape(fread(bid, nps, 'real*4'),ne,ne,6);
737     offm1 = (iz - 2)*nps*4;
738     fseek(bm1id, offm1, 'bof');
739     bm1 =reshape(fread(bm1id, nps, 'real*4'),ne,ne,6);
740    
741     % ( v'B' )/( dB/dz )
742     dbdz = (b - bm1)/(R(iz) - R(izm1));
743     ind0 = find(dbdz==0.0);
744     dbdz(ind0) = 1.0;
745     rdbdz = 1.0 / (dbdz);
746     rdbdz(ind0) = 0.0;
747     vpbpdbdz = 0.5*(llvpbp + llvpbpm1) .* rdbdz;
748    
749     % Write the results
750     mid = fopen('vpbpdbdz-1992-2003', 'a', 'ieee-be');
751     fwrite(mid, vpbpdbdz, 'real*4');
752     fclose(mid);
753    
754     % Plot results
755     clow = [ -10 ]; % -20;
756     chigh = [ 10 ]; % 20;
757     ll = zeros(6*510,510);
758     for i = 1:6
759     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
760     yb = 1; ye = 510;
761     ll(xb:xe,yb:ye) = vpbpdbdz(:,:,i);
762     end
763     shift=-1;
764     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
765     title([ '(v''B'')/(dB/dz) at ' ...
766     sprintf('%g',Rmid(iz)) ...
767     'm depth on the 510x510x6 cubesphere for 1992--2003 ["cube5"]']);
768     % print('-painters', '-dpng', '-r650', ...
769     % ['vpTpdTdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
770     print('-painters', '-dpng', '-r150', ...
771     ['vpBpdBdz_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
772    
773    
774     % Stress: calc, write, and plot
775     stress = 1000 * (2*pi/(24*3600)*sin(pi*yc/180));
776     stress = stress .* vpbpdbdz;
777     sid = fopen('stress-b-1992-2003', 'a', 'ieee-be');
778     fwrite(sid, stress, 'real*4');
779     fclose(sid);
780     clow = []; % [ -1 ];
781     chigh = []; % [ 1 ];
782     ll = zeros(6*510,510);
783     for i = 1:6
784     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
785     yb = 1; ye = 510;
786     ll(xb:xe,yb:ye) = stress(:,:,i);
787     end
788     shift=-1;
789     grph_CS(sq(ll),xcs,ycs,xgs,ygs,[clow],[chigh],shift);
790     title([ 'Stress at ' ...
791     sprintf('%g',Rmid(iz)) ...
792     'm depth on the 510x510x6 cubesphere for 1994--2003 ["cube5"]']);
793     % print('-painters', '-dpng', '-r650', ...
794     % ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_650.png'])
795     print('-painters', '-dpng', '-r150', ...
796     ['stress_' sprintf('%02d_%07.1f',iz,Rmid(iz)) 'm_150.png'])
797    
798    
799     % Next level
800     llvpbpm1 = llvpbp;
801    
802     end
803    
804    
805     %=======================================================
806     % Zonally average vpbpdbdz and stress
807    
808     clear all ; close all
809     tx = 85;
810     ty = 85;
811     nt = 216;
812     cx = 510;
813     cy = 510;
814     nz = 1;
815     ne = 510;
816     nps = ne * ne * 6;
817    
818     % XCS YCS XGS YGS
819     fnam = [ 'XC' ; 'YC'; 'XG'; 'YG' ];
820     for in = 1:4
821     uid = fopen([fnam(in,:) '.data'], 'r', 'ieee-be');
822     phi = unmangleJPL1( reshape(fread(uid, nps, 'real*4'), ...
823     tx*nt,ty), ne, tx );
824     fclose(uid);
825     eval([lower(fnam(in,:)) ' = phi;']);
826     a = zeros(6*510,510);
827     for i = 1:6
828     xb = (i-1)*510 + 1; xe = xb + 510 - 1;
829     yb = 1; ye = 510;
830     a(xb:xe,yb:ye) = phi(:,:,i);
831     end
832     eval([lower(fnam(in,:)) 's = a;'])
833     end
834     xcs = xcs + -360.0*(xcs > 180.0);
835     xgs = xgs + -360.0*(xgs > 180.0);
836     clear phi a
837    
838     nz = 1;
839     nr = 50;
840     nrm1 = nr - 1;
841     nlat = 181; nlatm1 = nlat - 1;
842    
843     hvals = linspace(-90,90,nlat);
844     % save indicies for zonal averages
845     i = 2;
846     for i = 2:nlat
847     inds = find(hvals(i-1)<yg & yg<hvals(i));
848     comm = sprintf('inds%04d = uint32(inds);',i-1);
849     eval(comm);
850     end
851    
852     % Zonally average vpbp
853     acc = zeros(nlatm1, nrm1);
854     num = zeros(size(acc));
855     ne = 510;
856     nps = ne * ne * 6;
857     zid = fopen('vpbpll-1992-2003', 'r', 'ieee-be');
858     iz = 1;
859     for iz = 1:50,
860     disp(sprintf('iz = %d',iz));
861     offset = (iz - 1)*nps*4;
862     fseek(zid, offset, 'bof');
863     vpbp = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
864    
865     for jj = 1:nlatm1
866     eval( sprintf('clear inds; inds = inds%04d;',jj) );
867     tmp = vpbp(inds);
868     nzinds = find(tmp ~= 0.0);
869     num(jj,iz) = length(nzinds);
870     acc(jj,iz) = sum(tmp(nzinds));
871     end
872     end
873     fclose(zid);
874     zvpbp = acc ./ num;
875     % save zvpbp zvpbp
876    
877     % zonally average vpbpdbdz
878     acc = zeros(nlatm1, nrm1);
879     num = zeros(size(acc));
880     ne = 510;
881     nps = ne * ne * 6;
882     zid = fopen('vpbpdbdz-1992-2003', 'r', 'ieee-be');
883     iz = 1;
884     for iz = 1:49,
885     disp(sprintf('iz = %d',iz));
886     offset = (iz - 1)*nps*4;
887     fseek(zid, offset, 'bof');
888     vpbpdbdz = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
889    
890     for jj = 1:nlatm1
891     eval( sprintf('clear inds; inds = inds%04d;',jj) );
892     tmp = vpbpdbdz(inds);
893     nzinds = find(tmp ~= 0.0);
894     num(jj,iz) = length(nzinds);
895     acc(jj,iz) = sum(tmp(nzinds));
896     end
897     end
898     fclose(zid);
899     zvpbpdbdz = acc ./ num;
900     % save zvpbpdbdz zvpbpdbdz
901    
902     % zonally average stress
903     acc = zeros(nlatm1, nrm1);
904     num = zeros(size(acc));
905     ne = 510;
906     nps = ne * ne * 6;
907     zid = fopen('stress-b-1992-2003', 'r', 'ieee-be');
908     iz = 1;
909     for iz = 1:49,
910     disp(sprintf('iz = %d',iz));
911     offset = (iz - 1)*nps*4;
912     fseek(zid, offset, 'bof');
913     stress = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
914     jj = 22;
915     for jj = 1:nlatm1
916     eval( sprintf('clear inds; inds = inds%04d;',jj) );
917     tmp = stress(inds);
918     nzinds = find(tmp ~= 0.0);
919     num(jj,iz) = length(nzinds);
920     acc(jj,iz) = sum(tmp(nzinds));
921     end
922     end
923     fclose(zid);
924     stress = acc ./ num;
925     % save stress stress
926    
927     %------- UNUSED -----------------------
928     % zonally average u
929     acc = zeros(nlatm1, 50);
930     num = zeros(size(acc));
931     ne = 510;
932     nps = ne * ne * 6;
933     zid = fopen('Ull_ave_1994--2003', 'r', 'ieee-be');
934     iz = 1;
935     for iz = 1:50,
936     disp(sprintf('iz = %d',iz));
937     offset = (iz - 1)*nps*4;
938     fseek(zid, offset, 'bof');
939     U = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
940     jj = 22;
941     for jj = 1:nlatm1
942     eval( sprintf('clear inds; inds = inds%04d;',jj) );
943     tmp = U(inds);
944     nzinds = find(tmp ~= 0.0);
945     num(jj,iz) = length(nzinds);
946     acc(jj,iz) = sum(tmp(nzinds));
947     end
948     end
949     fclose(zid);
950     zonalu = acc ./ num;
951     % save zonalu zonalu
952     %------- UNUSED -----------------------
953    
954     % zonally average B
955     acc = zeros(nlatm1, 50);
956     num = zeros(size(acc));
957     ne = 510;
958     nps = ne * ne * 6;
959     zid = fopen('b-ave-1992-2003', 'r', 'ieee-be');
960     iz = 1;
961     for iz = 1:50,
962     disp(sprintf('iz = %d',iz));
963     offset = (iz - 1)*nps*4;
964     fseek(zid, offset, 'bof');
965     B = reshape(fread(zid, nps, 'real*4'),ne,ne,6);
966     jj = 22;
967     for jj = 1:nlatm1
968     eval( sprintf('clear inds; inds = inds%04d;',jj) );
969     tmp = B(inds);
970     nzinds = find(tmp ~= 0.0);
971     num(jj,iz) = length(nzinds);
972     acc(jj,iz) = sum(tmp(nzinds));
973     end
974     end
975     fclose(zid);
976     zonalB = acc ./ num;
977     % save zonalB zonalB
978    
979     %------- UNUSED -----------------------
980     % Vertical gradient of zonally averaged B (dB/dz)
981     load zonalT
982     load allR
983     nz = size(zonalT,1);
984     dzaTdz = zeros(nz,size(zonalT,2)-1);
985     for iz = 2:50
986     dzaTdz(:,iz-1) = (zonalT(:,iz) - zonalT(:,iz-1)) ./ (R(iz) - R(iz-1));
987     end
988     % save dzaTdz dzaTdz
989     %------- UNUSED -----------------------
990    
991    
992     delR = [
993     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
994     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
995     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
996     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
997     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
998     341.50,364.50,387.50,410.50,433.50,456.50 ];
999     R = cumsum(delR) - 0.5*delR;
1000     Rmid = R(1:(length(R)-1)) + 0.5*diff(R);
1001     hmid = [ -89.5:1:89.5 ];
1002     % save allR delR R Rmid hmid
1003    
1004     %------- UNUSED -----------------------
1005     surf(hmid,-Rmid,num'), view(2), shading interp, colorbar
1006     caxis([ -1 1 ])
1007     caxis('manual')
1008     surf(hmid,-Rmid,mu'), view(2), shading interp, colorbar
1009     title('Zonally Averaged Stress for 1994--2003 [cube5]')
1010     xlabel('Latitude [deg]')
1011     zlabel('Depth [m]')
1012     axis([ -90 90 -5500 200 ])
1013    
1014     % print -dps -painters t_001.ps
1015     print('-painters', '-dpng', '-r150', 'za_str_94--03_r150.png')
1016     print('-painters', '-dpng', '-r650', 'za_str_94--03_r650.png')
1017    
1018     % plot and print bar(U)
1019     load zonalu
1020     surf(hmid,-R,zonalu'), view(2), shading interp, colorbar
1021     contour(hmid,-R,mu'), colorbar
1022     ac = 25;
1023     ac = [-0.3:.005:.3];
1024     [c,h] = contourf(hmid,-R,mu',ac); colorbar
1025     hold on, contour(hmid,-R,mu',ac), hold off
1026     title('Zonally Averaged U for 1994--2003 [cube5]')
1027     xlabel('Latitude [deg]')
1028     zlabel('Depth [m]')
1029     axis([ -90 90 -6000 200 ])
1030     print('-painters', '-dpng', '-r150', 'Ull_94--03_r150.png')
1031     print('-painters', '-dpng', '-r650', 'Ull_94--03_r650.png')
1032     %------- UNUSED -----------------------
1033    
1034    

  ViewVC Help
Powered by ViewVC 1.1.22