Back to home page

MITgcm

 
 

    


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