gcmfaces/ 0000755 0004354 0000144 00000000000 12375431052 012207 5 ustar gforget users gcmfaces/gcmfaces_misc/ 0000755 0004354 0000144 00000000000 12375431055 014775 5 ustar gforget users gcmfaces/gcmfaces_misc/ccaa.m 0000644 0004354 0000144 00000000076 11600706727 016046 0 ustar gforget users %abbreviation for clear all; close all;
clear all; close all;
gcmfaces/gcmfaces_misc/convertR8toR4.m 0000644 0004354 0000144 00000000464 11600706727 017623 0 ustar gforget users function []=convertR8toR4(fileIn,fileOut);
%object: read R8 fileIn, write R4 fileOut
%inputs: fileIn and fileOut
tmp1=dir(fileIn);
tmp1=tmp1.bytes/8;
fid=fopen(fileIn,'r','b'); tmp2=fread(fid,tmp1,'float64'); fclose(fid);
fid=fopen(fileOut,'w','b'); fwrite(fid,tmp2,'float32'); fclose(fid);
gcmfaces/gcmfaces_misc/depthStretch.m 0000644 0004354 0000144 00000002041 11600706727 017612 0 ustar gforget users function [z]=depthStretch(d,varargin)
%object: compute stretched depth coordinate
%inputs: d is the depth vector (>=0)
%optional: depthStretchLim are the stretching depth limits ([0 500 6000] by default)
%
%this routine maps
%[depthStretchLim(1) depthStretchLim(2)] to [1 0]
%[depthStretchLim(2) depthStretchLim(3)] to [0 -1]
if nargin>1; depthStretchLim=varargin{1}; else; depthStretchLim=[0 500 6000]; end;
zStretchLim=[1 0 -1];
%make sure depth is positive
d=abs(d);
%initialize stretched vertical coordinate
z=zeros(size(d));
%values between 0 and depthStretchLim(2) get half the range of z
tmp1=find(d<=depthStretchLim(2));
tmp2=(depthStretchLim(2)-d(tmp1))./(depthStretchLim(2)-depthStretchLim(1))...
.*(zStretchLim(1)-zStretchLim(2))+zStretchLim(2);
z(tmp1)=tmp2;
%values between depthStretchLim(2) and depthStretchLim(3) get the other half z range
tmp1=find(d>depthStretchLim(2));
tmp2=(d(tmp1)-depthStretchLim(2))./(depthStretchLim(3)-depthStretchLim(2))...
.*(zStretchLim(3)-zStretchLim(2))+zStretchLim(2);
z(tmp1)=tmp2;
gcmfaces/gcmfaces_misc/sym_g.m 0000644 0004354 0000144 00000003312 12266577513 016301 0 ustar gforget users function [field_out]=sym_g(field_in,sym_in,op_in);
%apply symmetry in the complex plane to field_in(z)
% where field_in and z are in the complex arrays.
%then, depending on op_in, we apply a compensating operation
%
%format of field_in:
% z(1,:) is z=-1+?i, z(end,:) is z=1+?i
% z(:,1) is z=?-i, z(:,end) is z=?+i
%
%sym_in: 1,2,3,4 are reflexions about x=0,y=0,x=y,x=-y
% 4+nRot is rotation by nRot*pi/2
% 0 just returns the input field as output
%op_in: 1 apply the opposite complex operation
field_out=field_in;
if isempty(whos('op_in')); op_in=0; end;
if sym_in==0;
%do nothing
elseif sym_in==1;
field_out=flipdim(field_in,1);
if op_in==1; field_out=-real(field_out)+i*imag(field_out); end;
elseif sym_in==2;
field_out=flipdim(field_in,2);
if op_in==1; field_out=real(field_out)-i*imag(field_out); end;
elseif sym_in==3;
field_out=permute(field_in,[2,1,3,4]);
if op_in==1; field_out=imag(field_out)+i*real(field_out); end;
elseif sym_in==4;
field_out=flipdim(permute(flipdim(field_in,1),[2,1,3,4]),1);
if op_in==1; field_out=-imag(field_out)-i*real(field_out); end;
elseif sym_in>=5&sym_in<=7;
for icur=1:sym_in-4;
field_out=flipdim(permute(field_out,[2,1,3,4]),1);
if op_in==1; field_out=i*field_out; end;
end
else;
fprintf('error in sym_g2\n'); return;
end;
%for test case:
%--------------
%xx=[1:10]'*ones(1,10); yy=xx';
%zz=zeros(10,10); zz(1:2,1:3)=1; zz(8,9)=2; zz(2,6:8)=-1; zz(7,3)=-2;
%sym_g(zz,1); %etc
%with the following uncommented:
%-------------------------------
%figure;
%subplot(2,2,1); imagesc(field_in);
%subplot(2,2,4); imagesc(field_out);
%xlabel('y'); xlabel('x');
gcmfaces/gcmfaces_misc/convertR4toR4nonan.m 0000644 0004354 0000144 00000000532 11600706727 020645 0 ustar gforget users function []=convertR4toR4nonan(fileIn,fileOut);
%object: read R4 fileIn, replace NaN with 0, write R4 fileOut
%inputs: fileIn and fileOut
tmp1=dir(fileIn);
tmp1=tmp1.bytes/4;
fid=fopen(fileIn,'r','b'); tmp2=fread(fid,tmp1,'float32'); fclose(fid);
tmp2(find(isnan(tmp2)))=0;
fid=fopen(fileOut,'w','b'); fwrite(fid,tmp2,'float32'); fclose(fid);
gcmfaces/gcmfaces_misc/diff_mat.m 0000644 0004354 0000144 00000006260 11625550015 016723 0 ustar gforget users function [report0]=diff_mat(file1,dir2,varargin);
%object: compute the %rms difference between fields in two mat files
%inputs: file1 is the file to compare with a reference
% dir2 is the reference file directory
%optional: file2 is the reference file (file1 by default)
%output: cell array containing the % rms difference display per field
if nargin>2; file2=varargin{1}; else; file2=file1; end;
new=open(file1); old=open([dir2 file2]);
%new list of fields:
%===================
IInew=fieldnames(new); NInew=length(IInew);
%extend to structure fields:
JJ={};
for ii=1:NInew;
tmp0=IInew{ii}; tmp1=getfield(new,tmp0);
if isstruct(tmp1);
tmp2=fieldnames(tmp1);
for jj=1:length(tmp2); JJ{length(JJ)+1}=[tmp0 '.' tmp2{jj}]; end;
else;
JJ{length(JJ)+1}=tmp0;
end;
end;
%overwrite list of fields:
IInew=JJ; NInew=length(IInew);
%old list of fields:
%===================
IIold=fieldnames(old); NIold=length(IIold);
%extend to structure fields:
JJ={};
for ii=1:NIold;
tmp0=IIold{ii}; tmp1=getfield(old,tmp0);
if isstruct(tmp1);
tmp2=fieldnames(tmp1);
for jj=1:length(tmp2); JJ{length(JJ)+1}=[tmp0 '.' tmp2{jj}]; end;
else;
JJ{length(JJ)+1}=tmp0;
end;
end;
%overwrite list of fields:
IIold=JJ; NIold=length(IIold);
%compare new to old fields:
%==========================
report0={};
report1={};
for ii=1:NInew;
test0=isempty(find(strcmp(IIold,IInew{ii})));
test1=strfind(IInew{ii},'gcm2faces');
if test0;
txt0=sprintf(['new contains ' IInew{ii} ' BUT ref does not']);
report1{length(report1)+1,1}=txt0;
elseif test1;
txt0=sprintf(['test of ' IInew{ii} ' was omitted']);
report1{length(report1)+1,1}=txt0;
else;
tmp0=IInew{ii}; tmp00=strfind(tmp0,'.');
if isempty(tmp00);
%get the field:
tmp1=getfield(new,IInew{ii}); tmp2=getfield(old,IInew{ii});
elseif length(tmp00)==1;
%get the sub field:
tmp11=IInew{ii}(1:tmp00-1);
tmp1=getfield(new,tmp11); tmp2=getfield(old,tmp11);
tmp11=IInew{ii}(tmp00+1:end);
tmp1=getfield(tmp1,tmp11); tmp2=getfield(tmp2,tmp11);
else;
error(sprintf('cannot compare %s',tmp0));
end;
%do the comparison:
if isa(tmp1,'double')|isa(tmp1,'gcmfaces');
%add blanks for display:
tmp11=ceil(length(tmp0)/10)*10-length(tmp0); tmp11=char(32*ones(1,tmp11));
%compute difference:
tmp22=nanstd(tmp1(:))^2; if tmp22==0; tmp22=nanmean( tmp1(:).^2 ); end;
tmp22=100*sqrt(nanmean( (tmp1(:)-tmp2(:)).^2 )./tmp22);
%full text for display:
txt0=[tmp0 tmp11];
%txt0=sprintf('%s differs by %i%% from %s',txt0,round(tmp22),[dir2 file2]);
txt0=sprintf('%s diff by %i%%',txt0,round(tmp22));
%do display:
if tmp22~=0; fprintf([txt0 '\n']); end;
%add to report:
report0{length(report0)+1,1}=txt0;
end;
end;
end;
%compare old field to new fields:
%================================
report2={};
for ii=1:NIold;
test0=~isempty(find(strcmp(IInew,IIold{ii})));
if ~test0;
txt0=sprintf(['ref contains ' IIold{ii} ' BUT new does not']);
report2{length(report2)+1,1}=txt0;
end;
end;
%combine the various reports:
%============================
report0={report0{:} report1{:}}';
report0={report0{:} report2{:}}';
gcmfaces/gcmfaces_misc/disp_budget_mean_zonal.m 0000644 0004354 0000144 00000003270 12276600407 021650 0 ustar gforget users function []=disp_budget_mean_zonal(lat,budg,uni,tit);
%object : plot a budget from budget_integr
%inputs: lat is the latitude vector
% budg is the result of budget_integr (in e.g. m^3/s)
% uni is the unit to be displayed
% tit is the variable description to be displayed
%display: cumulative sum of budg*median(dt)
gcmfaces_global; global myparms;
%scaled the budget terms, and integrate in time:
nyears=[myparms.yearInAve(2)-myparms.yearInAve(1)];
dt=365*86400*nyears;
cont=sum(dt*budg(1,:,:),3);
hdiv=sum(dt*budg(2,:,:),3);
zdiv=sum(dt*budg(3,:,:),3);
%do the plotting:
plot(lat,cont,'b','LineWidth',1); hold on;
plot(lat,zdiv,'r','LineWidth',1);
plot(lat,hdiv,'g','LineWidth',1);
plot(lat,cont-zdiv-hdiv,'k','LineWidth',1);
%finish the plot:
tmp1=max(abs(cont)); tmp2=max(abs(zdiv)); tmp3=max(tmp1,tmp2)*1.5;
if tmp3~=0&std(lat)~=0; axis([min(lat) max(lat) -tmp3 tmp3]); end;
grid on; legend('content','vert. div.','hor. div.','residual','Orientation','horizontal');
res=sqrt(mean((cont-zdiv-hdiv).^2));
title(sprintf('Budget %s scaled by %d years; in %s; residual : %0.1e',tit,nyears,uni,res));
%set the axis range : symmetric about 0 and rounded up
aa=axis;
tmp1=max(abs(hdiv));
tmp2=ceil(log10(tmp1));
tmp2=10^tmp2;
tmp3=2*tmp2*round(tmp1/tmp2*10)/10;
if ~isnan(tmp3); aa(3:4)=tmp3*[-1 1]; end;
axis(aa);
%print to screen:
if ~isempty(strfind(tit,'Mass')); tmp0='Mass';
elseif ~isempty(strfind(tit,'Heat')); tmp0='Heat';
elseif ~isempty(strfind(tit,'Salt')); tmp0='Salt';
else; tmp0='????';
end;
tmp1=round(log10(sqrt(mean((cont).^2))));
tmp2=round(log10(sqrt(mean((cont-zdiv-hdiv).^2))));
fprintf('%5s budget [log10(cont) | log10(res)] %3i | %3i\n',tmp0,tmp1,tmp2);
gcmfaces/gcmfaces_misc/figureL.m 0000644 0004354 0000144 00000000212 12023440611 016527 0 ustar gforget users
%figure; set(gcf,'Units','Normalized','Position',[0.1 0.2 0.8 0.8]);
figure; set(gcf,'Units','Normalized','Position',[0. 0.2 0.4 0.8]);
gcmfaces/gcmfaces_misc/runmean.m 0000644 0004354 0000144 00000005265 11626253126 016627 0 ustar gforget users function [fldOut]=runmean(fldIn,halfWindow,dim,varargin);
%object: compute running mean window ('rmw') over a dimension
%input: fldIn is the field to which the rmw will be applied
% halfWindow is the half width of the rmw
% dim is the dimension over which the rmw will be applied
%optional: doCycle states whether the boundary condition is cyclic (1) or
% not (0; default). If doCycle==0, the no. of averaged points
% decreases from 1+2*halfWindow to halfWindow at the edges,
% and we mask all of those edge points with NaNs.
%output: fldOut is the resulting field
%
%notes: - NaNs are discarded in the rmw, implying that an average is
% computed if the rmw contains at least 1 valid point.
% - setting halfWindow to 0 implies fldOut=fldIn.
%determine and check a couple things
fld_isa_gcmfaces=isa(fldIn,'gcmfaces');
if fld_isa_gcmfaces;
if dim<3;
error('for gcmfaces objects runmean excludes dim=1 or 2');
end;
end;
if nargin==3; doCycle=0; else; doCycle=varargin{1}; end;
%switch to array format if needed
if fld_isa_gcmfaces; fldIn=convert2array(fldIn); end;
%switch dim to first dimension
sizeIn=size(fldIn);
perm1to2=[1:length(sizeIn)];
perm1to2=[dim perm1to2(find(perm1to2~=dim))];
perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]];
sizeIn2=sizeIn(perm1to2);
fldIn=permute(fldIn,perm1to2);
sizeCur=size(fldIn);
if ~doCycle;
%add NaNs at both edges
sizeCur(1)=sizeCur(1)+2*halfWindow;
fldIn2=NaN*ones(sizeCur);
fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn;
fldIn=fldIn2; clear fldIn2;
end;
%create mask and remove NaNs:
fldMsk=~isnan(fldIn);
fldIn(isnan(fldIn))=0;
fldCnt=0*fldIn;
%apply the running mean
fldOut=zeros(sizeCur);
for tcur=-halfWindow:halfWindow
%To have halfWindow*2 coeffs rather than halfWindow*2+1, centered to the current
%point, it is convenient to reduce the weight of the farthest points to 1/2.
%This became necessary to get proper annual means, from monthly data, with halfWindow=6.
if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end;
tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]);
fldOut=fldOut+fac*tmp1;
tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]);
fldCnt=fldCnt+fac*tmp1;
end
fldCnt(fldCnt<2*halfWindow)=NaN;
fldOut=fldOut./fldCnt;
if ~doCycle;
fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:);
%consistent with old version bug (one point offset)
% fldOut=fldOut(halfWindow:end-halfWindow,:,:,:);
end;
%switch dimensions order back to original order
fldOut=permute(fldOut,perm2to1);
%switch back to gcmfaces format if needed
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end;
gcmfaces/gcmfaces_misc/input_list_check.m 0000644 0004354 0000144 00000004357 12357056601 020513 0 ustar gforget users function []=input_list_check(myFunction,myNargin);
%object: print warning/error messages related to
% changes in arguments lists (June 2011).
%inputs: myFunction: name of the calling routine
% myNargin: number of arguments passed to the calling routine
if strcmp(myFunction,'convert2gcmfaces')&myNargin~=1;
warning('inputCheck:convert2gcmfaces',...
['As of june 2011, convert2gcmfaces.m expects only \n' ...
' one input argument. The others were dismissed.\n' ...
' Type ''help convert2gcmfaces'' for details.']);
end;
if strcmp(myFunction,'convert2array')&myNargin~=1;
warning('inputCheck:convert2array',...
['As of june 2011, convert2array.m expects only \n' ...
' one input argument. The other ones are dismissed.\n' ...
' Type ''help convert2array'' for details.']);
end;
if strcmp(myFunction,'convert2pcol')&myNargin~=3;
warning('inputCheck:convert2pcol',...
['As of june 2011, convert2pcol.m expects only \n' ...
' three input arguments. The other ones are dismissed.\n' ...
' Type ''help convert2pcol'' for details.']);
end;
if strcmp(myFunction,'plot_one_field')&myNargin~=2;
error('inputCheck:plot_one_field',...
['As of june 2011, plot_one_field.m is a \n' ...
'function with two input argument. \n'...
'Type ''help plot_one_field'' for details.\n']);
end;
if strcmp(myFunction,'plot_std_field')&myNargin~=1;
error('inputCheck:plot_std_field',...
['As of june 2011, plot_std_field.m is a \n' ...
'function with one input argument. \n'...
'Type ''help plot_std_field'' for details.\n']);
end;
if strcmp(myFunction,'basic_diags_compute_v3_or_v4')&myNargin~=1;
error('inputCheck:basic_diags_compute_v3_or_v4',...
['As of june 2011, basic_diags_compute_v3_or_v4.m is a \n' ...
'function with one input argument. \n'...
'Type ''help basic_diags_compute_v3_or_v4'' for details.\n']);
end;
if strcmp(myFunction,'basic_diags_display_v3_or_v4')&myNargin~=1;
error('inputCheck:basic_diags_display_v3_or_v4',...
['As of june 2011, basic_diags_display_v3_or_v4.m is a \n' ...
'function with one input argument. \n'...
'Type ''help basic_diags_display_v3_or_v4'' for details.\n']);
end;
gcmfaces/gcmfaces_misc/disp_budget_mean_mask.m 0000644 0004354 0000144 00000003235 12072655270 021463 0 ustar gforget users function []=disp_budget_mean_mask(tim,budg,uni,tit);
%object : plot a budget from budget_integr
%inputs: tim is the time vector (in days)
% budg is the result of budget_integr (in e.g. m^3/s)
% uni is the unit to be displayed
% tit is the variable description to be displayed
%display: cumulative sum of budg*median(dt)
%scale the budget terms, and integrate in time:
dt=365.25*median(diff(tim));%in days
if isnan(dt); dt=1; end;%assume annual mean
dt=dt*86400;%in seconds
cont=dt*cumsum(budg(1,:));
hdiv=dt*cumsum(budg(2,:));
zdiv=dt*cumsum(budg(3,:));
%do the plotting:
plot(tim,cont,'b','LineWidth',1); hold on;
plot(tim,zdiv,'r','LineWidth',1);
plot(tim,hdiv,'g','LineWidth',1);
plot(tim,cont-zdiv-hdiv,'k','LineWidth',1);
%finish the plot:
tmp1=max(abs(cont)); tmp2=max(abs(zdiv)); tmp3=max(tmp1,tmp2)*1.5;
if tmp3~=0&std(tim)~=0; axis([min(tim) max(tim) -tmp3 tmp3]); end;
grid on; legend('content','vert. div.','hor. div.','residual','Orientation','horizontal');
res=sqrt(mean((cont-zdiv-hdiv).^2));
title(sprintf('Budget %s; in %s; residual : %0.1e',tit,uni,res));
%set the axis range : symmetric about 0 and rounded up
aa=axis;
tmp1=max(abs(cont));
tmp2=ceil(log10(tmp1));
tmp2=10^tmp2;
tmp3=2*tmp2*round(tmp1/tmp2*10)/10;
if ~isnan(tmp3); aa(3:4)=tmp3*[-1 1]; end;
axis(aa);
%print to screen:
if ~isempty(strfind(tit,'Mass')); tmp0='Mass';
elseif ~isempty(strfind(tit,'Heat')); tmp0='Heat';
elseif ~isempty(strfind(tit,'Salt')); tmp0='Salt';
else; tmp0='????';
end;
tmp1=round(log10(sqrt(mean((cont).^2))));
tmp2=round(log10(sqrt(mean((cont-zdiv-hdiv).^2))));
fprintf('%5s budget [log10(cont) | log10(res)] %3i | %3i\n',tmp0,tmp1,tmp2);
gcmfaces/gcmfaces_misc/density.m 0000644 0004354 0000144 00000006123 11625550015 016627 0 ustar gforget users function [RHOP,RHOIS,RHOR] = density(TN,SN,GDEPT,ZREF)
% [RHOP,RHOIS,RHOR] = density(Tpot,S,P,PREF)
%
% Adapted from OPA, subroutine EOS.F (bferron@ifremer.fr)
%
% ROUTINE EOS
% *******************
%
% PURPOSE :
% -------
% COMPUTE THE POTENTIAL DENSITY RHOP, IN SITU DENSITY RHOIS, AND, IF A
% REFERENCE LEVEL PREF IS SPECIFIED, THE DENSITY RHOR REFERENCED TO PREF
% FROM POTENTIAL TEMPERATURE AND SALINITY FIELDS.
%
% IMPORTANT NOTE: This version minimize the deviations of in-situ
% density compared to the UNESCO equation within
% [0-4500]dbars!!!!!! It is the version used in OPA
% (differences < 2e-4 kg/m3)
% Below 6000 dbar , differences are > 1e-3 kg/m-3
%
% METHOD :
% ------
% JACKETT AND McDOUGALL (1994) EQUATION OF STATE
% *********************
% THE IN SITU DENSITY IS COMPUTED DIRECTLY AS A FUNCTION OF
% POTENTIAL TEMPERATURE RELATIVE TO THE SURFACE, SALT AND PRESSURE
%
%
% RHO(T,S,P)
%
% WITH PRESSURE P DECIBARS
% POTENTIAL TEMPERATURE T DEG CELSIUS
% SALINITY S PSU
% DENSITY RHO KG/M**3
%
% CHECK VALUE: RHO = 1.04183326 KG/M**3 FOR P=10000 DBAR,
% T = 40 DEG CELCIUS, S=40 PSU
%
%
% JACKETT AND McDOUGALL (1994) FORMULATION
%--------------------------------------------
%
%... potential temperature, salinity and depth
ZT = TN;
ZS = SN;
if size(TN,2)>1 & size(TN,1)>1 & size(TN,1)==length(GDEPT)
ZH = GDEPT*ones(1,size(ZT,2)); % comented 13/04/99 for WHP Thorpe calc.
elseif size(TN,2)>1 & size(TN,1)>1 & size(TN,2)==length(GDEPT)
ZH = ones(size(ZT,1),1) *GDEPT';
elseif sum(size(TN))==2
ZH = ones(size(TN))*GDEPT;
else
ZH = GDEPT; % Added 13/04/99 for WHP Thorpe calc.
end
%... square root salinity
ZSR= sqrt(ZS);
%... compute density pure water at atm pressure
ZR1= ((((6.536332E-9*ZT-1.120083E-6).*ZT+1.001685E-4).*ZT ...
-9.095290E-3).*ZT+6.793952E-2).*ZT+999.842594;
%... seawater density atm pressure
ZR2= (((5.3875E-9*ZT-8.2467E-7).*ZT+7.6438E-5).*ZT ...
-4.0899E-3).*ZT+0.824493;
ZR3= (-1.6546E-6*ZT+1.0227E-4).*ZT-5.72466E-3;
ZR4= 4.8314E-4;
%
%... potential density (referenced to the surface)
RHOP= (ZR4*ZS + ZR3.*ZSR + ZR2).*ZS + ZR1;
%
%... add the compression terms
ZE = (-3.508914E-8*ZT-1.248266E-8).*ZT-2.595994E-6;
ZBW= ( 1.296821E-6*ZT-5.782165E-9).*ZT+1.045941E-4;
ZB = ZBW + ZE .* ZS;
%
ZD = -2.042967E-2;
ZC = (-7.267926E-5*ZT+2.598241E-3).*ZT+0.1571896;
ZAW= ((5.939910E-6*ZT+2.512549E-3).*ZT-0.1028859).*ZT-4.721788;
ZA = ( ZD*ZSR + ZC).*ZS + ZAW;
%
ZB1= (-0.1909078*ZT+7.390729).*ZT-55.87545;
ZA1= ((2.326469E-3*ZT+1.553190).*ZT-65.00517).*ZT+1044.077;
ZKW= (((-1.361629E-4*ZT-1.852732E-2).*ZT-30.41638).*ZT ...
+2098.925).*ZT+190925.6;
ZK0= (ZB1.*ZSR + ZA1).*ZS + ZKW;
%
%... in situ density
RHOIS = RHOP ./ (1.0-ZH./(ZK0-ZH.*(ZA-ZH.*ZB)));
%... density referenced to level ZREF
if nargin==4
RHOR = RHOP ./ (1.0-ZREF./(ZK0-ZREF.*(ZA-ZREF.*ZB)));
end
gcmfaces/gcmfaces_misc/depthStretchPlot.m 0000644 0004354 0000144 00000002741 11600706727 020460 0 ustar gforget users function [varargout]=depthStretchPlot(plotName,plotData,varargin);
%object: front end to plot/pcolor/etc. using a stretched depth coordinate (see depthStretch.m)
%inputs: plotName is the name of the plotting routine (e.g. 'pcolor')
% plotData is the list or arguments to pass over to the plotting routine
% in a cell array (e.g. {depth,time,temperature}), with the
% depth coordinate coming second (as usual for plot/pcolor/etc.)
%optional: depthTics is the vector of yTics depths
% depthStretchLim are the stretching depth limits ([0 500 6000] by default)
% to pass over to depthStretch (type 'help depthStretch' for details).
%
%notes: the depth coordinate must be first in plotData
%get depthStretchDef if provided
if nargin>2; depthTics=varargin{1}; else; depthTics=[0:100:500 750 1000:500:6000]; end;
if nargin>3; depthStretchDef=varargin{2}; else; depthStretchDef=[0 500 6000]; end;
%replace it with stretched coordinate:
plotData{2}=depthStretch(abs(plotData{2}),depthStretchDef);
%do the very plot:
eval(['h=' plotName '(plotData{:});']);
%take care of depth tics in stretched coordinate:
depthTics=sort(depthTics,'descend');
depthTicsLabel=[];
for kkk=1:length(depthTics)
depthTicsLabel=strvcat(depthTicsLabel,num2str(depthTics(kkk)));
end
depthTics=depthStretch(depthTics,depthStretchDef);
set(gca,'YTick',depthTics);
set(gca,'YTickLabel',depthTicsLabel);
%output plot handle:
if nargout>0; varargout={h}; end;
gcmfaces/gcmfaces_misc/annualmean.m 0000644 0004354 0000144 00000004164 12023434027 017267 0 ustar gforget users function [myMean,myAnom]=annualmean(myTimes,myFld,myYear);
%object: compute an annual mean for myYear (= 'first','last',n (int) or 'all')
% of a gcmfaces or array object (myFld) provided at myTimes (in years)
%input: myTimes,myFld,myYear
%output: myMean if the corresponding time mean
% (optional: myAnom is myFld-myMean)
%
%algorithm :
% - determine the number of records in one year (lYear)
% - determine the number of full years (nYears)
% - prepare list of records to average (depending on myYear, lYear and nYears)
% - in case when lYear<2 we use records as years
%
%by assumption:
% - the time dimension is last in myFld, and matches the length of myTimes
%determine the number of records in one year (lYear)
tmp1=mean(myTimes(2:end)-myTimes(1:end-1));
lYear=round(1/tmp1);
%in case when lYear<2 we use records as years
if ~(lYear>=2); lYear=1; myTimes=[1:length(myTimes)]; end;
%determine the number of full years (nYears)
nYears=floor(length(myTimes)/lYear);
%determine records that correspond to myYear, which
% may be 'first','last',n (an integer) or 'all'
if ischar(myYear);
if strcmp(myYear,'first');
recInAve=[1:lYear];
elseif strcmp(myYear,'last');
recInAve=[1:lYear]+(nYears-1)*lYear;
elseif strcmp(myYear,'all');
recInAve=[1:nYears*lYear];
else;
error('inconsistent specification of myYear');
end;
elseif (myYear>=1)&(myYear<=nYears);
recInAve=[1:lYear]+(myYear-1)*lYear;
else;
error('inconsistent specification of myYear');
end;
nRecs=length(recInAve);
%determine last dimension, which need to match the length myTimes
if strcmp(class(myFld),'gcmfaces'); nDim=size(myFld{1}); else; nDim=size(myFld); end;
if length(myTimes)>1;
if nDim(end)~=length(myTimes);
error('last dimension should match the length of myTimes');
end;
nDim=length(nDim); tt=''; for jj=1:nDim-1; tt=[tt ':,']; end;
else;
nDim=length(nDim)+1; tt=''; for jj=1:nDim-1; tt=[tt ':,']; end;
end;
%compute time mean:
eval(['myMean=mean(myFld(' tt 'recInAve),nDim);']);
%compute anomaly:
if nargout>1;
myAnom=myFld-repmat(myMean,[ones(1:nDim-1) length(myTimes)]);
end;
gcmfaces/gcmfaces_misc/write2tex.m 0000644 0004354 0000144 00000012157 12365534166 017124 0 ustar gforget users function []=write2tex(myFile,myStep,varargin);
%object: create/increment/complete/compile a tex file from within matlab
%input: myFile is the file name
% myStep is the operation to perform on the tex file
% 0 create file starting with title page (see myText)
% 1 add section or subsection (see myLev)
% 2 add a figure plus caption (see myFig)
% 3 add a paragraph
% 4 finish file
% 5 compile (latex x 2 then dvipdf)
% 6 remove temporary files (incl. *fig*.ps)
%optional myText is a cell array of text lines (for myStep=1 to 2)
% myLev is the title level (for myStep=1)
% 1=section, 2=subsection (not yet implemented)
% myFig is a figure handle (for myStep=2)
if isempty(myFile);
%if no file name, then do nothing
return;
end;
if iscell(myFile);
%extecute alternative command (myFile{1}) passing rest as arguments
eval([myFile{1} '(myFile{2:end});']);
return;
end;
myText=[]; myLev=[]; myFig=[]; myRdm=[];
if myStep<4; myText=varargin{1}; end;
if myStep==0; myRdm=varargin{2};
elseif myStep==1; myLev=varargin{2};
elseif myStep==2; myFig=varargin{2};
end;
%create file starting with write2tex.header
if myStep==0;
test0=dir(myFile);
if ~isempty(test0);
test0=input(['you are about to overwrite ' myFile ' !!! \n type 1 to proceed, 0 to stop \n']);
else;
test0=1;
end;
if ~test0;
return;
else;
fid=fopen(myFile,'w');
fprintf(fid,'\\documentclass[12pt]{beamer}\n');
fprintf(fid,'%%a nice series of examples for the beamer class:\n');
fprintf(fid,'%%http://www.informatik.uni-freiburg.de/~frank/ENG/beamer/example/beamer-class-example-en-5.html\n');
fprintf(fid,'\\usepackage{multicol}\n');
fprintf(fid,'\n');
fprintf(fid,'\\newcommand\\Fontvi{\\fontsize{6}{7.2}\\selectfont}\n');
fprintf(fid,'\n');
fprintf(fid,'\\begin{document} \n\n');
fprintf(fid,'\\title{\n');
for ii=1:length(myText); fprintf(fid,[myText{ii} '\\\\ \n']); end;
fprintf(fid,'}\n\n');
fprintf(fid,'\\date{\\today}\n\n');
fprintf(fid,'\\frame{\\titlepage}\n\n');
fprintf(fid,'\\frame{\n');
fprintf(fid,'\\frametitle{Table of contents}\n');
fprintf(fid,'\\begin{multicols}{2}\n');
fprintf(fid,'\\Fontvi\n');
fprintf(fid,'\\tableofcontents\n');
fprintf(fid,'\\end{multicols}\n');
fprintf(fid,'} \n\n');
if ~isempty(myRdm);
fprintf(fid,'\\frame{\n');
fprintf(fid,'\\section{README}\n');
fprintf(fid,'\\Fontvi\n');
for pp=1:length(myRdm);
fprintf(fid,[myRdm{pp} '\n\n']);
end;
fprintf(fid,'} \n\n');
end;
fclose(fid);
end;
myFigNumTex=0;
mySection='';
eval(['save ' myFile(1:end-4) '.mat myFigNumTex mySection;']);
end;
%open file and go to last line
fid=fopen(myFile,'a+');
eval(['load ' myFile(1:end-4) '.mat;']);
%add title or section page (see myLev)
if myStep==1;
mySection=myText;
if myLev==1; fprintf(fid,'\\section{\n');
else; fprintf(fid,'\\subsection{\n');
end;
fprintf(fid,mySection);
fprintf(fid,'} \n\n');
end;
%add a figure plus caption (see myFig)
if myStep==2;
figure(myFig);
drawnow;
myFigNumTex=myFigNumTex+1;
nn=strfind(myFile,filesep);
if ~isempty(nn);
dirTex=myFile(1:nn(end)); fileTex=myFile(nn(end)+1:end-4);
else;
dirTex='./'; fileTex=myFile(1:end-4)
end;
%print the very figure
print(myFig,'-depsc',[dirTex fileTex '.fig' num2str(myFigNumTex)]);
close;
%add figure to text file
fprintf(fid,'\\frame{ \n');
fprintf(fid,['\\frametitle{' mySection '} \n']);
fprintf(fid,'\\begin{figure}[tbh] \\centering \n');
% fprintf(fid,'\\includegraphics[width=\\textwidth,height=0.9\\textheight]');
fprintf(fid,'\\includegraphics[width=0.75\\textwidth]');
fprintf(fid,['{' fileTex '.fig' num2str(myFigNumTex) '}\n']);
fprintf(fid,'\\caption{');
for ii=1:length(myText); fprintf(fid,[myText{ii} '\n']); end;
fprintf(fid,'} \\end{figure} \n');
fprintf(fid,'} \n\n');
end;
%add a paragraph
if myStep==3;
for ii=1:length(myText);
fprintf(fid,[myText{ii} '\n']);
end;
end;
%finish file
if myStep==4; fprintf(fid,'\n\n \\end{document} \n\n'); end;
%close file
fprintf(fid,'\n\n');
fclose(fid);
eval(['save ' myFile(1:end-4) '.mat myFigNumTex mySection;']);
%compile
if myStep==5;
dirOrig=pwd;
nn=strfind(myFile,filesep);
if ~isempty(nn);
cd(myFile(1:nn(end))); fileTex=myFile(nn(end)+1:end-4);
else;
fileTex=myFile(1:end-4);
end;
system(['latex ' fileTex]);
system(['latex ' fileTex]);
system(['dvipdf ' fileTex]);
cd(dirOrig);
end;
%compile
if myStep==6&ispc;
fprintf('warning : compiling tex to pdf is bypassed on PCs\n');
end;
if myStep==6&~ispc;
dirOrig=pwd;
nn=strfind(myFile,filesep);
if ~isempty(nn);
cd(myFile(1:nn(end))); fileTex=myFile(nn(end)+1:end-4);
else;
fileTex=myFile(1:end-4);
end;
delete([fileTex '.fig*']);
delete([fileTex '.aux']);
delete([fileTex '.log']);
delete([fileTex '.out']);
delete([fileTex '.dvi']);
cd(dirOrig);
end;
gcmfaces/gcmfaces_misc/gcmfaces_msg.m 0000644 0004354 0000144 00000006754 11650105265 017601 0 ustar gforget users function [msgOut]=gcmfaces_msg(msgIn,headerOut,widthOut);
%object: print formatted message to screen
%input: msgIn is a char or a cell array of char
%optional: headerOut is the header(default: no header if
% myenv.verbose<2; calling tree if myenv.verbose>=2)
% widthOut is the line width (75 by default)
%output: msgOut is a char where line breaks were set
% so that all lines have approx. the same width.
%note: could also introduce form feeds (page break for print) or
%test case:
% msgIn={'xxxx \ryy\tzz\n' ['1 ' char(9) '2' char(10) '3' char(11) '4' char(12) '5']}
% widthOut=3;
% gcmfaces_msg(msgIn,widthOut);
gcmfaces_global;
if isempty(who('headerOut'))&myenv.verbose<2;
headerOut='';
elseif isempty(who('headerOut'))&myenv.verbose>=2;
%pause before each event
if myenv.verbose>=3; fprintf('\n > hit return to continue <\n'); pause; end;
%
[ST,I]=dbstack;
headerOut={'gcmfaces message from '};
if isempty(ST);
headerOut={headerOut{:} 'main workspace'};
else;
for ii=length(ST):-1:2;
tmp1=[char(45)*ones(1,length(ST)-ii+1)];
headerOut={headerOut{:} [tmp1 ST(ii).name ' at line ' num2str(ST(ii).line)]};
end;
end;
m=0; for ii=1:length(headerOut); m=max(m,length(headerOut{ii})); end;
tmpOut=45*ones(1,m+2);
for ii=1:length(headerOut);
n=length(headerOut{ii});
tmpOut=[tmpOut char(10) '|' headerOut{ii} 32*ones(1,m-n) '|'];
end;
tmpOut=[tmpOut char(10) 45*ones(1,m+2) char(10)];
headerOut=tmpOut;
end;
if isempty(who('widthOut')); widthOut=75; end;
%step 1: reformat msgIn to one char line, if needed
%=======
%1.1) make cell array
if ischar(msgIn)&size(msgIn,2)>1;
nLinesIn=size(msgIn,1); tmpIn={};
for ii=1:nLinesIn; tmpIn{ii}=msgIn(ii,:); end;
msgIn=tmpIn;
end;
if iscell(msgIn)&size(msgIn,2)>1; msgIn=msgIn'; end;
%2.2) cat to one char line
msgIn=strcat(cell2mat(msgIn'));
%step 2: deformat text -- rm/replace control characters from caller
%=======
%2.1) matlab control characters
ii=strfind(msgIn,'\n'); for jj=ii; msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)]; end;
ii=strfind(msgIn,'\r'); for jj=ii; msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)]; end;
ii=strfind(msgIn,'\t'); for jj=ii; msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)]; end;
ii=strfind(msgIn,'\b'); for jj=ii; msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)]; end;
ii=strfind(msgIn,'\f'); for jj=ii; msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)]; end;
%2.2) ANSI C control characters
tmp1=msgIn; tmp2=double(tmp1);
%substitute some with a space
jj=find(tmp2>8&tmp2<13); tmp1(jj)=32;
%remove the others
jj=find(tmp2>=32); tmp1=tmp1(jj);
msgIn=char(tmp1);
%2.3) double spaces
ii=strfind(msgIn,' ');
while ~isempty(ii);
jj=ii(1);
msgIn=[msgIn(1:jj-1) ' ' msgIn(jj+2:end)];
ii=strfind(msgIn,' ');
end;
%step 3: reformat text -- according to to gcmfaces_msg standards
%=======
ii=strfind(msgIn,' ');
if isempty(ii); ii=[ii length(msgIn)+1]; end;
if ii(end)~=length(msgIn); ii=[ii length(msgIn)+1]; end;
if ii(1)~=1; ii=[0 ii]; end;
msgOut=headerOut;
nn=0;
for jj=1:length(ii)-1;
kk=[ii(jj)+1 :ii(jj+1)-1];
tmp1=msgIn(kk);
if nn+length(tmp1)+1<=widthOut;
msgOut=[msgOut tmp1 ' '];
nn=nn+length(tmp1)+1;
else;
msgOut=[msgOut char(10)];
nn=0;
msgOut=[msgOut ' ' tmp1 ' '];%start new line with two spaces
nn=nn+length(tmp1)+1;
end;
end;
msgOut=[msgOut char(10)];
nn=0;
%step 4: print to screen
%=======
fprintf(msgOut);
gcmfaces/gcmfaces_misc/imagescnan.m 0000644 0004354 0000144 00000034277 11334543422 017271 0 ustar gforget users function [H,HNAN] = imagescnan(varargin)
%IMAGESCNAN Scale data and display as image with uncolored NaNs.
%
% SYNTAX:
% imagescnan(U)
% imagescnan(U,...,'NanColor',CNAN)
% imagescnan(U,...,'NanMask',MNAN)
% imagescnan(U,...,IOPT)
% imagescnan(X,Y,U,...)
% [H,HNAN] = imagescnan(...);
%
% INPUT:
% U - 2 dimensional N-by-M image or N-by-M-by-3 RGB image.
% X - 2 extrema X-axis data; or the M values; or the N-by-M values
% as obtained from MESHGRID (see DESCRIPTION below).
% DEFAULT: [1 N]
% Y - 2 extrema X-axis data; or the N values; or the N-by-M values
% as obtained from MESHGRID (see DESCRIPTION below).
% DEFAULT: [1 M]
% CNAN - Color for the NaNs elements. May be a char specifier or an [R
% G B] triplet specifying the color.
% DEFAULT: invisible (axes background color)
% MNAN - Elements to be ignored besides not finite values. May be an
% scalar or a logical M-by-N matrix indicating the elements to
% be ignored.
% DEFAULT: []
% IOPT - IMAGE function normal optional pair arguments like
% ('Parent',H) or/and CLIM like optional last argument as in
% IMAGESC.
% DEFAULT: none
%
% OUTPUT (all optional):
% H - Image handle
% HNAN - Handle of every ignored (NaN) value colored patch.
%
% DESCRIPTION:
% MATLAB function IMAGESC does not work properly with NaNs. This
% programs deals with this problem by including colored patches over
% this elements and maybe others specyfied by the user with MNAN.
%
% Besides, those functions does not work properly with X,Y values
% variable interval, but this functions does it by generating a whole
% new image of several rectangular patches, but whose centers may not
% lay in the specified coordinate (see NOTE below). This functionality
% is experimental and not recommended (see ADDITIONAL NOTES inside this
% program).
%
% In previous release, 2-dim input images were transformed into a
% 3-dim RGB image. This is not used anymore (see ADDITIONAL NOTES
% inside this file).
%
% NOTE:
% * Optional inputs use its DEFAULT value when not given or [].
% * Optional outputs may or not be called.
% * If X is a two element vector, min(X) will be the coordinate of the
% first column and max(X) of the last column.
% * If Y is a two element vector, min(Y) will be the coordinate of the
% first row and max(Y) of the last row.
% * If vector X-axis is decreasing U=fliplr(U) will be used.
% * If vector Y-axis is decreasing U=flipud(U) will be used.
% * When X or Y do not have a constant increasing/decreasing step, the
% vertices of the color rectangules are set in the middle of each
% pair of coordinates. For this reason its center may not lay on the
% specified coordinate, except on the coordinates at the edges where
% it always lays on the center.
% * To get a non-scaled image (IMAGE instead of IMAGESC) use:
% >> H = imagescnan(...);
% >> set(H,'CDataMapping','direct')
% * ADDITIONAL NOTES are included inside this file.
%
% EXAMPLE:
% % Compares with normal IMAGESC:
% N = 100;
% PNaNs = 0.10;
% U = peaks(N);
% U(round(1 + (N^2-1).*rand(N^2*PNaNs,1))) = NaN; % Adds NaNs
% subplot(221), imagesc(U)
% title('With IMAGESC: ugly NaNs')
% subplot(222), imagescnan(U)
% title('With IMAGESCNAN: uncolored NaNs')
% % Compares with SPY:
% subplot(223), spy(isnan(U))
% title('SPY(isnan(U))')
% subplot(224), imagescnan(isnan(U),'NaNMask',0), axis equal tight
% title('SPY with IMAGESCNAN')
%
% SEE ALSO:
% IMAGE, IMAGESC, COLORBAR, IMREAD, IMWRITE
% and
% CMAPPING, CBFREEZE by Carlos Vargas
% at http://www.mathworks.com/matlabcentral/fileexchange
%
%
% ---
% MFILE: imagescnan.m
% VERSION: 2.1 (Aug 20, 2009) (download)
% MATLAB: 7.7.0.471 (R2008b)
% AUTHOR: Carlos Adrian Vargas Aguilera (MEXICO)
% CONTACT: nubeobscura@hotmail.com
% ADDITIONAL NOTES:
% * I keep getting a kind of BUG with the edges of the patched NaNs. I
% added two NOTE inside this program that may fix this problem.
% Another way is to convert the intensity matrix U into RGB colors by
% using the CMAPPING function, as used by the first version of this
% program.
% * Besides, if the matrix is too large, sometimes there is an
% undocumented failure while drawing the patch NaNs. Is recommended
% to use U = cmapping(U,[],'k','discrete') instead, and change the
% CLIM to [min(U(:)) max(U(:))].
% * The use of not homogeneous step interval X,Y axes is not
% recommended because the program tries to put its value in the
% middle of the colored rectangule (as IMAGESC does) and soetimes the
% result may not be what the user wants. So this is for experimental
% use only.
% REVISIONS:
% 1.0 Released. (Jun 30, 2008)
% 1.1 Fixed bug when CAXIS used. Colorbar freezed colormap. Fixed
% bug in color vector input (Found by Greg King) and now
% accets RGB image as input. (Jul 14, 2008)
% 2.0 Totally rewritten code. Do not converts to RGB anymore. Do not
% freezes the colormap anymore. Do not output any colorbar. New
% X and Y variable steps accepted input. Now uses patches. (Jun
% 08, 2009)
% 2.1 Fixed bug with RGB input. Added a NOTE about the use of
% CMAPPING. (Aug 20, 2009)
% DISCLAIMER:
% imagescnan.m is provided "as is" without warranty of any kind, under
% the revised BSD license.
% Copyright (c) 2008,2009 Carlos Adrian Vargas Aguilera
% INPUTS CHECK-IN
% -------------------------------------------------------------------------
% Initializes:
X = [];
Y = [];
CNAN = [];
MNAN = [];
ha = [];
% Checks number of inputs:
if nargin<1
error('CVARGAS:imagescnan:notEnoughInputs',...
'At least 1 input is required.')
elseif nargout>2
error('CVARGAS:imagescnan:tooManyOutputs',...
'At most 2 outputs are allowed.')
end
% Gets X,Y,U:
if ((nargin==1) || (nargin==2))
U = varargin{1};
varargin(1) = [];
else
if (isnumeric(varargin{1}) && isnumeric(varargin{2}) && ...
isnumeric(varargin{3}))
X = varargin{1};
Y = varargin{2};
U = varargin{3};
varargin(1:3) = [];
else
U = varargin{1};
varargin(1) = [];
end
end
% Check U:
ndim = ndims(U);
if (ndim==2)
[M,N] = size(U);
O = 1;
elseif (ndim==3)
[M,N,O] = size(U);
if (O~=3)
error('CVARGAS:imagescnan:incorrectRgbImage',...
'RGB image must be of size M-by-N-by-3.')
end
else
error('CVARGAS:imagescnan:incorrectImageSize',...
'Image must be 2-dimensional or a 3-dim RGB image.')
end
% Check X:
aequal = true; % Equal intervals on x-axis?
dX = [];
if isempty(X)
X = [1 N];
else
if (ndims(X)>2)
error('CVARGAS:imagescnan:incorrectXDims',...
'X must be a vector or a matrix as a result of MESHGRID.')
end
if any(~isfinite(X(:)))
error('CVARGAS:imagescnan:incorrectXValue',...
'X elements must be numeric and finite.')
end
[Mx,Nx] = size(X);
if ((Mx*Nx)==2)
if X(2)(eps*max(abs(dX(:)))*1000))
error('CVARGAS:imagescnan:incorrectXMatrix',...
'X matrix must be as generated by MESHGRID.')
end
X = X(1,:);
elseif (~any([Mx Nx]==1) && ~((Mx*Nx)==N))
error('CVARGAS:imagescnan:incorrectXSize',...
'X must be an scalar or a matrix.')
end
% Forces ascending x-axis:
[X,I] = sort(X(:).');
for k = 1:O % Fixed bug Aug 2009
U(:,:,k) = U(:,I,k);
end
clear I
% Checks equal intervals:
dX = diff(X);
if any(abs(dX(1)-dX(2:end))>(eps*max(dX)*1000))
if aequal
aequal = false;
end
else
X = [X(1) X(end)];
dX = [];
end
end
end
% Check Y:
dY = [];
if isempty(Y)
Y = [1 M];
else
if (ndims(Y)>2)
error('CVARGAS:imagescnan:incorrectYDims',...
'Y must be a vector or a matrix as a result of MESHGRID.')
end
if any(~isfinite(Y(:)))
error('CVARGAS:imagescnan:incorrectYValue',...
'Y elements must be numeric and finite.')
end
[My,Ny] = size(Y);
if ((My*Ny)==2)
if Y(2)(eps*max(abs(dY(:)))*1000))
error('CVARGAS:imagescnan:incorrectYMatrix',...
'Y matrix must be as generated by MESHGRID.')
end
Y = Y(:,1);
elseif (~any([My Ny]==1) && ~((My*Ny)==M))
error('CVARGAS:imagescnan:incorrectYSize',...
'Y must be an scalar or a matrix.')
end
% Forces ascending y-axis:
[Y,I] = sort(Y(:).');
for k = 1:O % Fixed bug Aug 2009
U(:,:,k) = U(I,:,k);
end
clear I
% Checks equal intervals:
dY = diff(Y);
if any(abs(dY(1)-dY(2:end))>(eps*max(dY)*1000))
if aequal
aequal = false;
end
else
Y = [Y(1) Y(end)];
dY = [];
end
end
end
% Checks varargin:
ind = [];
Nopt = length(varargin);
for k = 1:Nopt-1
if (~isempty(varargin{k}) && ischar(varargin{k}))
if strncmpi(varargin{k},'NanColor',4)
CNAN = varargin{k+1};
ind = [ind k k+1];
elseif strncmpi(varargin{k},'NanMask',4)
MNAN = varargin{k+1};
ind = [ind k k+1];
elseif (strncmpi(varargin{k},'Parent',2) && isempty(CNAN))
try
CNAN = get(varargin{k+1},'Color');
ha = varargin{k+1};
catch
error('CVARGAS:imagescnan:incorrectParentHandle',...
'''Parent'' must be a valid axes handle.')
end
end
end
end
varargin(ind) = [];
Nargin = length(varargin);
% Check ha:
if isempty(ha)
ha = gca;
end
% Check CNAN:
if isempty(CNAN)
CNAN = get(ha,'Color');
elseif ischar(CNAN)
switch lower(CNAN)
case 'y', CNAN = [1 1 0];
case 'm', CNAN = [1 0 0];
case 'c', CNAN = [0 1 1];
case 'r', CNAN = [1 0 0];
case 'g', CNAN = [0 1 0];
case 'b', CNAN = [0 0 1];
case 'w', CNAN = [1 1 1];
case 'k', CNAN = [0 0 0];
otherwise
error('CVARGAS:imagescnan:incorrectNancString',...
'Color string must be a valid color identifier. One of ''ymcrgbwk''.')
end
elseif isnumeric(CNAN) && (length(CNAN)==3)
CNAN = CNAN(:).'; % Forces row vector.
else
error('CVARGAS:imagescnan:incorrectNancInput',...
'Not recognized CNAN input.')
end
% Check MNAN:
if isempty(MNAN)
MNAN = any(~isfinite(U),3);
else
if (ndims(MNAN)==2)
[Mm,Nm] = size(MNAN);
if ((Mm*Nm)==1)
MNAN = (any(~isfinite(U),3) | any(U==MNAN,3));
elseif ((Mm==M) && (Nm==N) && islogical(MNAN))
MNAN = (any(~isfinite(U),3) | MNAN);
else
error('CVARGAS:imagescnan:incorrectNanmSize',...
'MNAN must be an scalar or a logical matrix of size M-by-N.')
end
else
error('CVARGAS:imagescnan:incorrectNanmDims',...
'MNAN must be an scalar or a matrix.')
end
end
% -------------------------------------------------------------------------
% MAIN
% -------------------------------------------------------------------------
% Generates the image:
if aequal
% IMAGESC way.
H = imagesc(X,Y,U,varargin{:});
else
% PATCH way.
% Check clim:
if (rem(Nargin,2)==1)
clim = varargin{end};
varargin(end) = [];
if ((length(clim)~=2) || (clim(1)>clim(2)))
error('CVARGAS:imagescnan:incorrectClimInput',...
'clim must be a 2 element increasing vector.')
end
else
clim = [];
end
% Generates vertices between coordinates (coordinates may not be at the
% center of these vertices):
if (length(X)~=N)
X = (0:N-1)*((X(2)-X(1))/(N-1)) + X(1);
end
if (length(Y)~=M)
Y = (0:M-1)*((Y(2)-Y(1))/(M-1)) + Y(1);
end
if isempty(dX)
dX = diff(X);
end
if isempty(dY)
dY = diff(Y);
end
[X,Y] = meshgrid([X(1)-dX(1)/2 X+dX([1:N-1 N-1])/2],...
[Y(1)-dY(1)/2 Y+dY([1:M-1 M-1])/2]);
% Generates faces:
ind = (1:(M+1)*N)';
ind(M+1:M+1:end) = [];
% Generates patches:
H = patch(...
'Vertices' ,[X(:) Y(:)],...
'Faces' ,[ind ind+1 ind+M+2 ind+M+1],...
'FaceVertexCData',U(:),...
'FaceColor' ,'flat',...
'EdgeColor' ,'none',... % NOTE: Sometimes this is not required.
varargin{:});
set(ha,...
'YDir' ,'reverse',...
'View' ,[0 90],...
'Box' ,'on',...
'Layer','top')
axis(ha,'tight')
% Sets clim:
if ~isempty(clim)
set(ha,'CLim',clim)
else
set(ha,'CLimMode','auto')
end
end
% Adds NaNs patches:
if any(MNAN(:))
if aequal
% dX and dY is constant:
[MNAN,NNAN] = ind2sub([M,N],find(MNAN));
Nnan = length(MNAN);
dX = (X(2)-X(1))/(N-1)/2;
dY = (Y(2)-Y(1))/(M-1)/2;
HNAN = patch(repmat((X(1)+(NNAN(:)'-1)*(2*dX)),4,1) + ...
(repmat([-1 1 1 -1]'*dX,1,Nnan)),...
repmat((Y(1)+(MNAN(:)'-1)*(2*dY)),4,1) + ...
(repmat([1 1 -1 -1]'*dY,1,Nnan)),...
CNAN,...
'EdgeColor',CNAN,... 'EdgeColor','none',...
varargin{1:Nargin-rem(Nargin,2)});
else
% dX and/or dY is not constant:
MNAN = find(MNAN);
HNAN = patch(...
'Vertices' ,[X(:) Y(:)],...
'Faces' ,[ind(MNAN) ind(MNAN)+1 ind(MNAN)+M+2 ind(MNAN)+M+1],...
'FaceColor' ,CNAN,...
'EdgeColor' ,'none',... 'EdgeColor',CNAN,... % NOTE: may be better?
varargin{:});
end
else
HNAN = [];
end
% OUTPUTS CHECK-OUT
% -------------------------------------------------------------------------
% Clears outputs?:
if (nargout==0)
clear H
end
% [EOF] imagescnan.m gcmfaces/gcmfaces_legacy/ 0000755 0004354 0000144 00000000000 12375431055 015306 5 ustar gforget users gcmfaces/gcmfaces_legacy/line_greatC_TUV_MASKS_test.m 0000644 0004354 0000144 00000004011 11600767631 022431 0 ustar gforget users
%test:
%lons=[-81 -77]; lats=[28 27]; name='Florida Strait';
%lons=[-81 -79]; lats=[28 22]; name='Florida Strait W1';
%lons=[-76 -76]; lats=[21 8]; name='Florida Strait S1';
%lons=[-77 -77]; lats=[27 25]; name='Florida Strait E1';
%lons=[-77 -77]; lats=[25 22]; name='Florida Strait E2';
%lons=[-65 -50]; lats=[66 66]; name='Davis Strait';
%lons=[-35 -20]; lats=[67 65]; name='Denmark Strait';
%lons=[-16 -7]; lats=[65 62.5]; name='Iceland Feroe';
%lons=[-6.5 -4]; lats=[62.5 57]; name='Feroe England';
%lons=[-4 8]; lats=[57 62]; name='England Norway';
%lons=[-68 -63]; lats=[-54 -66]; name='Drake Passage';
%lons=[103 103]; lats=[-1 4]; name='Indonesia W1';
%lons=[104 109]; lats=[-3 -8]; name='Indonesia W2';
%lons=[113 118]; lats=[-8.5 -8.5]; name='Indonesia W3';
%lons=[118 127 ]; lats=[-8.5 -15]; name='Indonesia W4';
%lons=[127 127]; lats=[-25 -68]; name='Australia Antarctica';
%lons=[38 46]; lats=[-10 -22]; name='Madagascar Channel';
%lons=[46 46]; lats=[-22 -69]; name='Madagascar Antarctica';
%lons=[20 20]; lats=[-30 -69.5]; name='South Africa Antarctica';
%compute:
line_cur=line_greatC_TUV_mask(lons,lats);
line_cur.name=name;
%check:
ii=1; close all;
tmp1=mygrid.hFacW{ii}(:,:,1); tmp2=line_cur.mmuIn{ii}; tmp1(~isnan(tmp2))=-1; figure; imagesc(tmp1);
%axis([40 60 50 100]);
%axis([50 190 40 90]);
axis([230 270 0 50]);
drawnow; refresh; pause;
tmp1=mygrid.hFacS{ii}(:,:,1); tmp2=line_cur.mmvIn{ii}; tmp1(~isnan(tmp2))=-1; figure; imagesc(tmp1);
%axis([40 60 50 100]);
%axis([50 190 40 90]);
axis([230 270 0 50]);
drawnow; refresh; pause;
[X,Y,MSK]=convert_llc2pcol(mygrid.XC,mygrid.YC,mygrid.Depth,1,1); MSK(MSK==0)=NaN;
figure; pcolor(X,Y,MSK); axis([-180 180 -90 90]); shading flat; caxis([0 6000]);
hold on; plot(lons(1),lats(1),'r.','MarkerSize',32);
plot(lons(2),lats(2),'k.','MarkerSize',32);
[X,Y,mmtIn]=convert_llc2pcol(mygrid.XC,mygrid.YC,line_cur.mmtIn,1,1);
kk=find(mmtIn==1); plot(X(kk),Y(kk),'m.');
%axis([-90 -70 10 30]);
%axis([90 140 -70 30]);
%axis([15 30 -72 -30]);
axis([-80 20 50 70]);
drawnow; refresh; pause;
gcmfaces/gcmfaces_legacy/grid_load_occa.m 0000644 0004354 0000144 00000002067 11600376700 020375 0 ustar gforget users
global mygrid;
if isempty(mygrid);
dirGrid='/net/ross/raid2/gforget/1x1_50levels/GRID/';
list0={'XC','XG','YC','YG','RC','RF','RAC','DRC','DRF',...
'DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
%list0={'AngleCS','AngleSN','Depth','DRC','DRF','DXC','DXG','DYC','DYG',...
% 'hFacC','hFacS','hFacW','maskCtrlC','maskCtrlS','maskCtrlW',...
% 'PHrefC','PHrefF','RAC','RAS','RAW','RAZ','RC','RF','XC','XG','YC','YG'};
list0={'XC','XG','YC','YG','RAC','DXC','DXG','DYC','DYG',...
'hFacC','hFacS','hFacW','Depth'};
for iFld=1:length(list0);
tmp1=gcmfaces(1,'ll'); eval(['tmp2=rdmds([dirGrid ''' list0{iFld} '*'']);']);
tmp1.f1=tmp2; eval(['mygrid.' list0{iFld} '=tmp1;']);
end;
mygrid.AngleCS=mygrid.XC; mygrid.AngleCS(:,:)=1;
mygrid.AngleSN=mygrid.XC; mygrid.AngleSN(:,:)=0;
list0={'RC','RF','DRC','DRF'};
for iFld=1:length(list0);
eval(['mygrid.' list0{iFld} '=rdmds([dirGrid ''' list0{iFld} '*'']);']);
end;
mygrid.hFacCsurf=mygrid.hFacC;
for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
end;
gcmfaces/gcmfaces_legacy/convert_UV2UeVn.m 0000644 0004354 0000144 00000000142 11334544171 020431 0 ustar gforget users function [fldUe,fldVn]=convert_UV2UeVn(fldU,fldV);
[fldUe,fldVn]=calc_UEVNfromUXVY(fldU,fldV);
gcmfaces/gcmfaces_legacy/nanrunmean.m 0000644 0004354 0000144 00000001713 11600706724 017626 0 ustar gforget users %this function does a running mean along dimension dim_cur,
%averaging over indices [-nb_cur:nb_cur]
%
function [field_out]=nanrunmean(field_in,nb_cur,dim_cur);
if nb_cur~=0
size_cur=size(field_in);
perm1to2=[1:length(size_cur)];
perm1to2=[dim_cur perm1to2(find(perm1to2~=dim_cur))];
perm2to1=[[1:dim_cur-1]+1 1 [dim_cur+1:length(size_cur)]];
size_cur2=size_cur(perm1to2);
field_in2=permute(field_in,perm1to2);
field_mask2=1*~isnan(field_in2);
field_in2(isnan(field_in2))=0;
field_out2=zeros(size_cur2);
count_out2=zeros(size_cur2);
for tcur=-nb_cur:nb_cur
tmp1=[tcur:tcur-1+size_cur2(1)];
tmp2=find(tmp1>=1&tmp1<=size_cur2(1));
field_out2(tmp2,:)=field_out2(tmp2,:)+field_in2(tmp1(tmp2),:);
count_out2(tmp2,:)=count_out2(tmp2,:)+field_mask2(tmp1(tmp2),:);
end
tmp1=find(count_out2>0);
field_out2(tmp1)=field_out2(tmp1)./count_out2(tmp1);
field_out2(count_out2==0)=NaN;
field_out=permute(field_out2,perm2to1);
else
field_out=field_in;
end%if nb_cur~=0
gcmfaces/gcmfaces_legacy/grid_load_96x23.m 0000644 0004354 0000144 00000001707 11600376700 020263 0 ustar gforget users
global mygrid;
if isempty(mygrid);
dirGrid='/net/altix3700/raid4/gforget/mysetups/ecco_v4/RUNS/GRIDmds/';
list0={'XC','XG','YC','YG','RC','RF','RAC','DRC','DRF',...
'DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
%list0={'AngleCS','AngleSN','Depth','DRC','DRF','DXC','DXG','DYC','DYG',...
% 'hFacC','hFacS','hFacW','maskCtrlC','maskCtrlS','maskCtrlW',...
% 'PHrefC','PHrefF','RAC','RAS','RAW','RAZ','RC','RF','XC','XG','YC','YG'};
list0={'XC','XG','YC','YG','RAC','DXC','DXG','DYC','DYG',...
'hFacC','hFacS','hFacW','Depth','AngleCS','AngleSN'};
for iFld=1:length(list0);
eval(['mygrid.' list0{iFld} '=rdmds_compact2llc([dirGrid ''' list0{iFld} '*'']);']);
end;
list0={'RC','RF','DRC','DRF'};
for iFld=1:length(list0);
eval(['mygrid.' list0{iFld} '=rdmds([dirGrid ''' list0{iFld} '*'']);']);
end;
mygrid.hFacCsurf=mygrid.hFacC;
for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
end;
gcmfaces/gcmfaces_legacy/diffsmooth2Dstraight.m 0000644 0004354 0000144 00000000124 11334544171 021555 0 ustar gforget users function [fld]=diffsmooth2Dstraight(varargin);
[fld]=diffsmooth2D(varargin{:});
gcmfaces/gcmfaces_legacy/runmean_cycle.m 0000644 0004354 0000144 00000001404 11600706724 020305 0 ustar gforget users %this function does a running mean along dimension dim_cur,
%averaging over indices [-nb_cur:nb_cur]
%
function [field_out]=runmean_cycle(field_in,nb_cur,dim_cur);
if nb_cur~=0
size_cur=size(field_in);
perm1to2=[1:length(size_cur)];
perm1to2=[dim_cur perm1to2(find(perm1to2~=dim_cur))];
perm2to1=[[1:dim_cur-1]+1 1 [dim_cur+1:length(size_cur)]];
size_cur2=size_cur(perm1to2);
field_in2=permute(field_in,perm1to2);
field_out2=zeros(size_cur2);
count_out2=0;
for tcur=-nb_cur:nb_cur
field_out2=field_out2+circshift(field_in2,[tcur zeros(1,length(size_cur2)-1)]);
count_out2=count_out2+1;
end
field_out2=field_out2/count_out2;
%field_out=field_in-permute(field_out2,perm2to1);
field_out=permute(field_out2,perm2to1);
else
field_out=field_in;
end%if nb_cur~=0
gcmfaces/gcmfaces_legacy/line_greatC_TUV_MASKS_v4_until_evol4.m 0000644 0004354 0000144 00000004541 11600767631 024337 0 ustar gforget users function [LINES_out]=line_greatC_TUV_MASKS(varargin);
if nargin==1; doDisplay=varargin{1}; else; doDisplay=0; end;
global mygrid;
for iy=1:21;
switch iy;
case 1; lons=[-173 -164]; lats=[65.5 65.5]; name='Bering Strait';
case 2; lons=[-5 -5]; lats=[34 40]; name='Gibraltar';
case 3; lons=[-81 -77]; lats=[28 27]; name='Florida Strait';
case 4; lons=[-81 -79]; lats=[28 22]; name='Florida Strait W1';
case 5; lons=[-76 -76]; lats=[21 8]; name='Florida Strait S1';
case 6; lons=[-77 -77]; lats=[27 25]; name='Florida Strait E1';
case 7; lons=[-77 -77]; lats=[25 22]; name='Florida Strait E2';
case 8; lons=[-65 -50]; lats=[66 66]; name='Davis Strait';
case 9; lons=[-35 -20]; lats=[67 65]; name='Denmark Strait';
case 10; lons=[-16 -7]; lats=[65 62.5]; name='Iceland Feroe';
case 11; lons=[-6.5 -4]; lats=[62.5 57]; name='Feroe England';
case 12; lons=[-4 8]; lats=[57 62]; name='England Norway';
case 13; lons=[-68 -63]; lats=[-54 -66]; name='Drake Passage';
case 14; lons=[103 103]; lats=[-1 4]; name='Indonesia W1';
case 15; lons=[104 109]; lats=[-3 -8]; name='Indonesia W2';
case 16; lons=[113 118]; lats=[-8.5 -8.5]; name='Indonesia W3';
case 17; lons=[118 127 ]; lats=[-8.5 -15]; name='Indonesia W4';
case 18; lons=[127 127]; lats=[-25 -68]; name='Australia Antarctica';
case 19; lons=[38 46]; lats=[-10 -22]; name='Madagascar Channel';
case 20; lons=[46 46]; lats=[-22 -69]; name='Madagascar Antarctica';
case 21; lons=[20 20]; lats=[-30 -69.5]; name='South Africa Antarctica';
end;
%compute:
line_cur=line_greatC_TUV_mask(lons,lats);
line_cur.name=name;
if doDisplay;
ii=5;
tmp1=mygrid.hFacW{ii}(:,:,1); tmp2=line_cur.mmuIn{ii}; tmp1(~isnan(tmp2))=-1; figure; imagesc(tmp1);
tmp1=mygrid.hFacS{ii}(:,:,1); tmp2=line_cur.mmvIn{ii}; tmp1(~isnan(tmp2))=-1; figure; imagesc(tmp1);
[X,Y,MSK]=convert_llc2pcol(mygrid.XC,mygrid.YC,mygrid.Depth,1,1); MSK(MSK==0)=NaN;
figure; pcolor(X,Y,MSK); axis([-180 180 -90 90]); shading flat; caxis([0 6000]);
hold on; plot(lons(1),lats(1),'r.','MarkerSize',32);
plot(lons(2),lats(2),'k.','MarkerSize',32);
[X,Y,mmtIn]=convert_llc2pcol(mygrid.XC,mygrid.YC,line_cur.mmtIn,1,1);
kk=find(mmtIn==1); plot(X(kk),Y(kk),'m.');
%[X,Y,mmtOut]=convert_llc2pcol(mygrid.XC,mygrid.YC,line_cur.mmtOut,1,1);
%kk=find(mmtOut==1); plot(X(kk),Y(kk),'cx');
end;
if iy==1;
LINES_out=line_cur;
else;
LINES_out(iy)=line_cur;
end;
end;
gcmfaces/gcmfaces_legacy/line_zonal_TUV_mask.m 0000644 0004354 0000144 00000002151 11601141417 021355 0 ustar gforget users function [line_out]=line_zonal_TUV_mask(lat_in);
global mygrid;
y=mygrid.YC;
yP1=exch_T_N(y);
%T part:
yy=lat_in;
mm=yP1<=yy; mm(isnan(mm))=0;
for iF=1:y.nFaces;
eval(['tmp1=mm.f' num2str(iF) ';']);
tmp2=circshift(tmp1,[-1 0])+circshift(tmp1,[1 0])+...
circshift(tmp1,[0 -1])+circshift(tmp1,[0 1]);
tmp2(tmp1>0)=0; tmp2(tmp2~=0)=1; tmp2(tmp2==0)=NaN;
eval(['mm.f' num2str(iF) '=tmp2(2:end-1,2:end-1);']);
end;
mmt=mm;
%UV part:
nn=exch_T_N(mmt);
mm=yP1<=yy; mm(isnan(mm))=0;
mmu=y; mmv=y;
for iFace=1:y.nFaces;
iF=num2str(iFace);
eval(['tmp_u=mm.f' iF '(1:end-1,:)==1&nn.f' iF '(2:end,:)==1;']);
eval(['tmp_u2=mm.f' iF '(2:end,:)==1&nn.f' iF '(1:end-1,:)==1;']);
tmp_u=1*(tmp_u(1:end-1,2:end-1)-tmp_u2(1:end-1,2:end-1)); eval(['mmu.f' iF '=tmp_u;']);
eval(['tmp_v=mm.f' iF '(:,1:end-1)==1&nn.f' iF '(:,2:end)==1;']);
eval(['tmp_v2=mm.f' iF '(:,2:end)==1&nn.f' iF '(:,1:end-1)==1;']);
tmp_v=1*(tmp_v(2:end-1,1:end-1)-tmp_v2(2:end-1,1:end-1)); eval(['mmv.f' iF '=tmp_v;']);
end;
mmu(mmu==0)=NaN; mmv(mmv==0)=NaN;
%store:
line_out=struct('mmt',mmt,'mmu',mmu,'mmv',mmv,'lat',lat_in);
gcmfaces/gcmfaces_legacy/convert_compact2llc.m 0000644 0004354 0000144 00000000125 11334544171 021423 0 ustar gforget users function [v1]=convert_compact2llc(varargin);
[v1]=convert2gcmfaces(varargin{:});
gcmfaces/gcmfaces_legacy/diffsmooth2Dstraight_extrapolate.m 0000644 0004354 0000144 00000000153 11334544171 024167 0 ustar gforget users function [fld]=diffsmooth2Dstraight_extrapolate(varargin);
[fld]=diffsmooth2D_extrap_fwd(varargin{:});
gcmfaces/gcmfaces_legacy/myverticaldeform.m 0000644 0004354 0000144 00000001111 11600706724 021030 0 ustar gforget users function [varout]=myverticaldeform(varin)
%varin : profondeur en m (positif)
%varout : z correspondant (negatif), avec
% deformation selon le niveau
%var de sortie :
varout=zeros(size(varin));
%valeurs limite :
%val_lim=500;fprintf('deform v:500m\n');
val_lim=500;
%val_lim=250; fprintf('deform v:250m\n');
%val_lim=-0.0001;
val_lim2=2000;
%valeurs inferieures a 500m :
tmp1=find(varin<=val_lim);
tmp2=-100*varin(tmp1)/val_lim;
varout(tmp1)=tmp2;
%valeurs superieurs a 500m :
tmp1=find(varin>val_lim);
tmp2=-100-200*(varin(tmp1)-val_lim)./(val_lim2-val_lim);
varout(tmp1)=tmp2;
gcmfaces/gcmfaces_legacy/line_zonal_TUV_MASKS.m 0000644 0004354 0000144 00000000322 11601141417 021276 0 ustar gforget users function [LINES_out]=line_zonal_TUV_MASKS(Y_in);
global mygrid;
for iy=1:length(Y_in);
line_cur=line_zonal_TUV_mask(Y_in(iy));
if iy==1;
LINES_out=line_cur;
else;
LINES_out(iy)=line_cur;
end;
end;
gcmfaces/gcmfaces_legacy/rdmds_compact2llc.m 0000644 0004354 0000144 00000000131 11334544171 021051 0 ustar gforget users function [fldOut]=rdmds_compact2llc(varargin);
[fldOut]=rdmds2gcmfaces(varargin{:});
gcmfaces/gcmfaces_legacy/nanrunmean_cycle.m 0000644 0004354 0000144 00000001640 11600706724 021004 0 ustar gforget users %this function does a running mean along dimension dim_cur,
%averaging over indices [-nb_cur:nb_cur]
%
function [field_out]=nanrunmean_cycle(field_in,nb_cur,dim_cur);
if nb_cur~=0
size_cur=size(field_in);
perm1to2=[1:length(size_cur)];
perm1to2=[dim_cur perm1to2(find(perm1to2~=dim_cur))];
perm2to1=[[1:dim_cur-1]+1 1 [dim_cur+1:length(size_cur)]];
size_cur2=size_cur(perm1to2);
field_in2=permute(field_in,perm1to2);
field_mask2=1*~isnan(field_in2);
field_out2=zeros(size_cur2);
count_out2=zeros(size_cur2);
for tcur=-nb_cur:nb_cur
tmp1=circshift(field_in2,[tcur zeros(1,length(size_cur2)-1)]);
tmp2=find(~isnan(tmp1.*field_in2));
field_out2(tmp2)=field_out2(tmp2)+tmp1(tmp2);
count_out2(tmp2)=count_out2(tmp2)+1;
end
tmp1=find(count_out2>0);
field_out2(tmp1)=field_out2(tmp1)./count_out2(tmp1);
field_out2(count_out2==0)=NaN;
field_out=permute(field_out2,perm2to1);
else
field_out=field_in;
end%if nb_cur~=0
gcmfaces/gcmfaces_legacy/grid_load_90x50.m 0000644 0004354 0000144 00000001715 11600376700 020254 0 ustar gforget users
global mygrid;
if isempty(mygrid);
dirGrid='/net/altix3700/raid4/gforget/mysetups/ecco_v4/RUNS/GRIDmds_90x50/';
list0={'XC','XG','YC','YG','RC','RF','RAC','DRC','DRF',...
'DXC','DXG','DYC','DYG','hFacC','hFacS','hFacW','Depth'};
%list0={'AngleCS','AngleSN','Depth','DRC','DRF','DXC','DXG','DYC','DYG',...
% 'hFacC','hFacS','hFacW','maskCtrlC','maskCtrlS','maskCtrlW',...
% 'PHrefC','PHrefF','RAC','RAS','RAW','RAZ','RC','RF','XC','XG','YC','YG'};
list0={'XC','XG','YC','YG','RAC','DXC','DXG','DYC','DYG',...
'hFacC','hFacS','hFacW','Depth','AngleCS','AngleSN'};
for iFld=1:length(list0);
eval(['mygrid.' list0{iFld} '=rdmds_compact2llc([dirGrid ''' list0{iFld} '*'']);']);
end;
list0={'RC','RF','DRC','DRF'};
for iFld=1:length(list0);
eval(['mygrid.' list0{iFld} '=rdmds([dirGrid ''' list0{iFld} '*'']);']);
end;
mygrid.hFacCsurf=mygrid.hFacC;
for ff=1:mygrid.hFacC.nFaces; mygrid.hFacCsurf{ff}=mygrid.hFacC{ff}(:,:,1); end;
end;
gcmfaces/gcmfaces_legacy/verticalplot.m 0000644 0004354 0000144 00000004051 11600706724 020172 0 ustar gforget users %22/11/2003
%auteur : Gael Forget
%version : beta
%objet : a la realsiation d'un plot xz ou yz, on utilise
%des echelles differentes au dessus/en dessous de 500m
%remarque : old version etait spacifique a chaque type
% de plot, cette nouvelle se vaut generique...
function [varargout]=verticalplot(typeplot,cur_x,cur_z,varargin);
%variables d'entree :
%1) typeplot est le nom du type de plot (px : pcolor)
%2) cur_x et cur_z sont les chamsp x et z
%3) varargin sont les autres options a passer a la routine de plot
%nouvelles ordonnees :
xplot2=cur_x; yplot2=myverticaldeform(cur_z);
%************************************
%gestion des vars de sortie (partie 1)
tmptxt='';
if nargout==1
tmptxt='[tmp1]=';
elseif nargout>1
tmptxt='[tmp1';
for kkk=2:nargout
tmptxt=[tmptxt ',tmp' num2str(kkk)];
end
tmptxt=[tmptxt ']='];
end
%************************************
%le plot en lui meme :
eval([tmptxt typeplot '(xplot2,yplot2,varargin{:});']);
%************************************
%gestion des vars de sortie (partie 2)
if nargout>=1
for kkk=1:nargout
eval(['varargout(kkk) = {tmp' num2str(kkk) '};'])
end
end
%************************************
%reglage des ytics :
newtick=[floor(min(min(cur_z))) ceil(max(max(cur_z)))];
if (newtick(1)<=500)&(newtick(2)>500)
tmp1=100*floor(newtick(1)/100); tmp2=ceil(newtick(2)/500)*500;
%newtick=[0:1000:6000];
%newtick=[tmp1:50:250]; newtick=[ newtick(1:end) [500:500:tmp2]];
%newtick=[tmp1:100:500]; newtick=[ newtick(1:end-1) [500:500:tmp2]];
newtick=[tmp1:100:1000]; newtick=[ newtick(1:end-1) [1000:500:6000]];
elseif (newtick(1)<=500)
tmp1=floor(newtick(1)/50)*50; tmp2=ceil(newtick(2)/50)*50;
newtick=[tmp1:50:tmp2];
else
tmp1=floor(newtick(1)/500)*500; tmp2=ceil(newtick(2)/500)*500;
newtick=[tmp1:500:tmp2];
end
%on les range dans l'ordre approprie :
newtick=flipdim(newtick,2);
%format texte :
newticklabel=[];
for kkk=1:length(newtick)
newticklabel=strvcat(newticklabel,num2str(newtick(kkk)));
end
%mise a jour des yticks :
newtick=myverticaldeform(newtick);
set(gca,'YTick',newtick);
set(gca,'YTickLabel',newticklabel);
gcmfaces/gcmfaces_legacy/convert_llc2pcol.m 0000644 0004354 0000144 00000000130 11334544171 020726 0 ustar gforget users function [X,Y,fld]=convert_llc2pcol(varargin);
[X,Y,fld]=convert2pcol(varargin{:});
gcmfaces/gcmfaces_legacy/diffsmooth2Dstraight_extrapolate_inverse.m 0000644 0004354 0000144 00000000163 11334544171 025723 0 ustar gforget users function [fld]=diffsmooth2Dstraight_extrapolate_inverse(varargin);
[fld]=diffsmooth2D_extrap_inv(varargin{:});
gcmfaces/gcmfaces_legacy/line_greatC_TUV_mask.m 0000644 0004354 0000144 00000007117 11601141417 021446 0 ustar gforget users function [line_out]=line_greatC_TUV_mask(lonPair_in,latPair_in);
global mygrid;
%get carthesian coordinates:
%... of grid
lon=mygrid.XC; lat=mygrid.YC;
x=cos(lat*pi/180).*cos(lon*pi/180);
y=cos(lat*pi/180).*sin(lon*pi/180);
z=sin(lat*pi/180);
%... and of end points
x0=cos(latPair_in*pi/180).*cos(lonPair_in*pi/180);
y0=cos(latPair_in*pi/180).*sin(lonPair_in*pi/180);
z0=sin(latPair_in*pi/180);
%get the rotation matrix:
%1) rotate around x axis to put first point at z=0
theta=atan2(-z0(1),y0(1));
R1=[[1;0;0] [0;cos(theta);sin(theta)] [0;-sin(theta);cos(theta)]];
tmp0=[x0;y0;z0]; tmp1=R1*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
x0=x1; y0=y1; z0=z1;
%2) rotate around z axis to put first point at y=0
theta=atan2(x0(1),y0(1));
R2=[[cos(theta);sin(theta);0] [-sin(theta);cos(theta);0] [0;0;1]];
tmp0=[x0;y0;z0]; tmp1=R2*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
x0=x1; y0=y1; z0=z1;
%3) rotate around y axis to put second point at z=0
theta=atan2(-z0(2),-x0(2));
R3=[[cos(theta);0;-sin(theta)] [0;1;0] [sin(theta);0;cos(theta)]];
tmp0=[x0;y0;z0]; tmp1=R3*tmp0; x1=tmp1(1,:); y1=tmp1(2,:); z1=tmp1(3,:);
x0=x1; y0=y1; z0=z1;
%apply rotation to grid:
tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
tmp1=find(~isnan(tmpx));
tmpx2=tmpx(tmp1); tmpy2=tmpy(tmp1); tmpz2=tmpz(tmp1);
tmp2=[tmpx2';tmpy2';tmpz2'];
tmp3=R3*R2*R1*tmp2;
tmpx2=tmp3(1,:); tmpy2=tmp3(2,:); tmpz2=tmp3(3,:);
tmpx(tmp1)=tmpx2; tmpy(tmp1)=tmpy2; tmpz(tmp1)=tmpz2;
x=convert2array(tmpx); y=convert2array(tmpy); z=convert2array(tmpz);
%compute the great circle mask:
zz=exch_T_N(z);
mm=zz<=0; mm(isnan(mm))=0;
for iF=1:y.nFaces;
eval(['tmp1=mm.f' num2str(iF) ';']);
tmp2=circshift(tmp1,[-1 0])+circshift(tmp1,[1 0])+...
circshift(tmp1,[0 -1])+circshift(tmp1,[0 1]);
tmp2(tmp1>0)=0; tmp2(tmp2~=0)=1; tmp2(tmp2==0)=NaN;
eval(['mm.f' num2str(iF) '=tmp2(2:end-1,2:end-1);']);
end;
%UV part:
mmt=mm;
nn=exch_T_N(mmt);
mm=zz<=0; mm(isnan(mm))=0;
mmu=mmt; mmv=mmt;
for iFace=1:y.nFaces;
iF=num2str(iFace);
eval(['tmp_u=mm.f' iF '(1:end-1,:)==1&nn.f' iF '(2:end,:)==1;']);
eval(['tmp_u2=mm.f' iF '(2:end,:)==1&nn.f' iF '(1:end-1,:)==1;']);
tmp_u=1*(tmp_u(1:end-1,2:end-1)-tmp_u2(1:end-1,2:end-1)); eval(['mmu.f' iF '=tmp_u;']);
eval(['tmp_v=mm.f' iF '(:,1:end-1)==1&nn.f' iF '(:,2:end)==1;']);
eval(['tmp_v2=mm.f' iF '(:,2:end)==1&nn.f' iF '(:,1:end-1)==1;']);
tmp_v=1*(tmp_v(2:end-1,1:end-1)-tmp_v2(2:end-1,1:end-1)); eval(['mmv.f' iF '=tmp_v;']);
end;
mmu(mmu==0)=NaN; mmv(mmv==0)=NaN;
for kk=1:3;
%select field to treat:
switch kk;
case 1; mm=mmt;
case 2; mm=mmu;
case 3; mm=mmv;
end;
%split in two contours:
theta=[];
theta(1)=atan2(y0(1),x0(1));
theta(2)=atan2(y0(2),x0(2));
tmpx=convert2array(x); tmpy=convert2array(y); tmpz=convert2array(z);
tmptheta=atan2(tmpy,tmpx);
tmpm=convert2array(mm); tmpmIn=tmpm; tmpmOut=tmpmIn;
if theta(2)<0;
tmp00=find(tmptheta<=theta(2)); tmptheta(tmp00)=tmptheta(tmp00)+2*pi;
theta(2)=theta(2)+2*pi;
end;
tmpmIn(find(tmptheta>theta(2)|tmptheta=theta(1)))=NaN;
mmIn=convert2array(tmpmIn); mmOut=convert2array(tmpmOut);
%ensure that we select the shorther segment as mmIn:
if theta(2)-theta(1)>pi; tmp1=mmIn; mmIn=mmOut; mmOut=tmp1; end;
%store result:
switch kk;
case 1; mmtIn=mmIn; mmtOut=mmOut;
case 2; mmuIn=mmIn; mmuOut=mmOut;
case 3; mmvIn=mmIn; mmvOut=mmOut;
end;
end;
%output:
line_out=struct('mmt',mmt,'mmtIn',mmtIn,'mmtOut',mmtOut,...
'mmu',mmu,'mmuIn',mmuIn,'mmuOut',mmuOut,...
'mmv',mmv,'mmvIn',mmvIn,'mmvOut',mmvOut,...
'lonPair',lonPair_in,'latPair',latPair_in);
gcmfaces/ecco_v4/ 0000755 0004354 0000144 00000000000 12375431055 013534 5 ustar gforget users gcmfaces/ecco_v4/cost_seaicearea.m 0000644 0004354 0000144 00000017073 12365532071 017033 0 ustar gforget users function []=cost_seaicearea(dirModel,dirMat,doComp,dirTex,nameTex);
%object: compute cost function term for sea ice data
%inputs: dimodel is the model directory
% dirMat is the directory where diagnozed .mat files will be saved
% -> set it to '' to use the default [dirModel 'mat/']
% doComp is a switch (1->compute; 0->display)
%optional: dirTex is the directory where tex and figures files are created
% (if not specified then display all results to screen instead)
% nameTex is the tex file name (default : 'myPlots')
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); mkdir([dirMat]); end;
%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
if isempty(who('dirTex'));
addToTex=0;
else;
if ~ischar(dirTex); error('mis-specified dirTex'); end;
addToTex=1;
if isempty(who('nameTex')); nameTex='myPlots'; end;
fileTex=[dirTex nameTex '.tex'];
end;
%grid, params and inputs
gcmfaces_global; global myparms;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
if isfield(myparms,'yearsta'); yearsta=myparms.yearsta; yearend=myparms.yearend;
else; yearsta=1992; yearend=2011;
end;
[lonPairs,latPairs,names] = line_greatC_TUV_MASKS_core2_antarctic;
lonLims=[lonPairs(1:5,1);lonPairs(1,1)];
if doComp;
%grid, params and inputs
fld_err=ones(90,1170);
fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;
dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
fileModel=dir([dirModel 'barfiles/gbar_area*data']);
fileModel=['barfiles/' fileModel.name];
nyears=yearend-yearsta+1;
nmonths=12*nyears;
%misfits :
fld_dif=convert2gcmfaces(NaN*ones(90,90*13,nmonths));
fld_nsidc=convert2gcmfaces(NaN*ones(90,90*13,nmonths));
%monthly mean climatology :
climMod=convert2gcmfaces(zeros(90,90*13,12));
climObs=convert2gcmfaces(zeros(90,90*13,12));
climMsk=convert2gcmfaces(zeros(90,90*13,12));
climNb=convert2gcmfaces(zeros(90,90*13,12));
%monthly integrals :
IceAreaNorthMod=NaN*zeros(1,nmonths);
IceAreaNorthObs=NaN*zeros(1,nmonths);
IceAreaSouthMod=NaN*zeros(6,nmonths);
IceAreaSouthObs=NaN*zeros(6,nmonths);
%computational loop :
for ycur=yearsta:yearend;
tic;
for mcur=1:12;
mm=(ycur-yearsta)*12+mcur;
fld_dat=read_bin([dirData 'input_nsidc_all/nsidc79_monthly_' num2str(ycur)],mcur,0);
fld_dat=fld_dat.*mygrid.mskC(:,:,1);%land mask
fld_dat(find(fld_dat<-99))=NaN;%missing data
msk=1+0*fld_dat;%combined mask
fld_mod=read_bin([dirModel fileModel],mm,0);
fld_mod=fld_mod.*msk;%mask consistent with fld_dat
%misfits :
fld_dif(:,:,mm)=fld_mod-fld_dat;
fld_nsidc(:,:,mm)=fld_dat;
%climatology :
tmp1=msk; tmp1(isnan(tmp1))=0; climMsk(:,:,mcur)=climMsk(:,:,mcur)+tmp1;
tmp1=fld_mod; tmp1(isnan(tmp1))=0; climMod(:,:,mcur)=climMod(:,:,mcur)+tmp1;
tmp1=fld_dat; tmp1(isnan(tmp1))=0; climObs(:,:,mcur)=climObs(:,:,mcur)+tmp1;
%integrals :
fld=fld_mod.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthMod(mm)=nansum(fld);
fld=fld_dat.*mygrid.RAC.*(mygrid.YC>0); IceAreaNorthObs(mm)=nansum(fld);
fld=fld_mod.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthMod(1,mm)=nansum(fld);
fld=fld_dat.*mygrid.RAC.*(mygrid.YC<0); IceAreaSouthObs(1,mm)=nansum(fld);
for kk=1:5;
tmpmsk=0.*mygrid.XC;
if lonLims(kk+1) > lonLims(kk)
tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC < lonLims(kk+1)))=1.;
else
tmpmsk(find(mygrid.XC >= lonLims(kk) & mygrid.XC <= 180.))=1.;
tmpmsk(find(mygrid.XC >= -180. & mygrid.XC < lonLims(kk+1)))=1.;
end
tmpmsk=tmpmsk.*(mygrid.YC<0);
%
fld=fld_mod.*mygrid.RAC.*tmpmsk;
IceAreaSouthMod(kk+1,mm)=nansum(fld);
%
fld=fld_dat.*mygrid.RAC.*tmpmsk;
IceAreaSouthObs(kk+1,mm)=nansum(fld);
end
end;
toc;
end;
%misfits :
mis_rms=sqrt(nanmean(fld_dif.^2,3));
obs_std=nanstd(fld_nsidc,[],3);
mod_std=nanstd(fld_nsidc+fld_dif,[],3);
%climatology :
for mcur=1:12;
tmp1=climMsk(:,:,mcur); tmp1(tmp1==0)=NaN;
climNb(:,:,mcur)=tmp1;
climMod(:,:,mcur)=climMod(:,:,mcur)./tmp1;
climObs(:,:,mcur)=climObs(:,:,mcur)./tmp1;
end;
clear climMsk;
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_seaicearea.mat fld_err mis_rms obs_std mod_std IceArea* clim*;']);
else;%display previously computed results
if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
eval(['load ' dirMat '/cost_seaicearea.mat;']);
%variance maps:
figure; m_map_gcmfaces(mis_rms,0,{'myCaxis',[0:0.1:1.]});
myCaption={'modeled-observed rms -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(obs_std,0,{'myCaxis',[0:0.1:1.]});
myCaption={'observed std -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(mod_std,0,{'myCaxis',[0:0.1:1.]});
myCaption={'modelled std -- sea ice concentration'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
%arctic maps
figureL;
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),2,{'myCaxis',[0:0.1:1.]});
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),2,{'myCaxis',[0:0.1:1.]});
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
%southern ocean maps
figureL;
subplot(2,2,1); m_map_gcmfaces(climMod(:,:,3),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,2); m_map_gcmfaces(climObs(:,:,3),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,3); m_map_gcmfaces(climMod(:,:,9),3,{'myCaxis',[0:0.1:1.]});
subplot(2,2,4); m_map_gcmfaces(climObs(:,:,9),3,{'myCaxis',[0:0.1:1.]});
myCaption={'ECCO (left) and NSIDC (right, gsfc bootstrap) ice concentration in March (top) and September (bottom).'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
%northern/southern integrals
ny=length(IceAreaNorthMod)/12;
yy=yearsta+[0:ny-1];
figureL;
subplot(1,2,1);
plot(yy,IceAreaNorthMod(3:12:end),'LineWidth',2); hold on;
plot(yy,IceAreaNorthObs(3:12:end),'r','LineWidth',2);
plot(yy,IceAreaNorthMod(9:12:end),'LineWidth',2);
plot(yy,IceAreaNorthObs(9:12:end),'r','LineWidth',2);
axis([yearsta yearsta+ny-1 0 20e12]);
ylabel('m^2'); title('Northern Hemisphere');
subplot(1,2,2);
plot(yy,IceAreaSouthMod(1,3:12:end),'LineWidth',2); hold on;
plot(yy,IceAreaSouthObs(1,3:12:end),'r','LineWidth',2);
plot(yy,IceAreaSouthMod(1,9:12:end),'LineWidth',2);
plot(yy,IceAreaSouthObs(1,9:12:end),'r','LineWidth',2);
axis([yearsta yearsta+ny-1 0 20e12]);
ylabel('m^2'); title('Southern Hemisphere');
myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in March and September',...
'in Northern Hemisphere (left) and Southern Hemisphere (right)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
%southern basin integrals
for mm=[3 9];
figureL;
for ii=1:6;
subplot(3,2,ii);
if ii==1; iiTxt='Entire Southern Ocean';
else; iiTxt=sprintf('%dE to %dE',lonLims(ii-1),lonLims(ii));
end;
plot(yy,IceAreaSouthMod(ii,mm:12:end),'LineWidth',2); hold on;
plot(yy,IceAreaSouthObs(ii,mm:12:end),'r','LineWidth',2);
aa=axis; aa(1:2)=[yearsta yearsta+ny-1]; axis(aa);
ylabel('m^2'); title(iiTxt);
end;
if mm==3; mmTxt='March'; elseif mm==9; mmTxt='September'; else; '???'; end;
myCaption={'ECCO (blue) and NSIDC (red, gsfc bootstrap) ice concentration in ',mmTxt,' per Southern Ocean sector'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
end;
gcmfaces/ecco_v4/cost_altimeter_disp.m 0000644 0004354 0000144 00000012665 12365532071 017760 0 ustar gforget users function []=cost_altimeter_disp(dirMat,choicePlot,suf,dirTex,nameTex);
%object: plot the various sea level statistics
% (std model-obs, model, obs, leading to cost function terms)
%inputs: dirMat is the model run directory
% choicePlot is 1 (rms) 2 (prior uncertainty) or 3 (cost)
% suf is 'modMobs', 'obs' or 'mod'
%optional: dirTex is the directory where tex and figures files are created
% (if not specified then display all results to screen instead)
% nameTex is the tex file name (default : 'myPlots')
gcmfaces_global;
%backward compatibility test
test1=~isempty(dir([dirMat 'basic_diags_ecco_mygrid.mat']));
test2=~isempty(dir([dirMat 'diags_grid_parms.mat']));
if ~test1&~test2;
error('missing diags_grid_parms.mat')
elseif test2;
nameGrid='diags_grid_parms.mat';
else;
nameGrid='basic_diags_ecco_mygrid.mat';
end;
%here we always reload the grid from dirMat to make sure the same one is used throughout
eval(['load ' dirMat nameGrid ';']);
%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
if isempty(who('dirTex'));
addToTex=0;
else;
if ~ischar(dirTex); error('mis-specified dirTex'); end;
addToTex=1;
if isempty(who('nameTex')); nameTex='myPlots'; end;
fileTex=[dirTex nameTex '.tex'];
end;
%%%%%%%%%%%%%%%
%define pathes:
%%%%%%%%%%%%%%%
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);
%%%%%%%%%%%%%%%%%
%do computations:
%%%%%%%%%%%%%%%%%
if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
eval(['load ' dirMat 'cost_altimeter_' suf '.mat myflds;']);
%mask for plotting
if isfield(myflds,'msk_100pts');
MSK=myflds.msk_100pts;
else;
MSK=mygrid.mskC(:,:,1);
end;
%hack to omit any masking : MSK(:)=1;
if strcmp(suf,'modMobs'); tit='modeled-observed';
elseif strcmp(suf,'obs'); tit='observed';
elseif strcmp(suf,'mod'); tit='modeled';
else; error('unknown field');
end
if choicePlot==1; tit=[tit ' log(variance)']; uni='(m$^2$)';
elseif choicePlot==2; tit=' log(prior error variance)'; uni='(m$^2$)';
else; tit=[tit ' cost']; uni='';
end;
if choicePlot==1;%rms
if strcmp(suf,'modMobs');
cc=[-0.4:0.05:-0.25 -0.2:0.03:-0.05 -0.03:0.01:0.03 0.05:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(100*MSK.*myflds.dif_mdt,0,{'myCaxis',100*cc}); drawnow;
myCaption={'mean dynamic topography misfit (cm)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
cc=[-4:0.2:-1]-0.5;
figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
cc=[-4:0.2:-1];
figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
if 0;
figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_point35d.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- large time scales'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_35dMsmooth.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- point35d minus lsc'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(log10(MSK.*(myflds.rms_sladiff_pointMpoint35d.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- point minus point35d'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
elseif choicePlot==2;%uncertainty fields
cc=[0:0.005:0.02 0.03:0.01:0.05 0.06:0.02:0.1 0.14:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(100*myflds.sig_mdt,0,{'myCaxis',100*cc}); drawnow;
myCaption={'mean dynamic topography prior uncertainty (cm)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
cc=[-4:0.2:-1]-0.5;
figure; m_map_gcmfaces(log10(MSK.*(myflds.sig_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
cc=[-4:0.2:-1];
figure; m_map_gcmfaces(log10(MSK.*(myflds.sig_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
else;%cost
cc=[0:0.005:0.02 0.03:0.01:0.05 0.06:0.02:0.1 0.14:0.03:0.2 0.25:0.05:0.4]*100;
figure; m_map_gcmfaces(MSK.*((myflds.dif_mdt.^2)./(myflds.sig_mdt.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- mean dynamic topography ',uni};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(MSK.*((myflds.rms_sladiff_smooth.^2)./(myflds.sig_sladiff_smooth.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- large space/time scales'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure; m_map_gcmfaces(MSK.*((myflds.rms_sladiff_point.^2)./(myflds.sig_sladiff_point.^2)),0,{'myCaxis',cc}); drawnow;
myCaption={tit,'-- sea level anomaly ',uni,' -- pointwise'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
gcmfaces/ecco_v4/cost_xx.m 0000644 0004354 0000144 00000007066 12365532071 015411 0 ustar gforget users function []=cost_xx(dirModel,dirMat,doComp,dirTex,nameTex);
%object: compute cost function term for atmospheric controls
%inputs: dirModel is the model directory
% dirMat is the directory where diagnozed .mat files will be saved
% -> set it to '' to use the default [dirModel 'mat/']
% doComp is a switch (1->compute; 0->display)
%optional: dirTex is the directory where tex and figures files are created
% (if not specified then display all results to screen instead)
% nameTex is the tex file name (default : 'myPlots')
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); mkdir([dirMat]); end;
%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
if isempty(who('dirTex'));
addToTex=0;
else;
if ~ischar(dirTex); error('mis-specified dirTex'); end;
addToTex=1;
if isempty(who('nameTex')); nameTex='myPlots'; end;
fileTex=[dirTex nameTex '.tex'];
end;
for ii=1:6;
switch ii;
case 1; xxName='atemp'; sigName='cap_sigma_tmp2m_degC_eccollc.bin'; cc=2; uni='K';
case 2; xxName='aqh'; sigName='cap_sigma_spfh2m_eccollc.bin'; cc=2; uni='g/kg';
case 3; xxName='tauu'; sigName='cap_sigma_ustr_eccollc.bin'; cc=0.04; uni='N/m2';
case 4; xxName='tauv'; sigName='cap_sigma_vstr_eccollc.bin'; cc=0.04; uni='N/m2';
case 5; xxName='lwdown'; sigName='cap_sigma_dlw_eccollc.bin'; cc=20; uni='W/m2';
case 6; xxName='swdown'; sigName='cap_sigma_dsw_eccollc.bin'; cc=40; uni='W/m2';
end;
if doComp;
%load grid
gcmfaces_global;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
dirSig='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_all_from_pleiades/';
%read model cost output
tmp1=dir([dirModel 'ADXXfiles/xx_' xxName '.effective.*data']);
tmp2=size(convert2gcmfaces(mygrid.XC));
fld_xx=read2memory([dirModel 'ADXXfiles/' tmp1.name],[tmp2 tmp1.bytes/tmp2(1)/tmp2(2)/4]);
fld_xx=convert2gcmfaces(fld_xx);
%does not work when adjoint overwrites xx*effective*meta fld_xx=rdmds2gcmfaces([dirModel 'ADXXfiles/xx_' xxName '.effective.*']);
fld_sig=read_bin([dirSig sigName],1,0);
if strcmp(xxName,'aqh'); fld_xx=fld_xx*1000; fld_sig=fld_sig*1000; end;
%compute xx stats
fld_rms=sqrt(mean(fld_xx.^2,3));
fld_mean=mean(fld_xx,3);
fld_std=std(fld_xx,[],3);
%mask
fld_rms=fld_rms.*mygrid.mskC(:,:,1);
fld_sig=fld_sig.*mygrid.mskC(:,:,1);
fld_mean=fld_mean.*mygrid.mskC(:,:,1);
fld_std=fld_std.*mygrid.mskC(:,:,1);
clear fld_xx;
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_xx_' xxName '.mat fld_* cc uni;']);
else;%display previously computed results
global mygrid;
if ~isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
eval(['load ' dirMat '/cost_xx_' xxName '.mat;']);
figure;
m_map_gcmfaces(fld_sig,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['prior uncertainty -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure;
m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['rms adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure;
m_map_gcmfaces(fld_std,0,{'myCaxis',[0:0.05:0.5 0.6:0.1:1 1.25]*cc});
myCaption={['std adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figure;
m_map_gcmfaces(fld_mean,0,{'myCaxis',[-0.5:0.05:0.5]*cc});
myCaption={['mean adjustment -- ' xxName ' (' uni ')']};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
end;%for ii=1:6;
gcmfaces/ecco_v4/cost_read.m 0000644 0004354 0000144 00000004504 12103041413 015640 0 ustar gforget users function [cost]=cost_read(dir0,file0,list0);
%object: read a cost function file such as costfunction0000
%input: dir0 is the model run directory name
% file0 (optional) is the cost function file name (costfunction0000 or
% or costfunction0001 etc ). If file0 is not specified,
% of specified as '' then it will be automatically be set
% by inspecting dir0 (assuming there is only one in dir0)
% list0 (optional) is a cell array of cost term names
%output : if list0 is NOT specified, then cost is a structure vector
% of cost(kk).fc, cost(kk).no, cost(kk).fc
% if list0 is specified, then cost is the vector of cost terms
if isempty(who('file0')); file0=''; end;
if isempty(who('list0')); list0=''; end;
%take care of file name:
dir0=[dir0 '/'];
if isempty(file0);
file0=dir([dir0 'costfunction0*']);
if length(file0)>1&nargin>1;
fprintf('several costfunction0??? files were found:\n');
{file0(:).name}'
error('please be more specific');
end;
file0=file0.name;
end;
file0=[dir0 file0];
%read costfunction0??? file:
fid=fopen(file0); tmp2='';
while 1;
tline = fgetl(fid);
if ~ischar(tline), break, end
if isempty(tmp2); tmp2=[tline ' ; ']; else; tmp2=[tmp2 ' ' tline ' ; ']; end;
end
fclose(fid);
%get cost from file (already in text form in memory):
cost=[]; cost2=[]; kk=0;
tmp1=[-3 strfind(tmp2,' ; ')];
for ii=1:length(tmp1)-1;
tmp3=tmp2(tmp1(ii)+4:tmp1(ii+1));
jj=strfind(tmp3,'='); jj=jj(end)+1;
tmp_name=tmp3(1:jj-2); tmp_val=tmp3(jj:end-1);
%
tmp_val(strfind(tmp_val,'D'))='e'; tmp_val(strfind(tmp_val,'+'))='';
eval(['tmp_val=[ ' tmp_val ' ];']); tmp_val(tmp_val==0)=NaN;
%
if ~isnan(tmp_val(1));
kk=kk+1;
if kk==1;
cost.name=strtrim(tmp_name); cost.fc=tmp_val(1); cost.no=tmp_val(2);
else;
cost(kk).name=strtrim(tmp_name); cost(kk).fc=tmp_val(1); cost(kk).no=tmp_val(2);
end;
end;
end;
%output vector of cost?
if ~isempty(list0);
cost_bak=cost;
nn=length(list0);
cost=NaN*zeros(nn,1);
list1={cost_bak(:).name};
for kk=1:nn;
jj=find(strcmp(list1,list0{kk}));
if length(jj)>1;
error([list0{kk} ' is not specific enough']);
elseif length(jj)==1;
%cost(kk,:)=[cost_bak(jj).fc cost_bak(jj).no];
cost(kk)=cost_bak(jj).fc;
end;
end;
end;
gcmfaces/ecco_v4/v4_basin.m 0000644 0004354 0000144 00000005534 12365534162 015430 0 ustar gforget users function [varargout]=v4_basin(nameBasin,varargin);
%object: obtain the mask of an ocean basin
%
%input: nameBasin name of the basin of interest (atl, pac, ind, arctic, etc.)
%optional: msk0 value for the masked region (0 by default)
%
%output: mskC mask for tracer points (msk0=outside basin; 1=inside basin)
%optional: mskW,mskS mask for velocity points (0=oustide/land; 1/2=edge; 1=inside)
if nargin==2; msk0=varargin{1}; else; msk0=0; end;
gcmfaces_global;
dir0=which('gcmfaces_demo');
tmp1=strfind(dir0,filesep);
dir0=[dir0(1:tmp1(end)) 'sample_input/OCCAetcONv4GRID/'];
msk=read_bin([dir0 'basin_masks_eccollc_90x50.bin'],0,1);
%list of available basins:
list0={ 'pac','atl','ind','arct','bering',...
'southChina','mexico','okhotsk','hudson','med',...
'java','north','japan','timor','eastChina','red','gulf',...
'baffin','gin','barents'};
%list of selected basins:
if ischar(nameBasin); nameBasin={nameBasin}; end;
list1={};
for ii=1:length(nameBasin);
if strcmp(nameBasin{ii},'atlExt');
list1={list1{:},'atl','mexico','hudson','med','north','baffin','gin'};
elseif strcmp(nameBasin{ii},'pacExt');
list1={list1{:},'pac','bering','okhotsk','japan','eastChina'};
elseif strcmp(nameBasin{ii},'indExt');
list1={list1{:},'ind','southChina','java','timor','red','gulf'};
else;
list1={list1{:},nameBasin{ii}};
end;
end
%derive tracer points mask:
mskC=0*msk;
for ii=1:length(list1);
jj=find(strcmp(list1{ii},list0));
if ~isempty(jj); mskC(find(msk==jj))=1; end;
end;
%determine velocity points masks, if needed:
if nargout>1;
%flag velocity points according to neighboring pair:
fld=3*mskC+1*(~isnan(mygrid.mskC(:,:,1)));
FLD=exch_T_N(fld);
fldW=fld; fldS=fld;
for iF=1:FLD.nFaces;
tmpA=FLD{iF}(2:end-1,2:end-1);
tmpB=FLD{iF}(1:end-2,2:end-1);
fldW{iF}=(tmpA+tmpB)/2;
tmpA=FLD{iF}(2:end-1,2:end-1);
tmpB=FLD{iF}(2:end-1,1:end-2);
fldS{iF}=(tmpA+tmpB)/2;
end;
%compute corresponding masks:
mskW=0*mskC;
mskW(find(fldW==4))=1;%inside points
mskW(find(fldW==2.5))=0.5;%basin edge points
mskS=0*mskC;
mskS(find(fldS==4))=1;%inside points
mskS(find(fldS==2.5))=0.5;%basin edge points
%for checking:
if 0;
mskWout=0*mskC;
mskWout(find(fldW==1))=1;%outside points
mskWout(find(fldW==2.5))=0.5;%basin edge points
mskSout=0*mskC;
mskSout(find(fldS==1))=1;%outside points
mskSout(find(fldS==2.5))=0.5;%basin edge points
end;
end;
%replace 0 with msk0:
mskC(find(mskC==0))=msk0;
if nargout>1; mskW(find(mskW==0))=msk0; mskS(find(mskS==0))=msk0; end;
%output(s):
if nargout==1; varargout={mskC}; else; varargout={mskC,mskW,mskS}; end;
%for checking:
if 0;
figure;
msk0=1*(msk0>0); msk0(find(msk0==0))=NaN;
subplot(2,1,1); imagescnan(convert2array(msk0)'); axis xy; caxis([-1 2]);
subplot(2,1,2); imagescnan(convert2array(mskC.*msk0)'); axis xy; caxis([-1 2]);
drawnow;
end;
gcmfaces/ecco_v4/read_grace_fld.m 0000644 0004354 0000144 00000000730 12161145421 016603 0 ustar gforget users function [fld,fldsm]=read_grace_fld(fname);
%[fld]=read_grace_fld('GRACE_CSR_Error.asc');
fid=fopen(fname,'rt');
fgetl(fid);
fgetl(fid);
fld4columns=zeros(1e5,4);
ii=1;
while ~feof(fid);
fld4columns(ii,:)=str2num(fgetl(fid));
ii=ii+1;
end;
fclose(fid);
ii=find(fld4columns(:,2)>0&fld4columns(:,2)<360);
jj=(fld4columns(ii,1)+89.5)*360+fld4columns(ii,2)+0.5;
fld=NaN*zeros(360,180);
fld(jj)=fld4columns(ii,3);
fldsm=NaN*zeros(360,180);
fldsm(jj)=fld4columns(ii,4);
gcmfaces/ecco_v4/cost_summary.m 0000644 0004354 0000144 00000012050 11625547244 016442 0 ustar gforget users function []=cost_summary(dir0,file0,cost0);
%object: display summary of the various cost terms from costfunction0??? file
%input: dir0 is the model run directory name
% file0 is the costfunction0??? file name (if '' then it will
% be detected -- assuming there is only one in dir0)
% cost0 is a reference value that will be used for normalization (1e8 by default)
%note: hard coded ... mainly because multipliers are not in costfunction0??? file
if isempty(cost0); cost0=1e8; end;
%take care of file name:
dir0=[dir0 '/'];
if isempty(file0);
file0=dir([dir0 'costfunction0*']);
if length(file0)>1&nargin>1;
fprintf('several costfunction0??? files were found:\n');
{file0(:).name}'
error('please be more specific');
end;
file0=file0.name;
end;
file0=[dir0 file0];
%read costfunction0??? file:
fid=fopen(file0); tmp2='';
while 1;
tline = fgetl(fid);
if ~ischar(tline), break, end
if isempty(tmp2); tmp2=[tline ' ; ']; else; tmp2=[tmp2 ' ' tline ' ; ']; end;
end
fclose(fid);
%list of cost function terms and multipliers: (after apr1alpha/part2)
mylist={'fc','f_temp','f_salt','f_sst','f_tmi',...
'argo_pacific_MITprof_latest_4mygridprof_T','argo_pacific_MITprof_latest_4mygridprof_S',...
'argo_atlantic_MITprof_latest_4mygridprof_T','argo_atlantic_MITprof_latest_4mygridprof_S',...
'argo_indian_MITprof_latest_4mygridprof_T','argo_indian_MITprof_latest_4mygridprof_S',...
'seals_MITprof_latest_4mygridprof_T','seals_MITprof_latest_4mygridprof_S',...
'WOD09_CTD_4mygridprof_T','WOD09_CTD_4mygridprof_S','WOD09_XBT_4mygridprof_T',...
'area','area2sst','sshv4-mdt','sshv4-lsc',...
'sshv4-tp','sshv4-ers','sshv4-gfo',...
'sstv4-amsre-lsc','sstv4-amsre'};
mymult=[[0 0.15 0.15 2 2] [2 2 2 2 2 2 1 1 2 2 2] [0 50 200 40] 0.25*[1 1 1] 25 1];
%before apr1alpha part 2
%mylist={'fc','f_temp','f_salt','f_sst','f_tmi',...
% 'argo_pacific_MITprof_latest_4mygridprof_T','argo_pacific_MITprof_latest_4mygridprof_S',...
% 'argo_atlantic_MITprof_latest_4mygridprof_T','argo_atlantic_MITprof_latest_4mygridprof_S',...
% 'argo_indian_MITprof_latest_4mygridprof_T','argo_indian_MITprof_latest_4mygridprof_S',...
% 'seals_MITprof_latest_4mygridprof_T','seals_MITprof_latest_4mygridprof_S','XBT_v5_4mygridprof_T',...
% 'area','area2sst','sshv4-mdt','sshv4-lsc','sstv4-amsre-lsc',...
% 'sshv4-tp','sshv4-ers','sshv4-gfo','sstv4-amsre'};
%mymult0=[[0 0.5 0.5 1 1] [2 2 2 2 2 2 1 1 2] [2 50 10 10 10] [1 1 1 1]];%during no1beta part 4
%mymult1=[[0 0.15 0.15 1 1] [2 2 2 2 2 2 1 1 2] [1 50 100 20 50] 0.25*[1 1 1 1]];%during nov1beta part 5
%mymult=[[0 0.15 0.15 1 1] [2 2 2 2 2 2 1 1 2] [1 50 100 20 50] 0.25*[1 1 1 1]];%during apr1alpha/part1
%get cost from file (already in text form in memory):
mycost=[]; mycost2=[]; kk=0;
tmp1=[-3 strfind(tmp2,' ; ')];
for ii=1:length(tmp1)-1;
tmp3=tmp2(tmp1(ii)+4:tmp1(ii+1));
jj=strfind(tmp3,'='); jj=jj(end)+1;
tmp_name=tmp3(1:jj-2); tmp_val=tmp3(jj:end-1);
%
jj=strfind(tmp_name,'gencost'); if ~isempty(jj); tmp_name=deblank(tmp_name(1:jj-2)); end;
jj=find(~isspace(tmp_name)); tmp_name=tmp_name(jj);
%if isempty(strfind(tmp_name,' prof_')); tmp_name=deblank(tmp_name); end;
%
tmp_val(strfind(tmp_val,'D'))='e'; tmp_val(strfind(tmp_val,'+'))='';
eval(['tmp_val=[ ' tmp_val ' ];']); tmp_val(tmp_val==0)=NaN;
%
if ~isnan(tmp_val(1))&sum(strcmp(tmp_name,mylist))>0;
tmp_mult=mymult(find(strcmp(tmp_name,mylist)));
kk=kk+1;
if kk==1;
mycost.name={tmp_name}; mycost.fc=tmp_val(1); mycost.no=tmp_val(2); mycost.mult=tmp_mult;
mycost2.name=tmp_name; mycost2.fc=tmp_val(1); mycost2.no=tmp_val(2); mycost2.mult=tmp_mult;
else;
mycost.name(kk)={tmp_name}; mycost.fc(kk)=tmp_val(1); mycost.no(kk)=tmp_val(2); mycost.mult(kk)=tmp_mult;
mycost2(kk).name=tmp_name; mycost2(kk).fc=tmp_val(1); mycost2(kk).no=tmp_val(2); mycost2(kk).mult=tmp_mult;
end;
end;
end;
%check that I recover the total:
tmp1=100*sum(mycost.fc.*mycost.mult/mycost.fc(1));
fprintf('offline/online cost ratio is %3.3f%%\n',tmp1);
%share for each cost term:
myshare=mycost.fc.*mycost.mult; myshare=100*myshare/sum(myshare);
%groups of cost terms:
%myshare2=[sum(myshare(2:3)) sum(myshare([4 5 21 22])) sum(myshare(14:15)) sum(myshare([16:20])) sum(myshare(6:13))]';
fprintf('%15s contributes %3.2f%% \n','T/S atlas',sum(myshare(2:3)));
fprintf('%15s contributes %3.2f%% \n','SST',sum(myshare([4 5 24 25])));
fprintf('%15s contributes %3.2f%% \n','ice conc',sum(myshare(17:18)));
fprintf('%15s contributes %3.2f%% \n','SSH',sum(myshare(19:23)));
fprintf('%15s contributes %3.2f%% \n','in situ',sum(myshare(6:16)));
%plot the various contributions:
% myplot=myshare; aa=[0 1.5 0 length(myplot)+1];
% figure; barh(myshare); axis(aa); grid on; hold on; title('cost terms contribution in %');
myplot=mycost.fc.*mycost.mult/cost0; aa=[0 10 0 length(myplot)+1];
myplot(1)=mycost.fc(1)/cost0; %to plot fc despite mult=0
figure; barh(myplot); axis(aa); grid on; hold on; title(sprintf('cost terms divided by %0.3g',cost0));
for ii=1:length(myplot); xx=myplot(ii)+0.01; yy=ii-0.4; text(xx,yy,mycost.name(ii),'Interpreter','none'); end;
gcmfaces/ecco_v4/gcmfaces_remap.m 0000644 0004354 0000144 00000004331 11747566713 016663 0 ustar gforget users function []=gcmfaces_remap(dirIn,fileIn,gridIn,dirOut,fileOut);
%object: use bin average to remap a lat-lon grid product to a gcmfaces grid
%inputs: dirIn is the input directory
% fileIn is the input file name without the final four year characters
% (e.g. 'SST_monthly_r2_' to process 'SST_monthly_r2_1992' etc.)
% gridIn states the originating grid. It can be set to
% 1 implying that the grid is [0.5:359.5] & [-89.5:89.5]
% 2 implying that the grid is [0.5:359.5] & [-79.5:79.5]
% 3 implying that the grid is [.125:.25:360-.125] & [-90+.125:.25:90-.125]
% {x,y} where x and y are the position arrays
% dirOut and filOut are the corresponding output names
%
%assumption: mygrid has been read using grid_load
% for the input grid, lon must be 0-360 (see gcmfaces_remap_2d)
global mygrid mytri;
%create triangulation
gcmfaces_bindata;
%list files to be processed
listFiles=dir([dirIn fileIn '*']);
%case of user defined grid
if iscell(gridIn); x=gridIn{1}; y=gridIn{2}; [nx,ny]=size(x); mis=0;
else;
%standard cases
if gridIn==1;%gloabl 1 degree grid
x=[0.5:359.5]'*ones(1,180);
y=ones(360,1)*[-89.5:89.5];
[nx,ny]=size(x);
mis=0;
elseif gridIn==2;%ECCO 1 degree grid
x=[0.5:359.5]'*ones(1,160);
y=ones(360,1)*[-79.5:79.5];
[nx,ny]=size(x);
mis=0;
elseif gridIn==3;%1/4 degree grid (e.g. REMSS)
x=[.125:.25:360-.125]'*ones(1,180*4);
y=ones(360*4,1)*[-90+.125:.25:90-.125];
[nx,ny]=size(x);
mis=-9999;
else;
error('unknown grid');
end;
end;
%process one file after the other
for ii=1:length(listFiles);
yy=listFiles(ii).name(end-3:end); fprintf(['processing ' fileIn ' for year ' yy '\n']);
nt=listFiles(ii).bytes/nx/ny/4; if round(nt)~=nt; error('inconsistent sizes'); end;
%read data
fld=reshape(read2memory([dirIn listFiles(ii).name],[nx*ny*nt 1]),[nx ny nt]);
%mask land
fld(fld==mis)=NaN;
%map to v4 grid
FLD=convert2array(zeros(360,360,nt));
for tt=1:nt; FLD(:,:,tt)=gcmfaces_remap_2d(x,y,fld(:,:,tt),3); end;
%set missing value
FLD(find(isnan(FLD)))=mis;
%write data
write2file([dirOut fileOut yy],convert2gcmfaces(FLD));
end;
gcmfaces/ecco_v4/rads_noice_mad_plot.m 0000644 0004354 0000144 00000007261 12151726162 017703 0 ustar gforget users function [zBefore,zAfter]=rads_noice_mad_plot(choiceData);
if isempty(who('doTesting')); doTesting=0; end;
doPlot=0;
%directories
dirMad='/net/nares/raid10/gforget/2013mayInputs/rads/rads_ice_not_flagged/'
dirRads='/bay/scratch2/RADS_alongtrk_v4_ice_not_flagged_2013/'
dirOutput='./';
topexfile = 'tj_daily_ssh_v4_r1';
ersfile = 'en_daily_ssh_v4_r1';
gfofile = 'g1_daily_ssh_v4_r1';
dirIce='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_nsidc_all/';
fileIce='nsidc79_daily';
nday=zeros(20,1);
for yy=1992:2012;
tmp1=dir([dirRads topexfile '_' num2str(yy)]);
nday(yy-1991)=tmp1.bytes/90/1170/4;
end;
years=[1992:2012];
dayMax=datenum(years(end)+1,1,1)-datenum(years(1),1,1);
gcmfaces_global;
eval(['fileData=' choiceData ';']);
slaBefore=zeros(90*1170,dayMax);
slaAfter=slaBefore;
for yy=years;
fprintf(['reading ' fileData '_' num2str(yy) '\n']);
day1=1+datenum(yy,1,1)-datenum(years(1),1,1);
dayN=datenum(yy+1,1,1)-datenum(years(1),1,1);
%
fld=read2memory([dirRads fileData '_' num2str(yy)],[90*1170 1+dayN-day1]);
fld(fld<-999)=NaN;
slaBefore(:,day1:dayN)=fld;
%
fld=read2memory([dirMad fileData '_mad_' num2str(yy)],[90*1170 1+dayN-day1]);
fld(fld<-999)=NaN;
slaAfter(:,day1:dayN)=fld;
end;
%global+time mean stats :
cntNbs=[sum(~isnan(slaBefore(:))) sum(~isnan(slaAfter(:)))];
cntNbs(3)=100*(1-cntNbs(2)/cntNbs(1));
stdNbs=[nanstd(slaBefore(:)) nanstd(slaAfter(:))];
stdNbs(3)=100*(1-stdNbs(2)/stdNbs(1));
fid = fopen([dirMad fileData '_mad_stats.txt'],'w');
fprintf(fid,' instr : %s \n',fileData);
fprintf(fid,' # samples : %d \n',cntNbs(1));
fprintf(fid,' %% mad reject : %3.2f \n',cntNbs(3));
fclose(fid);
%global mean stats :
stdBefore=nanstd(slaBefore,[],1);
stdAfter=nanstd(slaAfter,[],1);
cntBefore=sum(~isnan(slaBefore),1);
cntAfter=sum(~isnan(slaAfter),1);
stdBefore=runmean(stdBefore,17,2);
stdAfter=runmean(stdAfter,17,2);
cntBefore=runmean(cntBefore,17,2);
cntAfter=runmean(cntAfter,17,2);
%zonal mean stats :
msk=mygrid.mskC(:,:,1);
z=nanstd(slaBefore,[],2); z=convert2gcmfaces(reshape(z,[90 1170]));
zBefore=z;%output
slaBeforeStd=calc_zonmean_T(msk.*z);
z=nanstd(slaAfter,[],2); z=convert2gcmfaces(reshape(z,[90 1170]));
zAfter=z;%output
slaAfterStd=calc_zonmean_T(msk.*z);
z=sum(~isnan(slaBefore),2); z=convert2gcmfaces(reshape(z,[90 1170]));
slaBeforeCnt=calc_zonmean_T(msk.*z);
z=sum(~isnan(slaAfter),2); z=convert2gcmfaces(reshape(z,[90 1170]));
slaAfterCnt=calc_zonmean_T(msk.*z);
%plotting:
figureL;
%
x=mygrid.LATS;
subplot(3,2,1);
plot(x,slaBeforeStd); hold on; plot(x,slaAfterStd,'r'); axis([-90 90 0 30]);
title('std before (b) and after (r)'); ylabel('cm'); xlabel('lat');
subplot(3,2,3);
plot(x,slaBeforeCnt); hold on; plot(x,slaAfterCnt,'r'); axis([-90 90 0 1e3])
title('count before (b) and after (r)'); ylabel('# obs'); xlabel('lat');
subplot(3,2,5);
plot(x,100*(1-slaAfterStd./slaBeforeStd),'k'); hold on;
plot(x,100*(1-slaAfterCnt./slaBeforeCnt),'m'); axis([-90 90 0 20])
title('std (k) and count (m) reduction : 100*(1-after/before)'); ylabel('%'); xlabel('lat');
%
x=[1:dayMax];
subplot(3,2,2);
plot(x,stdBefore); hold on; plot(x,stdAfter,'r'); axis([0 8e3 0 20]);
title('std before (b) and after (r)'); ylabel('cm'); xlabel('lat');
subplot(3,2,4);
plot(x,cntBefore); hold on; plot(x,cntAfter,'r'); axis([0 8e3 0 2e4]);
title('count before (b) and after (r)'); ylabel('# obs'); xlabel('lat');
subplot(3,2,6);
plot(x,100*(1-stdAfter./stdBefore),'k'); hold on;
plot(x,100*(1-cntAfter./cntBefore),'m'); axis([0 8e3 0 20]);
title('std (k) and count (m) reduction : 100*(1-after/before)'); ylabel('%'); xlabel('lat');
%
saveas(gcf,[dirMad fileData '_mad_stats'],'fig');
eval(['print -djpeg90 ' dirMad fileData '_mad_stats.jpg']);
gcmfaces/ecco_v4/v4_read_data.m 0000644 0004354 0000144 00000001113 12365534162 016225 0 ustar gforget users function [fld]=v4_read_data(fileName,irec);
%usage: fld=v4_read_data(fileName,irec); 'fast' read of 2D fields no irec in all [fileName '*.data']
gcmfaces_global;
dir0=strfind(fileName,filesep); if isempty(dir0); dir0='./'; else; dir0=fileName(1:dir0(end)); end;
fileList=dir([fileName '.data']);
nn=length(fileList);
fld=zeros(90,1170,nn);
for ii=1:nn;
fid_cur=fopen([dir0 fileList(ii).name],'r','b');
recl=90*1170*4; position0=recl*(irec-1);
status=fseek(fid_cur,position0,'bof');
fld(:,:,ii)=fread(fid_cur,[90 1170],'float32');
fclose(fid_cur);
end;
fld=convert2gcmfaces(fld);
gcmfaces/ecco_v4/cost_altimeter.m 0000644 0004354 0000144 00000025232 12365532071 016733 0 ustar gforget users function []=cost_altimeter(dirModel,dirMat);
%object: compute or plot the various sea level statistics
% (std model-obs, model, obs, leading to cost function terms)
%inputs: dirModel is the model run directory
% dirMat is the directory where diagnozed .mat files will be saved
% -> set it to '' to use the default [dirModel 'mat/']
% doComp is the switch from computation to plot
doComp=1;
if doComp==1;
doSave=1;
doPlot=0;
else;
doSave=0;
doPlot=1;
end;
doPrint=0;
for doDifObsOrMod=1:3;
if doDifObsOrMod==1; suf='modMobs'; elseif doDifObsOrMod==2; suf='obs'; else; suf='mod'; end;
%%%%%%%%%%%
%load grid:
%%%%%%%%%%%
gcmfaces_global;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
%%%%%%%%%%%%%%%
%define pathes:
%%%%%%%%%%%%%%%
dirSigma='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_err/';
%nameSigma={'sigma_MDT_glob_eccollc.bin','sigma_SLA_smooth_eccollc.bin','sigma_SLA_PWS07r2_glob_eccollc.bin'}; maxLatObs=66
nameSigma={'sigma_MDT_glob_eccollc.bin','slaerr_largescale_r4.err','slaerr_gridscale_r4.err'}; maxLatObs=90
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); mkdir([dirMat]); end;
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
runName=pwd; tmp1=strfind(runName,filesep); runName=runName(tmp1(end)+1:end);
%%%%%%%%%%%%%%%%%
%do computations:
%%%%%%%%%%%%%%%%%
if doComp==0;
eval(['load ' dirMat 'cost/cost_altimeter_' suf '.mat myflds;']);
else;
tic;
%mdt cost function term (misfit plot)
dif_mdt=rdmds2gcmfaces([dirModel 'barfiles/mdtdiff_smooth']);
sig_mdt=read_bin([dirSigma nameSigma{1}],1,0);
sig_mdt(find(sig_mdt==0))=NaN;
%store:
myflds.dif_mdt=dif_mdt;
myflds.sig_mdt=sig_mdt;
%skip blanks:
tmp1=dir([dirModel 'barfiles/sladiff_smooth*data']);
nrec=tmp1.bytes/90/1170/4;
ttShift=17;
%skip 1992
listRecs=[1+365-ttShift:17+365*18-ttShift];
nRecs=length(listRecs); TT=1993+[0:nRecs-1]/365.25;
toc; tic;
%pointwise/1day terms:
sladiff_point=repmat(0*mygrid.RAC,[1 1 nRecs]); count_point=sladiff_point;
for ii=1:3;
if ii==1; myset='tp'; elseif ii==2; myset='gfo'; else; myset='ers'; end;
%topex pointwise misfits:
if doDifObsOrMod==1;
sladiff_tmp=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],ttShift+listRecs);
elseif doDifObsOrMod==2;
sladiff_tmp=cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],ttShift+listRecs);
else;
sladiff_tmp=cost_altimeter_read([dirModel 'barfiles/sladiff_' myset '_raw'],ttShift+listRecs);
sladiff_tmp=sladiff_tmp+cost_altimeter_read([dirModel 'barfiles/slaobs_' myset '_raw'],ttShift+listRecs);
end;
%add to overall data set:
sladiff_point=sladiff_point+sladiff_tmp; count_point=count_point+(sladiff_tmp~=0);
%finalize data set:
sladiff_tmp(sladiff_tmp==0)=NaN;
%compute rms:
count_tmp=sum(~isnan(sladiff_tmp),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
eval(['myflds.rms_' myset '=sqrt(nansum(sladiff_tmp.^2,3)./count_tmp);']);
eval(['myflds.std_' myset '=nanstd(sladiff_tmp,0,3).*msk_tmp;']);
end;
%finalize overall data set:
count_point(count_point==0)=NaN;
sladiff_point=sladiff_point./count_point;
%compute overall rms,std:
count_tmp=sum(~isnan(sladiff_point),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
myflds.rms_sladiff_point=sqrt(nansum(sladiff_point.^2,3)./count_tmp);
myflds.std_sladiff_point=nanstd(sladiff_point,0,3).*msk_tmp;
%
msk_point=1*(count_point>0);
%fill blanks:
warning('off','MATLAB:divideByZero');
msk=mygrid.mskC(:,:,1); msk(find(abs(mygrid.YC)>maxLatObs))=NaN;
myflds.xtrp_rms_sladiff_point=diffsmooth2D_extrap_inv(myflds.rms_sladiff_point,msk);
myflds.xtrp_std_sladiff_point=diffsmooth2D_extrap_inv(myflds.std_sladiff_point,msk);
warning('on','MATLAB:divideByZero');
%get weight:
sig_sladiff_point=read_bin([dirSigma nameSigma{3}],1,0);
sig_sladiff_point(find(sig_sladiff_point==0))=NaN;
myflds.sig_sladiff_point=sig_sladiff_point;
%computational mask : only points that were actually observed, and in 35d average
msk_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
msk_point35d=1*(msk_point35d~=0);
tmp1=sum(msk_point35d==0&msk_point~=0);
tmp2=sum(msk_point35d~=0&msk_point~=0);
fprintf('after masking : %d omitted, %d retained \n',tmp1,tmp2);
msk_point=msk_point.*msk_point35d;
msk_point(msk_point==0)=NaN;
%plotting mask : regions of less than 100 observations are omitted
msk_100pts=1*(nansum(msk_point,3)>=100);
msk_100pts(msk_100pts==0)=NaN;
%store
myflds.msk_100pts=msk_100pts;
toc; tic;
%lsc cost function term:
if doDifObsOrMod==1;
sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
elseif doDifObsOrMod==2;
sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
else;
sladiff_smooth=cost_altimeter_read([dirModel 'barfiles/sladiff_smooth'],listRecs);
sladiff_smooth=sladiff_smooth+cost_altimeter_read([dirModel 'barfiles/slaobs_smooth'],listRecs);
end;
%mask missing points:
sladiff_smooth(sladiff_smooth==0)=NaN;
sladiff_smooth=sladiff_smooth.*msk_point;
count_tmp=sum(~isnan(sladiff_smooth),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
%compute rms:
rms_sladiff_smooth=msk_tmp.*sqrt(nanmean(sladiff_smooth.^2,3));
std_sladiff_smooth=msk_tmp.*nanstd(sladiff_smooth,0,3);
%get weight:
sig_sladiff_smooth=read_bin([dirSigma nameSigma{2}],1,0);
sig_sladiff_smooth(find(sig_sladiff_smooth==0))=NaN;
%store:
myflds.rms_sladiff_smooth=rms_sladiff_smooth;
myflds.std_sladiff_smooth=std_sladiff_smooth;
myflds.sig_sladiff_smooth=sig_sladiff_smooth;
toc; tic;
%pointwise/point35days cost function term:
if doDifObsOrMod==1;
sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
elseif doDifObsOrMod==2;
sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
else;
sladiff_point35d=cost_altimeter_read([dirModel 'barfiles/sladiff_raw'],listRecs);
sladiff_point35d=sladiff_point35d+cost_altimeter_read([dirModel 'barfiles/slaobs_raw'],listRecs);
end;
%mask missing points:
sladiff_point35d(sladiff_point35d==0)=NaN;
sladiff_point35d=sladiff_point35d.*msk_point;
count_tmp=sum(~isnan(sladiff_point35d),3); count_tmp(find(count_tmp<10))=NaN; msk_tmp=1+0*count_tmp;
%compute rms:
rms_sladiff_point35d=msk_tmp.*sqrt(nanmean(sladiff_point35d.^2,3));
std_sladiff_point35d=msk_tmp.*nanstd(sladiff_point35d,0,3);
%store:
myflds.rms_sladiff_point35d=rms_sladiff_point35d;
myflds.std_sladiff_point35d=std_sladiff_point35d;
if 0;
%difference between scales
rms_sladiff_35dMsmooth=msk_tmp.*sqrt(nanmean((sladiff_point35d-sladiff_smooth).^2,3));
std_sladiff_35dMsmooth=msk_tmp.*nanstd((sladiff_point35d-sladiff_smooth),0,3);
%store:
myflds.rms_sladiff_35dMsmooth=rms_sladiff_35dMsmooth;
myflds.std_sladiff_35dMsmooth=std_sladiff_35dMsmooth;
%difference between scales
rms_sladiff_pointMpoint35d=msk_tmp.*sqrt(nanmean((sladiff_point-sladiff_point35d).^2,3));
std_sladiff_pointMpoint35d=msk_tmp.*nanstd((sladiff_point-sladiff_point35d),0,3);
%store:
myflds.rms_sladiff_pointMpoint35d=rms_sladiff_pointMpoint35d;
myflds.std_sladiff_pointMpoint35d=std_sladiff_pointMpoint35d;
end;
%compute zonal mean/median:
for ii=1:4;
switch ii;
case 1; tmp1='mdt'; cost_fld=(mygrid.mskC(:,:,1).*myflds.dif_mdt./myflds.sig_mdt).^2;
case 2; tmp1='lsc'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_smooth./myflds.sig_sladiff_smooth).^2;
case 3; tmp1='point35d'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point35d./myflds.sig_sladiff_point).^2;
case 4; tmp1='point'; cost_fld=(mygrid.mskC(:,:,1).*myflds.rms_sladiff_point./myflds.sig_sladiff_point).^2;
end;
cost_zmean=calc_zonmean_T(cost_fld); eval(['mycosts_mean.' tmp1 '=cost_zmean;']);
cost_zmedian=calc_zonmedian_T(cost_fld); eval(['mycosts_median.' tmp1 '=cost_zmedian;']);
end;
toc; %write to disk:
if doSave; eval(['save ' dirMat 'cost/cost_altimeter_' suf '.mat myflds mycosts_mean mycosts_median;']); end;
end;%if doComp
if doPlot;
cc=[-0.4:0.05:-0.25 -0.2:0.03:-0.05 -0.03:0.01:0.03 0.05:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(myflds.dif_mdt,0,{'myCaxis',cc}); drawnow;
cc=[0:0.005:0.02 0.03:0.01:0.05 0.06:0.02:0.1 0.14:0.03:0.2 0.25:0.05:0.4];
figure; m_map_gcmfaces(myflds.rms_sladiff_smooth,0,{'myCaxis',cc}); drawnow;
figure; m_map_gcmfaces(myflds.rms_sladiff_point35d,0,{'myCaxis',cc}); drawnow;
figure; m_map_gcmfaces(myflds.rms_sladiff_point,0,{'myCaxis',cc}); drawnow;
end;
if doPlot&doPrint;
dirFig='../figs/altimeter/'; ff0=gcf-4;
for ff=1:4;
figure(ff+ff0); saveas(gcf,[dirFig runName '_' suf num2str(ff)],'fig');
eval(['print -depsc ' dirFig runName '_' suf num2str(ff) '.eps;']);
eval(['print -djpeg90 ' dirFig runName '_' suf num2str(ff) '.jpg;']);
end;
end;
end;%for doDifObsOrMod=1:3;
function [fldOut]=cost_altimeter_read(fileIn,recIn);
nrec=length(recIn);
global mygrid; siz=[size(convert2gcmfaces(mygrid.XC)) nrec];
lrec=siz(1)*siz(2)*4;
myprec='float32';
fldOut=zeros(siz);
fid=fopen([fileIn '.data'],'r','b');
for irec=1:nrec;
status=fseek(fid,(recIn(irec)-1)*lrec,'bof');
fldOut(:,:,irec)=fread(fid,siz(1:2),myprec);
end;
fldOut=convert2gcmfaces(fldOut);
gcmfaces/ecco_v4/cost_sst.m 0000644 0004354 0000144 00000015213 12365532071 015554 0 ustar gforget users function []=cost_sst(dirModel,dirMat,doComp,dirTex,nameTex);
%object: compute cost function term for sst data
%inputs: dimodel is the model directory
% dirMat is the directory where diagnozed .mat files will be saved
% -> set it to '' to use the default [dirModel 'mat/']
% doComp is a switch (1->compute; 0->display)
%optional: dirTex is the directory where tex and figures files are created
% (if not specified then display all results to screen instead)
% nameTex is the tex file name (default : 'myPlots')
if isempty(dirMat); dirMat=[dirModel 'mat/']; else; dirMat=[dirMat '/']; end;
if isempty(dir(dirMat)); mkdir([dirMat]); end;
%determine if and where to create tex and figures files
dirMat=[dirMat '/'];
if isempty(who('dirTex'));
addToTex=0;
else;
if ~ischar(dirTex); error('mis-specified dirTex'); end;
addToTex=1;
if isempty(who('nameTex')); nameTex='myPlots'; end;
fileTex=[dirTex nameTex '.tex'];
end;
if doComp;
%grid, params and inputs
gcmfaces_global; global myparms;
if ~isfield(mygrid,'XC'); grid_load('./GRID/',5,'compact'); end;
if ~isfield(mygrid,'LATS_MASKS'); gcmfaces_lines_zonal; end;
if isfield(myparms,'yearFirst'); yearFirst=myparms.yearFirst; yearLast=myparms.yearLast;
else; yearFirst=1992; yearLast=2011;
end;
dirData='/net/nares/raid11/ecco-shared/ecco-version-4/input/';
subdirReynolds='input_sst_reynolds_v4/'; nameReynolds='reynolds_oiv2_r1';
subdirRemss='input_sst_remss_v4/'; nameRemss='tmi_amsre_oisst_r1';
subdirL2p='input_sst_amsre_daily/'; nameL2p='tmi_amsre_l2p_r1';
if 0;%old versions
nameReynolds='v4_reynolds_monthly';
nameRemss='v4_TMI_AMSRE_sst_monav';
nameL2p='amsre_r2';
yearLast=2010;
end;
dirErr='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_all_from_pleiades/';
fld_err=read2memory([dirErr '/sigma_half.bin'],[90 1170]);
fld_err=convert2gcmfaces(fld_err);
fld_w=fld_err.^-2;
fileModel=dir([dirModel 'barfiles/tbar*data']); fileModel=['barfiles/' fileModel.name];
%computational loop
mod_m_rey=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
mod_m_remss=repmat(NaN*mygrid.mskC(:,:,1),[1 1 12*(yearLast-yearFirst+1)]);
%
zm_mod=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_mod_m_rey=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
zm_mod_m_remss=repmat(ones(length(mygrid.LATS),1),[1 12*(yearLast-yearFirst+1)]);
%
for ycur=yearFirst:yearLast;
fprintf(['starting ' num2str(ycur) '\n']);
tic;
for mcur=1:12;
%load Reynolds SST
file0=[dirData subdirReynolds nameReynolds '_' num2str(ycur)];
if ~isempty(dir(file0)); fld_rey=read_bin(file0,mcur,0); else; fld_rey=NaN*mygrid.mskC(:,:,1); end;
fld_rey(find(fld_rey==0))=NaN;
fld_rey(find(fld_rey<-99))=NaN;
%load Remss SST
file0=[dirData subdirRemss nameRemss '_' num2str(ycur)];
if ~isempty(dir(file0)); fld_remss=read_bin(file0,mcur,0); else; fld_remss=NaN*mygrid.mskC(:,:,1); end;
fld_remss(find(fld_remss==0))=NaN;
fld_remss(find(fld_remss<-99))=NaN;
%load model SST
mm=(ycur-yearFirst)*12+mcur;
fld_mod=read_bin([dirModel fileModel],mm,1).*mygrid.mskC(:,:,1);
%store misfit maps
mod_m_rey(:,:,mm)=fld_mod-fld_rey;
mod_m_remss(:,:,mm)=fld_mod-fld_remss;
%comp zonal means
zm_mod(:,mm)=calc_zonmean_T(fld_mod.*mygrid.mskC(:,:,1));
msk=mygrid.mskC(:,:,1); msk(fld_rey==0)=NaN;
zm_rey(:,mm)=calc_zonmean_T(fld_rey.*msk);
zm_mod_m_rey(:,mm)=calc_zonmean_T(fld_mod-fld_rey.*msk);
msk=mygrid.mskC(:,:,1); msk(fld_remss==0)=NaN;
zm_remss(:,mm)=calc_zonmean_T(fld_remss.*msk);
zm_mod_m_remss(:,mm)=calc_zonmean_T(fld_mod-fld_remss.*msk);
end;
toc;
end;
%compute rms maps
rms_to_rey=sqrt(nanmean(mod_m_rey.^2,3));
rms_to_remss=sqrt(nanmean(mod_m_remss.^2,3));
if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
eval(['save ' dirMat '/cost/cost_sst.mat fld_err rms_* zm_*;']);
else;%display previously computed results
global mygrid;
if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
eval(['load ' dirMat '/cost_sst.mat;']);
if ~isempty(who('fld_rms'));
figure; m_map_gcmfaces(fld_rms,0,{'myCaxis',[0:0.2:1.2 1.5:0.5:3 4:1:6 8 10]/2});
myCaption={'modeled-observed rms -- sea surface temperature (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
else;
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});
myCaption={'modeled-Reynolds rms -- sea surface temperature (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
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});
myCaption={'modeled-REMSS rms -- sea surface temperature (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
ny=size(zm_rey,2)/12;
[xx,yy]=meshgrid(1992+([1:ny*12]-0.5)/12,mygrid.LATS);
figureL;
obs=zm_rey; obsCycle=sst_cycle(obs);
mod=zm_rey+zm_mod_m_rey; modCycle=sst_cycle(mod);
mis=zm_mod_m_rey; misCycle=sst_cycle(mis);
subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
title('Reynolds sst anomaly');
subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
title('ECCO sst anomaly');
subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-Reynolds sst misfit');
myCaption={'ECCO and Reynolds zonal mean sst anomalies (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
figureL;
obs=zm_remss; obsCycle=sst_cycle(obs);
mod=zm_remss+zm_mod_m_remss; modCycle=sst_cycle(mod);
mis=zm_mod_m_remss; misCycle=sst_cycle(mis);
subplot(3,1,1); pcolor(xx,yy,obs-obsCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
title('REMSS sst anomaly');
subplot(3,1,2); pcolor(xx,yy,mod-modCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); set(gca,'XTick',[]); ylabel('latitude');
title('ECCO sst anomaly');
subplot(3,1,3); pcolor(xx,yy,mis-misCycle); shading flat; caxis([-1 1]*1); colorbar;
set(gca,'FontSize',14); ylabel('latitude'); title('ECCO-REMSS sst misfit');
myCaption={'ECCO and REMSS zonal mean sst anomalies (K)'};
if addToTex; write2tex(fileTex,2,myCaption,gcf); end;
end;
end;
function [zmCycle]=sst_cycle(zmIn);
ny=size(zmIn,2)/12;
zmCycle=NaN*zeros(179,12);
for mm=1:12;
zmCycle(:,mm)=nanmean(zmIn(:,mm:12:ny*12),2);
end;
zmCycle=repmat(zmCycle,[1 ny]);
gcmfaces/ecco_v4/map_grace.m 0000644 0004354 0000144 00000003345 12161145421 015625 0 ustar gforget users
xC1=[0.5:359.5]'*ones(1,180);
yC1=ones(360,1)*[-89.5:89.5];
v4_reclen=convert2gcmfaces(mygrid.RAC); v4_reclen=size(v4_reclen); v4_reclen=v4_reclen(1)*v4_reclen(2);
dir_out='./'; %output directory
dir_in='./'; %input directory
file_pref='GRACE_CSR_withland_';
%load error field
[fieldErr_nosmooth,fieldErr]=read_grace_fld([dir_in 'GRACE_CSR_Error.asc']);
%process grace fields
for yy=1992:2012;
for mm=1:12;
tt=sprintf('%4i%02i',yy,mm);
file0=dir([dir_in file_pref tt '.asc']);
if ~isempty(file0);
%read
file0=file0.name;
field0=read_grace_fld([dir_in file0]);
%map
x=[xC1-360;xC1]; y=[yC1;yC1]; z=[field0;field0];
x=[x x x]; y=[y-180 y y+180]; z=[flipdim(z,2) z flipdim(z,2)];
field1=0*mygrid.XC;
for ii=1:5;
xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
zi = interp2(x',y',z',xi,yi);
%zi = griddata(x,y,z,xi,yi);
field1{ii}=zi;
end;
%mask the model poles
tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;
else;
fprintf(['did not find : ' file_pref tt '.asc \n']);
field1=NaN*mygrid.RAC;
end;
%use -999 for mask
field1(isnan(field1))=-999;
%write to file
if mm==1; fid=fopen([dir_out file_pref num2str(yy)],'w','b'); end;
fwrite(fid,convert2gcmfaces(field1),'float32');
if mm==12; fclose(fid); end;
end;
end;
%process error estimate
z=[fieldErr;fieldErr]; z=[flipdim(z,2) z flipdim(z,2)];
field1=0*mygrid.XC;
for ii=1:5;
xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
zi = interp2(x',y',z',xi,yi);
%zi = griddata(x,y,z,xi,yi);
field1{ii}=zi;
end;
%mask the model poles
tmp1=find(mygrid.RAC<8e8&mygrid.YC>0); field1(tmp1)=NaN;
tmp1=find(mygrid.RAC<2e8&mygrid.YC<0); field1(tmp1)=NaN;
%use 0 for mask
field1(isnan(field1))=0;
%write to file
write2file([dir_out file_pref 'err'],convert2gcmfaces(field1));
gcmfaces/ecco_v4/rads_noice_mad.m 0000644 0004354 0000144 00000016014 12145162164 016640 0 ustar gforget users function []=rads_noice_mad(choiceData,l0,L0,doTesting,doNoice);
%object : outlier (& ice if doNoice) data flagging for altimetry
%inputs : choiceData is 'topexfile','ersfile' or 'gfofile'
% l0 is the box longitude (l0+-10 will be treated)
% L0 is the box latitude (L0+-10 will be treated)
% doTesting (0 by default) is a testing option
% doNoice (0 by default) activate the flagging
% of all ice covered points (based on nsidc)
%
if isempty(who('doTesting')); doTesting=0; end;
if isempty(who('doNoice')); doNoice=0; end;
%directories
dirData='/bay/scratch2/RADS_alongtrk_v4_ice_flagged_2013/';
topexfile = 'tj_daily_ssh_v4_r1';
ersfile = 'en_daily_ssh_v4_r1';
gfofile = 'g1_daily_ssh_v4_r1';
dirIce='/net/nares/raid11/ecco-shared/ecco-version-4/input/input_nsidc_all/';
fileIce='nsidc79_daily';
dirOutput='./';
%testing switches
if doTesting; l0=-170; L0=-80; choiceData='gfofile'; end;
%select region
gcmfaces_global;
XC=convert2gcmfaces(mygrid.XC.*mygrid.mskC(:,:,1)); XC=XC(:);
YC=convert2gcmfaces(mygrid.YC.*mygrid.mskC(:,:,1)); YC=YC(:);
ii=find(XC>=l0-10&XC<=l0+10&YC>=L0-10&YC<=L0+10);
ni=length(ii);
nt=7671;%corresponds to 1992-2012
%load data
eval(['fileData=' choiceData ';']);
allData=zeros(ni,nt);
allIce=zeros(ni,nt);
ndata=0;
for yy=1992:2012;
fprintf(['reading ' fileData '_' num2str(yy) '\n']);
tmp1=dir([dirData fileData '_' num2str(yy)]); tmp0=tmp1.bytes/90/1170/4;
tmp1=read2memory([dirData fileData '_' num2str(yy)],[90*1170 tmp0]);
allData(:,ndata+[1:tmp0])=tmp1(ii,:);
if doNoice;
tmp1=read2memory([dirIce fileIce '_' num2str(yy)],[90*1170 tmp0]);
allIce(:,ndata+[1:tmp0])=tmp1(ii,:);
end;
ndata=ndata+tmp0;
end;
allData(allData<-999)=NaN; allData=allData/100;
allIce(allIce<-999)=NaN; allIce=allIce;
if ndata~=nt; error('nt was mispecified'); end;
%if doTesting; eval(['save ' dirOutput 'rads_noice_mad.demo.mat allData allIce ndata;']); end;
%if doTesting; eval(['load ' dirOutput 'rads_noice_mad.demo.mat allData allIce ndata;']); end;
%only keep ice free data points
noiceData=allData; noiceData(allIce>0)=NaN;
iceRejectPerc=100*(1-nansum(~isnan(noiceData(:)))/nansum(~isnan(allData(:))));
%the outliers filtering
medianData=nanmedian(noiceData,2)*ones(1,nt);%center of distribution
madData=1.4826*mad(noiceData,1,2);%width of the distribution
cutoff=4*madData*ones(1,nt);%a factor 4 corresponds to a typical cost of 16
goodData=noiceData; goodData(abs(noiceData-medianData)>cutoff)=NaN;
madRejectPerc=100*(1-nansum(~isnan(goodData(:)))/nansum(~isnan(noiceData(:))));
%output result
[iceRejectPerc madRejectPerc]
%note : original testing was done with gfofile, l0=-170; L0=-80;
% including ice points madRejectPerc is ~ 3.9%
% excluding ice points madRejectPerc drops to 0.33%
if ~doTesting;
eval(['save ' dirOutput fileData 'l0is' num2str(l0) 'L0is' num2str(L0) ...
' ii goodData madData medianData iceRejectPerc madRejectPerc choiceData l0 L0 doNoice']);
end;
%------- testing case -----%
if doTesting;
%note : original testing was done with gfofile, l0=-170; L0=-80;
% the last point (see next line) happened to work as an example
jj=find(~isnan(madData));
%do example of mad detection of ice covered points as outliers at one grid point
kk=jj(end); xc_kk=XC(ii(kk)); yc_kk=YC(ii(kk));
all_kk=allData(kk,:);
median_kk=nanmedian(all_kk);%center of distribution
mad_kk=1.4826*mad(all_kk,1,2);%width of distribution
cutoff_kk=4*mad_kk*ones(1,nt);%a factor 4 corresponds to a typical cost of 16
good_kk=all_kk; good_kk(abs(good_kk-median_kk)>cutoff_kk)=NaN;
reject_kk=100*(1-nansum(~isnan(good_kk))/nansum(~isnan(all_kk)));
%compute the various statistics
msk=1+0*allIce; msk(isnan(allData))=NaN;
%in the testing case we dont filter out the ice points : msk(allIce>0)=NaN;
stdOut=nanstd(msk.*allData,0,2);%sample standard deviation
iqrOut=0.7413*iqr(msk.*allData,2);%intequartile range estimate of std
madOut=1.4826*mad(msk.*allData,1,2);%median absolute difference estimate of std
jj=find(~isnan(stdOut));
%plot region
if L0>60; myproj=2; elseif L0<-60; myproj=3; else; myproj=1; end;
figureL; m_map_gcmfaces(mygrid.Depth,myproj,{'myCaxis',[0 6e3]});
colormap('gray'); colorbar off;
m_map_gcmfaces({'plot',XC(ii),YC(ii),'m.'},myproj,{'doHold',1});
m_map_gcmfaces({'plot',xc_kk,yc_kk,'k.','MarkerSize',24},myproj,{'doHold',1});
saveas(gcf,[dirOutput 'outliers_0'],'fig');
eval(['print -djpeg90 ' dirOutput 'outliers_0.jpg']);
eval(['print -depsc ' dirOutput 'outliers_0.eps']);
%plot sensitivity of std estimates to outliers
figureL; set(gca,'FontSize',16);
plot(stdOut(jj),'.-'); hold on; plot(madOut(jj),'r.-'); plot(iqrOut(jj),'g.-');
legend('SSTD','MAD*1.4826','IQR*0.7413');
xlabel('grid point index'); ylabel('estimated standard deviation');
saveas(gcf,[dirOutput 'outliers_1'],'fig');
eval(['print -djpeg90 ' dirOutput 'outliers_1.jpg']);
eval(['print -depsc ' dirOutput 'outliers_1.eps']);
%plot distributions for the overall region
figureL;
subplot(3,1,1); set(gca,'FontSize',16);
%all data points
tmp1=allData(jj,:); hist(tmp1(:),[-1.5:0.02:1.5]);
aa=axis; aa(1:2)=[-1 1]; axis(aa); title('all data points');
%and the corresponding normal distribution
[xx,yy]=normal_distribution([-1.5:0.02:1.5],tmp1(:));
hold on; plot(xx,yy,'m','LineWidth',0.5);
%
subplot(3,1,2); set(gca,'FontSize',16);
tmp1=noiceData(jj,:); hist(tmp1(:),[-1.5:0.02:1.5]);
aa=axis; aa(1:2)=[-1 1]; axis(aa); title('ice free data points');
[xx,yy]=normal_distribution([-1.5:0.02:1.5],tmp1(:));
hold on; plot(xx,yy,'m','LineWidth',0.5);
subplot(3,1,3); set(gca,'FontSize',16);
tmp1=allData(jj,:); tmp1(~isnan(noiceData(jj,:)))=NaN; hist(tmp1(:),[-1.5:0.02:1.5]);
aa=axis; aa(1:2)=[-1 1]; axis(aa); title('icy data points');
[xx,yy]=normal_distribution([-1.5:0.02:1.5],tmp1(:));
hold on; plot(xx,yy,'m','LineWidth',0.5);
saveas(gcf,[dirOutput 'outliers_2'],'fig');
eval(['print -djpeg90 ' dirOutput 'outliers_2.jpg']);
eval(['print -depsc ' dirOutput 'outliers_2.eps']);
figureL;
subplot(2,1,1); set(gca,'FontSize',16);
hist(all_kk,[-1.5:0.02:1.5]);
aa=axis; aa(1:2)=[-1 1]; axis(aa); title('before MAD detection');
%and the corresponding normal distribution
[xx,yy]=normal_distribution([-1.5:0.02:1.5],all_kk(:));
hold on; plot(xx,yy,'m','LineWidth',2);
subplot(2,1,2); set(gca,'FontSize',16);
hist(good_kk,[-1.5:0.02:1.5]);
aa=axis; aa(1:2)=[-1 1]; axis(aa); title('after MAD detection');
%and the corresponding normal distribution
[xx,yy]=normal_distribution([-1.5:0.02:1.5],all_kk(:));
hold on; plot(xx,yy,'m','LineWidth',2);
saveas(gcf,[dirOutput 'outliers_3'],'fig');
eval(['print -djpeg90 ' dirOutput 'outliers_3.jpg']);
eval(['print -depsc ' dirOutput 'outliers_3.eps']);
end;
function [xx,yy]=normal_distribution(xx,data);
dx=diff(xx);
dx=median(dx);
%dx=unique(dx); if length(dx)>1; error('plot_normal_distribution expect regular spacing\n'); end;
mm=nanmedian(data);
ss=1.4826*mad(data,1);
nn=sum(~isnan(data));
%for testing : mm=0.2; ss=0.1; nn=100;
xx=xx+mm;%reset xx to center of distribution
yy=nn*dx/ss/sqrt(2*pi)*exp(-0.5*xx.*xx/ss/ss);%compute distribution
%for testing : sum(yy)/nn
gcmfaces/ecco_v4/extend_xx.m 0000644 0004354 0000144 00000004467 12365534162 015735 0 ustar gforget users function []=extend_xx(dirIn,dirOut,ndayIn,nyearOut);
%object: extend atmospheric controls to additional years
%inputs: dirIn is the input directory
% dirOut is the output directory
% ndayIn is the control vector frecuency, in days
% nyearOut is the extended number of years
%
%note: we will take the time mean seasonal cycle, and transition
% to it linearly over the last year that was already covered
%example:
%dirIn='/net/weddell/raid3/gforget/ecco_v4/output_files2/apr1alpha_it9/';
%dirOut='/net/weddell/raid3/gforget/ecco_v4/input_files2/apr1alpha_it9_extended_xx/';
%extend_xx(dirIn,dirOut,14,20);
gcmfaces_global;
%first copy time independent controls
listCp={'diffkr','kapgm','kapredi','salt','theta'};
for ii=1:length(listCp);
fileIn=dir([dirIn 'ADXXfiles/xx_' listCp{ii} '.00*data']);
copyfile([dirIn 'ADXXfiles/' fileIn.name] , [dirOut fileIn.name]);
end;
%then extend the time dependent ones
for ii=1:7;
switch ii;
case 1; xxName='atemp';
case 2; xxName='aqh';
case 3; xxName='tauu';
case 4; xxName='tauv';
case 5; xxName='lwdown';
case 6; xxName='swdown';
case 7; xxName='precip';
end;
%read model cost output
fld_xx=rdmds2gcmfaces([dirIn 'ADXXfiles/xx_' xxName '.00*']);
%determine already covered time period
nrec=size(fld_xx{1},3);
nyearIn=nrec*ndayIn/365;
nrecInOneYear=round(365/ndayIn);
nrecOut=nrecInOneYear*nyearOut+1;
%determine seasonal cycle
season_xx=0*fld_xx(:,:,1:nrecInOneYear);
for tt=1:nrecInOneYear;
season_xx(:,:,tt)=mean(fld_xx(:,:,tt:nrecInOneYear:nrec),3);
end;
%determine transition factor (fld_xx -> season_xx)
nrec0=nrecInOneYear*(floor(nyearIn)-1);
fac=([1:nrecOut]-nrec0)/(nrec-nrec0);
fac=max(min(fac,1),0);
%build extended time series
more_xx=zeros(fld_xx(:,:,1),nrecOut);
more_xx(:,:,1:nrec0)=fld_xx(:,:,1:nrec0);
for tt=nrec0+1:nrecOut;
ttt=mod(tt,nrecInOneYear);
if ttt==0; ttt=nrecInOneYear; end;
if tt