1 |
mlosch |
1.1 |
% This is a script that executes a series of commands for |
2 |
|
|
% extracting Levitus data, interpolating, and gap-filling |
3 |
|
|
% according to the mask-file generated by "topo.m". |
4 |
|
|
% |
5 |
|
|
% Uses data files: bathymetry.bin pmask.bin |
6 |
|
|
% Creates data files: lev_clim_*.bin lev_monthly_*.bin |
7 |
|
|
% Creates arrays: n/a |
8 |
|
|
% Uses m-scripts: setgrid grid_sphere finegrid gen_hxhy xyrecur_* |
9 |
|
|
% hfac |
10 |
|
|
% |
11 |
|
|
% Created 08/15/99 by adcroft@mit.edu |
12 |
|
|
% Modified 11/09/99 by adcroft@mit.edu |
13 |
|
|
% Maintained by adcroft@mit.edu, abiastoch@ucsd.edu |
14 |
|
|
|
15 |
|
|
skip_interp=0 |
16 |
|
|
|
17 |
|
|
if ~skip_interp |
18 |
|
|
%clear all |
19 |
|
|
format compact |
20 |
|
|
more off |
21 |
|
|
|
22 |
|
|
load FMT |
23 |
|
|
load HGRID.mat lon_lo lon_hi lat_lo lat_hi xc yc GRID |
24 |
|
|
dd=0; |
25 |
|
|
load VGRID.mat zc |
26 |
|
|
load MASKS |
27 |
|
|
GRID.mskc=msk; |
28 |
|
|
clear msk |
29 |
|
|
|
30 |
|
|
disp('Reading Levitus climatological T'); |
31 |
|
|
[TC,xl,yl,zlc]=... |
32 |
|
|
extract_levitus_clim('temp',lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd); |
33 |
|
|
if size(xc,3)==1 |
34 |
|
|
tc=interp3ddata(xl,yl,zlc,TC,xc,yc,zc,GRID.mskc); |
35 |
|
|
tc2=interp3ddata_int2(xl,yl,zlc,TC,xc,yc,zc,GRID.mskc); |
36 |
|
|
else |
37 |
|
|
for t=1:size(xc,3); |
38 |
|
|
tc(:,:,t,:)=interp3ddata(xl,yl,zlc,TC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:))); |
39 |
|
|
tc2(:,:,t,:)=interp3ddata_int2(xl,yl,zlc,TC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:))); |
40 |
|
|
end |
41 |
|
|
end |
42 |
|
|
check3dmask(tc); |
43 |
|
|
check3dmask(tc2); |
44 |
|
|
save tc tc tc2 |
45 |
|
|
%wrda('temp.bin',tc2,1,fmt,Ieee); |
46 |
|
|
|
47 |
|
|
disp('Reading Levitus climatological S'); |
48 |
|
|
[SC,xl,yl,zlc]=... |
49 |
|
|
extract_levitus_clim('salt',lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd); |
50 |
|
|
if size(xc,3)==1 |
51 |
|
|
sc=interp3ddata(xl,yl,zlc,SC,xc,yc,zc,GRID.mskc); |
52 |
|
|
sc2=interp3ddata_int2(xl,yl,zlc,SC,xc,yc,zc,GRID.mskc); |
53 |
|
|
else |
54 |
|
|
for t=1:size(xc,3); |
55 |
|
|
sc(:,:,t,:)=interp3ddata(xl,yl,zlc,SC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:))); |
56 |
|
|
sc2(:,:,t,:)=interp3ddata_int2(xl,yl,zlc,SC,xc(:,:,t),yc(:,:,t),zc,squeeze(GRID.mskc(:,:,t,:))); |
57 |
|
|
end |
58 |
|
|
end |
59 |
|
|
check3dmask(sc); |
60 |
|
|
check3dmask(sc2); |
61 |
|
|
save sc sc sc2 |
62 |
|
|
%wrda('salt.bin',sc2,1,fmt,Ieee); |
63 |
|
|
|
64 |
|
|
stats( mask(tc(:))-mask(sc(:)) ) |
65 |
|
|
%clear tc sc tc2 sc2 |
66 |
|
|
|
67 |
|
|
|
68 |
|
|
% Levitus climatology grid |
69 |
|
|
[LCGRID.xc,LCGRID.yc,LCGRID.xg,LCGRID.yg,LCGRID.rac]=grid_sphere( ... |
70 |
|
|
(xl(2)-xl(1))*ones(size(xl)),(yl(2)-yl(1))*ones(size(yl)), ... |
71 |
|
|
1.5*xl(1)-0.5*xl(2),1.5*yl(1)-0.5*yl(2)); |
72 |
|
|
LCGRID.mskc=mask(TC); |
73 |
|
|
[TCm,TCsd] = meanprofile(LCGRID,TC); |
74 |
|
|
[SCm,SCsd] = meanprofile(LCGRID,SC); |
75 |
|
|
stats( mask(TC(:))-mask(SC(:)) ) |
76 |
|
|
clear TC SC |
77 |
|
|
|
78 |
|
|
disp('Reading Levitus monthly T'); |
79 |
|
|
[TM,xl,yl,zlm]=... |
80 |
|
|
extract_levitus_monthly('temp',1:12,lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd); |
81 |
|
|
kk=find(zc>=min(zlm)); |
82 |
|
|
for m=1:12; |
83 |
|
|
if size(xc,3)==1 |
84 |
|
|
tm(:,:,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
85 |
|
|
tm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
86 |
|
|
else |
87 |
|
|
for t=1:size(xc,3); |
88 |
|
|
tm(:,:,t,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t))); |
89 |
|
|
tm2(:,:,t,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t))); |
90 |
|
|
%tm(:,:,:,m)=interp3ddata(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
91 |
|
|
%tm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,TM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
92 |
|
|
end |
93 |
|
|
end |
94 |
|
|
end |
95 |
|
|
%wrda('sst.bin',tm2(:,:,1,:),1,fmt,Ieee); |
96 |
|
|
save tm tm tm2 |
97 |
|
|
|
98 |
|
|
disp('Reading Levitus monthly S'); |
99 |
|
|
[SM,xl,yl,zlm]=... |
100 |
|
|
extract_levitus_monthly('salt',1:12,lon_lo-dd,lon_hi+dd,lat_lo-dd,lat_hi+dd); |
101 |
|
|
for m=1:12; |
102 |
|
|
if size(xc,3)==1 |
103 |
|
|
sm(:,:,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
104 |
|
|
sm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
105 |
|
|
else |
106 |
|
|
for t=1:size(xc,3); |
107 |
|
|
sm(:,:,t,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t))); |
108 |
|
|
sm2(:,:,t,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc(:,:,t),yc(:,:,t),zc(kk),squeeze(GRID.mskc(:,:,kk,t))); |
109 |
|
|
%sm(:,:,:,m)=interp3ddata(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
110 |
|
|
%sm2(:,:,:,m)=interp3ddata_int2(xl,yl,zlm,SM(:,:,:,m),xc,yc,zc(kk),GRID.mskc(:,:,kk)); |
111 |
|
|
end |
112 |
|
|
end |
113 |
|
|
end |
114 |
|
|
%wrda('sss.bin',sm2(:,:,1,:),1,fmt,Ieee); |
115 |
|
|
save sm sm sm2 |
116 |
|
|
|
117 |
|
|
% Levitus monthly grid |
118 |
|
|
LMGRID=LCGRID; |
119 |
|
|
LMGRID.mskc=mask(TM(:,:,:,1)); |
120 |
|
|
MM=mean(TM,4); |
121 |
|
|
[TMm,TMsd] = meanprofile(LMGRID,MM); |
122 |
|
|
clear TM MM |
123 |
|
|
MGRID=GRID; |
124 |
|
|
MGRID.mskc=MGRID.mskc(:,:,1:size(tm,3)); |
125 |
|
|
[tmm,tmsd] = meanprofile(MGRID,mean(tm,4)); |
126 |
|
|
[tm2m,tm2sd] = meanprofile(MGRID,mean(tm2,4)); |
127 |
|
|
|
128 |
|
|
MM=mean(SM,4); |
129 |
|
|
[SMm,SMsd] = meanprofile(LMGRID,MM); |
130 |
|
|
clear SM MM |
131 |
|
|
[smm,smsd] = meanprofile(MGRID,mean(sm,4)); |
132 |
|
|
[sm2m,sm2sd] = meanprofile(MGRID,mean(sm2,4)); |
133 |
|
|
|
134 |
|
|
km=size(tmm,2); |
135 |
|
|
end |
136 |
|
|
|
137 |
|
|
load tc |
138 |
|
|
[tcm,tcsd] = meanprofile(GRID,tc); |
139 |
|
|
[tc2m,tc2sd] = meanprofile(GRID,tc2); |
140 |
|
|
load sc |
141 |
|
|
[scm,scsd] = meanprofile(GRID,sc); |
142 |
|
|
[sc2m,sc2sd] = meanprofile(GRID,sc2); |
143 |
|
|
|
144 |
|
|
% $$$ figure(1) |
145 |
|
|
% $$$ subplot(121) |
146 |
|
|
% $$$ plot(TCm,zlc,'k',... |
147 |
|
|
% $$$ TMm,zlm,'r',... |
148 |
|
|
% $$$ tcm,zc,'g',... |
149 |
|
|
% $$$ tc2m,zc,'b',... |
150 |
|
|
% $$$ tmm,zc(1:km),'c',... |
151 |
|
|
% $$$ tm2m,zc(1:km),'m',... |
152 |
|
|
% $$$ [TCm'-TCsd' TCm'+TCsd'],zlc,'k--',... |
153 |
|
|
% $$$ [TMm'-TMsd' TMm'+TMsd'],zlm,'r--',... |
154 |
|
|
% $$$ [tcm'-tcsd' tcm'+tcsd'],zc,'g--',... |
155 |
|
|
% $$$ [tc2m'-tc2sd' tc2m'+tc2sd'],zc,'b--') |
156 |
|
|
% $$$ xlabel('\theta (^oC)');ylabel('Depth (m)'); |
157 |
|
|
% $$$ legend('Levitus Clim.','Levitus Monthly','Interpolated','Integrated',4) |
158 |
|
|
% $$$ |
159 |
|
|
% $$$ subplot(122) |
160 |
|
|
% $$$ plot(SCm,zlc,'k',... |
161 |
|
|
% $$$ SMm,zlm,'r',... |
162 |
|
|
% $$$ scm,zc,'g',... |
163 |
|
|
% $$$ sc2m,zc,'b',... |
164 |
|
|
% $$$ smm,zc(1:km),'c',... |
165 |
|
|
% $$$ sm2m,zc(1:km),'m',... |
166 |
|
|
% $$$ [SCm'-SCsd' SCm'+SCsd'],zlc,'k--',... |
167 |
|
|
% $$$ [SMm'-SMsd' SMm'+SMsd'],zlm,'r--',... |
168 |
|
|
% $$$ [scm'-scsd' scm'+scsd'],zc,'g--',... |
169 |
|
|
% $$$ [sc2m'-sc2sd' sc2m'+sc2sd'],zc,'b--') |
170 |
|
|
% $$$ xlabel('S (psu)');ylabel('Depth (m)'); |
171 |
|
|
% $$$ legend('Levitus Clim.','Levitus Monthly','Interpolated','Integrated',4) |
172 |
|
|
% $$$ |
173 |
|
|
% $$$ supertext(.5,1,... |
174 |
|
|
% $$$ 'Levitus climatology, monthly and interpolated datasets',... |
175 |
|
|
% $$$ 'FontSize',14,... |
176 |
|
|
% $$$ 'VerticalAlignment','Bottom',... |
177 |
|
|
% $$$ 'HorizontalAlignment','Center') |
178 |
|
|
% $$$ print -depsc2 comparison_profiles.eps |
179 |
|
|
% $$$ print -djpeg100 -r90 comparison_profiles.jpg |
180 |
|
|
% $$$ |
181 |
|
|
% $$$ figure(2) |
182 |
|
|
% $$$ subplot(121) |
183 |
|
|
% $$$ plot(TMm-TCm(1:19),zlm,'k',tc2m-tcm,zc,'r') |
184 |
|
|
% $$$ xlabel('\Delta \theta (^oC)');ylabel('Depth (m)'); |
185 |
|
|
% $$$ legend('Clim. - Monthly','Integr. - Intep.',4); |
186 |
|
|
% $$$ subplot(122) |
187 |
|
|
% $$$ plot(SMm-SCm(1:19),zlm,'k',sc2m-scm,zc,'r') |
188 |
|
|
% $$$ xlabel('\Delta S (psu)');ylabel('Depth (m)'); |
189 |
|
|
% $$$ legend('Clim. - Monthly','Integr. - Intep.',4); |
190 |
|
|
% $$$ |
191 |
|
|
% $$$ supertext(.5,1,... |
192 |
|
|
% $$$ 'Difference between datasets',... |
193 |
|
|
% $$$ 'FontSize',14,... |
194 |
|
|
% $$$ 'VerticalAlignment','Bottom',... |
195 |
|
|
% $$$ 'HorizontalAlignment','Center') |
196 |
|
|
% $$$ print -depsc2 levitus_problems.eps |
197 |
|
|
% $$$ print -djpeg100 -r90 levitus_problems.jpg |
198 |
|
|
% $$$ |
199 |
|
|
% $$$ |
200 |
|
|
% $$$ % Test EOS |
201 |
|
|
% $$$ P3=ini_poly3; |
202 |
|
|
% $$$ load tc tc2;T=tc2;clear tc2 |
203 |
|
|
% $$$ %T=rdda('temp.bin',[90 40 15],1,fmt,Ieee); |
204 |
|
|
% $$$ load sc sc2;S=sc2;clear sc2 |
205 |
|
|
% $$$ %S=rdda('salt.bin',[90 40 15],1,fmt,Ieee); |
206 |
|
|
% $$$ D=dens_poly3(P3,T,S); |
207 |
|
|
% $$$ |
208 |
|
|
% $$$ z=ones(90*40,1)*zc; |
209 |
|
|
% $$$ y=ones(90,1)*yc;y=y(:); |
210 |
|
|
% $$$ p=sw_pres(-z',y')'; |
211 |
|
|
% $$$ p=reshape(p,[90 40 15]); |
212 |
|
|
% $$$ clear z y |
213 |
|
|
% $$$ d=(sw_dens(S,T,p)-1000).*mask(T); |
214 |
|
|
% $$$ Dd=D-d; |
215 |
|
|
% $$$ |
216 |
|
|
% $$$ [Dm,Dsd]=meanprofile(GRID,D); |
217 |
|
|
% $$$ [dm,dsd]=meanprofile(GRID,d); |
218 |
|
|
% $$$ [Ddm,Ddsd]=meanprofile(GRID,D,d); |
219 |
|
|
% $$$ |
220 |
|
|
% $$$ figure(3) |
221 |
|
|
% $$$ subplot(121) |
222 |
|
|
% $$$ plot(dm,zc,'k',... |
223 |
|
|
% $$$ Dm,zc,'g',... |
224 |
|
|
% $$$ [dm'-dsd' dm'+dsd'],zc,'k--',... |
225 |
|
|
% $$$ [Dm'-Dsd' Dm'+Dsd'],zc,'g--') |
226 |
|
|
% $$$ axis([22 51 -5000 0]) |
227 |
|
|
% $$$ xlabel('\rho-1000 (kg/m^3)');ylabel('Depth (m)'); |
228 |
|
|
% $$$ legend('Seawater EOS','POLY3',3) |
229 |
|
|
% $$$ |
230 |
|
|
% $$$ subplot(122) |
231 |
|
|
% $$$ plot(Ddm,zc,'k',... |
232 |
|
|
% $$$ [Ddm'-dsd' Ddm'+dsd'],zc,'g',... |
233 |
|
|
% $$$ [Ddm'-Ddsd' Ddm'+Ddsd'],zc,'k--') |
234 |
|
|
% $$$ axis([[-1 1]*.5 -5000 0]) |
235 |
|
|
% $$$ xlabel('\Delta \rho (kg/m^3)');ylabel('Depth (m)'); |
236 |
|
|
% $$$ title('POLY3-EOS80') |
237 |
|
|
% $$$ legend('Mean difference','Observed SD of \rho','SD of P3-EOS80',4) |
238 |
|
|
% $$$ print -depsc2 eos_problems.eps |
239 |
|
|
% $$$ print -djpeg100 -r90 eos_problems.jpg |
240 |
|
|
% $$$ |
241 |
|
|
% $$$ figure(4) |
242 |
|
|
% $$$ subplot(221) |
243 |
|
|
% $$$ pcol(xc,yc, sq(D(:,:,8))' ); |
244 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
245 |
|
|
% $$$ colorbar('h'); |
246 |
|
|
% $$$ title('\rho @ z=-1250m'); |
247 |
|
|
% $$$ |
248 |
|
|
% $$$ subplot(222) |
249 |
|
|
% $$$ pcol(xc,yc, sq(Dd(:,:,8))' ); |
250 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
251 |
|
|
% $$$ %cax=caxis;caxis([-1 1]*max(abs(cax))); |
252 |
|
|
% $$$ colorbar('h'); |
253 |
|
|
% $$$ title('\Delta \rho (POLY3-EOS80) @ z=-1250m'); |
254 |
|
|
% $$$ |
255 |
|
|
% $$$ subplot(223) |
256 |
|
|
% $$$ pcol(xc,yc, sq(D(:,:,12))' ); |
257 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
258 |
|
|
% $$$ colorbar('h'); |
259 |
|
|
% $$$ title('\rho @ z=-3010m'); |
260 |
|
|
% $$$ |
261 |
|
|
% $$$ subplot(224) |
262 |
|
|
% $$$ pcol(xc,yc, sq(Dd(:,:,12))' ); |
263 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
264 |
|
|
% $$$ %cax=caxis;caxis([-1 1]*max(abs(cax))); |
265 |
|
|
% $$$ colorbar('h'); |
266 |
|
|
% $$$ title('\Delta \rho (POLY3-ESO80) @ z=-3010m'); |
267 |
|
|
% $$$ print -depsc2 eos_problems2.eps |
268 |
|
|
% $$$ print -djpeg100 -r90 eos_problems2.jpg |
269 |
|
|
% $$$ |
270 |
|
|
% $$$ figure(5) |
271 |
|
|
% $$$ subplot(221) |
272 |
|
|
% $$$ pcol(xc,yc, sq(mean(tm2(:,:,1,:)-tm(:,:,1,:),4))') |
273 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
274 |
|
|
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h') |
275 |
|
|
% $$$ title('\Delta\theta (Integrated-Interpolated) k=1') |
276 |
|
|
% $$$ |
277 |
|
|
% $$$ subplot(222) |
278 |
|
|
% $$$ pcol(xc,yc, sq(mean(sm2(:,:,1,:)-sm(:,:,1,:),4))') |
279 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
280 |
|
|
% $$$ caxis([-1 1]*0.25) |
281 |
|
|
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h') |
282 |
|
|
% $$$ title('{\Delta}S (Integrated-Interpolated) k=1') |
283 |
|
|
% $$$ |
284 |
|
|
% $$$ subplot(223) |
285 |
|
|
% $$$ pcol(xc,yc, sq(mean(tc2(:,:,10,:)-tc(:,:,10,:),4))') |
286 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
287 |
|
|
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h') |
288 |
|
|
% $$$ title('\Delta\theta (Integrated-Interpolated) k=10') |
289 |
|
|
% $$$ |
290 |
|
|
% $$$ subplot(224) |
291 |
|
|
% $$$ pcol(xc,yc, sq(mean(sc2(:,:,10,:)-sc(:,:,10,:),4))') |
292 |
|
|
% $$$ ylabel('Latitude (^oN)'); |
293 |
|
|
% $$$ caxis([-1 1]*max(abs(caxis)));colorbar('h') |
294 |
|
|
% $$$ title('{\Delta}S (Integrated-Interpolated) k=10') |
295 |
|
|
% $$$ |
296 |
|
|
% $$$ print -depsc2 intgr-interp.eps |
297 |
|
|
% $$$ print -djpeg100 -r90 intgr-interp.jpg |