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

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

  ViewVC Help
Powered by ViewVC 1.1.22