/[MITgcm]/MITgcm_contrib/ESMF/global_ocean.128x60x15/diags_matlab/mit_loadglobal.m
ViewVC logotype

Annotation of /MITgcm_contrib/ESMF/global_ocean.128x60x15/diags_matlab/mit_loadglobal.m

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


Revision 1.1 - (hide annotations) (download)
Tue Mar 30 03:58:59 2004 UTC (21 years, 4 months ago) by cnh
Branch: MAIN
CVS Tags: adoption_1_0_pre_A, HEAD
New test with different size

1 cnh 1.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