| 1 |
function []=cost_sst(dirModel,dirMat,doComp,dirTex,nameTex); |
| 2 |
%object: compute cost function term for sst 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 |
if doComp; |
| 26 |
|
| 27 |
%grid, params and inputs |
| 28 |
|
| 29 |
gcmfaces_global; global myparms; |
| 30 |
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end; |
| 31 |
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end; |
| 32 |
if isfield(myparms,'yearFirst'); yearFirst=myparms.yearFirst; yearLast=myparms.yearLast; |
| 33 |
else; yearFirst=1992; yearLast=2011; |
| 34 |
end; |
| 35 |
|
| 36 |
%search for nctiles files |
| 37 |
useNctiles=0; |
| 38 |
dirNctiles=[dirModel 'nctiles_remotesensing/sst/']; |
| 39 |
if ~isempty(dir(dirNctiles)); useNctiles=1; end; |
| 40 |
|
| 41 |
nameReynolds='reynolds_oiv2_r1'; |
| 42 |
nameRemss='tmi_amsre_oisst_r1'; |
| 43 |
nameL2p='tmi_amsre_l2p_r1'; |
| 44 |
|
| 45 |
dirData=dirModel; subdirReynolds=''; subdirRemss=''; subdirL2p=''; |
| 46 |
if ~isempty(dir([dirModel 'inputfiles/' nameReynolds '*'])); |
| 47 |
dirData=[dirModel 'inputfiles/']; subdirReynolds=''; subdirRemss=''; subdirL2p=''; |
| 48 |
end; |
| 49 |
|
| 50 |
nameSigma='sigma_surf_0p5.bin'; |
| 51 |
dirSigma=dirModel; |
| 52 |
if ~isempty(dir([dirModel 'inputfiles/' nameSigma])); |
| 53 |
dirSigma=[dirModel 'inputfiles/']; |
| 54 |
end; |
| 55 |
|
| 56 |
fld_err=NaN*mygrid.RAC; |
| 57 |
file0=[dirSigma nameSigma]; |
| 58 |
if ~isempty(dir(file0)); fld_err=read_bin(file0,1,0); end; |
| 59 |
if useNctiles; fld_err=read_nctiles([dirNctiles 'sst'],'sst_sigma'); end; |
| 60 |
fld_w=fld_err.^-2; |
| 61 |
|
| 62 |
nameTbar='m_sst_mon'; |
| 63 |
if ~isempty(dir([dirModel nameTbar '*'])); |
| 64 |
dirTbar=dirModel; |
| 65 |
else; |
| 66 |
dirTbar=[dirModel 'barfiles' filesep]; |
| 67 |
end; |
| 68 |
fileModel='unknown'; |
| 69 |
file0=[dirTbar nameTbar '*data']; |
| 70 |
if ~isempty(dir(file0)); |
| 71 |
fileModel=dir(file0); |
| 72 |
fileModel=fileModel.name; |
| 73 |
end; |
| 74 |
|
| 75 |
if ~useNctiles; |
| 76 |
%computational loop |
| 77 |
sst_mod=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]); |
| 78 |
sst_rey=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]); |
| 79 |
sst_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]); |
| 80 |
% |
| 81 |
for ycur=yearFirst:yearLast; |
| 82 |
fprintf(['starting ' num2str(ycur) '\n']); |
| 83 |
tic; |
| 84 |
for mcur=1:12; |
| 85 |
|
| 86 |
mm=(ycur-yearFirst)*12+mcur; |
| 87 |
|
| 88 |
%load Reynolds SST |
| 89 |
fld_rey=NaN*mygrid.RAC; |
| 90 |
file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)]; |
| 91 |
if ~isempty(dir(file0)); fld_rey=read_bin(file0,mcur,0); end; |
| 92 |
fld_rey(find(fld_rey==0))=NaN; |
| 93 |
fld_rey(find(fld_rey<-99))=NaN; |
| 94 |
sst_rey(:,:,mm)=fld_rey; |
| 95 |
|
| 96 |
%load Remss SST |
| 97 |
fld_remss=NaN*mygrid.RAC; |
| 98 |
file0=[dirData subdirRemss nameRemss '_' num2str(ycur)]; |
| 99 |
if ~isempty(dir(file0)); fld_remss=read_bin(file0,mcur,0); end; |
| 100 |
fld_remss(find(fld_remss==0))=NaN; |
| 101 |
fld_remss(find(fld_remss<-99))=NaN; |
| 102 |
sst_remss(:,:,mm)=fld_remss; |
| 103 |
|
| 104 |
%load model SST |
| 105 |
fld_mod=NaN*mygrid.RAC; |
| 106 |
file0=[dirTbar fileModel]; |
| 107 |
if ~isempty(dir(file0)); fld_mod=read_bin(file0,mm,0); end; |
| 108 |
fld_mod=fld_mod.*mygrid.mskC(:,:,1); |
| 109 |
sst_mod(:,:,mm)=fld_mod; |
| 110 |
end; |
| 111 |
toc; |
| 112 |
end; |
| 113 |
end;%if ~useNctiles; |
| 114 |
|
| 115 |
if useNctiles; |
| 116 |
sst_rey=read_nctiles([dirNctiles 'sst'],'sst_REYNOLDS'); |
| 117 |
sst_mod=read_nctiles([dirNctiles 'sst'],'sst_ECCOv4r2'); |
| 118 |
sst_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]); |
| 119 |
end; |
| 120 |
|
| 121 |
%store misfit maps |
| 122 |
mod_m_rey=sst_mod-sst_rey; |
| 123 |
mod_m_remss=sst_mod-sst_remss; |
| 124 |
|
| 125 |
%comp zonal means |
| 126 |
tic; |
| 127 |
zm_mod=calc_zonmean_T(sst_mod); |
| 128 |
zm_rey=calc_zonmean_T(sst_rey); |
| 129 |
zm_remss=calc_zonmean_T(sst_remss); |
| 130 |
zm_mod_m_rey=calc_zonmean_T(mod_m_rey); |
| 131 |
zm_mod_m_remss=calc_zonmean_T(mod_m_remss); |
| 132 |
toc; |
| 133 |
|
| 134 |
%compute rms maps |
| 135 |
rms_to_rey=sqrt(nanmean(mod_m_rey.^2,3)); |
| 136 |
rms_to_remss=sqrt(nanmean(mod_m_remss.^2,3)); |
| 137 |
|
| 138 |
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end; |
| 139 |
eval(['save ' dirMat '/cost/cost_sst.mat fld_err rms_* zm_*;']); |
| 140 |
|
| 141 |
else;%display previously computed results |
| 142 |
|
| 143 |
global mygrid; |
| 144 |
|
| 145 |
if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end; |
| 146 |
|
| 147 |
eval(['load ' dirMat '/cost_sst.mat;']); |
| 148 |
|
| 149 |
if ~isempty(who('fld_rms')); |
| 150 |
figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2}); |
| 151 |
myCaption={'modeled-observed rms -- sea surface temperature (K)'}; |
| 152 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 153 |
else; |
| 154 |
|
| 155 |
testREMSS=sum(~isnan(rms_to_remss))>0; |
| 156 |
|
| 157 |
figure; m_map_gcmfaces(rms_to_rey,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2}); |
| 158 |
myCaption={'modeled-Reynolds rms -- sea surface temperature (K)'}; |
| 159 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 160 |
|
| 161 |
if testREMSS; |
| 162 |
figure; m_map_gcmfaces(rms_to_remss,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2}); |
| 163 |
myCaption={'modeled-REMSS rms -- sea surface temperature (K)'}; |
| 164 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 165 |
end; |
| 166 |
|
| 167 |
ny=size(zm_rey,2)/12; |
| 168 |
[xx,yy]=meshgrid(1992+([1:ny*12]-0.5)/12,mygrid.LATS); |
| 169 |
|
| 170 |
figureL; |
| 171 |
obs=zm_rey; obsCycle=sst_cycle(obs); |
| 172 |
mod=zm_rey+zm_mod_m_rey; modCycle=sst_cycle(mod); |
| 173 |
mis=zm_mod_m_rey; misCycle=sst_cycle(mis); |
| 174 |
subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 175 |
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude'); |
| 176 |
title('Reynolds sst anomaly'); |
| 177 |
subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 178 |
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude'); |
| 179 |
title('ECCO sst anomaly'); |
| 180 |
subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 181 |
set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-Reynolds sst misfit'); |
| 182 |
myCaption={'ECCO and Reynolds zonal mean sst anomalies (K)'}; |
| 183 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 184 |
|
| 185 |
if testREMSS; |
| 186 |
figureL; |
| 187 |
obs=zm_remss; obsCycle=sst_cycle(obs); |
| 188 |
mod=zm_remss+zm_mod_m_remss; modCycle=sst_cycle(mod); |
| 189 |
mis=zm_mod_m_remss; misCycle=sst_cycle(mis); |
| 190 |
subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 191 |
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude'); |
| 192 |
title('REMSS sst anomaly'); |
| 193 |
subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 194 |
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude'); |
| 195 |
title('ECCO sst anomaly'); |
| 196 |
subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar; |
| 197 |
set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-REMSS sst misfit'); |
| 198 |
myCaption={'ECCO and REMSS zonal mean sst anomalies (K)'}; |
| 199 |
if addToTex; write2tex(fileTex,2,myCaption,gcf); end; |
| 200 |
end; |
| 201 |
|
| 202 |
end; |
| 203 |
|
| 204 |
end; |
| 205 |
|
| 206 |
function [zmCycle]=sst_cycle(zmIn); |
| 207 |
|
| 208 |
ny=size(zmIn,2)/12; |
| 209 |
zmCycle=NaN*zeros(179,12); |
| 210 |
for mm=1:12; |
| 211 |
zmCycle(:,mm)=nanmean(zmIn(:,mm:12:ny*12),2); |
| 212 |
end; |
| 213 |
zmCycle=repmat(zmCycle,[1 ny]); |
| 214 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph