/[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.13 - (hide annotations) (download)
Wed Apr 27 22:00:14 2005 UTC (20 years, 3 months ago) by edhill
Branch: MAIN
Changes since 1.12: +56 -42 lines
 o more quantities

1 edhill 1.1 %=======================================================
2     %
3 edhill 1.13 % $Id: plot_c22.m,v 1.12 2005/04/27 17:43:54 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 edhill 1.13 nc{'X'}(:) = lon - 180;
115 edhill 1.1 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     id = fields{ifld};
122     if ismember(fields{ifld},fields_3d)
123     ir = [ 1 2 3 5 10 15 20 25 30 35 40 50 ];
124     nc{ id } = { 'Z' 'Y' 'X' };
125     else
126     ir = [ 1 ];
127     nc{ id } = { 'Y' 'X' };
128     end
129     nc{ id }.missing_value = ncdouble(NaN);
130     nc{ id }.FillValue_ = ncdouble(0.0);
131     disp([ ' ' fields{ifld} ' :' ]);
132     fname = sprintf('%s/%s.ave',adir,fields{ifld});
133     fid = fopen(fname,'r','ieee-be');
134     ii = 1;
135     for ii = 1:length(ir)
136    
137     iz = ir(ii);
138     disp(sprintf(' iz = %g',iz));
139     fseek(fid,nslab*4*(iz-1),'bof');
140     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
141     tr = permute(reshape(tmp,[ 510 6 510 ]),[1 3 2]);
142     % surf(tr(:,:,1)), view(2), shading interp
143     xc360 = XC + 180;
144     trn = tr;
145     trn(find(tr == 0.0)) = NaN;
146     clear tmp tr
147     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
148     v = ll_regrid(xc360,YC,trn,lon,lat);
149     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
150     if length(ir) == 1
151     nc{ id }(:,:) = permute(v,[2 1]);
152     else
153     nc{ id }(ii,:,:) = permute(v,[2 1]);
154     end
155     end
156     fclose(fid);
157     end
158     nc = close(nc);
159     % ! ncdump cube_22_at1deg.nc | more
160     % ! scp cube_22_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
161    
162    
163    
164     %=======================================================
165     % Compute [uvw]'[tsb]'
166    
167     clear all
168     close all
169    
170     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
171     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
172     ne = 510;
173     nep1 = ne + 1;
174     iface = 1;
175     for iface = 1:6
176     fname = sprintf('grid/tile%03d.mitgrid', iface);
177     gid = fopen(fname, 'r', 'ieee-be');
178     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
179     fclose(gid);
180     % surf(tmp(:,:,1)), view(2), shading interp
181     % for jj = 1:length(gvars)
182     for jj = 1:7
183     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
184     [gvars{jj}], iface, jj);
185     eval(comm);
186     end
187     end
188     % surf(XC(:,:,1)), view(2), shading interp
189     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
190     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
191     % surf(YC(:,:,1)), view(2), shading interp
192     % surf(XG(:,:,1)), view(2), shading interp
193     % surf(YG(:,:,1)), view(2), shading interp
194     is = [1:ne];
195     vs = { 'XC','YC','DXF','DYF','RA' };
196     for i = 1:length(vs)
197     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
198     end
199    
200     delR = [ ...
201     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
202     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
203     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
204     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
205     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
206     341.50,364.50,387.50,410.50,433.50,456.50 ];
207     R = cumsum(delR) - 0.5*delR;
208    
209     n1 = ne - 1;
210     dux = zeros(size(XC));
211     duy = zeros(size(XC));
212     dvx = zeros(size(XC));
213     dvy = zeros(size(XC));
214     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
215     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
216     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
217     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
218     dux = dux + 360*double(dux < 180);
219     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
220     duy = duy + 360*double(duy < 180);
221     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
222     dvx = dvx + 360*double(dvx < 180);
223     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
224     dvy = dvy + 360*double(dvy < 180);
225     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
226 edhill 1.10 dut = sqrt(dux.^2 + duy.^2);
227     dvt = sqrt(dvx.^2 + dvy.^2);
228     llux = dux ./ dut;
229     lluy = duy ./ dut;
230     llvx = dvx ./ dvt;
231     llvy = dvy ./ dvt;
232     clear XC XG YG DXF DYF RA tmp
233     clear dux duy dvx dvy
234 edhill 1.1
235 edhill 1.2 lpath = 'ave_1992-2004/';
236 edhill 1.1 u__id = fopen( [lpath 'UVEL.ave'], 'r', 'ieee-be'); % 1
237     v__id = fopen( [lpath 'VVEL.ave'], 'r', 'ieee-be'); % 2
238     u2_id = fopen( [lpath 'UVELSQ.ave'], 'r', 'ieee-be'); % 3
239     v2_id = fopen( [lpath 'VVELSQ.ave'], 'r', 'ieee-be'); % 4
240     w2_id = fopen( [lpath 'WVELSQ.ave'], 'r', 'ieee-be'); % 5
241     um_id = fopen( [lpath 'UVELMASS.ave'], 'r', 'ieee-be'); % 6
242     vm_id = fopen( [lpath 'VVELMASS.ave'], 'r', 'ieee-be'); % 7
243     wm_id = fopen( [lpath 'WVELMASS.ave'], 'r', 'ieee-be'); % 8
244     t__id = fopen( [lpath 'THETA.ave'], 'r', 'ieee-be'); % 9
245     t2_id = fopen( [lpath 'THETASQ.ave'], 'r', 'ieee-be'); % 10
246     s__id = fopen( [lpath 'SALT.ave'], 'r', 'ieee-be'); % 11
247     s2_id = fopen( [lpath 'SALTSQ.ave'], 'r', 'ieee-be'); % 12
248     b__id = fopen( [lpath 'RHOAnoma.ave'], 'r', 'ieee-be'); % 13
249     b2_id = fopen( [lpath 'RHOANOSQ.ave'], 'r', 'ieee-be'); % 14
250     ut_id = fopen( [lpath 'UTHMASS.ave'], 'r', 'ieee-be'); % 15
251     vt_id = fopen( [lpath 'VTHMASS.ave'], 'r', 'ieee-be'); % 16
252     wt_id = fopen( [lpath 'WTHMASS.ave'], 'r', 'ieee-be'); % 17
253     us_id = fopen( [lpath 'USLTMASS.ave'], 'r', 'ieee-be'); % 18
254     vs_id = fopen( [lpath 'VSLTMASS.ave'], 'r', 'ieee-be'); % 19
255     ws_id = fopen( [lpath 'WSLTMASS.ave'], 'r', 'ieee-be'); % 20
256     ub_id = fopen( [lpath 'URHOMASS.ave'], 'r', 'ieee-be'); % 21
257     vb_id = fopen( [lpath 'VRHOMASS.ave'], 'r', 'ieee-be'); % 22
258     wb_id = fopen( [lpath 'WRHOMASS.ave'], 'r', 'ieee-be'); % 23
259     dr_id = fopen( [lpath 'DRHODR.ave'], 'r', 'ieee-be'); % 24
260     iids = [ u__id v__id u2_id v2_id w2_id um_id vm_id wm_id ...
261     t__id t2_id s__id s2_id b__id b2_id ...
262     ut_id vt_id wt_id us_id vs_id ws_id ...
263     ub_id vb_id wb_id dr_id ];
264    
265     % ! rm -rf primes_92_04 ; mkdir primes_92_04
266     opath = 'primes_92_04/';
267     up2___id = fopen([opath 'up2'], 'wb', 'ieee-be');
268     vp2___id = fopen([opath 'vp2'], 'wb', 'ieee-be');
269     wp2___id = fopen([opath 'wp2'], 'wb', 'ieee-be');
270     tp2___id = fopen([opath 'tp2'], 'wb', 'ieee-be');
271     sp2___id = fopen([opath 'sp2'], 'wb', 'ieee-be');
272     bp2___id = fopen([opath 'bp2'], 'wb', 'ieee-be');
273     uptp__id = fopen([opath 'uptp'], 'wb', 'ieee-be');
274     vptp__id = fopen([opath 'vptp'], 'wb', 'ieee-be');
275     wptp__id = fopen([opath 'wptp'], 'wb', 'ieee-be');
276     upsp__id = fopen([opath 'upsp'], 'wb', 'ieee-be');
277     vpsp__id = fopen([opath 'vpsp'], 'wb', 'ieee-be');
278     wpsp__id = fopen([opath 'wpsp'], 'wb', 'ieee-be');
279     upbp__id = fopen([opath 'upbp'], 'wb', 'ieee-be');
280     vpbp__id = fopen([opath 'vpbp'], 'wb', 'ieee-be');
281     wpbp__id = fopen([opath 'wpbp'], 'wb', 'ieee-be');
282     vpbpdzid = fopen([opath 'vpbp_dbdz'], 'wb', 'ieee-be');
283     str___id = fopen([opath 'stress'], 'wb', 'ieee-be');
284 edhill 1.10 dbdy__id = fopen([opath 'dbdy'], 'wb', 'ieee-be');
285     K_____id = fopen([opath 'K'], 'wb', 'ieee-be');
286 edhill 1.1
287     comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
288     '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
289     readslab = inline(comm,'id','nslab','ne');
290    
291     ne = 510;
292     nslab = 6 * ne * ne;
293     nztot = 50;
294     iz = 1;
295     for iz = 1:nztot
296    
297     disp(sprintf(' iz = %d',iz));
298     offset = (iz - 1)*nslab*4;
299     for iid = 1:length(iids)
300     fseek(iids(iid), offset, 'bof');
301     end
302 edhill 1.11 t = readslab(t__id,nslab,ne); t2 = readslab(t2_id,nslab,ne);
303     nan_inds = find(t == 0.0);
304 edhill 1.1 u = readslab(u__id,nslab,ne); u2 = readslab(u2_id,nslab,ne);
305     v = readslab(v__id,nslab,ne); v2 = readslab(v2_id,nslab,ne);
306 edhill 1.11 um = readslab(um_id,nslab,ne); vm = readslab(vm_id,nslab,ne);
307     wm = readslab(wm_id,nslab,ne); w2 = readslab(w2_id,nslab,ne);
308     s = readslab(s__id,nslab,ne); s2 = readslab(s2_id,nslab,ne);
309     b = readslab(b__id,nslab,ne); b2 = readslab(b2_id,nslab,ne);
310 edhill 1.7 % surf(v2(:,:,1)), view(2), shading interp, colorbar
311 edhill 1.11 vars = { 'u','v','um','vm','wm','u2','v2','w2',...
312     's','s2','t','t2','b','b2' };
313     for i = 1:length(vars)
314     eval(sprintf('%s(nan_inds) = NaN;',vars{i}));
315     end
316 edhill 1.1 if (iz < nztot)
317     wmp1 = readslab(wm_id,nslab,ne);
318     else
319     wmp1 = zeros(size(wm));
320     end
321     wmp05 = (wm + wmp1)/2.0;
322    
323     % "simple squared" quantities
324     up2 = u2 - u.^2;
325     vp2 = v2 - v.^2;
326     wp2 = w2 - wm.^2;
327     fwrite(up2___id, up2, 'real*4');
328     fwrite(vp2___id, vp2, 'real*4');
329     fwrite(wp2___id, wp2, 'real*4');
330     clear up2 vp2 wp2 u2 v2 w2
331     tp2 = t2 - t.^2;
332     sp2 = s2 - s.^2;
333     bp2 = b2 - b.^2;
334     fwrite(tp2___id, tp2, 'real*4');
335     fwrite(sp2___id, sp2, 'real*4');
336     fwrite(bp2___id, bp2, 'real*4');
337     clear tp2 sp2 bp2 t2 s2 b2
338 edhill 1.11
339 edhill 1.1 ut = readslab(ut_id,nslab,ne); vt = readslab(vt_id,nslab,ne);
340     us = readslab(us_id,nslab,ne); vs = readslab(vs_id,nslab,ne);
341     ub = readslab(ub_id,nslab,ne); vb = readslab(vb_id,nslab,ne);
342     drhodr = readslab(dr_id,nslab,ne);
343 edhill 1.11 vars = { 'ut','vt','us','vs','ub','vb' };
344     for i = 1:length(vars)
345     eval(sprintf('%s(nan_inds) = NaN;',vars{i}));
346     end
347 edhill 1.1 [ tonu, tonv ] = mass_on_u_v(t);
348     [ sonu, sonv ] = mass_on_u_v(s);
349     [ bonu, bonv ] = mass_on_u_v(b);
350 edhill 1.8 uptp = ut - um .* tonu; vptp = vt - vm .* tonv;
351     upsp = us - um .* sonu; vpsp = vs - vm .* sonv;
352     upbp = ub - um .* bonu; vpbp = vb - vm .* bonv;
353 edhill 1.7 % surf(upbp(:,:,1)), shading interp, view(2)
354     % caxis([-0.1 0.1]), colorbar
355 edhill 1.1
356     % llupbp = upbp .* llux + vpbp .* llvx;
357     llvpbp = upbp .* lluy + vpbp .* llvy;
358     if iz > 1
359     % ave_llupbp = (llupbp + old_llupbp)/2.0;
360     ave_llvpbp = (llvpbp + old_llvpbp)/2.0;
361     tmp = drhodr;
362 edhill 1.6 % tmp(find(tmp == 0.0)) = 1.0;
363 edhill 1.1 vpbpdbdz = ave_llvpbp ./ tmp;
364     fwrite(vpbpdzid, vpbpdbdz, 'real*4');
365 edhill 1.6 % \tau_x = 2\Omega sin(\phi) = 4\pi*sin(\phi)/(24*3600)
366     fac = 1000 * 4*pi/(24*3600);
367 edhill 1.1 stress = fac * sin(pi*YC/180) .* vpbpdbdz;
368     fwrite(str___id, stress, 'real*4');
369     end
370     % old_llupbp = llupbp;
371     old_llvpbp = llvpbp;
372 edhill 1.10
373     % determine diffusivities
374     dbdy = calc_dbdy(b, dut,dvt, lluy,llvy);
375     diffus = - vpbp ./ dbdy;
376 edhill 1.1
377     fwrite(uptp__id, uptp, 'real*4');
378     fwrite(vptp__id, vptp, 'real*4');
379     fwrite(upsp__id, upsp, 'real*4');
380     fwrite(vpsp__id, vpsp, 'real*4');
381     fwrite(upbp__id, upbp, 'real*4');
382     fwrite(vpbp__id, vpbp, 'real*4');
383     clear uptp vptp upsp vpsp upbp vpbp vpbpdbdz
384 edhill 1.11 wt = readslab(wt_id,nslab,ne); wt(nan_inds) = NaN;
385     ws = readslab(ws_id,nslab,ne); ws(nan_inds) = NaN;
386     wb = readslab(wb_id,nslab,ne); wb(nan_inds) = NaN;
387     wptp = wt - wmp05 .* t;
388     wpsp = ws - wmp05 .* s;
389     wpbp = wb - wmp05 .* b;
390 edhill 1.1 fwrite(wptp__id, wptp, 'real*4');
391     fwrite(wpsp__id, wpsp, 'real*4');
392     fwrite(wpbp__id, wpbp, 'real*4');
393 edhill 1.10 fwrite(dbdy__id, dbdy, 'real*4');
394     fwrite(K_____id, diffus, 'real*4');
395 edhill 1.1 end
396 edhill 1.10
397     clear uptp vptp upsp vpsp ubbp vpbp
398     clear wptp wpsp wpbp dbdy diffus
399 edhill 1.11
400     % save current_state_20050422
401     % load current_state_20050422
402    
403     do_plots = 0;
404 edhill 1.1
405     ne = 510;
406     nz = 1;
407     nr = 50;
408     nrm1 = nr - 1;
409     nlat = 181; nlatm1 = nlat - 1;
410     % save indicies for zonal averages
411     hvals = linspace(-90,90,nlat);
412     i = 2;
413     for i = 2:nlat
414     inds = find(hvals(i-1)<YC & YC<hvals(i));
415     comm = sprintf('inds%04d = uint32(inds);',i-1);
416     eval(comm);
417     end
418    
419     comm = [ 'reshape(fread( id ,nslab,''real*4'',0,' ...
420     '''ieee-be''),[ne ne 6]);' ];
421     readcubelev = inline(comm,'id','nslab','ne');
422    
423 edhill 1.11 % Zonally average the ll_upbp,ll_vpbp
424 edhill 1.8 clear acc num
425 edhill 1.11 acc = zeros(nlatm1, nrm1); acc = NaN;
426 edhill 1.1 num = zeros(size(acc));
427 edhill 1.11 %um_id = fopen( [lpath 'UVELMASS.ave'], 'r', 'ieee-be'); % 6
428     %vm_id = fopen( [lpath 'VVELMASS.ave'], 'r', 'ieee-be'); % 7
429     zidu = fopen('primes_92_04/upbp', 'r', 'ieee-be');
430     zidv = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
431 edhill 1.1 for iz = 1:50,
432     disp(sprintf('iz = %d',iz));
433     fseek(zidu, (iz - 1)*(ne*ne*6)*4, 'bof');
434     fseek(zidv, (iz - 1)*(ne*ne*6)*4, 'bof');
435 edhill 1.11 %fseek(um_id, (iz - 1)*(ne*ne*6)*4, 'bof');
436     %fseek(vm_id, (iz - 1)*(ne*ne*6)*4, 'bof');
437 edhill 1.1 upbp = readcubelev(zidu,nslab,ne);
438 edhill 1.2 vpbp = readcubelev(zidv,nslab,ne);
439 edhill 1.11 %um = readslab(um_id,nslab,ne);
440     %vm = readslab(vm_id,nslab,ne);
441 edhill 1.2 llupbp = upbp .* llux + vpbp .* llvx;
442 edhill 1.1 llvpbp = upbp .* lluy + vpbp .* llvy;
443 edhill 1.11 %llum = um .* llux + vm .* llvx;
444     %llvm = um .* lluy + vm .* llvy;
445     if do_plots == 1
446     figure(1)
447     subplot(2,2,1),surf(llupbp(:,:,1))
448     view(2),shading interp,caxis([-.02 .02])
449     subplot(2,2,2),surf(llvpbp(:,:,1))
450     view(2),shading interp,caxis([-.02 .02])
451     figure(2)
452     subplot(2,2,1),surf(um(:,:,6))
453     view(2),shading interp
454     subplot(2,2,2),surf(vm(:,:,6))
455     view(2),shading interp
456     subplot(2,2,3),surf(llum(:,:,6))
457     view(2),shading interp
458     subplot(2,2,4),surf(llvm(:,:,6))
459     view(2),shading interp
460     pause(2)
461     end
462 edhill 1.1 for jj = 1:nlatm1
463     eval( sprintf('clear inds; inds = inds%04d;',jj) );
464     tmp = llvpbp(inds);
465 edhill 1.8 % nzinds = find(tmp ~= 0.0);
466     finds = find(isfinite(tmp));
467     num(jj,iz) = length(finds);
468     acc(jj,iz) = sum(tmp(finds));
469 edhill 1.1 end
470     end
471 edhill 1.11 fclose(zidu); fclose(zidv);
472     fclose(um_id); fclose(vm_id);
473 edhill 1.1 llzvpbp = acc ./ num;
474 edhill 1.8 % surf(flipud(llzvpbp')), view(2), shading interp, caxis([-.1 .1])
475 edhill 1.4 % ! rm -f primes_92_04/za_llvpbp.mat
476     % save primes_92_04/za_llvpbp.mat llzvpbp
477 edhill 1.1
478     % zonally average ll_vpbp_dbdz
479 edhill 1.9 clear acc num
480 edhill 1.11 acc = zeros(nlatm1, nrm1); acc = NaN;
481 edhill 1.1 num = zeros(size(acc));
482     zid = fopen('primes_92_04/vpbp_dbdz', '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     vpbpdbdz = readcubelev(zid,nslab,ne);
488     % surf(vpbpdbdz(:,:,1)), view(2), shading interp, caxis([-50 50])
489     for jj = 1:nlatm1
490     eval( sprintf('clear inds; inds = inds%04d;',jj) );
491     tmp = vpbpdbdz(inds);
492 edhill 1.9 nzinds = find(isfinite(tmp));
493 edhill 1.1 num(jj,iz) = length(nzinds);
494     acc(jj,iz) = sum(tmp(nzinds));
495     end
496     end
497     fclose(zid);
498     za_ll_vpbp_dbdz = acc ./ num;
499 edhill 1.4 % ! rm -f primes_92_04/za_ll_vpbp_dbdz.mat
500 edhill 1.1 % save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz
501    
502     % zonally average stress
503 edhill 1.9 clear acc num
504 edhill 1.11 acc = zeros(nlatm1, nrm1); acc = NaN;
505 edhill 1.1 num = zeros(size(acc));
506     zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
507     iz = 1;
508     for iz = 1:49,
509     disp(sprintf('iz = %d',iz));
510     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
511     stress = readcubelev(zid,nslab,ne);
512 edhill 1.11 stress(find(abs(stress) > 40.0)) = NaN;
513     if do_plots == 1
514     surf(stress(:,:,6)), view(2), shading flat
515     caxis([-2 2]), colorbar
516     pause(2)
517     end
518 edhill 1.1 for jj = 1:nlatm1
519     eval( sprintf('clear inds; inds = inds%04d;',jj) );
520     tmp = stress(inds);
521 edhill 1.9 nzinds = find(isfinite(tmp));
522 edhill 1.1 num(jj,iz) = length(nzinds);
523     acc(jj,iz) = sum(tmp(nzinds));
524     end
525     end
526     fclose(zid);
527 edhill 1.3 za_stress = acc ./ num;
528 edhill 1.11 % surf(za_stress'), view(2), shading flat, colorbar
529 edhill 1.4 % ! rm -f primes_92_04/stress.mat
530 edhill 1.3 % save primes_92_04/stress.mat za_stress
531 edhill 1.1
532 edhill 1.3 % Average the stress over SSH-contours (streamlines)
533 edhill 1.6 clear num acc
534 edhill 1.3 eta_id = fopen('ave_1992-2004/ETAN.ave', 'r', 'ieee-be');
535     comm = [ 'permute(reshape(fread( id ,nslab,''real*4'',0,' ...
536     '''ieee-be''),[ne 6 ne]),[1 3 2]);' ];
537     readslab = inline(comm,'id','nslab','ne');
538     etan = readslab(eta_id,nslab,ne);
539     fclose(eta_id);
540     % surf(YC(:,:,1)'), view(2),shading interp,colorbar
541     % surf(etan(:,:,1)'), view(2),shading interp,colorbar
542     nssh = 50;
543     sshvals = linspace(-2.3,1.2,nssh);
544     i = 2;
545 edhill 1.12 ssh_yc = zeros(1,49);
546 edhill 1.3 for i = 2:nssh
547 edhill 1.12 % ALL VALUES SOUTH OF *** 30S ***
548     inds = find(sshvals(i-1)<etan & etan<sshvals(i) & YC<-30);
549 edhill 1.3 comm = sprintf('sshinds%04d = uint32(inds);',i-1);
550     eval(comm);
551 edhill 1.12 % get the average YC at each contour
552     ssh_yc(i-1) = mean(YC(inds));
553 edhill 1.3 end
554     zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
555 edhill 1.11 acc = zeros((nssh-1),49); acc = NaN;
556     num = zeros(size(acc));
557 edhill 1.3 iz = 1;
558     for iz = 1:49,
559     disp(sprintf('iz = %d',iz));
560     fseek(zid, (iz - 1)*(ne*ne*6)*4, 'bof');
561     stress = readcubelev(zid,nslab,ne);
562 edhill 1.11 stress(find(abs(stress) > 40.0)) = NaN;
563 edhill 1.3 for jj = 1:(nssh-1)
564     eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
565     tmp = stress(inds);
566 edhill 1.9 nzinds = find(isfinite(tmp));
567 edhill 1.3 num(jj,iz) = length(nzinds);
568     acc(jj,iz) = sum(tmp(nzinds));
569     end
570     end
571     fclose(zid);
572     ssha_stress = acc ./ num;
573 edhill 1.12 % surf(flipud(ssha_stress')), view(2),shading interp,colorbar
574 edhill 1.4 % ! rm -f primes_92_04/ssha_stress.mat
575 edhill 1.12 % save primes_92_04/ssha_stress.mat ssha_stress sshvals ssh_yc
576 edhill 1.2
577 edhill 1.11 % Average the b,ull,vll, and vpbpll over SSH-contours (streamlines)
578     clear n_u n_v n_b n_vb a_u a_v a_b a_vb
579 edhill 1.12 clear n_t n_s a_t a_s
580 edhill 1.11 a_u = zeros((nssh-1),49); a_u = NaN;
581     a_v = zeros((nssh-1),49); a_v = NaN;
582     a_b = zeros((nssh-1),49); a_b = NaN;
583     a_vb = zeros((nssh-1),49); a_vb = NaN;
584 edhill 1.12 a_t = zeros((nssh-1),49); a_t = NaN;
585     a_s = zeros((nssh-1),49); a_s = NaN;
586 edhill 1.11 zid_t = fopen('ave_1992-2004/THETA.ave', 'r', 'ieee-be');
587 edhill 1.12 zid_s = fopen('ave_1992-2004/SALT.ave', 'r', 'ieee-be');
588 edhill 1.10 zid_u = fopen('ave_1992-2004/UVELMASS.ave', 'r', 'ieee-be');
589 edhill 1.11 zid_v = fopen('ave_1992-2004/VVELMASS.ave', 'r', 'ieee-be');
590 edhill 1.10 zid_b = fopen('ave_1992-2004/RHOAnoma.ave', 'r', 'ieee-be');
591     zid_ub = fopen('primes_92_04/upbp', 'r', 'ieee-be');
592     zid_vb = fopen('primes_92_04/vpbp', 'r', 'ieee-be');
593     iz = 1;
594     for iz = 1:49,
595     disp(sprintf('iz = %d',iz));
596 edhill 1.11 fseek(zid_t, (iz - 1)*(ne*ne*6)*4, 'bof');
597 edhill 1.12 fseek(zid_s, (iz - 1)*(ne*ne*6)*4, 'bof');
598 edhill 1.10 fseek(zid_u, (iz - 1)*(ne*ne*6)*4, 'bof');
599     fseek(zid_v, (iz - 1)*(ne*ne*6)*4, 'bof');
600 edhill 1.11 fseek(zid_b, (iz - 1)*(ne*ne*6)*4, 'bof');
601 edhill 1.10 fseek(zid_ub, (iz - 1)*(ne*ne*6)*4, 'bof');
602     fseek(zid_vb, (iz - 1)*(ne*ne*6)*4, 'bof');
603 edhill 1.11 t = readslab(zid_t, nslab,ne);
604 edhill 1.12 s = readslab(zid_s, nslab,ne);
605 edhill 1.10 u = readslab(zid_u, nslab,ne);
606 edhill 1.11 v = readslab(zid_v, nslab,ne);
607 edhill 1.10 b = readslab(zid_b, nslab,ne);
608 edhill 1.11 maski = find(t == 0.0);
609 edhill 1.10 upbp = readcubelev(zid_ub, nslab,ne);
610     vpbp = readcubelev(zid_vb, nslab,ne);
611     % llupbp = upbp .* llux + vpbp .* llvx;
612     llvpbp = upbp .* lluy + vpbp .* llvy;
613 edhill 1.11 llu = u .* llux + v .* llvx;
614     llv = u .* lluy + v .* llvy;
615 edhill 1.12 t(maski) = NaN; s(maski) = NaN; b(maski) = NaN;
616     llu(maski) = NaN; llv(maski) = NaN; llvpbp(maski) = NaN;
617 edhill 1.11 if do_plots == 1
618     figure(1), subplot(1,1,1)
619     subplot(1,2,1), surf(llu(:,:,1)), view(2), shading flat, colorbar
620     subplot(1,2,2), surf(llv(:,:,1)), view(2), shading flat, colorbar
621     figure(2), subplot(1,1,1)
622     subplot(1,2,1), surf(b(:,:,1)), view(2), shading flat, colorbar
623     subplot(1,2,2), surf(llvpbp(:,:,1)), view(2), shading flat, colorbar
624     pause(2)
625     end
626 edhill 1.10 for jj = 1:(nssh-1)
627     eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
628 edhill 1.12 t_t = t(inds);
629     t_s = s(inds);
630 edhill 1.11 t_u = llu(inds);
631     t_v = llv(inds);
632 edhill 1.10 t_b = b(inds);
633     t_vb = llvpbp(inds);
634 edhill 1.12 i_t = find(isfinite( t_t ));
635     i_s = find(isfinite( t_s ));
636 edhill 1.10 i_u = find(isfinite( t_u ));
637 edhill 1.11 i_v = find(isfinite( t_v ));
638 edhill 1.10 i_b = find(isfinite( t_b ));
639     i_vb = find(isfinite( t_vb ));
640 edhill 1.12 n_t(jj,iz) = length( i_t );
641     n_s(jj,iz) = length( i_s );
642 edhill 1.10 n_u(jj,iz) = length( i_u );
643 edhill 1.11 n_v(jj,iz) = length( i_v );
644 edhill 1.10 n_b(jj,iz) = length( i_b );
645     n_vb(jj,iz) = length( i_vb );
646 edhill 1.12 a_t(jj,iz) = sum( t_t( i_t ));
647     a_s(jj,iz) = sum( t_s( i_s ));
648 edhill 1.10 a_u(jj,iz) = sum( t_u( i_u ));
649 edhill 1.11 a_v(jj,iz) = sum( t_v( i_v ));
650 edhill 1.10 a_b(jj,iz) = sum( t_b( i_b ));
651     a_vb(jj,iz) = sum( t_vb( i_vb ));
652     end
653     end
654 edhill 1.12 fclose(zid_t); fclose(zid_s);
655 edhill 1.11 fclose(zid_u); fclose(zid_b);
656     fclose(zid_ub); fclose(zid_vb);
657 edhill 1.12 ssha_t = a_t ./ n_t;
658     ssha_s = a_s ./ n_s;
659 edhill 1.11 ssha_llu = a_u ./ n_u;
660     ssha_llv = a_v ./ n_v;
661     ssha_b = a_b ./ n_b;
662     ssha_llvpbp = a_vb ./ n_vb;
663     % ! rm -f primes_92_04/ssha_uv_b_vpbp.mat
664 edhill 1.12 % save primes_92_04/ssha_uv_b_vpbp.mat ssha_t ssha_s ssha_llu ssha_llv ssha_b ssha_llvpbp
665    
666     % Average the dtdy and vptpll over SSH-contours (streamlines)
667     clear n_vt a_vt n_t2 a_t2
668     a_vt = zeros((nssh-1),49); a_vt = NaN;
669     a_t2 = zeros((nssh-1),49); a_t2 = NaN;
670     zid_t = fopen('ave_1992-2004/THETA.ave', 'r', 'ieee-be');
671     zid_ut = fopen('primes_92_04/uptp', 'r', 'ieee-be');
672     zid_vt = fopen('primes_92_04/vptp', 'r', 'ieee-be');
673     zid_t2 = fopen('primes_92_04/tp2', 'r', 'ieee-be');
674     iz = 1;
675     for iz = 1:49,
676     disp(sprintf('iz = %d',iz));
677     fseek(zid_t, (iz - 1)*(ne*ne*6)*4, 'bof');
678     fseek(zid_ut, (iz - 1)*(ne*ne*6)*4, 'bof');
679     fseek(zid_vt, (iz - 1)*(ne*ne*6)*4, 'bof');
680     fseek(zid_t2, (iz - 1)*(ne*ne*6)*4, 'bof');
681     t = readslab(zid_t, nslab,ne);
682     t2 = readslab(zid_t2, nslab,ne);
683     maski = find(t == 0.0);
684     uptp = readcubelev(zid_ut, nslab,ne);
685     vptp = readcubelev(zid_vt, nslab,ne);
686     t2 = readcubelev(zid_t2, nslab,ne);
687     % llupbp = upbp .* llux + vpbp .* llvx;
688     llvptp = uptp .* lluy + vptp .* llvy;
689     llvptp(maski) = NaN;
690     t2(maski) = NaN;
691     for jj = 1:(nssh-1)
692     eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
693     t_vt = llvptp(inds);
694     t_t2 = t2(inds);
695     i_vt = find(isfinite( t_vt ));
696     i_t2 = find(isfinite( t_t2 ));
697     n_vt(jj,iz) = length( i_vt );
698     n_t2(jj,iz) = length( i_t2 );
699     a_vt(jj,iz) = sum( t_vt( i_vt ));
700     a_t2(jj,iz) = sum( t_t2( i_t2 ));
701     end
702     end
703     fclose(zid_t); fclose(zid_t2);
704     fclose(zid_ut); fclose(zid_vt);
705     ssha_llvptp = a_vt ./ n_vt;
706     ssha_tp2 = a_t2 ./ n_t2;
707     % surf(flipud(ssha_llvptp')), view(2), shading flat, colorbar
708     % ! rm -f primes_92_04/ssha_vptp.mat
709     % save primes_92_04/ssha_vptp.mat ssha_llvptp ssha_tp2
710 edhill 1.10
711 edhill 1.1
712     clear all
713     close all
714    
715 edhill 1.11 do_plots = 0;
716    
717 edhill 1.1 %==================================================================
718     % Read the tile00?.mitgrid files
719     gvars = { 'XC','YC','DXF','DYF','RA','XG','YG','DXV', ...
720     'DYU','RAZ','DXC','DYC','RAW','RAS','DXG','DYG' };
721    
722     ne = 510;
723     nep1 = ne + 1;
724     iface = 1;
725     for iface = 1:6
726     fname = sprintf('grid/tile%03d.mitgrid', iface);
727     gid = fopen(fname, 'r', 'ieee-be');
728     tmp = reshape(fread(gid,inf,'real*8',0,'ieee-be'),[nep1,nep1,16]);
729     fclose(gid);
730     % surf(tmp(:,:,1)), view(2), shading interp
731     % for jj = 1:length(gvars)
732     for jj = 1:7
733     comm = sprintf('%s(:,:,%d) = tmp(:,:,%d);', ...
734     [gvars{jj}], iface, jj);
735     eval(comm);
736     end
737     end
738     % surf(XC(:,:,1)), view(2), shading interp
739     % subplot(2,1,1), a = [1:10]; surf(XC(a,a,1)), view(2)
740     % subplot(2,1,2), a = [(nep1-10):nep1]; surf(XC(a,a,1)), view(2)
741     % surf(YC(:,:,1)), view(2), shading interp
742     % surf(XG(:,:,1)), view(2), shading interp
743     % surf(YG(:,:,1)), view(2), shading interp
744     is = [1:ne];
745     vs = { 'XC','YC','DXF','DYF','RA' };
746     for i = 1:length(vs)
747     eval(sprintf('%s = %s(is,is,:);',vs{i},vs{i}));
748     end
749    
750     delR = [ ...
751     10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.00, 10.01, ...
752     10.03, 10.11, 10.32, 10.80, 11.76, 13.42, 16.04 , 19.82, 24.85, ...
753     31.10, 38.42, 46.50, 55.00, 63.50, 71.58, 78.90, 85.15, 90.18, ...
754     93.96, 96.58, 98.25, 99.25,100.01,101.33,104.56,111.33,122.83, ...
755     139.09,158.94,180.83,203.55,226.50,249.50,272.50,295.50,318.50, ...
756     341.50,364.50,387.50,410.50,433.50,456.50 ];
757     R = cumsum(delR) - 0.5*delR;
758 edhill 1.4 Ri = R(1:(length(R)-1)) ...
759     + 0.25*delR(1:(length(R)-1)) + 0.25*delR(2:(length(R)));
760    
761 edhill 1.1 n1 = ne - 1;
762     dux = zeros(size(XC));
763     duy = zeros(size(XC));
764     dvx = zeros(size(XC));
765     dvy = zeros(size(XC));
766     dux(:,:,:) = diff(XG(:,1:ne,:),1,1);
767     dvx(:,:,:) = diff(XG(1:ne,:,:),1,2);
768     duy(:,:,:) = diff(YG(:,1:ne,:),1,1);
769     dvy(:,:,:) = diff(YG(1:ne,:,:),1,2);
770     dux = dux + 360*double(dux < 180);
771     dux = dux - 360*double(dux > 180); % [ min(min(dux)) max(max(dux)) ]
772     duy = duy + 360*double(duy < 180);
773     duy = duy - 360*double(duy > 180); % [ min(min(duy)) max(max(duy)) ]
774     dvx = dvx + 360*double(dvx < 180);
775     dvx = dvx - 360*double(dvx > 180); % [ min(min(dvx)) max(max(dvx)) ]
776     dvy = dvy + 360*double(dvy < 180);
777     dvy = dvy - 360*double(dvy > 180); % [ min(min(dvy)) max(max(dvy)) ]
778     llux = dux ./ sqrt(dux.^2 + duy.^2);
779     lluy = duy ./ sqrt(dux.^2 + duy.^2);
780     llvx = dvx ./ sqrt(dvx.^2 + dvy.^2);
781     llvy = dvy ./ sqrt(dvx.^2 + dvy.^2);
782    
783     %==================================================================
784     % Project fields to lower-res 1-degree Lat-Lon and write
785     % as NetCDF for viewing with Ingrid
786     %
787     ne = 510;
788     nf = 6;
789     nz = 50;
790     nslab = ne*ne*nf;
791     adir = 'primes_92_04';
792     lat = [-90:90];
793     lon = [0:360];
794 edhill 1.2 ir = [ 1 2 3 5 10 15 20 25 30 35 40 ];
795 edhill 1.1 xc360 = XC + 180;
796    
797 edhill 1.2 % ! rm -f cube_22_primes_at1deg.nc
798     nc = netcdf(['cube_22_primes_at1deg.nc'], 'clobber');
799     nc.reference = [ 'The cube_22 primes from Dimitris Menemenlis' ...
800     ' regridded to 1-deg Lat-Lon' ];
801 edhill 1.1 nc.author = 'Ed Hill <eh3@mit.edu>';
802 edhill 1.2 nc.date = 'March 27, 2005';
803 edhill 1.1 nc('X') = length(lon);
804     nc('Y') = length(lat);
805     nc('Z') = length(ir);
806 edhill 1.4 nc('Zc') = length(R);
807     nc('Zi') = length(Ri);
808 edhill 1.1 nc{'X'} = 'X';
809     nc{'Y'} = 'Y';
810     nc{'Z'} = 'Z';
811 edhill 1.4 nc{'Zc'} = 'Zc';
812     nc{'Zi'} = 'Zi';
813 edhill 1.1 nc{'X'}.uniquename = 'X';
814     nc{'X'}.long_name = 'longitude';
815     nc{'X'}.gridtype = ncint(1);
816     nc{'X'}.units = 'degree_east';
817     nc{'Y'}.uniquename = 'Y';
818     nc{'Y'}.long_name = 'latitude';
819     nc{'Y'}.gridtype = ncint(0);
820     nc{'Y'}.units = 'degree_north';
821     nc{'Z'}.uniquename = 'Z';
822     nc{'Z'}.long_name = 'depth';
823     nc{'Z'}.gridtype = ncint(0);
824     nc{'Z'}.units = 'm';
825 edhill 1.4 nc{'Zc'}.uniquename = 'Zc';
826     nc{'Zc'}.long_name = 'depth_at_center';
827     nc{'Zc'}.gridtype = ncint(0);
828     nc{'Zc'}.units = 'm';
829     nc{'Zi'}.uniquename = 'Zi';
830     nc{'Zi'}.long_name = 'depth_at_interface';
831     nc{'Zi'}.gridtype = ncint(0);
832     nc{'Zi'}.units = 'm';
833 edhill 1.13 nc{'X'}(:) = lon - 180;
834 edhill 1.1 nc{'Y'}(:) = lat;
835     nc{'Z'}(:) = R(ir);
836 edhill 1.4 nc{'Zc'}(:) = R;
837     nc{'Zi'}(:) = Ri;
838    
839 edhill 1.10 f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'}, ...
840     {'stress'}, {'dbdy'}, {'K'} };
841 edhill 1.1
842     ifg = 1;
843     for ifg = 1:length(f_s_3d)
844     acell = f_s_3d{ifg};
845     tname = acell{1};
846     disp([ ' ' tname ' :' ]);
847     fname = sprintf('%s/%s',adir,tname);
848     fid = fopen(fname,'r','ieee-be');
849     id = tname;
850     nc{ id } = { 'Z' 'Y' 'X' };
851     nc{ id }.missing_value = ncdouble(NaN);
852     nc{ id }.FillValue_ = ncdouble(0.0);
853     ii = 1;
854     for ii = 1:length(ir)
855     iz = ir(ii);
856     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
857    
858     fseek(fid,nslab*4*(iz-1),'bof');
859     tmp = fread(fid,nslab,'real*4',0,'ieee-be');
860     tr = reshape(tmp,[ 510 510 6 ]);
861     % surf(tr(:,:,1)), view(2), shading interp
862     trn = tr;
863     trn(find(tr == 0.0)) = NaN;
864     clear tmp tr
865     % v = sdac_regrid(xc360,YC,trn,lonm,latm);
866     v = ll_regrid(xc360,YC,trn,lon,lat);
867     % surf(lon,lat,v'), caxis([25 40]), view(2), shading interp, colorbar
868     nc{ id }(ii,:,:) = permute(v,[2 1]);
869     end
870     fclose(fid);
871     end
872    
873     id = 'sum_up2_vp2';
874     nc{ id } = { 'Z' 'Y' 'X' };
875     nc{ id }.missing_value = ncdouble(NaN);
876     nc{ id }.FillValue_ = ncdouble(0.0);
877     fidu = fopen(sprintf('%s/%s',adir,'up2'),'r','ieee-be');
878     fidv = fopen(sprintf('%s/%s',adir,'vp2'),'r','ieee-be');
879     for ii = 1:length(ir)
880     iz = ir(ii);
881     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
882     fseek(fidu,nslab*4*(iz-1),'bof');
883     fseek(fidv,nslab*4*(iz-1),'bof');
884     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
885     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
886     trnu = tru; trnu(find(tru == 0.0)) = NaN;
887     trnv = trv; trnv(find(trv == 0.0)) = NaN;
888     clear tmp tru trv
889     lluv = ll_regrid(xc360,YC,trnu+trnv,lon,lat);
890     nc{ id }(ii,:,:) = permute(lluv,[2 1]);
891     end
892     fclose(fidu);
893     fclose(fidv);
894    
895     f_v_3d = { {'up2','vp2'}, ...
896     {'uptp','vptp'}, {'upsp','vpsp'}, {'upbp','vpbp'} };
897     for ip = 1:length(f_v_3d)
898     cell = f_v_3d{ip};
899     idu = cell{1};
900     idv = cell{2};
901     disp([' ' idu ' ' idv]);
902     nc{ idu } = { 'Z' 'Y' 'X' };
903     nc{ idu }.missing_value = ncdouble(NaN);
904     nc{ idu }.FillValue_ = ncdouble(0.0);
905     nc{ idv } = { 'Z' 'Y' 'X' };
906     nc{ idv }.missing_value = ncdouble(NaN);
907     nc{ idv }.FillValue_ = ncdouble(0.0);
908     fidu = fopen(sprintf('%s/%s',adir,idu),'r','ieee-be');
909     fidv = fopen(sprintf('%s/%s',adir,idv),'r','ieee-be');
910     for ii = 1:length(ir)
911     iz = ir(ii);
912     disp(sprintf(' iz = %3d R(iz) = %g',iz,R(iz)));
913     fseek(fidu,nslab*4*(iz-1),'bof');
914     fseek(fidv,nslab*4*(iz-1),'bof');
915     tru = reshape(fread(fidu,nslab,'real*4',0,'ieee-be'),[510 510 6]);
916     trv = reshape(fread(fidv,nslab,'real*4',0,'ieee-be'),[510 510 6]);
917     trnu = tru; trnu(find(tru == 0.0)) = NaN;
918     trnv = trv; trnv(find(trv == 0.0)) = NaN;
919     clear tmp tru trv
920     llru = trnu .* llux + trnv .* llvx;
921     llrv = trnu .* lluy + trnv .* llvy;
922     llu = ll_regrid(xc360,YC,llru,lon,lat);
923     llv = ll_regrid(xc360,YC,llrv,lon,lat);
924     nc{ idu }(ii,:,:) = permute(llu,[2 1]);
925     nc{ idv }(ii,:,:) = permute(llv,[2 1]);
926     end
927     end
928     fclose(fidu);
929     fclose(fidv);
930     nc = close(nc);
931    
932 edhill 1.4 % === zonal and stream-wise averages ===
933     % save primes_92_04/za_llvpbp.mat llzvpbp
934     % save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz
935 edhill 1.11 % save primes_9b2_04/stress.mat za_stress
936 edhill 1.4 % save primes_92_04/ssha_stress.mat ssha_stress sshvals
937     load primes_92_04/za_llvpbp.mat
938     load primes_92_04/za_ll_vpbp_dbdz.mat
939     load primes_92_04/stress.mat
940     load primes_92_04/ssha_stress.mat
941 edhill 1.11 load primes_92_04/ssha_uv_b_vpbp.mat
942 edhill 1.12 load primes_92_04/ssha_vptp.mat
943    
944 edhill 1.5 ssh = sshvals(1:(length(sshvals)-1)) + 0.5*diff(sshvals);
945 edhill 1.11 if do_plots == 1
946     surf(ssh,Ri,ssha_b'), view(2), shading flat, colorbar
947     surf(ssh,-Ri,ssha_llu'), view(2), shading flat, colorbar
948     end
949 edhill 1.1
950 edhill 1.4 lat = [-90:90];
951     latm = lat(1:(length(lat)-1)) + 0.5*diff(lat);
952 edhill 1.1
953 edhill 1.13 % fit = [ 0 -70 ; 5 -66 ; 15 -55 ; ...
954     % 25 -48 ; 32 -42 ; 34 -35 ; 40 -26.5 ; 50 -22 ];
955     fit = [ 0 -70 ; 5 -66 ; 15 -55 ; ...
956     25 -48 ; 32 -45 ; 34 -38 ; 40 -35 ; 50 -30 ];
957     ssh_yc_sm = interp1(fit(:,1),fit(:,2),[1:length(ssh_yc)]);
958     if do_plots == 1
959     plot(ssh_yc)
960     hold on, plot(fit(:,1),fit(:,2),'ro-'), hold off
961     hold on, plot(ssh_yc_sm,'go-'), hold off
962     grid on
963     end
964    
965     % Compute dT/dy and dB/dy from stream-wise averages
966     nsshi = size(ssha_t,1);
967     ssha_dtdy = zeros(size(ssha_t));
968     ssha_dtdy(:) = NaN;
969     ssha_dtdy1 = ssha_dtdy; ssha_dtdy2 = ssha_dtdy;
970     ssha_dbdy = ssha_dtdy;
971     ssha_dbdy1 = ssha_dtdy; ssha_dbdy2 = ssha_dtdy;
972     rad = 6.37*10^6;
973     for iz = 1:49,
974     % dy = rad * pi * abs(ssh_yc_sm(jj)-ssh_yc_sm(jj-1))/180;
975     dy = 0.8 * rad * pi / 180;
976     for jj = 1:(nsshi-1)
977     ssha_dtdy1(jj,iz) = (ssha_t(jj+1,iz) - ssha_t(jj,iz))/dy;
978     ssha_dbdy1(jj,iz) = (ssha_b(jj+1,iz) - ssha_b(jj,iz))/dy;
979     end
980     for jj = 2:nsshi
981     ssha_dtdy2(jj,iz) = (ssha_t(jj,iz) - ssha_t(jj-1,iz))/dy;
982     ssha_dbdy2(jj,iz) = (ssha_b(jj,iz) - ssha_b(jj-1,iz))/dy;
983     end
984     end
985     ssha_dtdy = (ssha_dtdy1 + ssha_dtdy2)/2;
986     ssha_dbdy = (ssha_dbdy1 + ssha_dbdy2)/2;
987     if do_plots == 1
988     surf(flipud(ssha_dtdy')), view(2), shading flat, colorbar
989     surf(flipud(ssha_dbdy')), view(2), shading flat, colorbar
990     end
991    
992     % Compute K = (v't')/(dt/dy)
993     ssha_K_t = - ssha_llvptp ./ ssha_dtdy;
994     ssha_K_b = - ssha_llvpbp ./ ssha_dbdy;
995     if do_plots == 1
996     surf(flipud(ssha_K_t')), view(2), shading flat, colorbar
997     surf(flipud(ssha_K_b')), view(2), shading flat, colorbar
998     end
999    
1000    
1001 edhill 1.4 % ! rm -f cube_22_zsa.nc
1002     nc = netcdf(['cube_22_zsa.nc'], 'clobber');
1003     nc.reference = [ 'The cube_22 zonal and stream-wise averages.' ];
1004     nc.author = 'Ed Hill <eh3@mit.edu>';
1005     nc.date = 'March 27, 2005';
1006     nc('Y') = length(latm);
1007     nc('Zc') = length(R);
1008     nc('Zi') = length(Ri);
1009 edhill 1.5 nc('SSH') = length(ssh);
1010 edhill 1.4 nc{'Y'} = 'Y';
1011     nc{'Zc'} = 'Zc';
1012     nc{'Zi'} = 'Zi';
1013 edhill 1.5 nc{'SSH'} = 'SSH';
1014 edhill 1.4 nc{'Y'}.uniquename = 'Y';
1015     nc{'Y'}.long_name = 'latitude';
1016     nc{'Y'}.gridtype = ncint(0);
1017     nc{'Y'}.units = 'degree_north';
1018     nc{'Zc'}.uniquename = 'Zc';
1019     nc{'Zc'}.long_name = 'depth_at_center';
1020     nc{'Zc'}.gridtype = ncint(0);
1021     nc{'Zc'}.units = 'm';
1022     nc{'Zi'}.uniquename = 'Zi';
1023     nc{'Zi'}.long_name = 'depth_at_interface';
1024     nc{'Zi'}.gridtype = ncint(0);
1025     nc{'Zi'}.units = 'm';
1026 edhill 1.5 nc{'SSH'}.uniquename = 'SSH';
1027     nc{'SSH'}.long_name = 'sea_surface_height';
1028     nc{'SSH'}.gridtype = ncint(0);
1029     nc{'SSH'}.units = 'm';
1030 edhill 1.4 nc{'Y'}(:) = latm;
1031     nc{'Zc'}(:) = R;
1032     nc{'Zi'}(:) = Ri;
1033 edhill 1.5 nc{'SSH'}(:) = ssh;
1034 edhill 1.1
1035 edhill 1.4 id = 'llzvpbp';
1036     nc{ id } = { 'Zc' 'Y' };
1037     nc{ id }.missing_value = ncdouble(NaN);
1038     nc{ id }.FillValue_ = ncdouble(0.0);
1039     nc{ id }(:) = permute(llzvpbp,[2 1]);
1040 edhill 1.1
1041 edhill 1.5 f_i = { {'za_ll_vpbp_dbdz'}, {'za_stress'} };
1042 edhill 1.4 for ip = 1:length(f_i)
1043     cell = f_i{ip};
1044     id = cell{1};
1045     disp([' ' id]);
1046     nc{ id } = { 'Zi' 'Y' };
1047     nc{ id }.missing_value = ncdouble(NaN);
1048     nc{ id }.FillValue_ = ncdouble(0.0);
1049     eval(sprintf('tmp = %s;', id));
1050     nc{ id }(:) = permute(tmp,[2 1]);
1051 edhill 1.1 end
1052 edhill 1.5
1053 edhill 1.11 f_i = { {'ssha_stress'}, ...
1054 edhill 1.13 {'ssha_llvptp'}, {'ssha_tp2'}, {'ssha_dtdy'}, ...
1055     {'ssha_K_t'}, {'ssha_K_b'}, ...
1056 edhill 1.12 {'ssha_t'}, {'ssha_s'}, ...
1057 edhill 1.11 {'ssha_llu'}, {'ssha_llv'}, {'ssha_b'}, {'ssha_llvpbp'} };
1058     for ip = 1:length(f_i)
1059     cell = f_i{ip};
1060     id = cell{1};
1061     disp([' ' id]);
1062     nc{ id } = { 'Zi' 'SSH' };
1063     nc{ id }.missing_value = ncdouble(NaN);
1064     nc{ id }.FillValue_ = ncdouble(0.0);
1065     eval(sprintf('nc{ id }(:) = permute(%s,[2 1]);',id));
1066     end
1067 edhill 1.4 nc = close(nc);
1068 edhill 1.12
1069    
1070     % ! rm -f cube_22_zsa_elat.nc
1071     nc = netcdf(['cube_22_zsa_elat.nc'], 'clobber');
1072     nc.reference = [ 'The cube_22 zonal and stream-wise averages.' ];
1073     nc.author = 'Ed Hill <eh3@mit.edu>';
1074     nc.date = 'March 27, 2005';
1075     nc('Y') = length(latm);
1076     nc('Zc') = length(R);
1077     nc('Zi') = length(Ri);
1078     nc('elat') = length(ssh_yc);
1079     nc{'Y'} = 'Y';
1080     nc{'Zc'} = 'Zc';
1081     nc{'Zi'} = 'Zi';
1082     nc{'elat'} = 'elat';
1083     nc{'Y'}.uniquename = 'Y';
1084     nc{'Y'}.long_name = 'latitude';
1085     nc{'Y'}.gridtype = ncint(0);
1086     nc{'Y'}.units = 'degree_north';
1087     nc{'Zc'}.uniquename = 'Zc';
1088     nc{'Zc'}.long_name = 'depth_at_center';
1089     nc{'Zc'}.gridtype = ncint(0);
1090     nc{'Zc'}.units = 'm';
1091     nc{'Zi'}.uniquename = 'Zi';
1092     nc{'Zi'}.long_name = 'depth_at_interface';
1093     nc{'Zi'}.gridtype = ncint(0);
1094     nc{'Zi'}.units = 'm';
1095     nc{'elat'}.uniquename = 'elat';
1096     nc{'elat'}.long_name = 'equivalent_latitude';
1097     nc{'elat'}.gridtype = ncint(0);
1098     nc{'elat'}.units = 'degree_north';
1099     nc{'Y'}(:) = latm;
1100     nc{'Zc'}(:) = R;
1101     nc{'Zi'}(:) = Ri;
1102     nc{'elat'}(:) = ssh_yc_sm;
1103    
1104     id = 'llzvpbp';
1105     nc{ id } = { 'Zc' 'Y' };
1106     nc{ id }.missing_value = ncdouble(NaN);
1107     nc{ id }.FillValue_ = ncdouble(0.0);
1108     nc{ id }(:) = permute(llzvpbp,[2 1]);
1109    
1110     f_i = { {'za_ll_vpbp_dbdz'}, {'za_stress'} };
1111     for ip = 1:length(f_i)
1112     cell = f_i{ip};
1113     id = cell{1};
1114     disp([' ' id]);
1115     nc{ id } = { 'Zi' 'Y' };
1116     nc{ id }.missing_value = ncdouble(NaN);
1117     nc{ id }.FillValue_ = ncdouble(0.0);
1118     eval(sprintf('tmp = %s;', id));
1119     nc{ id }(:) = permute(tmp,[2 1]);
1120     end
1121    
1122     f_i = { {'ssha_stress'}, ...
1123 edhill 1.13 {'ssha_llvptp'}, {'ssha_tp2'}, {'ssha_dtdy'}, ...
1124     {'ssha_K_t'}, {'ssha_K_b'}, ...
1125 edhill 1.12 {'ssha_t'}, {'ssha_s'}, ...
1126     {'ssha_llu'}, {'ssha_llv'}, {'ssha_b'}, {'ssha_llvpbp'} };
1127     for ip = 1:length(f_i)
1128     cell = f_i{ip};
1129     id = cell{1};
1130     disp([' ' id]);
1131     nc{ id } = { 'Zi' 'elat' };
1132     nc{ id }.missing_value = ncdouble(NaN);
1133     nc{ id }.FillValue_ = ncdouble(0.0);
1134     eval(sprintf('nc{ id }(:) = permute(%s,[2 1]);',id));
1135     end
1136    
1137     nc = close(nc);
1138    
1139 edhill 1.6 % surf(,,log(ssha_stress')), view(2),shading interp,colorbar
1140    
1141 edhill 1.13 % ! scp cube_22_primes_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1142 edhill 1.5 % ! scp cube_22_zsa.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1143 edhill 1.12 % ! scp cube_22_zsa_elat.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1144 edhill 1.4 % ! mv cube_22_primes_at1deg.nc primes_92_04
1145 edhill 1.1

  ViewVC Help
Powered by ViewVC 1.1.22