/[MITgcm]/MITgcm/verification/global_ocean.90x40x15/diags_matlab/mit_loadglobal.m
ViewVC logotype

Contents of /MITgcm/verification/global_ocean.90x40x15/diags_matlab/mit_loadglobal.m

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


Revision 1.1.2.1 - (show annotations) (download)
Wed Oct 23 18:26:35 2002 UTC (21 years, 6 months ago) by mlosch
Branch: release1
CVS Tags: release1_p12, release1_p10, release1_p16, release1_p15, release1_p11, release1_p14, release1_p13, release1_p17, release1_p8, release1_p9, release1_p6, release1_p7, release1_p13_pre, release1_p12_pre
Branch point for: release1_50yr
Changes since 1.1: +0 -0 lines
o fixed the verification/global_ocean.90x40x15 experiment:
 - new bathymetry (the world according to A., JMC, and M.)
 - new initial fields and forcing fields (*.bin files)
 - new POLY3.COEFFS (for the next release one should switch to a full
   equation of state)
 - fixed several errors and redundancies in the data file
 - experiment uses looped cells
 - added matlab directory with diagnostic scripts for plotting of output

1 % m-file: mit_loadglobal.m
2 % plot various diagnostics of the global 4x4 degree runs
3
4 % uses
5 % mit_loadgrid
6 % mit_getparm
7 % mit_oceanmasks
8 % mit_timesteps
9 % mit_loadglobaldata
10 % mit_globalmovie
11 % mit_plotstreamfunctions
12 % mit_plotzonalvelocity
13 % mit_globalmean
14 % mit_zonalmean
15 % mit_meridflux
16 % mit_plotmeandrift
17
18
19 % read in all grid files, etc. This has to be done at the very beginning!
20 grd = mit_loadgrid('.');
21 dname = grd.dname;
22 % add some masks
23 grd = mit_oceanmasks(grd);
24 %
25 if ~exist('meanfields','var')
26 meanfields = 1; % read in mean fields instead of instantaneous ones,
27 % include some more diagnostics
28 % meanfields = 1 is the default
29 end
30 %
31 mit_timesteps
32
33 if meanfields
34 disp('reading averaged fields (*tave.*)')
35 %
36 mit_loadglobaldata
37 %
38 u=rdmds('uVeltave',timesteps(2:end)');
39 v=rdmds('vVeltave',timesteps(2:end)');
40 w=rdmds('wVeltave',timesteps(2:end)');
41 t=rdmds('Ttave',timesteps(2:end)');
42 s=rdmds('Stave',timesteps(2:end)');
43 eta=rdmds('ETAtave',timesteps(2:end)').*repmat(squeeze(grd.cmask(:,:,1)),[1 1 nt-1]);
44 k=1;
45 u=cat(4,zeros([nx ny nz 1]),u);
46 v=cat(4,zeros([nx ny nz 1]),v);
47 w=cat(4,zeros([nx ny nz 1]),w);
48 % hack
49 t=cat(4,rdmds('T',timesteps(1)'),t);
50 s=cat(4,rdmds('S',timesteps(1)'),s);
51 eta=cat(3,0.*grd.cmask(:,:,1),eta);
52 else
53 % load snap shots
54 disp('reading snap shot fields')
55 u=rdmds('U',timesteps');
56 v=rdmds('V',timesteps');
57 w=rdmds('W',timesteps');
58 t=rdmds('T',timesteps');
59 s=rdmds('S',timesteps');
60 eta=rdmds('Eta',timesteps').*repmat(squeeze(grd.cmask(:,:,1)),[1 1 nt]);
61 end
62 % orientation and masking
63 if ~strcmp(grd.buoyancy,'OCEANIC');
64 u = u(:,:,end:-1:1,:);
65 v = v(:,:,end:-1:1,:);
66 w = w(:,:,end:-1:1,:).*repmat(grd.cmask,[1 1 1 nt]);
67 t = t(:,:,end:-1:1,:).*repmat(grd.cmask,[1 1 1 nt]);
68 s = s(:,:,end:-1:1,:).*repmat(grd.cmask,[1 1 1 nt]);
69 else
70 % $$$ u = u;
71 % $$$ v = v;
72 w = w.*repmat(grd.cmask,[1 1 1 nt]);
73 t = t.*repmat(grd.cmask,[1 1 1 nt]);
74 s = s.*repmat(grd.cmask,[1 1 1 nt]);
75 end
76
77 % transport through Drake Passage
78 TDP = NaN*ones(nt,1);
79 kx = 73;
80 kyg = [5 6 7];
81 %kyg = [4 5 6 7];
82 da = grd.dz*grd.dyg(kx,kyg);
83 for k=kt;
84 TDP(k) = sum(nansum(squeeze(u(kx,kyg,:,k))'.*da))*1e-6; % in Sv
85 TDPz(:,k) = nansum(squeeze(u(kx,kyg,:,k)).*da')'*1e-6; % in Sv
86 end
87
88 % set some depth
89 if ~meanfields
90 mit_globalmovie
91 end
92
93 % $$$ [x y z]=meshgrid(grd.lonc,grd.grd.latc,-grd.zc);
94 % $$$ tk = permute(squeeze(t(:,:,:,k)),[2 1 3]);
95 % $$$ xplane = [0:30:360];
96 % $$$ yplane = 30;
97 % $$$ zplane = -[50 500:500:6000];
98 % $$$ slice(x,y,z,tk,xplane,yplane,zplane);
99 % $$$ shading flat, colorbar('h')
100
101 mit_plotstreamfunctions
102
103 mit_plotzonalvelocity
104
105 if meanfields
106 % compute actual heat and E-P fluxes from time averages
107 % this only makes sense when you load the time averaged fields
108 tauTheta = mit_getparm('data','tauThetaClimRelax');
109 tauSalt = mit_getparm('data','tauSaltClimRelax');
110 rhonil = mit_getparm('data','rhonil');
111 if isempty(rhonil)
112 rhonil = 1035;
113 end
114 Cp = 3994;
115 if isempty(tauTheta)
116 tauTheta = 0;
117 end
118 if tauTheta==0
119 recip_tauTheta = 0.;
120 else
121 recip_tauTheta = 1./tauTheta;
122 end
123 if isempty(tauSalt)
124 tauSalt = 0;
125 end
126 if tauSalt==0
127 recip_tauSalt = 0.;
128 else
129 recip_tauSalt = 1./tauSalt;
130 end
131 figure('PaperPosition',[0.31 0.25 10.5 7.88],'PaperOrientation','landscape')
132 k=kmax;
133 qnet_diag = (t(:,:,1,k)-tdmean)*grd.dz(1).*grd.hfacc(:,:,1)*Cp*rhonil*recip_tauTheta;
134 empr_diag = - (s(:,:,1,k)-sdmean)*grd.dz(1).*grd.hfacc(:,:,1)*recip_tauSalt/35;
135 clear sh;
136 sh(1) = subplot(2,2,1);imagesc(grd.lonc,grd.latc,qnet_diag');
137 axis xy; axis image; colorbar('h')
138 title(['Q_{diag} from T-T_{lev} at the surface'])
139 sh(2) =subplot(2,2,2);imagesc(grd.lonc,grd.latc,qnet'.*grd.cmask(:,:,1)');
140 axis xy; axis image; colorbar('h')
141 title('annual mean Q_{net} (data)')
142 sh(3) =subplot(2,2,3);imagesc(grd.lonc,grd.latc,empr_diag');
143 axis xy; axis image; colorbar('h')
144 title(['(E-P)_{diag} from S-S_{lev} at the surface'])
145 sh(4) =subplot(2,2,4);imagesc(grd.lonc,grd.latc,empr'.*grd.cmask(:,:,1)');
146 axis xy; axis image; colorbar('h')
147 title('annual mean E-P (data)')
148 suptitle(['experiment ' dname ', timestep = ' num2str(timesteps(k)) ...
149 ', ' tuname ' = ' num2str(tim(k))])
150 set(sh(1:2),'clim',[-1 1]*200);
151 set(gcf,'currentAxes',sh(1));colorbar('h')
152 set(gcf,'currentAxes',sh(2));colorbar('h')
153 set(sh(3:4),'clim',[-1 1]*1.e-7);
154 set(gcf,'currentAxes',sh(3));colorbar('h')
155 set(gcf,'currentAxes',sh(4));colorbar('h')
156
157 clear sh
158
159 if length(tim) > 1
160 rac3d = repmat(grd.rac,[1 1 nz]).*grd.hfacc;
161 % horizontal averages of t and s in a hovmoeller-type diagramme
162 mdt = repmat(NaN,grd.nz,length(kt));
163 mdtt= mdt;
164 mdtstd= mdt;
165 mdt_alt = mdt;
166 mdtt_alt= mdt;
167 mdtstd_alt = mdt;
168 mdt_pac = mdt;
169 mdtt_pac= mdt;
170 mdtstd_pac = mdt;
171 mdt_sou = mdt;
172 mdtt_sou= mdt;
173 mdtstd_sou = mdt;
174 mds = mdt;
175 mdss= mdt;
176 mdsstd = mdt;
177 mds_alt = mdt;
178 mdss_alt= mdt;
179 mdsstd_alt = mdt;
180 mds_pac = mdt;
181 mdss_pac= mdt;
182 mdsstd_pac = mdt;
183 mds_sou = mdt;
184 mdss_sou= mdt;
185 mdsstd_sou = mdt;
186 indopac_hfacc = ...
187 change(change(grd.pacific_hfacc,'==',NaN,0) ...
188 +change(grd.indic_hfacc,'==',NaN,0),'==',0,NaN);
189 indopac_hfacc(:,1:12,:) = NaN;
190 southern_hfacc = grd.hfacc;
191 southern_hfacc(:,13:end,:) = NaN;
192 atl_hfacc = grd.atlantic_hfacc;
193 atl_hfacc(:,1:12,:) = NaN;
194 for k=kt;
195 dt = t(:,:,:,k)-tdatamean; % anomalies
196 ds = s(:,:,:,k)-sdatamean; % anomalies
197 [mdt(:,k) mdtt(:,k) mdtstd(:,k)] = mit_globalmean(dt,rac3d);
198 [mds(:,k) mdss(:,k) mdsstd(:,k)] = mit_globalmean(ds,rac3d);
199 [mdt_atl(:,k) mdtt_atl(:,k) mdtstd_atl(:,k)] = mit_globalmean(dt,rac3d.*atl_hfacc);
200 [mds_atl(:,k) mdss_atl(:,k) mdsstd_atl(:,k)] = mit_globalmean(ds,rac3d.*atl_hfacc);
201 [mdt_pac(:,k) mdtt_pac(:,k) mdtstd_pac(:,k)] = mit_globalmean(dt,rac3d.*indopac_hfacc);
202 [mds_pac(:,k) mdss_pac(:,k) mdsstd_pac(:,k)] = mit_globalmean(ds,rac3d.*indopac_hfacc);
203 [mdt_sou(:,k) mdtt_sou(:,k) mdtstd_sou(:,k)] = mit_globalmean(dt,rac3d.*southern_hfacc);
204 [mds_sou(:,k) mdss_sou(:,k) mdsstd_sou(:,k)] = mit_globalmean(ds,rac3d.*southern_hfacc);
205 end
206
207 figure('PaperPosition',[0.31 0.25 10.5 7.88],'PaperOrientation','landscape')
208 if length(tim) > 2
209 clear sh
210 sh(1) = subplot(2,2,1);
211 contourf(tim(2:end),-grd.zc,mdt(:,2:end),20);
212 caxis([-1 1]*max(abs(caxis))); shading flat; colorbar('v')
213 title('T-T_{lev} horizontally averaged')
214 xlabel(['Time [' timeunit ']']); ylabel('Depth [m]')
215 sh(2) = subplot(2,2,2);
216 contourf(tim(2:end),-grd.zc,mds(:,2:end),20);
217 caxis([-1 1]*max(abs(caxis))); shading flat; colorbar('v')
218 title('S-S_{lev} horizontally averaged')
219 xlabel(['Time [' timeunit ']']); ylabel('Depth [m]')
220 %
221 sh(3) = subplot(2,2,3);
222 contourf(tim(2:end),-grd.zc,mdtt(:,2:end),20);
223 shading flat; colorbar('v')
224 title('|T-T_{lev}| horizontally averaged');
225 xlabel(['Time [' timeunit ']']); ylabel('Depth [m]')
226 sh(4) = subplot(2,2,4);
227 contourf(tim(2:end),-grd.zc,mdss(:,2:end),20);
228 shading flat; colorbar('v')
229 title('|S-S_{lev}| horizontally averaged')
230 xlabel(['Time [' timeunit ']']); ylabel('Depth [m]')
231 set(sh,'layer','top');
232 suptitle(['experiment ' dname])
233 end %if length(tim) > 2
234
235 k=kmax;
236 figure('PaperOrientation','landscape','PaperPosition',[0.25 0.3125 10.5 7.875])
237 clear sh
238 sh(1) = subplot(2,4,1);
239 plot(mdt(:,k),-grd.zc,'b-',mdt(:,k)+mdtstd(:,k),-grd.zc,'r--',mdt(:,k)-mdtstd(:,k),-grd.zc,'r--');
240 xlabel('T-T_{lev} [degC]'); ylabel('depth [m]'); title('global')
241 legend('horiz. mean','std',3);
242 sh(2) = subplot(2,4,2);
243 plot(mdt_atl(:,k),-grd.zc,'b-',mdt_atl(:,k)+mdtstd_atl(:,k),-grd.zc,'r--',mdt_atl(:,k)-mdtstd_atl(:,k),-grd.zc,'r--');
244 xlabel('T-T_{lev} [degC]'); %ylabel('depth [m]');
245 title('atlantic')
246 sh(3) = subplot(2,4,3);
247 plot(mdt_pac(:,k),-grd.zc,'b-',mdt_pac(:,k)+mdtstd_pac(:,k),-grd.zc,'r--',mdt_pac(:,k)-mdtstd_pac(:,k),-grd.zc,'r--');
248 xlabel('T-T_{lev} [degC]'); %ylabel('depth [m]');
249 title('indo-pacific')
250 sh(4) = subplot(2,4,4);
251 plot(mdt_sou(:,k),-grd.zc,'b-',mdt_sou(:,k)+mdtstd_sou(:,k),-grd.zc,'r--',mdt_sou(:,k)-mdtstd_sou(:,k),-grd.zc,'r--');
252 xlabel('T-T_{lev} [degC]'); ylabel('depth [m]'); title('southern')
253 %
254 sh(5) = subplot(2,4,5);
255 plot(mds(:,k),-grd.zc,'b-',mds(:,k)+mdsstd(:,k),-grd.zc,'r--',mds(:,k)-mdsstd(:,k),-grd.zc,'r--');
256 xlabel('S-S_{lev} [PSU]'); ylabel('depth [m]');
257 sh(6) = subplot(2,4,6);
258 plot(mds_atl(:,k),-grd.zc,'b-',mds_atl(:,k)+mdsstd_atl(:,k),-grd.zc,'r--',mds_atl(:,k)-mdsstd_atl(:,k),-grd.zc,'r--');
259 xlabel('S-S_{lev} [PSU]'); %ylabel('depth [m]');
260 sh(7) = subplot(2,4,7);
261 plot(mds_pac(:,k),-grd.zc,'b-',mds_pac(:,k)+mdsstd_pac(:,k),-grd.zc,'r--',mds_pac(:,k)-mdsstd_pac(:,k),-grd.zc,'r--');
262 xlabel('S-S_{lev} [PSU]'); %ylabel('depth [m]');
263 sh(8) = subplot(2,4,8);
264 plot(mds_sou(:,k),-grd.zc,'b-',mds_sou(:,k)+mdsstd_sou(:,k),-grd.zc,'r--',mds_sou(:,k)-mdsstd_sou(:,k),-grd.zc,'r--');
265 xlabel('S-S_{lev} [PSU]'); ylabel('depth [m]');
266 set(sh,'xgrid','on','ygrid','on')
267 set([sh(4); sh(8)],'YAxisLocation','right')
268 suptitle(['experiment ' dname ', timestep = ' num2str(timesteps(k)) ...
269 ', ' tuname ' = ' num2str(tim(k))])
270 end %if length(tim) > 1
271
272 % zonal mean of temperature and salinity for different ocean basins
273 k=kmax;
274 tdzm_glo = mit_zonalmean(tdatamean,grd.hfacc,grd.dxc);
275 tdzm_atl = mit_zonalmean(tdatamean,grd.atlantic_hfacc,grd.dxc);
276 tdzm_pac = mit_zonalmean(tdatamean,grd.pacific_hfacc,grd.dxc);
277 tdzm_ind = mit_zonalmean(tdatamean,grd.indic_hfacc,grd.dxc);
278 sdzm_glo = mit_zonalmean(sdatamean,grd.hfacc,grd.dxc);
279 sdzm_atl = mit_zonalmean(sdatamean,grd.atlantic_hfacc,grd.dxc);
280 sdzm_pac = mit_zonalmean(sdatamean,grd.pacific_hfacc,grd.dxc);
281 sdzm_ind = mit_zonalmean(sdatamean,grd.indic_hfacc,grd.dxc);
282
283 figure('PaperPosition',[0.25 0.368552 8 10.2629]);
284 tlev = -5:.5:5;
285 slev = -1:.1:1;
286 clear sh clh
287 clh = cell(8,1);
288 delay = 1;
289 for k = kmax
290 tzm_glo = mit_zonalmean(t(:,:,:,k),grd.hfacc,grd.dxc);
291 tzm_atl = mit_zonalmean(t(:,:,:,k),grd.atlantic_hfacc,grd.dxc);
292 tzm_pac = mit_zonalmean(t(:,:,:,k),grd.pacific_hfacc,grd.dxc);
293 tzm_ind = mit_zonalmean(t(:,:,:,k),grd.indic_hfacc,grd.dxc);
294 szm_glo = mit_zonalmean(s(:,:,:,k),grd.hfacc,grd.dxc);
295 szm_atl = mit_zonalmean(s(:,:,:,k),grd.atlantic_hfacc,grd.dxc);
296 szm_pac = mit_zonalmean(s(:,:,:,k),grd.pacific_hfacc,grd.dxc);
297 szm_ind = mit_zonalmean(s(:,:,:,k),grd.indic_hfacc,grd.dxc);
298 caxt = [min(tzm_glo(:)-tdzm_glo(:)) max(tzm_glo(:)-tdzm_glo(:))];
299 tlev = -5:.5:5;
300 if max(abs(caxt)) < 1;
301 tlev = -1:.1:1;
302 end
303 if max(abs(caxt)) < .1;
304 tlev = .1*tlev;
305 end
306 caxs = [min(szm_glo(:)-sdzm_glo(:)) max(szm_glo(:)-sdzm_glo(:))];
307 slev = -1.:.05:1;
308 if max(abs(caxs)) < 0.2;
309 slev = -.20:.01:.20;
310 end
311 if max(abs(caxs)) < .02;
312 slev = .1*slev;
313 end
314 sh(1) = subplot(4,2,1);
315 [cs h] = contourf(grd.latc,-grd.zc,(tzm_glo-tdzm_glo)',tlev);
316 if ~isempty(h);
317 clh{1} = clabel(cs,h,[0 0]);
318 end
319 title('\theta-\theta_{lev} [degC]: global ocean')
320 sh(3) = subplot(4,2,3);
321 [cs h] = contourf(grd.latc,-grd.zc,(tzm_atl-tdzm_atl)',tlev);
322 if ~isempty(h);
323 clh{3} = clabel(cs,h,[0 0]);
324 end
325 title('atlantic ocean')
326 sh(5) = subplot(4,2,5);
327 [cs h] = contourf(grd.latc,-grd.zc,(tzm_pac-tdzm_pac)',tlev);
328 if ~isempty(h);
329 clh{5} = clabel(cs,h,[0 0]);
330 end
331 title('pacific ocean')
332 sh(7) = subplot(4,2,7);
333 [cs h] = contourf(grd.latc,-grd.zc,(tzm_ind-tdzm_ind)',tlev);
334 if ~isempty(h);
335 clh{7} = clabel(cs,h,[0 0]);
336 end
337 title('indian ocean')
338 set(sh(1:2:end),'clim',[tlev(1) tlev(end)])
339 colorbar('h')
340 sh(2) = subplot(4,2,2);
341 [cs h] = contourf(grd.latc,-grd.zc,(szm_glo-sdzm_glo)',slev);
342 if ~isempty(h);
343 clh{2} = clabel(cs,h,[0 0]);
344 end
345 title('S-S_{lev} [PSU]: global ocean')
346 sh(4) = subplot(4,2,4);
347 [cs h] = contourf(grd.latc,-grd.zc,(szm_atl-sdzm_atl)',slev);
348 if ~isempty(h);
349 clh{4} = clabel(cs,h,[0 0]);
350 end
351 title('atlantic ocean')
352 sh(6) = subplot(4,2,6);
353 [cs h] = contourf(grd.latc,-grd.zc,(szm_pac-sdzm_pac)',slev);
354 if ~isempty(h);
355 clh{6} = clabel(cs,h,[0 0]);
356 end
357 title('pacific ocean')
358 sh(8) = subplot(4,2,8);
359 [cs h] = contourf(grd.latc,-grd.zc,(szm_ind-sdzm_ind)',slev);
360 if ~isempty(h);
361 clh{8} = clabel(cs,h,[0 0]);
362 end
363 title('indian ocean')
364 set(sh(2:2:end),'clim',[slev(1) slev(end)])
365 colorbar('h')
366 set(sh,'layer','top')
367 if ~isempty(cat(1,clh{:}))
368 set(cat(1,clh{:}),'fontsize',8);
369 end
370
371 suptitle(['experiment ' dname ', timestep = ' num2str(timesteps(k)) ...
372 ', ' tuname ' = ' num2str(tim(k)) ', zonal averages'])
373 drawnow;
374 pause(delay)
375 end
376 clear clh sh
377
378 % meridional transport of heat (internal energy) and fresh water
379 % diagnosed from surface fluxes
380 petawatts = 1e-15;
381 sverdrups = 1e-6;
382 heat_data = -mit_meridflux(qnet,grd.dxc,grd.dyc)*petawatts;
383 heat_diag = -mit_meridflux(qnet_diag,grd.dxc,grd.dyc)*petawatts;
384 freshwater_data = -mit_meridflux(empr,grd.dxc,grd.dyc)*sverdrups;
385 freshwater_diag = -mit_meridflux(empr_diag,grd.dxc,grd.dyc)*sverdrups;
386 figure
387 latgf = [2*grd.latc-grd.latg];
388 clear sh
389 sh(1) = subplot(2,1,1);
390 hh = plot(latgf,heat_data,'-',latgf,heat_diag,'--', ...
391 latgf,heat_data+heat_diag,'-.');
392 title('heat (internal energy) flux [PW]')
393 suptitle(['experiment ' dname ', timestep = ' num2str(timesteps(k)) ...
394 ', ' tuname ' = ' num2str(tim(k))])
395 legend('surface flux','restoring','net',2)
396 sh(2) = subplot(2,1,2);
397 fwh = plot(latgf,freshwater_data,'-',latgf,freshwater_diag,'--', ...
398 latgf,freshwater_data+freshwater_diag,'-.');
399 title('freshwater flux [Sv]')
400 legend('surface flux','restoring','net',3)
401 set(sh,'xlim',[latgf(1) latgf(end)],'xgrid','on','ygrid','on');
402
403 end % meanfields
404
405 mit_plotmeandrift
406
407 fname = [grd.dname '_' sprintf('%u',tim(kmax)) timeunit];
408 in = findstr(fname,'\');
409 fname(in) = [];
410 fn = gcf;
411 if meanfields
412 fname = [fname '.ps'];
413 f0 = fn-7;
414 else
415 fname = [fname '_snap.ps'];
416 f0 = fn-4;
417 end
418 print(f0,'-dpsc',fname);
419 for k = f0+1:fn
420 print(k,'-dpsc',fname,'-append')
421 end

  ViewVC Help
Powered by ViewVC 1.1.22