| 1 |
function []=cost_seaicearea(dirModel,dirMat,doComp,dirTex,nameTex); |
| 2 |
%object: compute cost function term for sea ice data |
| 3 |
%inputs: dimodel is the model directory |
| 4 |
% dirMat is the directory where diagnozed .mat files will be saved |
| 5 |
% -> set it to '' to use the default [dirModel 'mat/'] |
| 6 |
% doComp is a switch (1->compute; 0->display) |
| 7 |
%optional: dirTex is the directory where tex and figures files are created |
| 8 |
% (if not specified then display all results to screen instead) |
| 9 |
% nameTex is the tex file name (default : 'myPlots') |
| 10 |
|
| 11 |
if isempty(dirMat); dirMat=[dirModel 'mat' filesep]; else; dirMat=[dirMat filesep]; end; |
| 12 |
if isempty(dir(dirMat)); mkdir([dirMat]); end; |
| 13 |
|
| 14 |
%determine if and where to create tex and figures files |
| 15 |
dirMat=[dirMat filesep]; |
| 16 |
if isempty(who('dirTex')); |
| 17 |
addToTex=0; |
| 18 |
else; |
| 19 |
if ~ischar(dirTex); error('mis-specified dirTex'); end; |
| 20 |
addToTex=1; |
| 21 |
if isempty(who('nameTex')); nameTex='myPlots'; end; |
| 22 |
fileTex=[dirTex filesep nameTex '.tex']; |
| 23 |
end; |
| 24 |
|
| 25 |
%grid, params and inputs |
| 26 |
gcmfaces_global; global myparms; |
| 27 |
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end; |
| 28 |
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end; |
| 29 |
if isfield(myparms,'yearsta'); yearsta=myparms.yearsta; yearend=myparms.yearend; |
| 30 |
else; yearsta=1992; yearend=2011; |
| 31 |
end; |
| 32 |
|
| 33 |
[lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic; |
| 34 |
lonLims=[lonPairs(1:5,1);lonPairs(1,1)]; |
| 35 |
|
| 36 |
%integrals : |
| 37 |
nlon=(length(lonLims)-1); |
| 38 |
areamskN=NaN*repmat(mygrid.RAC,[1 1 nlon+1]); |
| 39 |
areamskS=NaN*repmat(mygrid.RAC,[1 1 nlon+1]); |
| 40 |
areamskN(:,:,1)=mygrid.RAC.*(mygrid.YC>0); |
| 41 |
areamskS(:,:,1)=mygrid.RAC.*(mygrid.YC<0); |
| 42 |
for kk=1:nlon; |
| 43 |
tmpmsk=0.*mygrid.XC; |
| 44 |
if lonLims(kk+1) > lonLims(kk) |
| 45 |
tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.; |
| 46 |
else |
| 47 |
tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.; |
| 48 |
tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.; |
| 49 |
end |
| 50 |
areamskN(:,:,1+kk)=mygrid.RAC.*tmpmsk.*(mygrid.YC>0); |
| 51 |
areamskS(:,:,1+kk)=mygrid.RAC.*tmpmsk.*(mygrid.YC<0); |
| 52 |
end; |
| 53 |
|
| 54 |
%search for nctiles files |
| 55 |
useNctiles=0; |
| 56 |
dirNctiles=[dirModel 'nctiles_remotesensing/sic/']; |
| 57 |
if ~isempty(dir(dirNctiles)); useNctiles=1; end; |
| 58 |
|
| 59 |
if doComp; |
| 60 |
%grid, params and inputs |
| 61 |
fld_err=0.5*ones(90,1170); fld_err=convert2gcmfaces(fld_err); |
| 62 |
if useNctiles; fld_err=read_nctiles([dirNctiles 'sic'],'sic_sigma'); end; |
| 63 |
fld_w=fld_err.^-2; |
| 64 |
|
| 65 |
dirEcco=dirModel; |
| 66 |
if ~isempty(dir([dirModel 'barfiles'])); |
| 67 |
dirEcco=[dirModel 'barfiles' filesep]; |
| 68 |
end; |
| 69 |
|
| 70 |
nameData='nsidc79_monthly_'; |
| 71 |
dirData=dirModel; |
| 72 |
if ~isempty(dir([dirModel 'inputfiles/' nameData '*'])); |
| 73 |
dirData=[dirModel 'inputfiles/']; |
| 74 |
end; |
| 75 |
|
| 76 |
fileModel='unknown'; |
| 77 |
file0=[dirEcco 'm_siv4_area*data']; |
| 78 |
if ~isempty(dir(file0)); fileModel=dir(file0); fileModel=fileModel.name; end; |
| 79 |
|
| 80 |
nyears=yearend-yearsta+1; |
| 81 |
nmonths=12*nyears; |
| 82 |
|
| 83 |
if ~useNctiles; |
| 84 |
%misfits : |
| 85 |
fld_estim=convert2gcmfaces(NaN*ones(90,90*13,nmonths)); |
| 86 |
fld_nsidc=convert2gcmfaces(NaN*ones(90,90*13,nmonths)); |
| 87 |
fld_dif=convert2gcmfaces(NaN*ones(90,90*13,nmonths)); |
| 88 |
|
| 89 |
%computational loop : |
| 90 |
for ycur=yearsta:yearend; |
| 91 |
tic; |
| 92 |
for mcur=1:12; |
| 93 |
mm=(ycur-yearsta)*12+mcur; |
| 94 |
|
| 95 |
fld_dat=NaN*mygrid.RAC; |
| 96 |
file0=[dirData nameData num2str(ycur)]; |
| 97 |
if ~isempty(dir(file0)); fld_dat=read_bin(file0,mcur,0); end; |
| 98 |
fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask |
| 99 |
fld_dat(find(fld_dat<-99))=NaN;%missing data |
| 100 |
msk=1+0*fld_dat;%combined mask |
| 101 |
|
| 102 |
fld_mod=NaN*mygrid.RAC; |
| 103 |
file0=[dirEcco fileModel]; |
| 104 |
if ~isempty(dir(file0)); fld_mod=read_bin(file0,mm,0); end; |
| 105 |
fld_mod=fld_mod.*msk;%mask consistent with fld_dat |
| 106 |
|
| 107 |
%misfits : |
| 108 |
fld_estim(:,:,mm)=fld_mod; |
| 109 |
fld_nsidc(:,:,mm)=fld_dat; |
| 110 |
fld_dif(:,:,mm)=fld_mod-fld_dat; |
| 111 |
end; |
| 112 |
toc; |
| 113 |
end; |
| 114 |
end;%if ~useNctiles; |
| 115 |
|
| 116 |
if useNctiles; |
| 117 |
fld_estim=read_nctiles([dirNctiles 'sic'],'sic_ECCOv4r2'); |
| 118 |
fld_nsidc=read_nctiles([dirNctiles 'sic'],'sic_NSIDC'); |
| 119 |
fld_estim(isnan(fld_nsidc))=NaN; |
| 120 |
fld_nsidc(isnan(fld_estim))=NaN; |
| 121 |
fld_dif=fld_estim-fld_nsidc; |
| 122 |
end; |
| 123 |
|
| 124 |
%monthly mean climatology : |
| 125 |
climMod=reshape(fld_estim,[1 1 12 nyears]); |
| 126 |
climMod=nanmean(climMod,4); |
| 127 |
climObs=reshape(fld_nsidc,[1 1 12 nyears]); |
| 128 |
climObs=nanmean(climObs,4); |
| 129 |
|
| 130 |
%monthly integrals : |
| 131 |
IceAreaNorthMod=NaN*zeros(1+nlon,nmonths); |
| 132 |
IceAreaNorthObs=NaN*zeros(1+nlon,nmonths); |
| 133 |
IceAreaSouthMod=NaN*zeros(1+nlon,nmonths); |
| 134 |
IceAreaSouthObs=NaN*zeros(1+nlon,nmonths); |
| 135 |
for mm=1:nlon+1; |
| 136 |
tmp1=fld_estim.*repmat(areamskN(:,:,mm),[1 1 nmonths]); |
| 137 |
IceAreaNorthMod(mm,:)=nansum(tmp1,0); |
| 138 |
tmp1=fld_nsidc.*repmat(areamskN(:,:,mm),[1 1 nmonths]); |
| 139 |
IceAreaNorthObs(mm,:)=nansum(tmp1,0); |
| 140 |
% |
| 141 |
tmp1=fld_estim.*repmat(areamskS(:,:,mm),[1 1 nmonths]); |
| 142 |
IceAreaSouthMod(mm,:)=nansum(tmp1,0); |
| 143 |
tmp1=fld_nsidc.*repmat(areamskS(:,:,mm),[1 1 nmonths]); |
| 144 |
IceAreaSouthObs(mm,:)=nansum(tmp1,0); |
| 145 |
end; |
| 146 |
|
| 147 |
%misfits : |
| 148 |
mis_rms=sqrt(nanmean(fld_dif.^2,3)); |
| 149 |
obs_std=nanstd(fld_nsidc,[],3); |
| 150 |
mod_std=nanstd(fld_nsidc+fld_dif,[],3); |
| 151 |
|
| 152 |
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end; |
| 153 |
eval(['save ' dirMat '/cost/cost_seaicearea.mat fld_err mis_rms obs_std mod_std IceArea* clim*;']); |
| 154 |
|
| 155 |
else;%display previously computed results |
| 156 |
|
| 157 |
if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end; |
| 158 |
|
| 159 |
eval(['load ' dirMat '/cost_seaicearea.mat;']); |
| 160 |
|
| 161 |
%variance maps: |
| 162 |
figure; m_map_gcmfaces(mis_rms,0,{'myCaxis',[0:0.1:1.]}); |
| 163 |
myCaption={'modeled-observed rms -- sea ice concentration'}; |
| 164 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 165 |
|
| 166 |
figure; m_map_gcmfaces(obs_std,0,{'myCaxis',[0:0.1:1.]}); |
| 167 |
myCaption={'observed std -- sea ice concentration'}; |
| 168 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 169 |
|
| 170 |
figure; m_map_gcmfaces(mod_std,0,{'myCaxis',[0:0.1:1.]}); |
| 171 |
myCaption={'modelled std -- sea ice concentration'}; |
| 172 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 173 |
|
| 174 |
%arctic maps |
| 175 |
figureL; |
| 176 |
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),2,{'myCaxis',[0:0.1:1.]}); |
| 177 |
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),2,{'myCaxis',[0:0.1:1.]}); |
| 178 |
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),2,{'myCaxis',[0:0.1:1.]}); |
| 179 |
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),2,{'myCaxis',[0:0.1:1.]}); |
| 180 |
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'}; |
| 181 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 182 |
|
| 183 |
%southern ocean maps |
| 184 |
figureL; |
| 185 |
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),3,{'myCaxis',[0:0.1:1.]}); |
| 186 |
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),3,{'myCaxis',[0:0.1:1.]}); |
| 187 |
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),3,{'myCaxis',[0:0.1:1.]}); |
| 188 |
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),3,{'myCaxis',[0:0.1:1.]}); |
| 189 |
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'}; |
| 190 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 191 |
|
| 192 |
%northern/southern integrals |
| 193 |
ny=length(IceAreaNorthMod)/12; |
| 194 |
yy=yearsta+[0:ny-1]; |
| 195 |
|
| 196 |
figureL; |
| 197 |
subplot(1,2,1); |
| 198 |
plot(yy,IceAreaNorthMod(1,3:12:end),'LineWidth',2); hold on; |
| 199 |
plot(yy,IceAreaNorthObs(1,3:12:end),'r','LineWidth',2); |
| 200 |
plot(yy,IceAreaNorthMod(1,9:12:end),'LineWidth',2); |
| 201 |
plot(yy,IceAreaNorthObs(1,9:12:end),'r','LineWidth',2); |
| 202 |
axis([yearsta yearsta+ny-1 0 20e12]); |
| 203 |
ylabel('m^2'); title('Northern Hemisphere'); |
| 204 |
subplot(1,2,2); |
| 205 |
plot(yy,IceAreaSouthMod(1,3:12:end),'LineWidth',2); hold on; |
| 206 |
plot(yy,IceAreaSouthObs(1,3:12:end),'r','LineWidth',2); |
| 207 |
plot(yy,IceAreaSouthMod(1,9:12:end),'LineWidth',2); |
| 208 |
plot(yy,IceAreaSouthObs(1,9:12:end),'r','LineWidth',2); |
| 209 |
axis([yearsta yearsta+ny-1 0 20e12]); |
| 210 |
ylabel('m^2'); title('Southern Hemisphere'); |
| 211 |
|
| 212 |
myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in March and September',... |
| 213 |
'in Northern Hemisphere (left) and Southern Hemisphere (right)'}; |
| 214 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 215 |
|
| 216 |
%southern basin integrals |
| 217 |
|
| 218 |
for mm=[3 9]; |
| 219 |
figureL; |
| 220 |
for ii=1:6; |
| 221 |
subplot(3,2,ii); |
| 222 |
if ii==1; iiTxt='Entire Southern Ocean'; |
| 223 |
else; iiTxt=sprintf('%dE to %dE',lonLims(ii-1),lonLims(ii)); |
| 224 |
end; |
| 225 |
plot(yy,IceAreaSouthMod(ii,mm:12:end),'LineWidth',2); hold on; |
| 226 |
plot(yy,IceAreaSouthObs(ii,mm:12:end),'r','LineWidth',2); |
| 227 |
aa=axis; aa(1:2)=[yearsta yearsta+ny-1]; axis(aa); |
| 228 |
ylabel('m^2'); title(iiTxt); |
| 229 |
end; |
| 230 |
if mm==3; mmTxt='March'; elseif mm==9; mmTxt='September'; else; '???'; end; |
| 231 |
myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in ',mmTxt,' per Southern Ocean sector'}; |
| 232 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 233 |
end; |
| 234 |
|
| 235 |
end; |
| 236 |
|
| 237 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph