Back to home page

MITgcm

 
 

    


Warning, /verification/offline_exf_seaice/input/grph_res.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 rDir1='./'; rDir2=rDir1;
                0002 %rDir2='res_c00/';
                0003 %rDir2='./';
                0004 
                0005 %kdif=0; %- do diff plot (kdif=1) or plot both exp./var/time (kdif=0)
                0006 kdif=0; nf=0; k=0; ccB=[0 0]; ccD=[0 0];
                0007 
                0008 %- iters to load
                0009 deltaT=3600; iter=[0:5]*24;
                0010 %iter=[0 5]*24;
                0011 %iter=iter(2:end); %- iter=0 missing for some fields
                0012 Nit=length(iter);
                0013 
                0014 %- select variables to plot
                0015 %namv1='ice_fract'; ccB=[-.1 1.1];
                0016 namv1='ice_iceH'; ccB=[.1 1.1]/2.1; ccD=[-1 1]*0.06; %ccD=[-1 10]/2;
                0017 %namv1='ice_Tsrf'; ccB=[-5 .1];
                0018 %namv1='T'; ccB=[-.1 0.8]-1.62;
                0019 %namv1='AREA'; ccB=[-.1 1.1];
                0020 %namv1='HEFF'; ccB=[.1 1.1]/1.1; ccD=[-1 1]*0.06; %ccD=[-1 10]/2;
                0021 %namv1='Qsw'; ccB=[-1.2 0.1]*10;
                0022 %namv1='UICE'; nf=1; ccB=[-.1 0.7];
                0023 %namv1='VICE'; nf=2; ccB=[-1 1]*.21;
                0024 namv2=namv1; cc2=ccB;
                0025 %namv1='U'; namv2='V'; nf=1; ccB=[-1 6]*.1; cc2=[-1 1]*.2;
                0026 %namv2='ice_iceH';
                0027 %--------------------------
                0028 %namv='UICE'; nf=3;
                0029 %namv='U'; k=1;
                0030 %namv='KPPdiffKzS'; k=2;
                0031 
                0032 %- select time index to plot
                0033 nt1=Nit;
                0034 nt2=nt1;
                0035 nt1=2; nt2=Nit;
                0036 
                0037 %- load grid-files:
                0038 G=load_grid(rDir1,0);
                0039 nx=G.dims(1); ny=G.dims(2);
                0040 msk1=ceil(G.hFacC(:,:,1));
                0041 
                0042 %- load fields:
                0043 v1=rdmds([rDir1,namv1],iter);
                0044 v2=rdmds([rDir2,namv2],iter);
                0045 
                0046 %- kconv=1 : to convert in-situ thickness to effective thickness
                0047 %  kconv=-1: or the other way
                0048 kconv=0;
                0049 if kconv==1,
                0050  namv3='ice_fract';
                0051  v3=rdmds([rDir1,namv3],iter);
                0052  v4=rdmds([rDir2,namv3],iter);
                0053  if strcmp(namv1,'ice_iceH'), v1=v1.*v3; namv1='Ice-Vol'; end
                0054  if strcmp(namv2,'ice_iceH'), v2=v2.*v4; namv2='Ice-Vol'; end
                0055 end
                0056 if kconv==-1,
                0057  namv3='AREA';
                0058  v3=rdmds([rDir1,namv3],iter);
                0059  v4=rdmds([rDir2,namv3],iter);
                0060  v3(find(v3==0))=-1; v3=1./v3; v3=max(v3,0);
                0061  v4(find(v4==0))=-1; v4=1./v4; v4=max(v4,0);
                0062  if strcmp(namv1,'HEFF'), v1=v1.*v3; namv1='Ice-H'; end
                0063  if strcmp(namv2,'HEFF'), v2=v2.*v4; namv2='Ice-H'; end
                0064 end
                0065 
                0066 var=abs(v2-v1); %var(1:5,:)=0;
                0067 nDim=length(size(v1)); if Nit > 1, nDim=nDim-1; end
                0068 if nDim == 2, ccM=max(max(var)); k=0;
                0069 else ccM=max(max(max(var))); end
                0070 if max(ccM) > 0, fprintf(' %e',ccM) ; fprintf('\n'); end
                0071 %return
                0072 
                0073 titex1=rDir1(1:end-1); titex1=strrep(titex1,'_','\_');
                0074 titex2=rDir2(1:end-1); titex2=strrep(titex2,'_','\_');
                0075 titva1=namv1; titva1=strrep(titva1,'_','\_');
                0076 titva2=namv2; titva2=strrep(titva2,'_','\_');
                0077 titim1=sprintf('t= %4.1f d',iter(nt1)*deltaT/86400);
                0078 titim2=sprintf('t= %4.1f d',iter(nt2)*deltaT/86400);
                0079 
                0080 xax=[1:nx]-.5 ; yax=[1:ny]-.5;
                0081 xtxt=xax(1); ytxt=yax(1)-5*(yax(2)-yax(1));
                0082 xydp=[50 20]; xyp0=[50 20]+nf*xydp; xysp=[500 700]; xydp=[100 40];
                0083 
                0084 nf=nf+1; xyp0=xyp0+xydp;
                0085 figure(nf); set(nf,'position',[xyp0 xysp]);clf;
                0086 subplot(211);
                0087 %for nt=1:Nit,
                0088  if k > 0, var=v1(:,:,k,nt1); else var=v1(:,:,nt1); end
                0089  if strcmp(namv1,'U') & strcmp(namv2,'V'),
                0090   var=var.*var; vtmp=var;
                0091     vtmp(1:nx-1,:)=vtmp(1:nx-1,:)+var(2:nx,:); vtmp(nx,:)=vtmp(nx,:)+var(1,:);
                0092   var=v2(:,:,nt2); var=var.*var; vtmp=vtmp+var;
                0093     vtmp(:,1:ny-1)=vtmp(:,1:ny-1)+var(:,2:ny);
                0094   var=sqrt(vtmp*0.5); titva1='|uVel|';
                0095  end
                0096  var(find(msk1==0))=NaN;
                0097  mnV=min(var(:)); MxV=max(var(:));
                0098  imagesc(xax,yax,var'); set(gca,'YDir','normal');
                0099  if ccB(2) > ccB(1), caxis(ccB/1); end
                0100  change_colmap(-1);
                0101  colorbar;
                0102  grid
                0103  text(xtxt,ytxt,sprintf('min,Max= %9.5g  , %9.5g', mnV, MxV))
                0104  titplot=[titex1,' : ',titva1,' ; ',titim1];
                0105 title(titplot);
                0106 
                0107  if k > 0, var=v2(:,:,k,nt2)-kdif*v1(:,:,k,nt1); else var=v2(:,:,nt2)-kdif*v1(:,:,nt1); end
                0108  var(find(msk1==0))=NaN;
                0109  mnV=min(var(:)); MxV=max(var(:));
                0110 subplot(212);
                0111  imagesc(xax,yax,var'); set(gca,'YDir','normal');
                0112  if cc2(2) > cc2(1) & kdif == 0, caxis(cc2); end
                0113  if ccD(2) > ccD(1) & kdif == 1, caxis(ccD); end
                0114  change_colmap(-1);
                0115  colorbar;
                0116  grid
                0117  text(xtxt,ytxt,sprintf('min,Max= %9.5g  , %9.5g', mnV, MxV))
                0118  if kdif == 1,
                0119    %titplot=['Diff: ',titex2,' - ',titex1];
                0120    tt1=' -'; tt2='';
                0121    if ~strcmp(rDir1,rDir2), tt1=[tt1,' ',titex1];  tt2=[tt2,' ',titex2]; end
                0122    if ~strcmp(namv1,namv2), tt1=[tt1,' ',titva1];  tt2=[tt2,' ',titva2]; end
                0123    if nt1 ~= nt2, tt1=[tt1,' ',titim1];  tt2=[tt2,' ',titim2]; end
                0124    titplot=['Diff:',tt2,tt1];
                0125  else titplot=[titex2,' : ',titva2,' ; ',titim2]; end
                0126 title(titplot);
                0127 put_date;
                0128 %return
                0129 
                0130 clin='kbcmrgy';
                0131 nf=nf+1; xyp0=xyp0+[xysp(1) 0]; xysp=[400 500];
                0132 figure(nf); set(nf,'position',[xyp0 xysp]);clf;
                0133  is=61; is2=21;
                0134 %is=67; is=38;
                0135 subplot(211);
                0136  if k > 0, var=v1(is,:,k,:); else var=v1(is,:,:); end
                0137  var=reshape(var,[ny Nit]);
                0138  m1d=squeeze(msk1(is,:)); [J]=find(m1d==0);
                0139  for nt=1:Nit,
                0140    jt=1+rem(nt-1,length(clin));
                0141    if nt > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
                0142    var(J,nt)=NaN;
                0143    plot(yax,var(:,nt),ccln);
                0144    if nt==1, hold on; end
                0145  end
                0146  hold off;
                0147  AA=axis; axis([0 ny AA(3:4)]);
                0148  grid
                0149 titplot=[titex1,' : ',titva1,' ; is=',int2str(is)];
                0150 title(titplot);
                0151 
                0152 subplot(212);
                0153  if kdif == 1, is2=is; end
                0154  if k > 0, var=v2(is2,:,k,:)-kdif*v1(is,:,k,:); else var=v2(is2,:,:)-kdif*v1(is,:,:); end
                0155  var=reshape(var,[ny Nit]);
                0156  m2d=squeeze(msk1(is2,:)); [J]=find(m2d==0);
                0157  for nt=1:Nit,
                0158    jt=1+rem(nt-1,length(clin));
                0159    if nt > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
                0160    var(J,nt)=NaN;
                0161    plot(yax,var(:,nt),ccln);
                0162    if nt==1, hold on; end
                0163  end
                0164  hold off;
                0165  AA=axis; axis([0 ny AA(3:4)]);
                0166  grid
                0167  if kdif == 1, titplot=['Diff: ',titex2,' - ',titex1];
                0168  else titplot=[titex2,' : ',titva2,' ; is=',int2str(is2)]; end
                0169 title(titplot);
                0170 put_date;
                0171 %return
                0172 
                0173 clin='kbcmrgy';
                0174 nf=nf+1; xyp0=xyp0+[0 xysp(2)];
                0175 figure(nf); set(nf,'position',[xyp0 xysp]);clf;
                0176  js=41; js2=6;
                0177 subplot(211);
                0178  if k > 0, var=v1(:,js,k,:); else var=v1(:,js,:); end
                0179  var=reshape(var,[nx Nit]);
                0180  m1d=squeeze(msk1(:,js)); [J]=find(m1d==0);
                0181  for nt=1:Nit,
                0182    jt=1+rem(nt-1,length(clin));
                0183    if nt > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
                0184    var(J,nt)=NaN;
                0185    plot(xax,var(:,nt),ccln);
                0186    if nt==1, hold on; end
                0187  end
                0188  hold off;
                0189  grid
                0190 titplot=[titex1,' : ',titva1,' ; js=',int2str(js)];
                0191 title(titplot);
                0192 
                0193 subplot(212);
                0194  if kdif == 1, js2=js; end
                0195  if k > 0, var=v2(:,js2,k,:)-kdif*v1(:,js,k,:); else var=v2(:,js2,:)-kdif*v1(:,js,:); end
                0196  var=reshape(var,[nx Nit]);
                0197  m2d=squeeze(msk1(:,js2)); [J]=find(m2d==0);
                0198  for nt=1:Nit,
                0199    jt=1+rem(nt-1,length(clin));
                0200    if nt > length(clin); ccln=[clin(jt:jt),'--']; else ccln=[clin(jt:jt),'-'];end
                0201    var(J,nt)=NaN;
                0202    plot(xax,var(:,nt),ccln);
                0203    if nt==1, hold on; end
                0204  end
                0205  hold off;
                0206  grid
                0207  if kdif == 1, titplot=['Diff: ',titex2,' - ',titex1];
                0208  else titplot=[titex2,' : ',titva2,' ; js=',int2str(js2)]; end
                0209 title(titplot);
                0210 put_date;
                0211 return
                0212 
                0213 for nG=1:2,
                0214  if nG==1, var=v1; titplot=[titex1,' : ',titvar,' ; it=',int2str(nit)];
                0215 %  else var=v2; titplot=[titex2,' : ',titvar,' ; it=',int2str(nit)]; end
                0216    else var=v2-v1; titplot=['Diff: ',titex2,' - ',titex1]; end
                0217  subplot(210+nG);
                0218 
                0219  ccM=max(abs(var));
                0220  if ccM > 0, ccE=floor(log10(ccM)); ccM=ceil(ccM*10^-ccE)*10^ccE; else ccM = 1; end
                0221 %ccM=0;
                0222 %if nG == 2, ccM=ccM/1000; end
                0223  ccB=[-1 1]*ccM;
                0224 plot(xax,var','k');
                0225 
                0226 %mnV=min(min(var)); MxV=max(max(var));
                0227 %fprintf('min,Max= %9.5g  , %9.5g \n', mnV, MxV)
                0228 %text(xtxt,ytxt,sprintf('min,Max= %9.5g  , %9.5g', mnV, MxV))
                0229 
                0230  title(titplot);
                0231 end
                0232