Warning, /verification/offline_exf_seaice/input/grph_diag.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
94e7a5d378 Jean*0001 if size(who('kpr'),1) > 0,
0002 fprintf('kpr is defined and = %i \n',kpr);
0003 else
0004 fprintf('kpr undefined ; set to 1 \n'); kpr=1 ;
0005 end
0006
0007 Nexp=1; %- plot results from 1 (in rDir1) or from 2 experiments (in rDir1 & rDir2)
0008 kdif=0; %- if Nexp=2, do diff plot (kdif=1) or plot both exp. field (kdif=0)
0009 deltaT=3600;
0010 rDir1='./'; iter=[1:5]*24;
0011 rDir2=rDir1; ite2=iter;
0012 %rDir2='res_n05/'; %ite2=24;
0013 gDir=rDir1;
0014
0015 namF='iceDiag';
0016 %namF='snapshot'; iter=iter-1; ite2=iter;
0017
0018 ii=strfind(rDir1,'/'); if length(ii) > 1, ii=1+ii(end-1); else ii=1; end
0019 titexp1=rDir1(ii:end-1); titexp1=strrep(titexp1,'_','\_');
0020 ii=strfind(rDir2,'/'); if length(ii) > 1, ii=1+ii(end-1); else ii=1; end
0021 titexp2=rDir2(ii:end-1); titexp2=strrep(titexp2,'_','\_');
0022
0023 G=load_grid(gDir,0);
0024 nx=G.dims(1); ny=G.dims(2);
0025 xc=[1:nx]; xc=xc-mean(xc); yc=[1:ny]-.5;
0026
0027 msk1=squeeze(G.hFacC); msk1=ceil(msk1); msk1=min(msk1,1);
0028
0029 Nit=length(iter);
0030 if kpr > 0,
3fc7124653 Jean*0031 clear missingValue ;
94e7a5d378 Jean*0032 [v4d1,its,M]=rdmds([rDir1,namF],iter); %v4d1=squeeze(v4d1);
0033 eval(M); namV=char(fldList) ; nV=size(namV,1);
0034 if Nexp > 1, [v4d2,its,M]=rdmds([rDir2,namF],ite2); end
0035 jA=0; [J]=find(strcmp(fldList,'SI_Fract')); if length(J) == 1 & jA == 0, jA=J; end
0036 [J]=find(strcmp(fldList,'SIarea ')); if length(J) == 1 & jA == 0, jA=J; end
0037 jH=0; [J]=find(strcmp(fldList,'SI_Thick')); if length(J) == 1 & jH == 0, jH=J; end
0038 jE=0; [J]=find(strcmp(fldList,'SIheff ')); if length(J) == 1 & jE == 0, jE=J; end
3fc7124653 Jean*0039 if size(who('missingValue'),1) > 0,
0040 fprintf('take missingValue from meta file:');
0041 if strcmp(dataprec,'float32'), misVal=single(missingValue); else misVal=missingValue; end
0042 else
0043 fprintf('no missingValue defined ; set'); misVal=-999.;
0044 end
0045 fprintf(' misVal= %f\n',misVal);
94e7a5d378 Jean*0046 end
0047
0048 % add/replace Effective thickness
0049 nVo=nV;
0050 %- to disable addition/replacement:
0051 %jA=-1;
0052 if jE == 0 & jA*jH > 0,
0053 jE=12;
0054 if jE > nV, jE=nV+1;
0055 var=v4d1; v4d1=zeros([dimList jE Nit]); v4d1(:,:,1:nV,:)=var;
0056 if Nexp > 1, var=v4d2; v4d2=zeros([dimList jE Nit]); v4d2(:,:,1:nV,:)=var; end
0057 nV=jE;
0058 end
0059 namV(jE,:)='SI_Heff ';
0060 v4d1(:,:,jE,:)=v4d1(:,:,jA,:).*v4d1(:,:,jH,:);
0061 if Nexp > 1, v4d2(:,:,jE,:)=v4d2(:,:,jA,:).*v4d2(:,:,jH,:); end
0062 end
0063 if jH == 0 & jA*jE > 0,
0064 jH=12;
0065 if jH > nV, jH=nV+1;
0066 var=v4d1; v4d1=zeros([dimList jE Nit]); v4d1(:,:,1:nV,:)=var;
0067 if Nexp > 1, var=v4d2; v4d2=zeros([dimList jE Nit]); v4d2(:,:,1:nV,:)=var; end
0068 nV=jH;
0069 end
0070 namV(jH,:)='SI_thick';
0071 var=v4d1(:,:,jA,:); var(find(var==0))=-1; var=1./var; var=max(var,0);
0072 v4d1(:,:,jH,:)=v4d1(:,:,jE,:).*var;
0073 if Nexp > 1,
0074 var=v4d2(:,:,jA,:); var(find(var==0))=-1; var=1./var; var=max(var,0);
0075 v4d2(:,:,jH,:)=v4d2(:,:,jE,:).*var;
0076 end
0077 end
0078
0079 % nt : select which time index to plot in 2-d
0080 % is,js : select which section to plot
0081 nf=0; nt=Nit; kplot=ones(1,nV);
0082 %kplot(11:nV)=0;
0083 is=61; js=41;
0084
0085 xtxt=xc(1); ytxt=yc(1)-5*(yc(2)-yc(1)); clin='kbcmrgy';
0086 xydp=[50 20]; xyp0=[50 20]+nf*xydp; xysp=[500 700]; xydp=[100 40];
0087
0088 for j=1:nV,
0089 facM=1.; ccB=[0 0];
0090 if j <= nVo,
0091 %- convert Fresh-water flux (from kg/m^2/s & m/s) to [mm/d]:
0092 if strcmp(fldList(j),'SIfrw2oc') | strcmp(fldList(j),'SIfrwAtm'), facM=86400.; end
0093 if strcmp(fldList(j),'EXFempmr'), facM=1000*86400.; end
0094 end
0095 if kplot(j) == 1,
0096 nf=nf+1;
0097 %if kplot(nf) == 1,
0098 xyp0=xyp0+xydp;
0099 figure(nf); set(nf,'position',[xyp0 xysp]);clf;
0100 ns=210;
0101 ns=ns+1; subplot(ns);
0102 jg=j;
0103 titv=namV(jg,:); titv=strrep(titv,'_','\_');
0104 titime=sprintf('t= %4.1f d',iter(nt)*deltaT/86400);
0105 var=squeeze(v4d1(:,:,jg,nt)); var(find(var==misVal))=NaN;
0106 var(find(msk1==0))=NaN;
0107 var=facM*var; mnV=min(var(:)); MxV=max(var(:));
0108 if MxV > mnV,
0109 imagesc(xc,yc,var'); set(gca,'YDir','normal');
0110 ccB=[mnV MxV] + [-1 1]*(MxV-mnV)/10; caxis(ccB); change_colmap(-1);
0111 if strcmp(namV(jg,:),'SI_thick'), ccB=[-.1 1.1]; end
0112 if ccB(2) <= ccB(1), ccB=[mnV MxV] + [-1 1]*(MxV-mnV)/10; end
0113 caxis(ccB); change_colmap(-1);
0114 %contourf(xc,yc,var',14);
0115 BB=colorbar('EastOutside');
0116 %[cs,h]=contour(yc,yc,var',cc);clabel(cs); isoline0(h);
0117 else fprintf(' %s , uniform = %e\n',namV(jg,:),MxV); end
0118 title([titexp1,' ; ',titv,' ; ',titime]);
0119 text(xtxt,ytxt,sprintf('min,Max= %9.5g , %9.5g', mnV, MxV))
0120 if Nexp > 1, ns=ns+1; subplot(ns);
0121 vv2=squeeze(v4d2(:,:,jg,nt)); vv2(find(vv2==misVal))=NaN;
0122 var=facM*vv2-kdif*var;
0123 var(find(msk1==0))=NaN;
0124 mnV=min(var(:)); MxV=max(var(:));
0125 if MxV > mnV,
0126 imagesc(xc,yc,var'); set(gca,'YDir','normal');
0127 %- note: 2nd plot has same color-range (except if diff plot)
0128 if (ccB(2)-ccB(1))*(kdif-1) >= 0, ccB=[mnV MxV] + [-1 1]*(MxV-mnV)/10; end
0129 caxis(ccB); change_colmap(-1);
0130 %contourf(xc,yc,var',12);
0131 BB=colorbar('EastOutside');
0132 %[cs,h]=contour(yc,yc,var',cc);clabel(cs); isoline0(h);
0133 else fprintf(' %s , uniform = %e\n',namV(jg,:),MxV); end
0134 if kdif == 1, titplot=['Diff: ',titexp2,' - ',titexp1];
0135 else titplot=[titexp2,' ; ',titv]; end ; title(titplot);
0136 text(xtxt,ytxt,sprintf('min,Max= %9.5g , %9.5g', mnV, MxV))
0137 else
0138 subplot(223);
0139 var=squeeze(v4d1(is,:,jg,:)); var(find(var==misVal))=NaN;
0140 var=facM*var; mnV=min(var(:)); MxV=max(var(:));
0141 m1d=squeeze(msk1(is,:)); [J]=find(m1d==0);
0142 for n=1:Nit,
0143 jt=1+rem(n-1,length(clin));
0144 if n > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
0145 var(J,n)=NaN;
0146 plot(yc,var(:,n),ccln);
0147 if n==1, hold on; end
0148 end
0149 hold off; AA=axis; axis([0 ny AA(3:4)]);
0150 grid
0151 titplot=[titexp1,' ; ',titv,' ; is=',int2str(is)];
0152 title(titplot);
0153 subplot(224);
0154 var=squeeze(v4d1(:,js,jg,:)); var(find(var==misVal))=NaN;
0155 var=facM*var; mnV=min(var(:)); MxV=max(var(:));
0156 m1d=squeeze(msk1(:,js)); [J]=find(m1d==0);
0157 for n=1:Nit,
0158 jt=1+rem(n-1,length(clin));
0159 if n > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
0160 var(J,n)=NaN;
0161 plot(xc,var(:,n),ccln);
0162 if n==1, hold on; end
0163 end
0164 hold off;
0165 grid
0166 titplot=[titexp1,' ; ',titv,' ; js=',int2str(js)];
0167 title(titplot);
0168 end
0169 put_date;
0170 end
0171 end
0172
0173 return
