%%
%% February 2012, /hao/aim3/emery/ctmt/rd_tie.m
%% Read TIEGCM output for CTMT tides 

clear all; close all

glon = -180:5:180;
glat = -87.5:5:87.5;
ctmt_glat = -90:5:90;
mo12 = 1:1:12;
ut = 0;

iflg = 1;
nskp1 = 568;
%infile=sprintf('ctmt_09_0ut_out');
iflg = 2;
nskp1 = 754;
infile=sprintf('ctmt2_09_0ut_out');

IN = fopen(infile,'r');
% Skip over 5 lines (mon and glat)
for n=1:nskp1
 linesk = fgetl(IN);
end
if (iflg==1)
% Read DE3 den at 390km 12 mo for 37 lats
linesk
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  de3amp(mo,jl) = datmo(mo+2);
 end
end
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  de3ph(mo,jl) = datmo(mo+2);
 end
end
% Read DE2 den at 390km 12 mo for 37 lats
linesk = fgetl(IN);
linesk = fgetl(IN)
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  de2amp(mo,jl) = datmo(mo+2);
 end
end
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  de2ph(mo,jl) = datmo(mo+2);
 end
end
% Skip to next (875-644-2*37=157
for n=1:157
 linesk = fgetl(IN);
end
% Read SE2 den at 390km 12 mo for 37 lats
linesk
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  se2amp(mo,jl) = datmo(mo+2);
 end
end
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  se2ph(mo,jl) = datmo(mo+2);
 end
end
% Skip to next (952-875-2*37=3
for n=1:3
 linesk = fgetl(IN);
end
% Read S0 den at 390km 12 mo for 37 lats
linesk
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  s0amp(mo,jl) = datmo(mo+2);
 end
end
for jl=1:37
 linemo = fgetl(IN);
 datmo = sscanf(linemo,'%d %s %f %f %f %f %f %f %f %f %f %f %f %f');
 for mo=1:12
  s0ph(mo,jl) = datmo(mo+2);
 end
end

cmin = 0;
for mt=1:4
figure;
clf;
for nt=1:2
subplot(2,1,nt)
hold on;
if (nt==2)
 cmax = 24
 if (mt>2)
  cmax = 12;
 end
end
if (mt==1)
if (nt==1)
cmax = 2.2;
 contourf(mo12,ctmt_glat,de3amp');
 tit = sprintf('CTMT 390km DE3 amps from TIEGCM');
else
 contourf(mo12,ctmt_glat,de3ph');
 tit = sprintf('CTMT 390km DE3 ph from TIEGCM');
end
end
if (mt==2)
if (nt==1)
cmax = 2.;
 contourf(mo12,ctmt_glat,de2amp');
 tit = sprintf('CTMT 390km DE2 amps from TIEGCM');
else
 contourf(mo12,ctmt_glat,de2ph');
 tit = sprintf('CTMT 390km DE2 ph from TIEGCM');
end
end
if (mt==3)
if (nt==1)
cmax = 3;
 contourf(mo12,ctmt_glat,se2amp');
 tit = sprintf('CTMT 390km SE2 amps from TIEGCM');
else
 contourf(mo12,ctmt_glat,se2ph');
 tit = sprintf('CTMT 390km SE2 ph from TIEGCM');
end
end
if (mt==4)
if (nt==1)
cmax = 4.4;
 contourf(mo12,ctmt_glat,s0amp');
 tit = sprintf('CTMT 390km S0 amps from TIEGCM');
else
 contourf(mo12,ctmt_glat,s0ph');
 tit = sprintf('CTMT 390km S0 ph from TIEGCM');
end
end
 h=title([tit]);
shading flat
caxis([cmin cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
xlabel('month','fontsize',12,'fontweight','bold');
ylabel('latitude','fontsize',12,'fontweight','bold');	
end

hold off;
end	% mt=1:4
	
% Skip to next (1098-952-2*37=72
for n=1:72
 linesk = fgetl(IN);
end
end		% if (iflg==1)

cmax = 57;
cmin = -cmax;

for il=1:72
 linelon = fgetl(IN);
 datlon = sscanf(linelon,'%d %s %f %s');
 gglon(il) = datlon(3)
 for k=1:3
 linelat = fgetl(IN);
 datlat = sscanf(linelat,'%d %s %f %f %f %f %f %f %f %f %f %f');
  for n=1:10
   jl = (k-1)*10 + n;
   u(il,jl) = datlat(n+2);
  end
 end
 linelat = fgetl(IN);  
 datlat = sscanf(linelat,'%d %s %f %f %f %f %f %f');
  for n=1:6
   jl = 30 + n;
   u(il,jl) = datlat(n+2);
  end
end
il = 73;
for jl=1:36
 u(il,jl) = u(1,jl);
end


figure;
clf;

hold on;
 contourf(glon,glat,u');
 tit = sprintf('CTMT 97km Sep zonal tides UT=0 from TIEGCM');
 h=title([tit]);
shading flat
caxis([cmin cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
xlabel('longitude','fontsize',12,'fontweight','bold');
ylabel('latitude','fontsize',12,'fontweight','bold');	


hold off;	




