%%
%% February 2013, /hao/aim3/emery/ctmt/rd_jens.m
%% Read Jens Oberheide's CTMT .nc files 
%% Aug 2014, add U,W,T,D 100-400km plots for Fig 7 S0

clear all; close all

% data(25,145,73)=(LT,glon,glat) where LT=0:1:24, glon=0:2.5:360, glat=-90:5:90
% am_##_v(month,alt,glat)=(12,161,37) where month=1:1:12, alt=0:2.5:400, glat=-90:5:90
mon = 1:1:12;
glat = -90:5:90;
alt=0:2.5:400;
ncd=netcdf('tira_diurnal_2002_2008.nc','nowrite');
ncsd=netcdf('tira_semidiurnal_2002_2008.nc','nowrite');

lat = ncd{'lat'}(:);
lev = ncd{'lev'}(:);
month = ncd{'monrh'}(:);
am_w2_t = ncd{'amp_w2_t'}(:);
ph_w2_t = ncd{'phase_w2_t'}(:);
am_w2_u = ncd{'amp_w2_u'}(:);
ph_w2_u = ncd{'phase_w2_u'}(:);
am_w2_v = ncd{'amp_w2_v'}(:);
ph_w2_v = ncd{'phase_w2_v'}(:);
am_w2_w = ncd{'amp_w2_w'}(:);
ph_w2_w = ncd{'phase_w2_w'}(:);
am_w2_d = ncd{'amp_w2_d'}(:);
ph_w2_d = ncd{'phase_w2_d'}(:);
am_w1_t = ncd{'amp_w1_t'}(:);
ph_w1_t = ncd{'phase_w1_t'}(:);
am_w1_u = ncd{'amp_w1_u'}(:);
ph_w1_u = ncd{'phase_w1_u'}(:);
am_w1_v = ncd{'amp_w1_v'}(:);
ph_w1_v = ncd{'phase_w1_v'}(:);
am_w1_w = ncd{'amp_w1_w'}(:);
ph_w1_w = ncd{'phase_w1_w'}(:);
am_w1_d = ncd{'amp_w1_d'}(:);
ph_w1_d = ncd{'phase_w1_d'}(:);
am_s0_t = ncd{'amp_s0_t'}(:);
ph_s0_t = ncd{'phase_s0_t'}(:);
am_s0_u = ncd{'amp_s0_u'}(:);
ph_s0_u = ncd{'phase_s0_u'}(:);
am_s0_v = ncd{'amp_s0_v'}(:);
ph_s0_v = ncd{'phase_s0_v'}(:);
am_s0_w = ncd{'amp_s0_w'}(:);
ph_s0_w = ncd{'phase_s0_w'}(:);
am_s0_d = ncd{'amp_s0_d'}(:);
ph_s0_d = ncd{'phase_s0_d'}(:);
am_e1_t = ncd{'amp_e1_t'}(:);
ph_e1_t = ncd{'phase_e1_t'}(:);
am_e1_u = ncd{'amp_e1_u'}(:);
ph_e1_u = ncd{'phase_e1_u'}(:);
am_e1_v = ncd{'amp_e1_v'}(:);
ph_e1_v = ncd{'phase_e1_v'}(:);
am_e1_w = ncd{'amp_e1_w'}(:);
ph_e1_w = ncd{'phase_e1_w'}(:);
am_e1_d = ncd{'amp_e1_d'}(:);
ph_e1_d = ncd{'phase_e1_d'}(:);
am_e2_t = ncd{'amp_e2_t'}(:);
ph_e2_t = ncd{'phase_e2_t'}(:);
am_e2_u = ncd{'amp_e2_u'}(:);
ph_e2_u = ncd{'phase_e2_u'}(:);
am_e2_v = ncd{'amp_e2_v'}(:);
ph_e2_v = ncd{'phase_e2_v'}(:);
am_e2_w = ncd{'amp_e2_w'}(:);
ph_e2_w = ncd{'phase_e2_w'}(:);
am_e2_d = ncd{'amp_e2_d'}(:);
ph_e2_d = ncd{'phase_e2_d'}(:);
am_e3_t = ncd{'amp_e3_t'}(:);
ph_e3_t = ncd{'phase_e3_t'}(:);
am_e3_u = ncd{'amp_e3_u'}(:);
ph_e3_u = ncd{'phase_e3_u'}(:);
am_e3_v = ncd{'amp_e3_v'}(:);
ph_e3_v = ncd{'phase_e3_v'}(:);
am_e3_w = ncd{'amp_e3_w'}(:);
ph_e3_w = ncd{'phase_e3_w'}(:);
am_e3_d = ncd{'amp_e3_d'}(:);
ph_e3_d = ncd{'phase_e3_d'}(:);
% 8 Semidiurnal tides (12,161,37)
am2_w4_t = ncsd{'amp_w4_t'}(:);
ph2_w4_t = ncsd{'phase_w4_t'}(:);
am2_w3_t = ncsd{'amp_w3_t'}(:);
ph2_w3_t = ncsd{'phase_w3_t'}(:);
am2_w2_t = ncsd{'amp_w2_t'}(:);
ph2_w2_t = ncsd{'phase_w2_t'}(:);
am2_w1_t = ncsd{'amp_w1_t'}(:);
ph2_w1_t = ncsd{'phase_w1_t'}(:);
am2_s0_t = ncsd{'amp_s0_t'}(:);
ph2_s0_t = ncsd{'phase_s0_t'}(:);
am2_e1_t = ncsd{'amp_e1_t'}(:);
ph2_e1_t = ncsd{'phase_e1_t'}(:);
am2_e2_t = ncsd{'amp_e2_t'}(:);
ph2_e2_t = ncsd{'phase_e2_t'}(:);
am2_e3_t = ncsd{'amp_e3_t'}(:);
ph2_e3_t = ncsd{'phase_e3_t'}(:);
am2_w4_u = ncsd{'amp_w4_u'}(:);
ph2_w4_u = ncsd{'phase_w4_u'}(:);
am2_w3_u = ncsd{'amp_w3_u'}(:);
ph2_w3_u = ncsd{'phase_w3_u'}(:);
am2_w2_u = ncsd{'amp_w2_u'}(:);
ph2_w2_u = ncsd{'phase_w2_u'}(:);
am2_w1_u = ncsd{'amp_w1_u'}(:);
ph2_w1_u = ncsd{'phase_w1_u'}(:);
am2_s0_u = ncsd{'amp_s0_u'}(:);
ph2_s0_u = ncsd{'phase_s0_u'}(:);
am2_e1_u = ncsd{'amp_e1_u'}(:);
ph2_e1_u = ncsd{'phase_e1_u'}(:);
am2_e2_u = ncsd{'amp_e2_u'}(:);
ph2_e2_u = ncsd{'phase_e2_u'}(:);
am2_e3_u = ncsd{'amp_e3_u'}(:);
ph2_e3_u = ncsd{'phase_e3_u'}(:);
am2_w4_v = ncsd{'amp_w4_v'}(:);
ph2_w4_v = ncsd{'phase_w4_v'}(:);
am2_w3_v = ncsd{'amp_w3_v'}(:);
ph2_w3_v = ncsd{'phase_w3_v'}(:);
am2_w2_v = ncsd{'amp_w2_v'}(:);
ph2_w2_v = ncsd{'phase_w2_v'}(:);
am2_w1_v = ncsd{'amp_w1_v'}(:);
ph2_w1_v = ncsd{'phase_w1_v'}(:);
am2_s0_v = ncsd{'amp_s0_v'}(:);
ph2_s0_v = ncsd{'phase_s0_v'}(:);
am2_e1_v = ncsd{'amp_e1_v'}(:);
ph2_e1_v = ncsd{'phase_e1_v'}(:);
am2_e2_v = ncsd{'amp_e2_v'}(:);
ph2_e2_v = ncsd{'phase_e2_v'}(:);
am2_e3_v = ncsd{'amp_e3_v'}(:);
ph2_e3_v = ncsd{'phase_e3_v'}(:);
am2_w4_w = ncsd{'amp_w4_w'}(:);
ph2_w4_w = ncsd{'phase_w4_w'}(:);
am2_w3_w = ncsd{'amp_w3_w'}(:);
ph2_w3_w = ncsd{'phase_w3_w'}(:);
am2_w2_w = ncsd{'amp_w2_w'}(:);
ph2_w2_w = ncsd{'phase_w2_w'}(:);
am2_w1_w = ncsd{'amp_w1_w'}(:);
ph2_w1_w = ncsd{'phase_w1_w'}(:);
am2_s0_w = ncsd{'amp_s0_w'}(:);
ph2_s0_w = ncsd{'phase_s0_w'}(:);
am2_e1_w = ncsd{'amp_e1_w'}(:);
ph2_e1_w = ncsd{'phase_e1_w'}(:);
am2_e2_w = ncsd{'amp_e2_w'}(:);
ph2_e2_w = ncsd{'phase_e2_w'}(:);
am2_e3_w = ncsd{'amp_e3_w'}(:);
ph2_e3_w = ncsd{'phase_e3_w'}(:);
am2_w4_d = ncsd{'amp_w4_d'}(:);
ph2_w4_d = ncsd{'phase_w4_d'}(:);
am2_w3_d = ncsd{'amp_w3_d'}(:);
ph2_w3_d = ncsd{'phase_w3_d'}(:);
am2_w2_d = ncsd{'amp_w2_d'}(:);
ph2_w2_d = ncsd{'phase_w2_d'}(:);
am2_w1_d = ncsd{'amp_w1_d'}(:);
ph2_w1_d = ncsd{'phase_w1_d'}(:);
am2_s0_d = ncsd{'amp_s0_d'}(:);
ph2_s0_d = ncsd{'phase_s0_d'}(:);
am2_e1_d = ncsd{'amp_e1_d'}(:);
ph2_e1_d = ncsd{'phase_e1_d'}(:);
am2_e2_d = ncsd{'amp_e2_d'}(:);
ph2_e2_d = ncsd{'phase_e2_d'}(:);
am2_e3_d = ncsd{'amp_e3_d'}(:);
ph2_e3_d = ncsd{'phase_e3_d'}(:);


% Plot amp and phase of Sep sd s0 from 85-105km=alt(35-43) or 100-400km to match Figs 6 and 7 of  Oberheide et al, JGR, 2011
for il=1:37
 for nh=35:43
  alt85_105(nh-34) = alt(nh);
  sep85_v_sd_amp_s0(il,nh-34) = am2_s0_v(9,nh,il);
  sep85_v_sd_phs_s0(il,nh-34) = ph2_s0_v(9,nh,il);
 end
 for nh=41:161
  alt100_400(nh-40) = alt(nh);
  sep100_v_sd_amp_s0(il,nh-40) = am2_s0_v(9,nh,il);
  sep100_v_sd_phs_s0(il,nh-40) = ph2_s0_v(9,nh,il);
  sep100_u_sd_amp_s0(il,nh-40) = am2_s0_u(9,nh,il);
  sep100_u_sd_phs_s0(il,nh-40) = ph2_s0_u(9,nh,il);
  sep100_w_sd_amp_s0(il,nh-40) = am2_s0_w(9,nh,il);
  sep100_w_sd_phs_s0(il,nh-40) = ph2_s0_w(9,nh,il);
  sep100_t_sd_amp_s0(il,nh-40) = am2_s0_t(9,nh,il);
  sep100_t_sd_phs_s0(il,nh-40) = ph2_s0_t(9,nh,il);
  sep100_d_sd_amp_s0(il,nh-40) = am2_s0_d(9,nh,il);
  sep100_d_sd_phs_s0(il,nh-40) = ph2_s0_d(9,nh,il);
 end
end

figure;
clf;

for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glat,alt85_105,sep85_v_sd_amp_s0');
 tit = sprintf('CTMT Sep sd v amp');
 cmax = 10;
end
if (nt==2)
 contourf(glat,alt85_105,sep85_v_sd_phs_s0');
 tit = sprintf('CTMT Sep sd v phase');
 cmax = 12;
end
if (nt==3)
 contourf(glat,alt100_400,sep100_v_sd_amp_s0');
 tit = sprintf('CTMT Sep sd v amp');
 cmax = 10;
end
if (nt==4)
 contourf(glat,alt100_400,sep100_v_sd_phs_s0');
 tit = sprintf('CTMT Sep sd v phase');
 cmax = 12;
end
 h=title([tit]);
shading flat
caxis([0 cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('altitude (km)','fontsize',12,'fontweight','bold');
xlabel('latitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:4
hold off;	

% Plot Amp and Ph for U,W,T,D 100-400km Fig 7
figure;
clf;

for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glat,alt100_400,sep100_u_sd_amp_s0');
 tit = sprintf('CTMT Sep sd u amp');
 cmax = 6;
end
if (nt==2)
 contourf(glat,alt100_400,sep100_w_sd_amp_s0');
 tit = sprintf('CTMT Sep sd w amp');
 cmax = 0.14;
end
if (nt==3)
 contourf(glat,alt100_400,sep100_t_sd_amp_s0');
 tit = sprintf('CTMT Sep sd t amp');
 cmax = 8;
end
if (nt==4)
 contourf(glat,alt100_400,sep100_d_sd_amp_s0');
 tit = sprintf('CTMT Sep sd %d amp');
 cmax = 0.04;
end
 h=title([tit]);
shading flat
caxis([0 cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('altitude (km)','fontsize',12,'fontweight','bold');
xlabel('latitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:4
hold off;	

figure;
clf;
cmax = 12;

for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glat,alt100_400,sep100_u_sd_phs_s0');
 tit = sprintf('CTMT Sep sd u phase');
end
if (nt==2)
 contourf(glat,alt100_400,sep100_w_sd_phs_s0');
 tit = sprintf('CTMT Sep sd w phase');
end
if (nt==3)
 contourf(glat,alt100_400,sep100_t_sd_phs_s0');
 tit = sprintf('CTMT Sep sd t phase');
end
if (nt==4)
 contourf(glat,alt100_400,sep100_d_sd_phs_s0');
 tit = sprintf('CTMT Sep sd %d phase');
end
 h=title([tit]);
shading flat
caxis([0 cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('altitude (km)','fontsize',12,'fontweight','bold');
xlabel('latitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:4
hold off;

% TEMP	
stop

% Plot amp and phase of monthly 390km (157) density for DE2, DE3, S0, SE2 to match Figs 8-11 of  Oberheide et al, JGR, 2011
nh=157;
for il=1:37
 for mo=1:12
  mo390_amp_de2(mo,il) = am_e2_d(mo,nh,il)*100.;
  mo390_phs_de2(mo,il) = ph_e2_d(mo,nh,il);
  mo390_amp_de3(mo,il) = am_e3_d(mo,nh,il)*100.;
  mo390_phs_de3(mo,il) = ph_e3_d(mo,nh,il);
  mo390_amp_s0(mo,il) = am2_s0_d(mo,nh,il)*100.;
  mo390_phs_s0(mo,il) = ph2_s0_d(mo,nh,il);
  mo390_amp_se2(mo,il) = am2_e2_d(mo,nh,il)*100.;
  mo390_phs_se2(mo,il) = ph2_e2_d(mo,nh,il);
 end
end

figure;
clf;
for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(month,glat,mo390_amp_de2');
 tit = sprintf('CTMT DE2 amp');
 cmax = 2;
end
if (nt==2)
 contourf(month,glat,mo390_phs_de2');
 tit = sprintf('CTMT DE2 phase');
 cmax = 24;
end
if (nt==3)
 contourf(month,glat,mo390_amp_de3');
 tit = sprintf('CTMT DE3 amp');
 cmax = 2.2;
end
if (nt==4)
 contourf(month,glat,mo390_phs_de3');
 tit = sprintf('CTMT DE3 phase');
 cmax = 24;
end
 h=title([tit]);
shading flat
caxis([0 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	% for nt=1:4
hold off;


figure;
clf;
for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(month,glat,mo390_amp_s0');
 tit = sprintf('CTMT S0 amp');
 cmax = 4.4;
end
if (nt==2)
 contourf(month,glat,mo390_phs_s0');
 tit = sprintf('CTMT S0 phase');
 cmax = 12;
end
if (nt==3)
 contourf(month,glat,mo390_amp_se2');
 tit = sprintf('CTMT SE2 amp');
 cmax = 3;
end
if (nt==4)
 contourf(month,glat,mo390_phs_se2');
 tit = sprintf('CTMT SE2 phase');
 cmax = 12;
end
 h=title([tit]);
shading flat
caxis([0 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	% for nt=1:4
hold off;	
	
% Sum Sep tides at 100(#41) and 250km (#101) for LT=0,6,12,18
% D=e3,e2,e1,s0,w1,w2; SD=e3,e2,e1,s0,w1,w2,w3,w4
%glon=-180:5:180;
glon=0:5:360;

mo = 9;
nh = 40;
lt0 = 0;
lt12 = 12;
for il=1:37
 for jlon=1:73
% for LT=0, lt0=0 for jlon=37, and lt=glon(jlon)/15.+lt0 otherwise
%  lt =lt0 + glon(jlon)/15.;	% for real glon lt changes 
   ut = lt0 - glon(jlon)/15.;	% for UT where LT=lT0
% ut = lt0 gives a plot for a constant UT, while ut=lt0-glon/15 is for a constant LT
%  BUT PROBLEM SINCE s0 is constant for const UT which is WRONG!!??
%  ut = lt0;
% cos(omega*(phase-UT) +/- s*glon) where omega=1 for D, 2 for SD, s=0-4; -=W,+=E
%  DEF is: cos[n*omega*UT+s*glon-phase] so + is for W (NOT E), and W1=migrating diurnal
%    and W2=migrating sd
% X Zhang diff definition of tides so tide=amp*cos[n*pi/12*(LT-ph)] where
%  ph is not the usual phase, but ph={phase*n*pi/12-(s-n)*glon}*12/(pi*n)
%  where n=1,2 for d,sd
%  Thus tide=amp*cos[n*omega*LT +(s-n)glon -phase] or amp*cos[n*omega*UT +s*glon -phase]
tide100_00lt_d(jlon,il)=100.*am_e3_d(mo,nh,il)*cos(2.*pi*(ph_e3_d(mo,nh,il)-ut)/24. + 3.*glon(jlon)*2*pi/360.) + ...
              	100.*am_e2_d(mo,nh,il)*cos(2.*pi*(ph_e2_d(mo,nh,il)-ut)/24. + 2.*glon(jlon)*2*pi/360.) + ...
		100.*am_e1_d(mo,nh,il)*cos(2.*pi*(ph_e1_d(mo,nh,il)-ut)/24. + glon(jlon)*2*pi/360.) + ...
		100.*am_s0_d(mo,nh,il)*cos(2.*pi*(ph_s0_d(mo,nh,il)-ut)/24.) + ...
		100.*am_w1_d(mo,nh,il)*cos(2.*pi*(ph_w1_d(mo,nh,il)-ut)/24. - glon(jlon)*2*pi/360.) + ...
		100.*am_w2_d(mo,nh,il)*cos(2.*pi*(ph_w2_d(mo,nh,il)-ut)/24. - 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e3_d(mo,nh,il)*cos(2.*pi*(ph2_e3_d(mo,nh,il)-ut)/12. + 3.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e2_d(mo,nh,il)*cos(2.*pi*(ph2_e2_d(mo,nh,il)-ut)/12. + 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e1_d(mo,nh,il)*cos(2.*pi*(ph2_e1_d(mo,nh,il)-ut)/12. + glon(jlon)*2*pi/360.) + ...
		100.*am2_s0_d(mo,nh,il)*cos(2.*pi*(ph2_s0_d(mo,nh,il)-ut)/12.) + ...
		100.*am2_w1_d(mo,nh,il)*cos(2.*pi*(ph2_w1_d(mo,nh,il)-ut)/12. - glon(jlon)*2*pi/360.) + ...
		100.*am2_w2_d(mo,nh,il)*cos(2.*pi*(ph2_w2_d(mo,nh,il)-ut)/12. - 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_w3_d(mo,nh,il)*cos(2.*pi*(ph2_w3_d(mo,nh,il)-ut)/12. - 3.*glon(jlon)*2*pi/360.) + ...
		100.*am2_w4_d(mo,nh,il)*cos(2.*pi*(ph2_w4_d(mo,nh,il)-ut)/12. - 4.*glon(jlon)*2*pi/360.); 
tide100_00lt_u(jlon,il)=am_e3_u(mo,nh,il)*cos(2.*pi*(ph_e3_u(mo,nh,il)-ut)/24. + 3.*glon(jlon)*2*pi/360.) + ...
              	am_e2_u(mo,nh,il)*cos(2.*pi*(ph_e2_u(mo,nh,il)-ut)/24. + 2.*glon(jlon)*2*pi/360. ) + ...
		am_e1_u(mo,nh,il)*cos(2.*pi*(ph_e1_u(mo,nh,il)-ut)/24. + glon(jlon)*2*pi/360.) + ...
		am_s0_u(mo,nh,il)*cos(2.*pi*(ph_s0_u(mo,nh,il)-ut)/24.) + ...
		am_w1_u(mo,nh,il)*cos(2.*pi*(ph_w1_u(mo,nh,il)-ut)/24. - glon(jlon)*2*pi/360.) + ...
		am_w2_u(mo,nh,il)*cos(2.*pi*(ph_w2_u(mo,nh,il)-ut)/24. - 2.*glon(jlon)*2*pi/360.) + ...
		am2_e3_u(mo,nh,il)*cos(2.*pi*(ph2_e3_u(mo,nh,il)-ut)/12. + 3.*glon(jlon)*2*pi/360.) + ...
		am2_e2_u(mo,nh,il)*cos(2.*pi*(ph2_e2_u(mo,nh,il)-ut)/12. + 2.*glon(jlon)*2*pi/360.) + ...
		am2_e1_u(mo,nh,il)*cos(2.*pi*(ph2_e1_u(mo,nh,il)-ut)/12. + glon(jlon)*2*pi/360.) + ...
		am2_s0_u(mo,nh,il)*cos(2.*pi*(ph2_s0_u(mo,nh,il)-ut)/12.) + ...
		am2_w1_u(mo,nh,il)*cos(2.*pi*(ph2_w1_u(mo,nh,il)-ut)/12. - glon(jlon)*2*pi/360.) + ...
		am2_w2_u(mo,nh,il)*cos(2.*pi*(ph2_w2_u(mo,nh,il)-ut)/12. - 2.*glon(jlon)*2*pi/360.) + ...
		am2_w3_u(mo,nh,il)*cos(2.*pi*(ph2_w3_u(mo,nh,il)-ut)/12. - 3.*glon(jlon)*2*pi/360.) + ...
		am2_w4_u(mo,nh,il)*cos(2.*pi*(ph2_w4_u(mo,nh,il)-ut)/12. - 4.*glon(jlon)*2*pi/360.); 
   ut = lt12 - glon(jlon)/15.;	% for UT where LT=lT12
tide100_12lt_d(jlon,il)=100.*am_e3_d(mo,nh,il)*cos(2.*pi*(ph_e3_d(mo,nh,il)-ut)/24. + 3.*glon(jlon)*2*pi/360.) + ...
              	100.*am_e2_d(mo,nh,il)*cos(2.*pi*(ph_e2_d(mo,nh,il)-ut)/24. + 2.*glon(jlon)*2*pi/360.) + ...
		100.*am_e1_d(mo,nh,il)*cos(2.*pi*(ph_e1_d(mo,nh,il)-ut)/24. + glon(jlon)*2*pi/360.) + ...
		100.*am_s0_d(mo,nh,il)*cos(2.*pi*(ph_s0_d(mo,nh,il)-ut)/24.) + ...
		100.*am_w1_d(mo,nh,il)*cos(2.*pi*(ph_w1_d(mo,nh,il)-ut)/24. - glon(jlon)*2*pi/360.) + ...
		100.*am_w2_d(mo,nh,il)*cos(2.*pi*(ph_w2_d(mo,nh,il)-ut)/24. - 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e3_d(mo,nh,il)*cos(2.*pi*(ph2_e3_d(mo,nh,il)-ut)/12. + 3.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e2_d(mo,nh,il)*cos(2.*pi*(ph2_e2_d(mo,nh,il)-ut)/12. + 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_e1_d(mo,nh,il)*cos(2.*pi*(ph2_e1_d(mo,nh,il)-ut)/12. + glon(jlon)*2*pi/360.) + ...
		100.*am2_s0_d(mo,nh,il)*cos(2.*pi*(ph2_s0_d(mo,nh,il)-ut)/12.) + ...
		100.*am2_w1_d(mo,nh,il)*cos(2.*pi*(ph2_w1_d(mo,nh,il)-ut)/12. - glon(jlon)*2*pi/360.) + ...
		100.*am2_w2_d(mo,nh,il)*cos(2.*pi*(ph2_w2_d(mo,nh,il)-ut)/12. - 2.*glon(jlon)*2*pi/360.) + ...
		100.*am2_w3_d(mo,nh,il)*cos(2.*pi*(ph2_w3_d(mo,nh,il)-ut)/12. - 3.*glon(jlon)*2*pi/360.) + ...
		100.*am2_w4_d(mo,nh,il)*cos(2.*pi*(ph2_w4_d(mo,nh,il)-ut)/12. - 4.*glon(jlon)*2*pi/360.); 
tide100_12lt_u(jlon,il)=am_e3_u(mo,nh,il)*cos(2.*pi*(ph_e3_u(mo,nh,il)-ut)/24. + 3.*glon(jlon)*2*pi/360.) + ...
              	am_e2_u(mo,nh,il)*cos(2.*pi*(ph_e2_u(mo,nh,il)-ut)/24. + 2.*glon(jlon)*2*pi/360. ) + ...
		am_e1_u(mo,nh,il)*cos(2.*pi*(ph_e1_u(mo,nh,il)-ut)/24. + glon(jlon)*2*pi/360.) + ...
		am_s0_u(mo,nh,il)*cos(2.*pi*(ph_s0_u(mo,nh,il)-ut)/24.) + ...
		am_w1_u(mo,nh,il)*cos(2.*pi*(ph_w1_u(mo,nh,il)-ut)/24. - glon(jlon)*2*pi/360.) + ...
		am_w2_u(mo,nh,il)*cos(2.*pi*(ph_w2_u(mo,nh,il)-ut)/24. - 2.*glon(jlon)*2*pi/360.) + ...
		am2_e3_u(mo,nh,il)*cos(2.*pi*(ph2_e3_u(mo,nh,il)-ut)/12. + 3.*glon(jlon)*2*pi/360.) + ...
		am2_e2_u(mo,nh,il)*cos(2.*pi*(ph2_e2_u(mo,nh,il)-ut)/12. + 2.*glon(jlon)*2*pi/360.) + ...
		am2_e1_u(mo,nh,il)*cos(2.*pi*(ph2_e1_u(mo,nh,il)-ut)/12. + glon(jlon)*2*pi/360.) + ...
		am2_s0_u(mo,nh,il)*cos(2.*pi*(ph2_s0_u(mo,nh,il)-ut)/12.) + ...
		am2_w1_u(mo,nh,il)*cos(2.*pi*(ph2_w1_u(mo,nh,il)-ut)/12. - glon(jlon)*2*pi/360.) + ...
		am2_w2_u(mo,nh,il)*cos(2.*pi*(ph2_w2_u(mo,nh,il)-ut)/12. - 2.*glon(jlon)*2*pi/360.) + ...
		am2_w3_u(mo,nh,il)*cos(2.*pi*(ph2_w3_u(mo,nh,il)-ut)/12. - 3.*glon(jlon)*2*pi/360.) + ...
		am2_w4_u(mo,nh,il)*cos(2.*pi*(ph2_w4_u(mo,nh,il)-ut)/12. - 4.*glon(jlon)*2*pi/360.); 
 end
end


figure;
clf;
for nt=1:4

subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glon,glat,tide100_00lt_d');
 tit = sprintf('0LT CTMT pc den 97.5km');
 cmax = 20;
end
if (nt==2)
 contourf(glon,glat,tide100_00lt_u');
 tit = sprintf('0LT CTMT u 97.5km');
 cmax = 60;
end
if (nt==3)
 contourf(glon,glat,tide100_12lt_d');
 tit = sprintf('12LT CTMT pc den 97.5km');
 cmax = 20;
end
if (nt==4)
 contourf(glon,glat,tide100_12lt_u');
 tit = sprintf('12LT CTMT u 97.5km');
 cmax = 60;
end
 h=title([tit]);
shading flat
caxis([-cmax cmax])	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
xlabel('longitude','fontsize',12,'fontweight','bold');
ylabel('latitude','fontsize',12,'fontweight','bold');	

end	% for nt=1:2
hold off;	
	

