/[MITgcm]/MITgcm/verification/natl_box/matlab/comp_output.m
ViewVC logotype

Annotation of /MITgcm/verification/natl_box/matlab/comp_output.m

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


Revision 1.1 - (hide annotations) (download)
Mon Nov 13 16:02:32 2000 UTC (23 years, 5 months ago) by heimbach
Branch: MAIN
Adding verification experiment for KPP.

1 heimbach 1.1 % Compare output of new and reference (c32) codes.
2    
3     p1='../output/'; % reference (c32)output location
4     p2='../../../exe/'; % new output file location
5    
6     lat=13:2:43; lon=322:2:360; % latitude, longitude
7    
8     % model depths
9     dZ=[10 10 15 20 20 25 35 50 75 100 150 200 275 ...
10     350 415 450 500 500 500 500 500 500 500];
11     dpt=dZ/2;
12     for i=2:length(dZ)
13     dpt(i)=dpt(i)+sum(dZ(1:(i-1)));
14     end
15    
16     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18    
19     % load and compare temperature
20    
21     tx=4; % time index
22     dp=1; % depth level
23     cx=[15 28]; cxd=[-1 1]*.005; % color axes
24    
25     T1=readbin([p1 'T.001.001.data'],[20 16 23 4],1);
26     T2=readbin([p2 'T.001.001.data'],[20 16 23 4],1);
27    
28     figure(1), clf reset
29     set(gcf,'PaperOrientation','portrait')
30     set(gcf,'PaperPosition',[0.5 0.5 7.5 10.])
31     subplot(311), pcolor(lon,lat,T1(:,:,dp,tx)');
32     shading flat, caxis(cx), colorbar
33     title(['c32 temperature (deg C), day 30' ...
34     ', ' int2str(dpt(dp)) ' m depth' ...
35     ', ' ' min=' num2str(min(min(T1(:,:,dp,tx))),4) ...
36     ', ' ' max=' num2str(max(max(T1(:,:,dp,tx))),4) ])
37     ylabel('Latitude North')
38    
39     subplot(312), pcolor(lon,lat,T2(:,:,dp,tx)');
40     shading flat, caxis(cx), colorbar
41     title(['new temperature (deg C), day 30' ...
42     ', ' ' min=' num2str(min(min(T1(:,:,dp,tx))),4) ...
43     ', ' ' max=' num2str(max(max(T1(:,:,dp,tx))),4) ])
44     ylabel('Latitude North')
45    
46     subplot(313), pcolor(lon,lat,T2(:,:,dp,tx)'-T1(:,:,dp,tx)');
47     shading flat, caxis(cxd), colorbar
48     title(['Difference (deg C)' ...
49     ', ' ' min=' num2str(min(min(T2(:,:,dp,tx)'-T1(:,:,dp,tx)')),4) ...
50     ', ' ' max=' num2str(max(max(T2(:,:,dp,tx)'-T1(:,:,dp,tx)')),4) ])
51     ylabel('Latitude North'), xlabel('Longitude East')
52     orient tall
53     filename = 'comp_temp.eps'
54     eval([ 'print -depsc ', filename ])
55    
56     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58    
59     % load and compare boundary layer depth
60    
61     tx=4; % time index
62     cx=[0 55]; cxd=[-1 1]*4; % color axes
63    
64     H1=readbin([p1 'KPPhbl.001.001.data'],[20 16 1 4],1);
65     H2=readbin([p2 'KPPhbl.001.001.data'],[20 16 1 4],1);
66    
67     figure(2), clf reset
68     set(gcf,'PaperOrientation','portrait')
69     set(gcf,'PaperPosition',[0.5 0.5 7.5 10.])
70     subplot(311), pcolor(lon,lat,H1(:,:,1,tx)');
71     shading flat, caxis(cx), colorbar
72     title(['c32 boundary layer depth (m), day 30' ...
73     ', ' ' min=' num2str(min(min(H1(:,:,1,tx))),4) ...
74     ', ' ' max=' num2str(max(max(H1(:,:,1,tx))),4) ])
75     ylabel('Latitude North')
76    
77     subplot(312), pcolor(lon,lat,H2(:,:,1,tx)');
78     shading flat, caxis(cx), colorbar
79     title([ 'new boundary layer depth (m), day 30' ...
80     ', ' ' min=' num2str(min(min(H2(:,:,1,tx))),4) ...
81     ', ' ' max=' num2str(max(max(H2(:,:,1,tx))),4)] )
82     ylabel('Latitude North')
83    
84     subplot(313), pcolor(lon,lat,H2(:,:,1,tx)'-H1(:,:,1,tx)');
85     shading flat, caxis(cxd), colorbar
86     title(['Difference (m)' ...
87     ', ' ' min=' num2str(min(min(H2(:,:,1,tx)'-H1(:,:,1,tx)')),4) ...
88     ', ' ' max=' num2str(max(max(H2(:,:,1,tx)'-H1(:,:,1,tx)')),4)] )
89     ylabel('Latitude North'), xlabel('Longitude East')
90     orient tall
91     filename = 'comp_bldepth.eps'
92     eval([ 'print -depsc ', filename ])
93    
94     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
95     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96    
97     % load and compare KPP diffusivity
98    
99     tx=4; % time index
100     id=2:23; % depth index
101     cx=[-5 -1.4]; cxd=[-1 1]*.003; % color axes
102     yticks=['.00001';'.0001 ';'.001 ';'.01 '];
103     dticks=[3000 1000 300 100 30];
104    
105     D1=readbin([p1 'KPPdiffKzT.001.001.data'],[20 16 23 4],1);
106     D2=readbin([p2 'KPPdiffKzT.001.001.data'],[20 16 23 4],1);
107    
108     figure(3), clf reset
109     set(gcf,'PaperOrientation','portrait')
110     set(gcf,'PaperPosition',[0.5 0.5 7.5 10.])
111    
112     tmp=squeeze(D1(10,:,id,tx))'; tmp(find(~tmp))=1e-10; tmp1=log10(tmp);
113     subplot(311), pcolor(lat,-log10(dpt(id)),tmp1);
114     shading interp, caxis(cx), h=colorbar;
115     set(h,'YTick',-5:-2,'YTickLabel',yticks)
116     set(gca,'YTick',-log10(dticks),'YTickLabel',dticks)
117     title(['c32 KPP diffusivity (m^2/s), day 30' ...
118     ', ' ' min=' num2str(min(min(D1(10,:,id,tx))),4) ...
119     ', ' ' max=' num2str(max(max(D1(10,:,id,tx))),4)])
120     ylabel('Depth (m)')
121    
122     tmp=squeeze(D2(10,:,id,tx))'; tmp(find(~tmp))=1e-10; tmp2=log10(tmp);
123     subplot(312), pcolor(lat,-log10(dpt(id)),tmp2);
124     shading interp, caxis(cx), h=colorbar;
125     set(h,'YTick',-5:-2,'YTickLabel',yticks)
126     set(gca,'YTick',-log10(dticks),'YTickLabel',dticks)
127     title([ 'new KPP diffusivity (m^2/s), day 30' ...
128     ', ' ' min=' num2str(min(min(D2(10,:,id,tx))),4) ...
129     ', ' ' max=' num2str(max(max(D2(10,:,id,tx))),4) ])
130     ylabel('Depth (m)')
131    
132     subplot(313), pcolor(lat,-log10(dpt(id)),10.^tmp2-10.^tmp1);
133     shading interp, caxis(cxd), colorbar
134     set(gca,'YTick',-log10(dticks),'YTickLabel',dticks)
135     title([ 'Difference (m^2/s)' ...
136     ', ' ' min=' num2str(min(min(10.^tmp2-10.^tmp1)),4) ...
137     ', ' ' max=' num2str(max(max(10.^tmp2-10.^tmp1)),4) ])
138     ylabel('Depth (m)'), xlabel('Latitude North')
139     orient tall
140     filename = 'comp_diffus.eps'
141     eval([ 'print -depsc ', filename ])

  ViewVC Help
Powered by ViewVC 1.1.22