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