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; |