| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph