/[MITgcm]/MITgcm_contrib/atnguyen/verification/lab_sea/matlab/lookat_adfields_ncgetvar.m
ViewVC logotype

Annotation of /MITgcm_contrib/atnguyen/verification/lab_sea/matlab/lookat_adfields_ncgetvar.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Mon Jun 1 23:25:15 2015 UTC (10 years, 2 months ago) by atn
Branch: MAIN
CVS Tags: HEAD
matlab routines to analyze latest seaice thermo adjoint using IFenty code

1 atn 1.1 %attempt to reproduce plots from
2     %http://mitgcm.org/viewvc/*checkout*/MITgcm/MITgcm_contrib/ifenty/Fenty_seaice_thermo_code_merge/documentation/
3     %Seaice_Groth_Forward_and_Adjoint_comparisons.pdf?revision=1.2
4     %same as lookat_adfields.m but for glaciers
5    
6     clear all
7     dircur=pwd;
8     dirProf='~/matlab_gcmfaces/MITprof';cd(dirProf);MITprof_path;cd(dircur);
9    
10     dirRoot='/net/barents/raid16/atnguyen/MITgcm_c65h/MITgcm/verification/lab_sea/';
11     RunStr=input(['dirRun: "run_ad_stripped_reborn2" ']);dirRun=[dirRoot RunStr '/'];
12     ext=input('ext: ["+crash01/" or empty ');%ext='8784_cost2_icedyn_crash_nomultsnow3/';%ext='8784_ad/';
13     if(length(ext)==0);ext='';else;if(ext(end)~='/');ext=[ext '/'];end;end
14     dirOut=[dirRun ext 'matlab/'];if(exist(dirOut)==0);mkdir(dirOut);end;
15     nx=20;ny=16;dt=3600;
16     gg=ncopen([dirRun ext 'mnc_test_0001/grid.t001.nc'],'nowrite');
17     c=permute(ncgetvar(gg,'HFacC'),[3 2 1]);c(find(c>0))=1;c(find(c==0))=nan;c1=c(:,:,1);ncclose(gg);
18     temp=dir([dirRun ext 'mnc_test_0001/Main_3D_avg_daily*.nc']);
19     idot=find(temp.name=='.');idot=idot(1)+1:idot(2)-1;
20     ts=temp.name(idot);%ts='0000000001';
21    
22     id1=[27 79 131 183 235 287 339];L=length(id1);
23     id2=[53 105 157 209 261 313 365];
24     ix=[10 9 9 8 10 10 10 11 13 14 15 7];Lx=length(ix);
25     iy=[ 8 10 11 12 10 11 12 12 10 12 9 7];
26     cc1='brgmbrgmbrgm';
27     for iloop=1:3;
28     %state
29     if(iloop==1);
30     varname={'SALT','THETA','UVELMASS','VVELMASS'};
31     varname1={'adS','adT','adU','adV'};
32     fInfw=[dirRun ext 'mnc_test_0001/Main_3D_avg_daily.' ts '.t001.nc'];
33     if(exist(fInfw)>0);f1=ncopen(fInfw,'nowrite');end;
34     fInad=[dirRun ext 'mnc_test_0001/adstate.' ts '.t001.nc'];
35     if(exist(fInad)>0);f2=ncopen(fInad,'nowrite');fname2=ncvars(f2);end;
36     strOut='state';k1=1;k2=1;%k2=5;
37     fname=ncvars(f1);
38     elseif(iloop==2);
39     varname={'SIarea','SIheff','SIuice','SIvice','SIhsnow','SItices'};
40     varname1={'adarea','adheff','aduice','advice','adhsnow','adtices'};
41     f1=ncopen([dirRun ext 'mnc_test_0001/ICE_avg_daily.' ts '.t001.nc'],'nowrite');
42     f1u=ncopen([dirRun ext 'mnc_test_0001/ICE_U_avg_daily.' ts '.t001.nc'],'nowrite');
43     f1v=ncopen([dirRun ext 'mnc_test_0001/ICE_V_avg_daily.' ts '.t001.nc'],'nowrite');
44     if(exist(fInad)>0);f2=ncopen([dirRun ext 'mnc_test_0001/adseaice.' ts '.t001.nc'],'nowrite');fname2=ncvars(f2);end;
45     strOut='seaice';
46     elseif(iloop==3);
47     varname={'EXFaqh','EXFatemp','EXFswdn','EXFlwdn','EXFuwind','EXFvwind','EXFpreci'};
48     varname1={'adaqh','adatemp','adswdown','adlwdown','aduwind','advwind','adprecip'};
49     f1=ncopen([dirRun ext 'mnc_test_0001/EXF_avg_daily.' ts '.t001.nc'],'nowrite');
50     if(exist(fInad)>0);f2=ncopen([dirRun ext 'mnc_test_0001/adexf.' ts '.t001.nc'],'nowrite');fname2=ncvars(f2);end;
51     strOut='exf';
52     end;
53     Lv=size(varname,2);
54     for k=1:Lv;
55     tfw=ncgetvar(f1,'T')./dt/24;Ld=length(tfw);
56     if(exist(fInad)>0);tad=(ncgetvar(f2,'T'))./dt/24;else;tad=tfw;end;
57     if(iloop==1);
58     temp=strfind(fname,varname{k});temp=cell2mat(temp);
59     if(temp==1);fw=permute(ncgetvar(f1,varname{k}),[4 3 2 1]);
60     fw0=permute(fw(:,:,k1:k2,:),[1 2 4 3]);fw=squeeze(fw(:,:,k1,:));else;fw=nan(nx,ny,Ld);end;
61     if(exist(fInad)>0);ad=permute(ncgetvar(f2,varname1{k}),[4 3 2 1]);
62     ad0=permute(ad(:,:,k1:k2,:),[1 2 4 3]);ad=squeeze(ad(:,:,k1,:));else;ad=nan.*fw;end;
63     else;
64     if(iloop==2&k==3);
65     fw=permute(squeeze(ncgetvar(f1u,varname{k})),[3 2 1]);
66     elseif(iloop==2&k==4);
67     fw=permute(squeeze(ncgetvar(f1v,varname{k})),[3 2 1]);
68     else
69     fw=permute(squeeze(ncgetvar(f1,varname{k})),[3 2 1]);
70     end;
71     if(exist(fInad)>0);
72     temp=strfind(fname2,varname1{k});temp=cell2mat(temp);
73     if(temp==1);ad=permute(squeeze(ncgetvar(f2,varname1{k})),[3 2 1]);
74     else;ad=nan.*fw;end;
75     else;ad=nan.*fw;end;
76     end;
77     fw1=[];ad1=[];
78     for i=length(id1):-1:1;
79     i1=[];i2=[];temp1=nan(nx,ny);temp2=temp1;
80     i1=find(tfw==id1(i));i2=find(tfw==id2(i));
81     if(length(i1)>0);temp1=fw(1:nx,1:ny,i1);end;if(length(i2)>0);temp2=fw(1:nx,1:ny,i2);end;
82     temp=[squeeze(temp1);squeeze(temp2)].*[c1;c1];fw1=[fw1 temp];
83     end;
84     for i=length(id1):-1:1;
85     i1=[];i2=[];temp1=nan(nx,ny);temp2=temp1;
86     i1=find(tad==id1(i));i2=find(tad==id2(i));
87     if(length(i1)>0);temp1=ad(1:nx,1:ny,i1);end;if(length(i2)>0);temp2=ad(1:nx,1:ny,i2);end;
88     temp=[squeeze(temp1);squeeze(temp2)].*[c1;c1];ad1=[ad1 temp];
89     end;
90     %plotting time-series:
91     if(iloop==1);kk=k2;else;kk=1;end;
92     for klev=1:kk
93     if(klev>1);ad=squeeze(ad0(:,:,:,klev));fw=squeeze(fw0(:,:,:,klev));end;
94     figure((klev-1)*Lv+k+2);%clf;
95     cc0=0;
96     for l=1:3
97     if(l==1);ii=1:4;elseif(l==2);ii=5:8;else;ii=9:12;end;
98     ll=[];
99     for m=1:length(ii);
100     cc0=cc0+1;
101     subplot(2,3,l);plot(tfw,squeeze(fw(ix(ii(m)),iy(ii(m)),1:Ld)),'.-','color',cc1(ii(m)));hold on;
102     %check nan
103     temp=squeeze(ad(ix(ii(m)),iy(ii(m)),:));
104     ij=[];ij=find(isnan(temp)==1);if(length(ij)>0);temp(ij(1):end)=nan;end;
105     subplot(2,3,l+3);plot(tad,temp,'.-','color',cc1(ii(m)));hold on;
106     ll=strvcat(ll,[num2str(cc0) ':[' num2str(ix(ii(m))) ',' num2str(iy(ii(m))) ']']);
107     end;
108     subplot(2,3,l);hold off;grid;ylabel([varname{k} '(k=' num2str(klev) ')']);
109     tenPercentAboveBelowLeftRight;legend(ll,'location','best');
110     subplot(2,3,l+3);hold off;grid;ylabel([varname1{k} '(k=' num2str(klev) ')']);tenPercentAboveBelowLeftRight;
111     end;
112     figure((klev-1)*Lv+k+2);set(gcf,'paperunits','inches','paperposition',[0 0 12 6]);
113     if(klev>1);
114     fpr=[dirOut varname{k} '_' varname1{k} 'klev' myint2str(klev,2) '_ts.png'];print(fpr,'-dpng');
115     else;
116     fpr=[dirOut varname{k} '_' varname1{k} '_ts.png'];print(fpr,'-dpng');
117     end;
118     end;
119    
120     figure(1);if(k==1);clf;colormap(seismic(21));end;
121     subplot(1,Lv,k);imagescnan(fw1','NanColor',[.7,.7,.7]);axis xy;grid;cc=caxis;title(varname{k});
122     if(iloop==1&k==2);caxis([cc(1),10]);elseif(iloop==1&k==1);caxis([30.5,34.5]);end;
123     thincolorbar;set(gca,'Xtick',[0:nx:2*nx],'Ytick',[0:ny:L*ny]);grid;
124     for i=1:length(id1);text(nx-4,i*ny-4,num2str(id1(L-i+1)));text(2*nx-4,i*ny-4,num2str(id2(L-i+1)));end;
125    
126     figure(2);if(k==1);clf;colormap(seismic(21));end;
127     subplot(1,Lv,k);imagescnan(ad1','NanColor',[.7,.7,.7]);axis xy;grid;cc=caxis;title(varname1{k});
128     thincolorbar;set(gca,'Xtick',[0:nx:2*nx],'Ytick',[0:ny:L*ny]);grid;
129     for i=1:length(id1);text(nx-4,i*ny-4,num2str(id1(L-i+1)));text(2*nx-4,i*ny-4,num2str(id2(L-i+1)));end;
130     end;
131     figure(1);set(gcf,'paperunits','inches','paperposition',[0 0 4*Lv 9]);
132     fpr=[dirOut strOut '.png'];print(fpr,'-dpng');fprintf('%s\n',fpr);
133    
134     figure(2);set(gcf,'paperunits','inches','paperposition',[0 0 4*Lv 9]);
135     fpr=[dirOut 'ad' strOut '.png'];print(fpr,'-dpng');fprintf('%s\n',fpr);
136    
137     if(iloop==2);
138     temp=permute(squeeze(ncgetvar(f1,'SIheff')),[3 2 1]);
139     figure(10);clf;imagescnan(temp(:,:,1)'.*c1(:,:,1)');grid;axis xy;
140     for k=1:length(ix);text(ix(k)-.2,iy(k)-.2,num2str(k));end;
141     figure(10);fpr=[dirOut 'Map_index.png'];print(fpr,'-dpng');
142     end;
143     ncclose(f1);ncclose(f2);if(iloop==2);ncclose(f1u);ncclose(f1v);end;
144     fprintf('type <return> to continue\n');
145     keyboard
146     end;%iloop
147    
148     lookat_movie=0;
149     if(lookat_movie==1);
150     f2=ncopen([dirRun 'mnc_test_0001/adstate.0000000001.t001.nc'],'nowrite');
151     vv='adS';
152     h=permute(ncgetvar(f2,vv),[4 3 2 1]);
153     for k=3279:3281;
154     figure(10);clf;imagescnan([squeeze(h(1:nx,1:ny,1,k-1))',squeeze(h(1:nx,1:ny,1,k))']);axis xy;thincolorbar;
155     title([vv ' L: ' num2str(k-1) ' R: ' num2str(k) ',' num2str(366-k/24)]);
156     pause;
157     end;
158     ncclose(f2);
159     end;

  ViewVC Help
Powered by ViewVC 1.1.22