/[MITgcm]/MITgcm_contrib/high_res_cube/eddy_flux/c22/plot_c22.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/eddy_flux/c22/plot_c22.m

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


Revision 1.10 - (hide annotations) (download)
Wed Apr 20 21:43:11 2005 UTC (20 years, 3 months ago) by edhill
Branch: MAIN
Changes since 1.9: +61 -15 lines
 o more stream-wise averages

1 edhill 1.1 %=======================================================
2     %
3 edhill 1.10 % $Id: plot_c22.m,v 1.9 2005/04/14 01:40:29 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_22" or "c22" integration.
12    
13     % ssh eddy
14 edhill 1.4 % cd /r/r0/edhill/eddy_stats/c22
15 edhill 1.1
16     % matlab -nojvm
17     % matlab -nojvm -nodisplay
18    
19     clear all
20     close all
21    
22    
23     %==================================================================
24     % Read the tile00?.mitgrid files
25     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
26     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
27    
28     ne = 510;
29     nep1 = ne + 1;
30     iface = 1;
31     for iface = 1:6
32     fname = sprintf('grid/tile%03d.mitgrid', iface);
33     gid = fopen(fname, 'r', 'ieee-be');
34     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
35     fclose(gid);
36     % surf(tmp(:,:,1)), view(2), shading interp
37     % for jj = 1:length(gvars)
38     for jj = 1:7
39     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
40     [gvars{jj}], iface, jj);
41     eval(comm);
42     end
43     end
44     % surf(XC(:,:,1)), view(2), shading interp
45     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
46     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
47     % surf(YC(:,:,1)), view(2), shading interp
48     % surf(XG(:,:,1)), view(2), shading interp
49     % surf(YG(:,:,1)), view(2), shading interp
50     is = [1:ne];
51     vs = { 'XC','YC','DXF','DYF','RA' };
52     for i = 1:length(vs)
53     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
54     end
55    
56     delR = [ ...
57     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
58     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
59     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
60     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
61     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
62     341.50,364.50,387.50,410.50,433.50,456.50 ];
63     R = cumsum(delR) - 0.5*delR;
64    
65     %==================================================================
66     % Project fields to lower-res 1-degree Lat-Lon and write
67     % as NetCDF for viewing with Ingrid
68     %
69     fields_3d = { ...
70     'DRHODR', 'RHOANOSQ', 'RHOAnoma', 'SALT', 'SALTSQ', ...
71     'THETA', 'THETASQ', 'URHOMASS', ...
72     'USLTMASS', 'UTHMASS', 'UVEL', 'UVELMASS', 'UVELSQ', ...
73     'UV_VEL_Z', 'VRHOMASS', 'VSLTMASS', 'VTHMASS', ...
74     'VVEL', 'VVELMASS', 'VVELSQ', 'WRHOMASS', 'WSLTMASS', ...
75     'WTHMASS', 'WU_VEL', 'WVELMASS', 'WVELSQ', 'WV_VEL' };
76     fields_2d = { ...
77     'ETAN', 'ETANSQ', 'EmPmRtave', 'PHIBOT', 'SFLUX', 'SRELAX', ...
78     'TAUX', 'TAUY', 'TFLUX', ...
79     'TICE', 'TRELAX', 'UICEtave', 'VICEtave' };
80    
81     ne = 510;
82     nf = 6;
83     nz = 50;
84     nslab = ne*ne*nf;
85     adir = 'ave_1992-2004';
86     lat = [-90:90];
87     lon = [0:360];
88     ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
89    
90     % ! rm -f cube_22_at1deg.nc
91     nc = netcdf(['cube_22_at1deg.nc'], 'clobber');
92     nc.reference = [ 'The cube_22 averages from Dimitris Menemenlis' ...
93     ' regridded to 1-deg Lat-Lon' ];
94     nc.author = 'Ed Hill <eh3@mit.edu>';
95     nc.date = 'Mar 27, 2005';
96     nc('X') = length(lon);
97     nc('Y') = length(lat);
98     nc('Z') = length(ir);
99     nc{'X'} = 'X';
100     nc{'Y'} = 'Y';
101     nc{'Z'} = 'Z';
102     nc{'X'}.uniquename = 'X';
103     nc{'X'}.long_name = 'longitude';
104     nc{'X'}.gridtype = ncint(1);
105     nc{'X'}.units = 'degree_east';
106     nc{'Y'}.uniquename = 'Y';
107     nc{'Y'}.long_name = 'latitude';
108     nc{'Y'}.gridtype = ncint(0);
109     nc{'Y'}.units = 'degree_north';
110     nc{'Z'}.uniquename = 'Z';
111     nc{'Z'}.long_name = 'depth';
112     nc{'Z'}.gridtype = ncint(0);
113     nc{'Z'}.units = 'm';
114     nc{'X'}(:) = lon;
115     nc{'Y'}(:) = lat;
116     nc{'Z'}(:) = R(ir);
117    
118     ifld = 5;
119     fields = union(fields_2d, fields_3d);
120     for ifld = 1:length(fields)
121    
122     id = fields{ifld};
123     if ismember(fields{ifld},fields_3d)
124     ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
125     nc{ id } = { 'Z' 'Y' 'X' };
126     else
127     ir = [ 1 ];
128     nc{ id } = { 'Y' 'X' };
129     end
130     nc{ id }.missing_value = ncdouble(NaN);
131     nc{ id }.FillValue_ = ncdouble(0.0);
132    
133     disp([ ' ' fields{ifld} ' :' ]);
134    
135     fname = sprintf('%s/%s.ave',adir,fields{ifld});
136     fid = fopen(fname,'r','ieee-be');
137    
138     ii = 1;
139     for ii = 1:length(ir)
140    
141     iz = ir(ii);
142     disp(sprintf(' iz = %g',iz));
143    
144     fseek(fid,nslab*4*(iz-1),'bof');
145     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
146     tr = permute(reshape(tmp,[ 510 6 510 ]),[1 3 2]);
147     % surf(tr(:,:,1)), view(2), shading interp
148     xc360 = XC + 180;
149     trn = tr;
150     trn(find(tr == 0.0)) = NaN;
151     clear tmp tr
152     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
153     v = ll_regrid(xc360,YC,trn,lon,lat);
154     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
155    
156     if length(ir) == 1
157     nc{ id }(:,:) = permute(v,[2 1]);
158     else
159     nc{ id }(ii,:,:) = permute(v,[2 1]);
160     end
161    
162     end
163    
164     fclose(fid);
165    
166     end
167    
168     nc = close(nc);
169    
170    
171     % ! ncdump cube_22_at1deg.nc | more
172     % ! scp cube_22_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
173    
174    
175    
176    
177     %=======================================================
178     % Compute [uvw]'[tsb]'
179    
180     clear all
181     close all
182    
183     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
184     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
185     ne = 510;
186     nep1 = ne + 1;
187     iface = 1;
188     for iface = 1:6
189     fname = sprintf('grid/tile%03d.mitgrid', iface);
190     gid = fopen(fname, 'r', 'ieee-be');
191     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
192     fclose(gid);
193     % surf(tmp(:,:,1)), view(2), shading interp
194     % for jj = 1:length(gvars)
195     for jj = 1:7
196     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
197     [gvars{jj}], iface, jj);
198     eval(comm);
199     end
200     end
201     % surf(XC(:,:,1)), view(2), shading interp
202     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
203     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
204     % surf(YC(:,:,1)), view(2), shading interp
205     % surf(XG(:,:,1)), view(2), shading interp
206     % surf(YG(:,:,1)), view(2), shading interp
207     is = [1:ne];
208     vs = { 'XC','YC','DXF','DYF','RA' };
209     for i = 1:length(vs)
210     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
211     end
212    
213     delR = [ ...
214     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
215     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
216     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
217     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
218     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
219     341.50,364.50,387.50,410.50,433.50,456.50 ];
220     R = cumsum(delR) - 0.5*delR;
221    
222     n1 = ne - 1;
223     dux = zeros(size(XC));
224     duy = zeros(size(XC));
225     dvx = zeros(size(XC));
226     dvy = zeros(size(XC));
227     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
228     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
229     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
230     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
231     dux = dux + 360*double(dux < 180);
232     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
233     duy = duy + 360*double(duy < 180);
234     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
235     dvx = dvx + 360*double(dvx < 180);
236     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
237     dvy = dvy + 360*double(dvy < 180);
238     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
239 edhill 1.10 dut = sqrt(dux.^2 + duy.^2);
240     dvt = sqrt(dvx.^2 + dvy.^2);
241     llux = dux ./ dut;
242     lluy = duy ./ dut;
243     llvx = dvx ./ dvt;
244     llvy = dvy ./ dvt;
245     clear XC XG YG DXF DYF RA tmp
246     clear dux duy dvx dvy
247 edhill 1.1
248 edhill 1.2 lpath = 'ave_1992-2004/';
249 edhill 1.1 u__id = fopen( [lpath 'UVEL.ave'], 'r', 'ieee-be'); % 1
250     v__id = fopen( [lpath 'VVEL.ave'], 'r', 'ieee-be'); % 2
251     %w__id = fopen( [lpath 'WVEL.ave'], 'r', 'ieee-be');
252     u2_id = fopen( [lpath 'UVELSQ.ave'], 'r', 'ieee-be'); % 3
253     v2_id = fopen( [lpath 'VVELSQ.ave'], 'r', 'ieee-be'); % 4
254     w2_id = fopen( [lpath 'WVELSQ.ave'], 'r', 'ieee-be'); % 5
255     um_id = fopen( [lpath 'UVELMASS.ave'], 'r', 'ieee-be'); % 6
256     vm_id = fopen( [lpath 'VVELMASS.ave'], 'r', 'ieee-be'); % 7
257     wm_id = fopen( [lpath 'WVELMASS.ave'], 'r', 'ieee-be'); % 8
258     t__id = fopen( [lpath 'THETA.ave'], 'r', 'ieee-be'); % 9
259     t2_id = fopen( [lpath 'THETASQ.ave'], 'r', 'ieee-be'); % 10
260     s__id = fopen( [lpath 'SALT.ave'], 'r', 'ieee-be'); % 11
261     s2_id = fopen( [lpath 'SALTSQ.ave'], 'r', 'ieee-be'); % 12
262     b__id = fopen( [lpath 'RHOAnoma.ave'], 'r', 'ieee-be'); % 13
263     b2_id = fopen( [lpath 'RHOANOSQ.ave'], 'r', 'ieee-be'); % 14
264     ut_id = fopen( [lpath 'UTHMASS.ave'], 'r', 'ieee-be'); % 15
265     vt_id = fopen( [lpath 'VTHMASS.ave'], 'r', 'ieee-be'); % 16
266     wt_id = fopen( [lpath 'WTHMASS.ave'], 'r', 'ieee-be'); % 17
267     us_id = fopen( [lpath 'USLTMASS.ave'], 'r', 'ieee-be'); % 18
268     vs_id = fopen( [lpath 'VSLTMASS.ave'], 'r', 'ieee-be'); % 19
269     ws_id = fopen( [lpath 'WSLTMASS.ave'], 'r', 'ieee-be'); % 20
270     ub_id = fopen( [lpath 'URHOMASS.ave'], 'r', 'ieee-be'); % 21
271     vb_id = fopen( [lpath 'VRHOMASS.ave'], 'r', 'ieee-be'); % 22
272     wb_id = fopen( [lpath 'WRHOMASS.ave'], 'r', 'ieee-be'); % 23
273     dr_id = fopen( [lpath 'DRHODR.ave'], 'r', 'ieee-be'); % 24
274     iids = [ u__id v__id u2_id v2_id w2_id um_id vm_id wm_id ...
275     t__id t2_id s__id s2_id b__id b2_id ...
276     ut_id vt_id wt_id us_id vs_id ws_id ...
277     ub_id vb_id wb_id dr_id ];
278    
279     % ! rm -rf primes_92_04 ; mkdir primes_92_04
280     opath = 'primes_92_04/';
281     up2___id = fopen([opath 'up2'], 'wb', 'ieee-be');
282     vp2___id = fopen([opath 'vp2'], 'wb', 'ieee-be');
283     wp2___id = fopen([opath 'wp2'], 'wb', 'ieee-be');
284     tp2___id = fopen([opath 'tp2'], 'wb', 'ieee-be');
285     sp2___id = fopen([opath 'sp2'], 'wb', 'ieee-be');
286     bp2___id = fopen([opath 'bp2'], 'wb', 'ieee-be');
287     uptp__id = fopen([opath 'uptp'], 'wb', 'ieee-be');
288     vptp__id = fopen([opath 'vptp'], 'wb', 'ieee-be');
289     wptp__id = fopen([opath 'wptp'], 'wb', 'ieee-be');
290     upsp__id = fopen([opath 'upsp'], 'wb', 'ieee-be');
291     vpsp__id = fopen([opath 'vpsp'], 'wb', 'ieee-be');
292     wpsp__id = fopen([opath 'wpsp'], 'wb', 'ieee-be');
293     upbp__id = fopen([opath 'upbp'], 'wb', 'ieee-be');
294     vpbp__id = fopen([opath 'vpbp'], 'wb', 'ieee-be');
295     wpbp__id = fopen([opath 'wpbp'], 'wb', 'ieee-be');
296     vpbpdzid = fopen([opath 'vpbp_dbdz'], 'wb', 'ieee-be');
297     str___id = fopen([opath 'stress'], 'wb', 'ieee-be');
298 edhill 1.10 dbdy__id = fopen([opath 'dbdy'], 'wb', 'ieee-be');
299     K_____id = fopen([opath 'K'], 'wb', 'ieee-be');
300 edhill 1.1
301     comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
302     '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
303     readslab = inline(comm,'id','nslab','ne');
304    
305     ne = 510;
306     nslab = 6 * ne * ne;
307     nztot = 50;
308     iz = 1;
309     for iz = 1:nztot
310    
311     disp(sprintf(' iz = %d',iz));
312     offset = (iz - 1)*nslab*4;
313     for iid = 1:length(iids)
314     fseek(iids(iid), offset, 'bof');
315     end
316     u = readslab(u__id,nslab,ne); u2 = readslab(u2_id,nslab,ne);
317     v = readslab(v__id,nslab,ne); v2 = readslab(v2_id,nslab,ne);
318     um = readslab(um_id,nslab,ne);
319     vm = readslab(vm_id,nslab,ne);
320     wm = readslab(wm_id,nslab,ne);
321 edhill 1.7 % surf(v2(:,:,1)), view(2), shading interp, colorbar
322 edhill 1.1 if (iz < nztot)
323     wmp1 = readslab(wm_id,nslab,ne);
324     else
325     wmp1 = zeros(size(wm));
326     end
327     wmp05 = (wm + wmp1)/2.0;
328     w2 = readslab(w2_id,nslab,ne);
329     t = readslab(t__id,nslab,ne); t2 = readslab(t2_id,nslab,ne);
330     s = readslab(s__id,nslab,ne); s2 = readslab(s2_id,nslab,ne);
331     b = readslab(b__id,nslab,ne); b2 = readslab(b2_id,nslab,ne);
332    
333     % "simple squared" quantities
334     up2 = u2 - u.^2;
335     vp2 = v2 - v.^2;
336     wp2 = w2 - wm.^2;
337     fwrite(up2___id, up2, 'real*4');
338     fwrite(vp2___id, vp2, 'real*4');
339     fwrite(wp2___id, wp2, 'real*4');
340     clear up2 vp2 wp2 u2 v2 w2
341     tp2 = t2 - t.^2;
342     sp2 = s2 - s.^2;
343     bp2 = b2 - b.^2;
344     fwrite(tp2___id, tp2, 'real*4');
345     fwrite(sp2___id, sp2, 'real*4');
346     fwrite(bp2___id, bp2, 'real*4');
347     clear tp2 sp2 bp2 t2 s2 b2
348     ut = readslab(ut_id,nslab,ne); vt = readslab(vt_id,nslab,ne);
349     us = readslab(us_id,nslab,ne); vs = readslab(vs_id,nslab,ne);
350     ub = readslab(ub_id,nslab,ne); vb = readslab(vb_id,nslab,ne);
351     drhodr = readslab(dr_id,nslab,ne);
352 edhill 1.7 nan_inds = find(t == 0.0);
353     t(nan_inds) = NaN; s(nan_inds) = NaN; b(nan_inds) = NaN;
354 edhill 1.1 [ tonu, tonv ] = mass_on_u_v(t);
355     [ sonu, sonv ] = mass_on_u_v(s);
356     [ bonu, bonv ] = mass_on_u_v(b);
357 edhill 1.8 uptp = ut - um .* tonu; vptp = vt - vm .* tonv;
358     upsp = us - um .* sonu; vpsp = vs - vm .* sonv;
359     upbp = ub - um .* bonu; vpbp = vb - vm .* bonv;
360 edhill 1.7 % surf(upbp(:,:,1)), shading interp, view(2)
361     % caxis([-0.1 0.1]), colorbar
362 edhill 1.1
363     % llupbp = upbp .* llux + vpbp .* llvx;
364     llvpbp = upbp .* lluy + vpbp .* llvy;
365     if iz > 1
366     % ave_llupbp = (llupbp + old_llupbp)/2.0;
367     ave_llvpbp = (llvpbp + old_llvpbp)/2.0;
368     tmp = drhodr;
369 edhill 1.6 % tmp(find(tmp == 0.0)) = 1.0;
370 edhill 1.1 vpbpdbdz = ave_llvpbp ./ tmp;
371     fwrite(vpbpdzid, vpbpdbdz, 'real*4');
372 edhill 1.6 % \tau_x = 2\Omega sin(\phi) = 4\pi*sin(\phi)/(24*3600)
373     fac = 1000 * 4*pi/(24*3600);
374 edhill 1.1 stress = fac * sin(pi*YC/180) .* vpbpdbdz;
375     fwrite(str___id, stress, 'real*4');
376     end
377     % old_llupbp = llupbp;
378     old_llvpbp = llvpbp;
379 edhill 1.10
380     % determine diffusivities
381     dbdy = calc_dbdy(b, dut,dvt, lluy,llvy);
382     diffus = - vpbp ./ dbdy;
383 edhill 1.1
384     fwrite(uptp__id, uptp, 'real*4');
385     fwrite(vptp__id, vptp, 'real*4');
386     fwrite(upsp__id, upsp, 'real*4');
387     fwrite(vpsp__id, vpsp, 'real*4');
388     fwrite(upbp__id, upbp, 'real*4');
389     fwrite(vpbp__id, vpbp, 'real*4');
390     clear uptp vptp upsp vpsp upbp vpbp vpbpdbdz
391     wt = readslab(wt_id,nslab,ne); wptp = wt - wmp05 .* t;
392     ws = readslab(ws_id,nslab,ne); wpsp = ws - wmp05 .* s;
393     wb = readslab(wb_id,nslab,ne); wpbp = wb - wmp05 .* b;
394     fwrite(wptp__id, wptp, 'real*4');
395     fwrite(wpsp__id, wpsp, 'real*4');
396     fwrite(wpbp__id, wpbp, 'real*4');
397 edhill 1.10 fwrite(dbdy__id, dbdy, 'real*4');
398     fwrite(K_____id, diffus, 'real*4');
399 edhill 1.1 end
400 edhill 1.10
401     clear uptp vptp upsp vpsp ubbp vpbp
402     clear wptp wpsp wpbp dbdy diffus
403 edhill 1.1
404     ne = 510;
405     nz = 1;
406     nr = 50;
407     nrm1 = nr - 1;
408     nlat = 181; nlatm1 = nlat - 1;
409     % save indicies for zonal averages
410     hvals = linspace(-90,90,nlat);
411     i = 2;
412     for i = 2:nlat
413     inds = find(hvals(i-1)<YC & YC<hvals(i));
414     comm = sprintf('inds%04d = uint32(inds);',i-1);
415     eval(comm);
416     end
417    
418     comm = [ 'reshape(fread( id ,nslab,''real*4'',0,' ...
419     '''ieee-be''),[ne ne 6]);' ];
420     readcubelev = inline(comm,'id','nslab','ne');
421    
422     % Zonally average the ll_vpbp
423 edhill 1.8 clear acc num
424 edhill 1.1 acc = zeros(nlatm1, nrm1);
425     num = zeros(size(acc));
426     zidu = fopen('primes_92_04/upbp', 'r', 'ieee-be');
427     zidv = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
428     for iz = 1:50,
429     disp(sprintf('iz = %d',iz));
430     fseek(zidu, (iz - 1)*(ne*ne*6)*4, 'bof');
431     fseek(zidv, (iz - 1)*(ne*ne*6)*4, 'bof');
432     upbp = readcubelev(zidu,nslab,ne);
433 edhill 1.2 vpbp = readcubelev(zidv,nslab,ne);
434     llupbp = upbp .* llux + vpbp .* llvx;
435 edhill 1.1 llvpbp = upbp .* lluy + vpbp .* llvy;
436 edhill 1.8 % surf(llupbp(:,:,1)), view(2), shading interp, caxis([-.02 .02])
437     % surf(llvpbp(:,:,1)), view(2), shading interp, caxis([-.02 .02])
438 edhill 1.1 for jj = 1:nlatm1
439     eval( sprintf('clear inds; inds = inds%04d;',jj) );
440     tmp = llvpbp(inds);
441 edhill 1.8 % nzinds = find(tmp ~= 0.0);
442     finds = find(isfinite(tmp));
443     num(jj,iz) = length(finds);
444     acc(jj,iz) = sum(tmp(finds));
445 edhill 1.1 end
446     end
447     fclose(zidu);
448     fclose(zidv);
449     llzvpbp = acc ./ num;
450 edhill 1.8 % surf(flipud(llzvpbp')), view(2), shading interp, caxis([-.1 .1])
451 edhill 1.4 % ! rm -f primes_92_04/za_llvpbp.mat
452     % save primes_92_04/za_llvpbp.mat llzvpbp
453 edhill 1.1
454     % zonally average ll_vpbp_dbdz
455 edhill 1.9 clear acc num
456 edhill 1.1 acc = zeros(nlatm1, nrm1);
457     num = zeros(size(acc));
458     zid = fopen('primes_92_04/vpbp_dbdz', 'r', 'ieee-be');
459     iz = 1;
460     for iz = 1:49,
461     disp(sprintf('iz = %d',iz));
462     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
463     vpbpdbdz = readcubelev(zid,nslab,ne);
464     % surf(vpbpdbdz(:,:,1)), view(2), shading interp, caxis([-50 50])
465     for jj = 1:nlatm1
466     eval( sprintf('clear inds; inds = inds%04d;',jj) );
467     tmp = vpbpdbdz(inds);
468 edhill 1.9 nzinds = find(isfinite(tmp));
469 edhill 1.1 num(jj,iz) = length(nzinds);
470     acc(jj,iz) = sum(tmp(nzinds));
471     end
472     end
473     fclose(zid);
474     za_ll_vpbp_dbdz = acc ./ num;
475 edhill 1.4 % ! rm -f primes_92_04/za_ll_vpbp_dbdz.mat
476 edhill 1.1 % save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz
477    
478     % zonally average stress
479 edhill 1.9 clear acc num
480 edhill 1.1 acc = zeros(nlatm1, nrm1);
481     num = zeros(size(acc));
482     zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
483     iz = 1;
484     for iz = 1:49,
485     disp(sprintf('iz = %d',iz));
486     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
487     stress = readcubelev(zid,nslab,ne);
488     for jj = 1:nlatm1
489     eval( sprintf('clear inds; inds = inds%04d;',jj) );
490     tmp = stress(inds);
491 edhill 1.9 nzinds = find(isfinite(tmp));
492 edhill 1.1 num(jj,iz) = length(nzinds);
493     acc(jj,iz) = sum(tmp(nzinds));
494     end
495     end
496     fclose(zid);
497 edhill 1.3 za_stress = acc ./ num;
498 edhill 1.4 % ! rm -f primes_92_04/stress.mat
499 edhill 1.3 % save primes_92_04/stress.mat za_stress
500 edhill 1.1
501 edhill 1.3 % Average the stress over SSH-contours (streamlines)
502 edhill 1.6 clear num acc
503 edhill 1.3 eta_id = fopen('ave_1992-2004/ETAN.ave', 'r', 'ieee-be');
504     comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
505     '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
506     readslab = inline(comm,'id','nslab','ne');
507     etan = readslab(eta_id,nslab,ne);
508     fclose(eta_id);
509     % surf(YC(:,:,1)'), view(2),shading interp,colorbar
510     % surf(etan(:,:,1)'), view(2),shading interp,colorbar
511     nssh = 50;
512     sshvals = linspace(-2.3,1.2,nssh);
513     i = 2;
514     for i = 2:nssh
515     inds = find(sshvals(i-1)<etan & etan<sshvals(i) & YC<-20);
516     comm = sprintf('sshinds%04d = uint32(inds);',i-1);
517     eval(comm);
518     end
519     zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
520     iz = 1;
521     for iz = 1:49,
522     disp(sprintf('iz = %d',iz));
523     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
524     stress = readcubelev(zid,nslab,ne);
525     for jj = 1:(nssh-1)
526     eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
527     tmp = stress(inds);
528 edhill 1.9 nzinds = find(isfinite(tmp));
529 edhill 1.3 num(jj,iz) = length(nzinds);
530     acc(jj,iz) = sum(tmp(nzinds));
531     end
532     end
533     fclose(zid);
534     ssha_stress = acc ./ num;
535     % surf(log(ssha_stress')), view(2),shading interp,colorbar
536 edhill 1.4 % ! rm -f primes_92_04/ssha_stress.mat
537     % save primes_92_04/ssha_stress.mat ssha_stress sshvals
538 edhill 1.2
539 edhill 1.10 % Average the b, ull, and vpbpll over SSH-contours (streamlines)
540     clear n_u n_b n_vb a_u a_b a_vb
541     zid_u = fopen('ave_1992-2004/UVELMASS.ave', 'r', 'ieee-be');
542     zid_b = fopen('ave_1992-2004/RHOAnoma.ave', 'r', 'ieee-be');
543     zid_ub = fopen('primes_92_04/upbp', 'r', 'ieee-be');
544     zid_vb = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
545     iz = 1;
546     for iz = 1:49,
547     disp(sprintf('iz = %d',iz));
548     fseek(zid_u, (iz - 1)*(ne*ne*6)*4, 'bof');
549     fseek(zid_v, (iz - 1)*(ne*ne*6)*4, 'bof');
550     fseek(zid_ub, (iz - 1)*(ne*ne*6)*4, 'bof');
551     fseek(zid_vb, (iz - 1)*(ne*ne*6)*4, 'bof');
552     u = readslab(zid_u, nslab,ne);
553     b = readslab(zid_b, nslab,ne);
554     upbp = readcubelev(zid_ub, nslab,ne);
555     vpbp = readcubelev(zid_vb, nslab,ne);
556     % llupbp = upbp .* llux + vpbp .* llvx;
557     llvpbp = upbp .* lluy + vpbp .* llvy;
558     for jj = 1:(nssh-1)
559     eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
560     t_u = u(inds);
561     t_b = b(inds);
562     t_vb = llvpbp(inds);
563     i_u = find(isfinite( t_u ));
564     i_b = find(isfinite( t_b ));
565     i_vb = find(isfinite( t_vb ));
566     n_u(jj,iz) = length( i_u );
567     n_b(jj,iz) = length( i_b );
568     n_vb(jj,iz) = length( i_vb );
569     a_u(jj,iz) = sum( t_u( i_u ));
570     a_b(jj,iz) = sum( t_b( i_b ));
571     a_vb(jj,iz) = sum( t_vb( i_vb ));
572     end
573     end
574     fclose(zid);
575     ssha_stress = acc ./ num;
576     % save primes_92_04/ssha_u_b_vpbp.mat ssha_ll_u ssha_b ssha_vpbp
577    
578 edhill 1.1
579     clear all
580     close all
581    
582     %==================================================================
583     % Read the tile00?.mitgrid files
584     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
585     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
586    
587     ne = 510;
588     nep1 = ne + 1;
589     iface = 1;
590     for iface = 1:6
591     fname = sprintf('grid/tile%03d.mitgrid', iface);
592     gid = fopen(fname, 'r', 'ieee-be');
593     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
594     fclose(gid);
595     % surf(tmp(:,:,1)), view(2), shading interp
596     % for jj = 1:length(gvars)
597     for jj = 1:7
598     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
599     [gvars{jj}], iface, jj);
600     eval(comm);
601     end
602     end
603     % surf(XC(:,:,1)), view(2), shading interp
604     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
605     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
606     % surf(YC(:,:,1)), view(2), shading interp
607     % surf(XG(:,:,1)), view(2), shading interp
608     % surf(YG(:,:,1)), view(2), shading interp
609     is = [1:ne];
610     vs = { 'XC','YC','DXF','DYF','RA' };
611     for i = 1:length(vs)
612     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
613     end
614    
615     delR = [ ...
616     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
617     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
618     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
619     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
620     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
621     341.50,364.50,387.50,410.50,433.50,456.50 ];
622     R = cumsum(delR) - 0.5*delR;
623 edhill 1.4 Ri = R(1:(length(R)-1)) ...
624     + 0.25*delR(1:(length(R)-1)) + 0.25*delR(2:(length(R)));
625    
626 edhill 1.1 n1 = ne - 1;
627     dux = zeros(size(XC));
628     duy = zeros(size(XC));
629     dvx = zeros(size(XC));
630     dvy = zeros(size(XC));
631     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
632     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
633     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
634     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
635     dux = dux + 360*double(dux < 180);
636     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
637     duy = duy + 360*double(duy < 180);
638     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
639     dvx = dvx + 360*double(dvx < 180);
640     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
641     dvy = dvy + 360*double(dvy < 180);
642     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
643     llux = dux ./ sqrt(dux.^2 + duy.^2);
644     lluy = duy ./ sqrt(dux.^2 + duy.^2);
645     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
646     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
647    
648     %==================================================================
649     % Project fields to lower-res 1-degree Lat-Lon and write
650     % as NetCDF for viewing with Ingrid
651     %
652     ne = 510;
653     nf = 6;
654     nz = 50;
655     nslab = ne*ne*nf;
656     adir = 'primes_92_04';
657     lat = [-90:90];
658     lon = [0:360];
659 edhill 1.2 ir = [ 1 2 3 5 10 15 20 25 30 35 40 ];
660 edhill 1.1 xc360 = XC + 180;
661    
662 edhill 1.2 % ! rm -f cube_22_primes_at1deg.nc
663     nc = netcdf(['cube_22_primes_at1deg.nc'], 'clobber');
664     nc.reference = [ 'The cube_22 primes from Dimitris Menemenlis' ...
665     ' regridded to 1-deg Lat-Lon' ];
666 edhill 1.1 nc.author = 'Ed Hill <eh3@mit.edu>';
667 edhill 1.2 nc.date = 'March 27, 2005';
668 edhill 1.1 nc('X') = length(lon);
669     nc('Y') = length(lat);
670     nc('Z') = length(ir);
671 edhill 1.4 nc('Zc') = length(R);
672     nc('Zi') = length(Ri);
673 edhill 1.1 nc{'X'} = 'X';
674     nc{'Y'} = 'Y';
675     nc{'Z'} = 'Z';
676 edhill 1.4 nc{'Zc'} = 'Zc';
677     nc{'Zi'} = 'Zi';
678 edhill 1.1 nc{'X'}.uniquename = 'X';
679     nc{'X'}.long_name = 'longitude';
680     nc{'X'}.gridtype = ncint(1);
681     nc{'X'}.units = 'degree_east';
682     nc{'Y'}.uniquename = 'Y';
683     nc{'Y'}.long_name = 'latitude';
684     nc{'Y'}.gridtype = ncint(0);
685     nc{'Y'}.units = 'degree_north';
686     nc{'Z'}.uniquename = 'Z';
687     nc{'Z'}.long_name = 'depth';
688     nc{'Z'}.gridtype = ncint(0);
689     nc{'Z'}.units = 'm';
690 edhill 1.4 nc{'Zc'}.uniquename = 'Zc';
691     nc{'Zc'}.long_name = 'depth_at_center';
692     nc{'Zc'}.gridtype = ncint(0);
693     nc{'Zc'}.units = 'm';
694     nc{'Zi'}.uniquename = 'Zi';
695     nc{'Zi'}.long_name = 'depth_at_interface';
696     nc{'Zi'}.gridtype = ncint(0);
697     nc{'Zi'}.units = 'm';
698 edhill 1.1 nc{'X'}(:) = lon;
699     nc{'Y'}(:) = lat;
700     nc{'Z'}(:) = R(ir);
701 edhill 1.4 nc{'Zc'}(:) = R;
702     nc{'Zi'}(:) = Ri;
703    
704 edhill 1.10 f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'}, ...
705     {'stress'}, {'dbdy'}, {'K'} };
706 edhill 1.1
707     ifg = 1;
708     for ifg = 1:length(f_s_3d)
709     acell = f_s_3d{ifg};
710     tname = acell{1};
711     disp([ ' ' tname ' :' ]);
712     fname = sprintf('%s/%s',adir,tname);
713     fid = fopen(fname,'r','ieee-be');
714     id = tname;
715     nc{ id } = { 'Z' 'Y' 'X' };
716     nc{ id }.missing_value = ncdouble(NaN);
717     nc{ id }.FillValue_ = ncdouble(0.0);
718     ii = 1;
719     for ii = 1:length(ir)
720     iz = ir(ii);
721     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
722    
723     fseek(fid,nslab*4*(iz-1),'bof');
724     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
725     tr = reshape(tmp,[ 510 510 6 ]);
726     % surf(tr(:,:,1)), view(2), shading interp
727     trn = tr;
728     trn(find(tr == 0.0)) = NaN;
729     clear tmp tr
730     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
731     v = ll_regrid(xc360,YC,trn,lon,lat);
732     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
733     nc{ id }(ii,:,:) = permute(v,[2 1]);
734     end
735     fclose(fid);
736     end
737    
738     id = 'sum_up2_vp2';
739     nc{ id } = { 'Z' 'Y' 'X' };
740     nc{ id }.missing_value = ncdouble(NaN);
741     nc{ id }.FillValue_ = ncdouble(0.0);
742     fidu = fopen(sprintf('%s/%s',adir,'up2'),'r','ieee-be');
743     fidv = fopen(sprintf('%s/%s',adir,'vp2'),'r','ieee-be');
744     for ii = 1:length(ir)
745     iz = ir(ii);
746     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
747     fseek(fidu,nslab*4*(iz-1),'bof');
748     fseek(fidv,nslab*4*(iz-1),'bof');
749     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
750     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
751     trnu = tru; trnu(find(tru == 0.0)) = NaN;
752     trnv = trv; trnv(find(trv == 0.0)) = NaN;
753     clear tmp tru trv
754     lluv = ll_regrid(xc360,YC,trnu+trnv,lon,lat);
755     nc{ id }(ii,:,:) = permute(lluv,[2 1]);
756     end
757     fclose(fidu);
758     fclose(fidv);
759    
760     f_v_3d = { {'up2','vp2'}, ...
761     {'uptp','vptp'}, {'upsp','vpsp'}, {'upbp','vpbp'} };
762     for ip = 1:length(f_v_3d)
763     cell = f_v_3d{ip};
764     idu = cell{1};
765     idv = cell{2};
766     disp([' ' idu ' ' idv]);
767     nc{ idu } = { 'Z' 'Y' 'X' };
768     nc{ idu }.missing_value = ncdouble(NaN);
769     nc{ idu }.FillValue_ = ncdouble(0.0);
770     nc{ idv } = { 'Z' 'Y' 'X' };
771     nc{ idv }.missing_value = ncdouble(NaN);
772     nc{ idv }.FillValue_ = ncdouble(0.0);
773     fidu = fopen(sprintf('%s/%s',adir,idu),'r','ieee-be');
774     fidv = fopen(sprintf('%s/%s',adir,idv),'r','ieee-be');
775     for ii = 1:length(ir)
776     iz = ir(ii);
777     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
778     fseek(fidu,nslab*4*(iz-1),'bof');
779     fseek(fidv,nslab*4*(iz-1),'bof');
780     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
781     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
782     trnu = tru; trnu(find(tru == 0.0)) = NaN;
783     trnv = trv; trnv(find(trv == 0.0)) = NaN;
784     clear tmp tru trv
785     llru = trnu .* llux + trnv .* llvx;
786     llrv = trnu .* lluy + trnv .* llvy;
787     llu = ll_regrid(xc360,YC,llru,lon,lat);
788     llv = ll_regrid(xc360,YC,llrv,lon,lat);
789     nc{ idu }(ii,:,:) = permute(llu,[2 1]);
790     nc{ idv }(ii,:,:) = permute(llv,[2 1]);
791     end
792     end
793     fclose(fidu);
794     fclose(fidv);
795     nc = close(nc);
796    
797 edhill 1.4 % === zonal and stream-wise averages ===
798     % save primes_92_04/za_llvpbp.mat llzvpbp
799     % save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz
800     % save primes_92_04/stress.mat za_stress
801     % save primes_92_04/ssha_stress.mat ssha_stress sshvals
802     load primes_92_04/za_llvpbp.mat
803     load primes_92_04/za_ll_vpbp_dbdz.mat
804     load primes_92_04/stress.mat
805     load primes_92_04/ssha_stress.mat
806 edhill 1.5
807     ssh = sshvals(1:(length(sshvals)-1)) + 0.5*diff(sshvals);
808 edhill 1.1
809 edhill 1.4 lat = [-90:90];
810     latm = lat(1:(length(lat)-1)) + 0.5*diff(lat);
811 edhill 1.1
812 edhill 1.4 % ! rm -f cube_22_zsa.nc
813     nc = netcdf(['cube_22_zsa.nc'], 'clobber');
814     nc.reference = [ 'The cube_22 zonal and stream-wise averages.' ];
815     nc.author = 'Ed Hill <eh3@mit.edu>';
816     nc.date = 'March 27, 2005';
817     nc('Y') = length(latm);
818     nc('Zc') = length(R);
819     nc('Zi') = length(Ri);
820 edhill 1.5 nc('SSH') = length(ssh);
821 edhill 1.4 nc{'Y'} = 'Y';
822     nc{'Zc'} = 'Zc';
823     nc{'Zi'} = 'Zi';
824 edhill 1.5 nc{'SSH'} = 'SSH';
825 edhill 1.4 nc{'Y'}.uniquename = 'Y';
826     nc{'Y'}.long_name = 'latitude';
827     nc{'Y'}.gridtype = ncint(0);
828     nc{'Y'}.units = 'degree_north';
829     nc{'Zc'}.uniquename = 'Zc';
830     nc{'Zc'}.long_name = 'depth_at_center';
831     nc{'Zc'}.gridtype = ncint(0);
832     nc{'Zc'}.units = 'm';
833     nc{'Zi'}.uniquename = 'Zi';
834     nc{'Zi'}.long_name = 'depth_at_interface';
835     nc{'Zi'}.gridtype = ncint(0);
836     nc{'Zi'}.units = 'm';
837 edhill 1.5 nc{'SSH'}.uniquename = 'SSH';
838     nc{'SSH'}.long_name = 'sea_surface_height';
839     nc{'SSH'}.gridtype = ncint(0);
840     nc{'SSH'}.units = 'm';
841 edhill 1.4 nc{'Y'}(:) = latm;
842     nc{'Zc'}(:) = R;
843     nc{'Zi'}(:) = Ri;
844 edhill 1.5 nc{'SSH'}(:) = ssh;
845 edhill 1.1
846 edhill 1.4 id = 'llzvpbp';
847     nc{ id } = { 'Zc' 'Y' };
848     nc{ id }.missing_value = ncdouble(NaN);
849     nc{ id }.FillValue_ = ncdouble(0.0);
850     nc{ id }(:) = permute(llzvpbp,[2 1]);
851 edhill 1.1
852 edhill 1.5 f_i = { {'za_ll_vpbp_dbdz'}, {'za_stress'} };
853 edhill 1.4 for ip = 1:length(f_i)
854     cell = f_i{ip};
855     id = cell{1};
856     disp([' ' id]);
857     nc{ id } = { 'Zi' 'Y' };
858     nc{ id }.missing_value = ncdouble(NaN);
859     nc{ id }.FillValue_ = ncdouble(0.0);
860     eval(sprintf('tmp = %s;', id));
861     nc{ id }(:) = permute(tmp,[2 1]);
862 edhill 1.1 end
863 edhill 1.5
864     id = 'ssha_stress';
865     nc{ id } = { 'Zi' 'SSH' };
866     nc{ id }.missing_value = ncdouble(NaN);
867     nc{ id }.FillValue_ = ncdouble(0.0);
868     nc{ id }(:) = permute(ssha_stress,[2 1]);
869    
870 edhill 1.4 nc = close(nc);
871 edhill 1.1
872 edhill 1.6 % surf(,,log(ssha_stress')), view(2),shading interp,colorbar
873    
874 edhill 1.5 % ! scp cube_22_primes_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
875     % ! scp cube_22_zsa.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
876 edhill 1.4 % ! mv cube_22_primes_at1deg.nc primes_92_04
877 edhill 1.1

  ViewVC Help
Powered by ViewVC 1.1.22