| 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; |
Parent Directory
| 
Revision Log
| 
Revision Graph