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

Contents 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.14 - (show annotations) (download)
Sun May 1 01:53:17 2005 UTC (20 years, 3 months ago) by edhill
Branch: MAIN
CVS Tags: HEAD
Changes since 1.13: +22 -2 lines
 o add stress from sw-averaged-first dbdz

1 %=======================================================
2 %
3 % $Id: plot_c22.m,v 1.13 2005/04/27 22:00:14 edhill Exp $
4 %
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 % cd /r/r0/edhill/eddy_stats/c22
15
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 - 180;
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 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 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
235 lpath = 'ave_1992-2004/';
236 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 dbdy__id = fopen([opath 'dbdy'], 'wb', 'ieee-be');
285 K_____id = fopen([opath 'K'], 'wb', 'ieee-be');
286
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 t = readslab(t__id,nslab,ne); t2 = readslab(t2_id,nslab,ne);
303 nan_inds = find(t == 0.0);
304 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 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 % surf(v2(:,:,1)), view(2), shading interp, colorbar
311 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 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
339 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 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 [ tonu, tonv ] = mass_on_u_v(t);
348 [ sonu, sonv ] = mass_on_u_v(s);
349 [ bonu, bonv ] = mass_on_u_v(b);
350 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 % surf(upbp(:,:,1)), shading interp, view(2)
354 % caxis([-0.1 0.1]), colorbar
355
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 % tmp(find(tmp == 0.0)) = 1.0;
363 vpbpdbdz = ave_llvpbp ./ tmp;
364 fwrite(vpbpdzid, vpbpdbdz, 'real*4');
365 % \tau_x = 2\Omega sin(\phi) = 4\pi*sin(\phi)/(24*3600)
366 fac = 1000 * 4*pi/(24*3600);
367 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
373 % determine diffusivities
374 dbdy = calc_dbdy(b, dut,dvt, lluy,llvy);
375 diffus = - vpbp ./ dbdy;
376
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 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 fwrite(wptp__id, wptp, 'real*4');
391 fwrite(wpsp__id, wpsp, 'real*4');
392 fwrite(wpbp__id, wpbp, 'real*4');
393 fwrite(dbdy__id, dbdy, 'real*4');
394 fwrite(K_____id, diffus, 'real*4');
395 end
396
397 clear uptp vptp upsp vpsp ubbp vpbp
398 clear wptp wpsp wpbp dbdy diffus
399
400 % save current_state_20050422
401 % load current_state_20050422
402
403 do_plots = 0;
404
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 % Zonally average the ll_upbp,ll_vpbp
424 clear acc num
425 acc = zeros(nlatm1, nrm1); acc = NaN;
426 num = zeros(size(acc));
427 %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 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 %fseek(um_id, (iz - 1)*(ne*ne*6)*4, 'bof');
436 %fseek(vm_id, (iz - 1)*(ne*ne*6)*4, 'bof');
437 upbp = readcubelev(zidu,nslab,ne);
438 vpbp = readcubelev(zidv,nslab,ne);
439 %um = readslab(um_id,nslab,ne);
440 %vm = readslab(vm_id,nslab,ne);
441 llupbp = upbp .* llux + vpbp .* llvx;
442 llvpbp = upbp .* lluy + vpbp .* llvy;
443 %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 for jj = 1:nlatm1
463 eval( sprintf('clear inds; inds = inds%04d;',jj) );
464 tmp = llvpbp(inds);
465 % nzinds = find(tmp ~= 0.0);
466 finds = find(isfinite(tmp));
467 num(jj,iz) = length(finds);
468 acc(jj,iz) = sum(tmp(finds));
469 end
470 end
471 fclose(zidu); fclose(zidv);
472 fclose(um_id); fclose(vm_id);
473 llzvpbp = acc ./ num;
474 % surf(flipud(llzvpbp')), view(2), shading interp, caxis([-.1 .1])
475 % ! rm -f primes_92_04/za_llvpbp.mat
476 % save primes_92_04/za_llvpbp.mat llzvpbp
477
478 % zonally average ll_vpbp_dbdz
479 clear acc num
480 acc = zeros(nlatm1, nrm1); acc = NaN;
481 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 nzinds = find(isfinite(tmp));
493 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 % ! rm -f primes_92_04/za_ll_vpbp_dbdz.mat
500 % save primes_92_04/za_ll_vpbp_dbdz.mat za_ll_vpbp_dbdz
501
502 % zonally average stress
503 clear acc num
504 acc = zeros(nlatm1, nrm1); acc = NaN;
505 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 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 for jj = 1:nlatm1
519 eval( sprintf('clear inds; inds = inds%04d;',jj) );
520 tmp = stress(inds);
521 nzinds = find(isfinite(tmp));
522 num(jj,iz) = length(nzinds);
523 acc(jj,iz) = sum(tmp(nzinds));
524 end
525 end
526 fclose(zid);
527 za_stress = acc ./ num;
528 % surf(za_stress'), view(2), shading flat, colorbar
529 % ! rm -f primes_92_04/stress.mat
530 % save primes_92_04/stress.mat za_stress
531
532 % Average the stress over SSH-contours (streamlines)
533 clear num acc
534 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 ssh_yc = zeros(1,49);
546 for i = 2:nssh
547 % ALL VALUES SOUTH OF *** 30S ***
548 inds = find(sshvals(i-1)<etan & etan<sshvals(i) & YC<-30);
549 comm = sprintf('sshinds%04d = uint32(inds);',i-1);
550 eval(comm);
551 % get the average YC at each contour
552 ssh_yc(i-1) = mean(YC(inds));
553 end
554 zid = fopen('primes_92_04/stress', 'r', 'ieee-be');
555 acc = zeros((nssh-1),49); acc = NaN;
556 num = zeros(size(acc));
557 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 stress(find(abs(stress) > 40.0)) = NaN;
563 for jj = 1:(nssh-1)
564 eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
565 tmp = stress(inds);
566 nzinds = find(isfinite(tmp));
567 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 % surf(flipud(ssha_stress')), view(2),shading interp,colorbar
574 % ! rm -f primes_92_04/ssha_stress.mat
575 % save primes_92_04/ssha_stress.mat ssha_stress sshvals ssh_yc
576
577 % 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 clear n_t n_s a_t a_s
580 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 a_t = zeros((nssh-1),49); a_t = NaN;
585 a_s = zeros((nssh-1),49); a_s = NaN;
586 zid_t = fopen('ave_1992-2004/THETA.ave', 'r', 'ieee-be');
587 zid_s = fopen('ave_1992-2004/SALT.ave', 'r', 'ieee-be');
588 zid_u = fopen('ave_1992-2004/UVELMASS.ave', 'r', 'ieee-be');
589 zid_v = fopen('ave_1992-2004/VVELMASS.ave', 'r', 'ieee-be');
590 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 fseek(zid_t, (iz - 1)*(ne*ne*6)*4, 'bof');
597 fseek(zid_s, (iz - 1)*(ne*ne*6)*4, 'bof');
598 fseek(zid_u, (iz - 1)*(ne*ne*6)*4, 'bof');
599 fseek(zid_v, (iz - 1)*(ne*ne*6)*4, 'bof');
600 fseek(zid_b, (iz - 1)*(ne*ne*6)*4, 'bof');
601 fseek(zid_ub, (iz - 1)*(ne*ne*6)*4, 'bof');
602 fseek(zid_vb, (iz - 1)*(ne*ne*6)*4, 'bof');
603 t = readslab(zid_t, nslab,ne);
604 s = readslab(zid_s, nslab,ne);
605 u = readslab(zid_u, nslab,ne);
606 v = readslab(zid_v, nslab,ne);
607 b = readslab(zid_b, nslab,ne);
608 maski = find(t == 0.0);
609 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 llu = u .* llux + v .* llvx;
614 llv = u .* lluy + v .* llvy;
615 t(maski) = NaN; s(maski) = NaN; b(maski) = NaN;
616 llu(maski) = NaN; llv(maski) = NaN; llvpbp(maski) = NaN;
617 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 for jj = 1:(nssh-1)
627 eval( sprintf('clear inds; inds = sshinds%04d;',jj) );
628 t_t = t(inds);
629 t_s = s(inds);
630 t_u = llu(inds);
631 t_v = llv(inds);
632 t_b = b(inds);
633 t_vb = llvpbp(inds);
634 i_t = find(isfinite( t_t ));
635 i_s = find(isfinite( t_s ));
636 i_u = find(isfinite( t_u ));
637 i_v = find(isfinite( t_v ));
638 i_b = find(isfinite( t_b ));
639 i_vb = find(isfinite( t_vb ));
640 n_t(jj,iz) = length( i_t );
641 n_s(jj,iz) = length( i_s );
642 n_u(jj,iz) = length( i_u );
643 n_v(jj,iz) = length( i_v );
644 n_b(jj,iz) = length( i_b );
645 n_vb(jj,iz) = length( i_vb );
646 a_t(jj,iz) = sum( t_t( i_t ));
647 a_s(jj,iz) = sum( t_s( i_s ));
648 a_u(jj,iz) = sum( t_u( i_u ));
649 a_v(jj,iz) = sum( t_v( i_v ));
650 a_b(jj,iz) = sum( t_b( i_b ));
651 a_vb(jj,iz) = sum( t_vb( i_vb ));
652 end
653 end
654 fclose(zid_t); fclose(zid_s);
655 fclose(zid_u); fclose(zid_b);
656 fclose(zid_ub); fclose(zid_vb);
657 ssha_t = a_t ./ n_t;
658 ssha_s = a_s ./ n_s;
659 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 % 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
711
712 clear all
713 close all
714
715 do_plots = 0;
716
717 %==================================================================
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 Ri = R(1:(length(R)-1)) ...
759 + 0.25*delR(1:(length(R)-1)) + 0.25*delR(2:(length(R)));
760
761 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 ir = [ 1 2 3 5 10 15 20 25 30 35 40 ];
795 xc360 = XC + 180;
796
797 % ! 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 nc.author = 'Ed Hill <eh3@mit.edu>';
802 nc.date = 'March 27, 2005';
803 nc('X') = length(lon);
804 nc('Y') = length(lat);
805 nc('Z') = length(ir);
806 nc('Zc') = length(R);
807 nc('Zi') = length(Ri);
808 nc{'X'} = 'X';
809 nc{'Y'} = 'Y';
810 nc{'Z'} = 'Z';
811 nc{'Zc'} = 'Zc';
812 nc{'Zi'} = 'Zi';
813 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 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 nc{'X'}(:) = lon - 180;
834 nc{'Y'}(:) = lat;
835 nc{'Z'}(:) = R(ir);
836 nc{'Zc'}(:) = R;
837 nc{'Zi'}(:) = Ri;
838
839 f_s_3d = { {'tp2'}, {'sp2'}, {'bp2'}, {'vpbp_dbdz'}, ...
840 {'stress'}, {'dbdy'}, {'K'} };
841
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 % === 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 % save primes_9b2_04/stress.mat za_stress
936 % 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 load primes_92_04/ssha_uv_b_vpbp.mat
942 load primes_92_04/ssha_vptp.mat
943
944 ssh = sshvals(1:(length(sshvals)-1)) + 0.5*diff(sshvals);
945 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
950 lat = [-90:90];
951 latm = lat(1:(length(lat)-1)) + 0.5*diff(lat);
952
953 % 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 % ssha_stress = (v'b')/(db/dz) from stream-wise averages
1001 ssha_dbdz = zeros(size(ssha_b));
1002 ssha_dbdz(:) = NaN;
1003 ssha_dbdz(:,1) = (ssha_b(:,2) - ssha_b(:,1))/(Ri(2)-Ri(1));
1004 for iz = 2:48,
1005 ssha_dbdz(:,iz) = ...
1006 0.5*(ssha_b(:,iz) - ssha_b(:,iz-1))/(Ri(iz)-Ri(iz-1)) ...
1007 + 0.5*(ssha_b(:,iz+1) - ssha_b(:,iz))/(Ri(iz+1)-Ri(iz));
1008 end
1009 ssha_dbdz(:,49) = (ssha_b(:,49) - ssha_b(:,48))/(Ri(49)-Ri(48));
1010 for jj = 1:nsshi
1011 fac = 1000 * 4*pi/(24*3600) * sin(pi*ssh_yc(jj)/180);
1012 ssha_str(jj,:) = fac * (ssha_llvpbp(jj,:) ./ ssha_dbdz(jj,:));
1013 end
1014 if do_plots == 1
1015 surf(flipud(ssha_str')), view(2), shading flat, colorbar
1016 end
1017
1018
1019 % ! rm -f cube_22_zsa.nc
1020 nc = netcdf(['cube_22_zsa.nc'], 'clobber');
1021 nc.reference = [ 'The cube_22 zonal and stream-wise averages.' ];
1022 nc.author = 'Ed Hill <eh3@mit.edu>';
1023 nc.date = 'March 27, 2005';
1024 nc('Y') = length(latm);
1025 nc('Zc') = length(R);
1026 nc('Zi') = length(Ri);
1027 nc('SSH') = length(ssh);
1028 nc{'Y'} = 'Y';
1029 nc{'Zc'} = 'Zc';
1030 nc{'Zi'} = 'Zi';
1031 nc{'SSH'} = 'SSH';
1032 nc{'Y'}.uniquename = 'Y';
1033 nc{'Y'}.long_name = 'latitude';
1034 nc{'Y'}.gridtype = ncint(0);
1035 nc{'Y'}.units = 'degree_north';
1036 nc{'Zc'}.uniquename = 'Zc';
1037 nc{'Zc'}.long_name = 'depth_at_center';
1038 nc{'Zc'}.gridtype = ncint(0);
1039 nc{'Zc'}.units = 'm';
1040 nc{'Zi'}.uniquename = 'Zi';
1041 nc{'Zi'}.long_name = 'depth_at_interface';
1042 nc{'Zi'}.gridtype = ncint(0);
1043 nc{'Zi'}.units = 'm';
1044 nc{'SSH'}.uniquename = 'SSH';
1045 nc{'SSH'}.long_name = 'sea_surface_height';
1046 nc{'SSH'}.gridtype = ncint(0);
1047 nc{'SSH'}.units = 'm';
1048 nc{'Y'}(:) = latm;
1049 nc{'Zc'}(:) = R;
1050 nc{'Zi'}(:) = Ri;
1051 nc{'SSH'}(:) = ssh;
1052
1053 id = 'llzvpbp';
1054 nc{ id } = { 'Zc' 'Y' };
1055 nc{ id }.missing_value = ncdouble(NaN);
1056 nc{ id }.FillValue_ = ncdouble(0.0);
1057 nc{ id }(:) = permute(llzvpbp,[2 1]);
1058
1059 f_i = { {'za_ll_vpbp_dbdz'}, {'za_stress'} };
1060 for ip = 1:length(f_i)
1061 cell = f_i{ip};
1062 id = cell{1};
1063 disp([' ' id]);
1064 nc{ id } = { 'Zi' 'Y' };
1065 nc{ id }.missing_value = ncdouble(NaN);
1066 nc{ id }.FillValue_ = ncdouble(0.0);
1067 eval(sprintf('tmp = %s;', id));
1068 nc{ id }(:) = permute(tmp,[2 1]);
1069 end
1070
1071 f_i = { {'ssha_stress'}, ...
1072 {'ssha_llvptp'}, {'ssha_tp2'}, {'ssha_dtdy'}, ...
1073 {'ssha_K_t'}, {'ssha_K_b'}, ...
1074 {'ssha_dbdz'}, {'ssha_str'}, ...
1075 {'ssha_t'}, {'ssha_s'}, ...
1076 {'ssha_llu'}, {'ssha_llv'}, {'ssha_b'}, {'ssha_llvpbp'} };
1077 for ip = 1:length(f_i)
1078 cell = f_i{ip};
1079 id = cell{1};
1080 disp([' ' id]);
1081 nc{ id } = { 'Zi' 'SSH' };
1082 nc{ id }.missing_value = ncdouble(NaN);
1083 nc{ id }.FillValue_ = ncdouble(0.0);
1084 eval(sprintf('nc{ id }(:) = permute(%s,[2 1]);',id));
1085 end
1086 nc = close(nc);
1087
1088
1089 % ! rm -f cube_22_zsa_elat.nc
1090 nc = netcdf(['cube_22_zsa_elat.nc'], 'clobber');
1091 nc.reference = [ 'The cube_22 zonal and stream-wise averages.' ];
1092 nc.author = 'Ed Hill <eh3@mit.edu>';
1093 nc.date = 'March 27, 2005';
1094 nc('Y') = length(latm);
1095 nc('Zc') = length(R);
1096 nc('Zi') = length(Ri);
1097 nc('elat') = length(ssh_yc);
1098 nc{'Y'} = 'Y';
1099 nc{'Zc'} = 'Zc';
1100 nc{'Zi'} = 'Zi';
1101 nc{'elat'} = 'elat';
1102 nc{'Y'}.uniquename = 'Y';
1103 nc{'Y'}.long_name = 'latitude';
1104 nc{'Y'}.gridtype = ncint(0);
1105 nc{'Y'}.units = 'degree_north';
1106 nc{'Zc'}.uniquename = 'Zc';
1107 nc{'Zc'}.long_name = 'depth_at_center';
1108 nc{'Zc'}.gridtype = ncint(0);
1109 nc{'Zc'}.units = 'm';
1110 nc{'Zi'}.uniquename = 'Zi';
1111 nc{'Zi'}.long_name = 'depth_at_interface';
1112 nc{'Zi'}.gridtype = ncint(0);
1113 nc{'Zi'}.units = 'm';
1114 nc{'elat'}.uniquename = 'elat';
1115 nc{'elat'}.long_name = 'equivalent_latitude';
1116 nc{'elat'}.gridtype = ncint(0);
1117 nc{'elat'}.units = 'degree_north';
1118 nc{'Y'}(:) = latm;
1119 nc{'Zc'}(:) = R;
1120 nc{'Zi'}(:) = Ri;
1121 nc{'elat'}(:) = ssh_yc_sm;
1122
1123 id = 'llzvpbp';
1124 nc{ id } = { 'Zc' 'Y' };
1125 nc{ id }.missing_value = ncdouble(NaN);
1126 nc{ id }.FillValue_ = ncdouble(0.0);
1127 nc{ id }(:) = permute(llzvpbp,[2 1]);
1128
1129 f_i = { {'za_ll_vpbp_dbdz'}, {'za_stress'} };
1130 for ip = 1:length(f_i)
1131 cell = f_i{ip};
1132 id = cell{1};
1133 disp([' ' id]);
1134 nc{ id } = { 'Zi' 'Y' };
1135 nc{ id }.missing_value = ncdouble(NaN);
1136 nc{ id }.FillValue_ = ncdouble(0.0);
1137 eval(sprintf('tmp = %s;', id));
1138 nc{ id }(:) = permute(tmp,[2 1]);
1139 end
1140
1141 f_i = { {'ssha_stress'}, ...
1142 {'ssha_llvptp'}, {'ssha_tp2'}, {'ssha_dtdy'}, ...
1143 {'ssha_K_t'}, {'ssha_K_b'}, ...
1144 {'ssha_dbdz'}, {'ssha_str'}, ...
1145 {'ssha_t'}, {'ssha_s'}, ...
1146 {'ssha_llu'}, {'ssha_llv'}, {'ssha_b'}, {'ssha_llvpbp'} };
1147 for ip = 1:length(f_i)
1148 cell = f_i{ip};
1149 id = cell{1};
1150 disp([' ' id]);
1151 nc{ id } = { 'Zi' 'elat' };
1152 nc{ id }.missing_value = ncdouble(NaN);
1153 nc{ id }.FillValue_ = ncdouble(0.0);
1154 eval(sprintf('nc{ id }(:) = permute(%s,[2 1]);',id));
1155 end
1156
1157 nc = close(nc);
1158
1159 % surf(,,log(ssha_stress')), view(2),shading interp,colorbar
1160
1161 % ! scp cube_22_primes_at1deg.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1162 % ! scp cube_22_zsa.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1163 % ! scp cube_22_zsa_elat.nc channel.mit.edu:/home/edhill/INGRID_PEOPLE/EH3/eddy_flux/cube_22/
1164 % ! mv cube_22_primes_at1deg.nc primes_92_04
1165

  ViewVC Help
Powered by ViewVC 1.1.22