%%
%% February 2012, /hao/aim3/emery/ctmt/rd_jens.m
%% Read Jens Oberheide's CTMT .nc files 

% Calculating the zonal mean ends up all zeros since tides have no zonal mean!

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, 
mon = 1:1:12;
glat = -90:5:90;
glon=-180:5:180;
%glon=0:5:360;

% Read zonal means of WINDII-HRDI winds at 95km
INhrdi=fopen('/home/emery/cnofs/mackjr/WINDII-HRDI_zmUV_95km.dat','r');
lathrdi=-70:10:70;
nhlat = 15;
for il=1:15
 linhrdi = fgetl(INhrdi);
 dat1hrdi = sscanf(linhrdi,'%f');
 linhrdi = fgetl(INhrdi);
 dat2hrdi = sscanf(linhrdi,'%f');
 for mo=1:6
  vn95_hrdi(mo,il) = dat1hrdi(mo);
  vn95_hrdi(mo+6,il) = dat2hrdi(mo);
 end
end
for il=1:15
 linhrdi = fgetl(INhrdi);
 dat1hrdi = sscanf(linhrdi,'%f');
 linhrdi = fgetl(INhrdi);
 dat2hrdi = sscanf(linhrdi,'%f');
 for mo=1:6
  un95_hrdi(mo,il) = dat1hrdi(mo);
  un95_hrdi(mo+6,il) = dat2hrdi(mo);
 end
end


% Read zonal means of TIDI winds at 95km (#7/11)
INtidi=fopen('/home/emery/tidi/tidi_un_vn_zav0_vs_alt','r');
lattidi=-87:6:87;
ntlat = 30;
  lintidi = fgetl(INtidi);
n95 = 9;
for il=1:30
 for nd=1:4	% (4 day numbers 80,172,263,355)
  lintidi = fgetl(INtidi);
  dat1tidi = sscanf(lintidi,'%d %d %f %f %f %f %f %f %f %f %f %f %f');
  un95_tidi(nd,il) = dat1tidi(n95);
  lintidi = fgetl(INtidi);
  dat1tidi = sscanf(lintidi,'%d %d %f %f %f %f %f %f %f %f %f %f %f');
  vn95_tidi(nd,il) = dat1tidi(n95);
 end	% for nd=1:4 	% (4 day numbers 80,172,263,355)
end

% Calculate linear interp for zonal means to glat

il = 1;
  un95t_mar(il) = un95_tidi(1,il);
  un95h_mar(il) = un95_tidi(1,il);
  vn95t_mar(il) = vn95_tidi(1,il);
  vn95h_mar(il) = vn95_tidi(1,il);
il=37;
  un95t_mar(il) = un95_tidi(1,ntlat);
  un95h_mar(il) = un95_tidi(1,ntlat);
  vn95t_mar(il) = vn95_tidi(1,ntlat);
  vn95h_mar(il) = vn95_tidi(1,ntlat);
for il=2:36
   jlt = 0;
 for jl=1:ntlat-1
  if (glat(il)>=lattidi(jl) & glat(il)<=lattidi(jl+1))
   frac1 = (lattidi(jl+1)-glat(il))/(lattidi(jl+1)-lattidi(jl));
   frac2 = 1.-frac1;
   jlt = jl;
  end
 end
  if (jlt==0)
stop
  end
  un95t_mar(il) = frac1*un95_tidi(1,jlt) + frac2*un95_tidi(1,jlt+1);
  vn95t_mar(il) = frac1*vn95_tidi(1,jlt) + frac2*vn95_tidi(1,jlt+1);
  if (glat(il)<lathrdi(1) | glat(il)>lathrdi(nhlat))
   un95h_mar(il) = un95t_mar(il);
   vn95h_mar(il) = vn95t_mar(il);
  else
 for jl=1:nhlat-1
  if (glat(il)>=lathrdi(jl) & glat(il)<=lathrdi(jl+1))
   frac1 = (lathrdi(jl+1)-glat(il))/(lathrdi(jl+1)-lathrdi(jl));
   frac2 = 1.-frac1;
   jlt = jl;
  end
 end
  if (jlt==0)
stop
  end
  un95h_mar(il) = frac1*un95_hrdi(3,jlt) + frac2*un95_hrdi(3,jlt+1);
  vn95h_mar(il) = frac1*vn95_hrdi(3,jlt) + frac2*vn95_hrdi(3,jlt+1);
  end
end	% for il=2:36

figure;
set (gcf,'position',[142 99 647 822])
H=gcf;
clf;

zero = [0,0];
totglat = [-90,90];

subplot(2,1,1);
hold on;
 plot (totglat,zero,':b');
 plot (glat,un95h_mar,'--r');
 plot (glat,un95t_mar,'-b');
 ylabel('Un 95km ZM (m/s)','FontSize',12);
 xlabel('glat','FontSize',12);
 title('Zonal means: WINDII-HRDI(--red), TIDI(blue)','FontSize',15);
hold off;
subplot(2,1,2);
hold on;
 plot (totglat,zero,':b');
 plot (glat,vn95h_mar,'--r');
 plot (glat,vn95t_mar,'-b');
 ylabel('Vn 95km ZM (m/s)','FontSize',12);
 xlabel('glat','FontSize',12);
hold off;



%OUT = fopen('2_ctmt_un_vn_ut0_105km_glon_vs_glat','w');
OUT = fopen('ctmt_un_vn_ut0_95km_glon_vs_glat','w');
% Calculate and print out UT=0 at 105km (#43) for Vn and Un for plot_alt_planes.m
%fprintf (OUT,'CTMT Un and Vn for 15 March at UT=0 105km for 73 glon -180:5:180 and 37 glat -90:5:90; glat; glon Un; glon Vn \n');
fprintf (OUT,'CTMT Un and Vn for 15 March at UT=0 95km for 73 glon -180:5:180 and 37 glat -90:5:90; glat; u,vn95_tidi_mar; u,vn95_hrdi_mar; glon Un; glon Vn \n');
fprintf (OUT,'%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',glat);

fprintf (OUT,'%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',un95t_mar);
fprintf (OUT,'%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',vn95t_mar);
fprintf (OUT,'%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',un95h_mar);
fprintf (OUT,'%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',vn95h_mar);


alt=0:2.5:400;

nvt = 2;
nvu = 3;
nvv = 4;
nvd = 1;

ncd=netcdf('tira_diurnal_2002_2008.nc','nowrite');
ncsd=netcdf('tira_semidiurnal_2002_2008.nc','nowrite');


% 8 Semidiurnal tides (12,161,37)
am(2,8,2,:,:,:) = ncsd{'amp_w4_t'}(:);
ph(2,8,2,:,:,:) = ncsd{'phase_w4_t'}(:);
am(2,8,3,:,:,:) = ncsd{'amp_w4_u'}(:);
ph(2,8,3,:,:,:) = ncsd{'phase_w4_u'}(:);
am(2,8,4,:,:,:) = ncsd{'amp_w4_v'}(:);
ph(2,8,4,:,:,:) = ncsd{'phase_w4_v'}(:);
am(2,8,5,:,:,:) = ncsd{'amp_w4_w'}(:);
ph(2,8,5,:,:,:) = ncsd{'phase_w4_w'}(:);
am(2,8,1,:,:,:) = ncsd{'amp_w4_d'}(:);
ph(2,8,1,:,:,:) = ncsd{'phase_w4_d'}(:);
am(2,7,2,:,:,:) = ncsd{'amp_w3_t'}(:);
ph(2,7,2,:,:,:) = ncsd{'phase_w3_t'}(:);
am(2,7,3,:,:,:) = ncsd{'amp_w3_u'}(:);
ph(2,7,3,:,:,:) = ncsd{'phase_w3_u'}(:);
am(2,7,4,:,:,:) = ncsd{'amp_w3_v'}(:);
ph(2,7,4,:,:,:) = ncsd{'phase_w3_v'}(:);
am(2,7,5,:,:,:) = ncsd{'amp_w3_w'}(:);
ph(2,7,5,:,:,:) = ncsd{'phase_w3_w'}(:);
am(2,7,1,:,:,:) = ncsd{'amp_w3_d'}(:);
ph(2,7,1,:,:,:) = ncsd{'phase_w3_d'}(:);
am(2,6,2,:,:,:) = ncsd{'amp_w2_t'}(:);
ph(2,6,2,:,:,:) = ncsd{'phase_w2_t'}(:);
am(2,6,3,:,:,:) = ncsd{'amp_w2_u'}(:);
ph(2,6,3,:,:,:) = ncsd{'phase_w2_u'}(:);
am(2,6,4,:,:,:) = ncsd{'amp_w2_v'}(:);
ph(2,6,4,:,:,:) = ncsd{'phase_w2_v'}(:);
am(2,6,5,:,:,:) = ncsd{'amp_w2_w'}(:);
ph(2,6,5,:,:,:) = ncsd{'phase_w2_w'}(:);
am(2,6,1,:,:,:) = ncsd{'amp_w2_d'}(:);
ph(2,6,1,:,:,:) = ncsd{'phase_w2_d'}(:);
am(2,5,2,:,:,:) = ncsd{'amp_w1_t'}(:);
ph(2,5,2,:,:,:) = ncsd{'phase_w1_t'}(:);
am(2,5,3,:,:,:) = ncsd{'amp_w1_u'}(:);
ph(2,5,3,:,:,:) = ncsd{'phase_w1_u'}(:);
am(2,5,4,:,:,:) = ncsd{'amp_w1_v'}(:);
ph(2,5,4,:,:,:) = ncsd{'phase_w1_v'}(:);
am(2,5,5,:,:,:) = ncsd{'amp_w1_w'}(:);
ph(2,5,5,:,:,:) = ncsd{'phase_w1_w'}(:);
am(2,5,1,:,:,:) = ncsd{'amp_w1_d'}(:);
ph(2,5,1,:,:,:) = ncsd{'phase_w1_d'}(:);
am(2,4,2,:,:,:) = ncsd{'amp_s0_t'}(:);
ph(2,4,2,:,:,:) = ncsd{'phase_s0_t'}(:);
am(2,4,3,:,:,:) = ncsd{'amp_s0_u'}(:);
ph(2,4,3,:,:,:) = ncsd{'phase_s0_u'}(:);
am(2,4,4,:,:,:) = ncsd{'amp_s0_v'}(:);
ph(2,4,4,:,:,:) = ncsd{'phase_s0_v'}(:);
am(2,4,5,:,:,:) = ncsd{'amp_s0_w'}(:);
ph(2,4,5,:,:,:) = ncsd{'phase_s0_w'}(:);
am(2,4,1,:,:,:) = ncsd{'amp_s0_d'}(:);
ph(2,4,1,:,:,:) = ncsd{'phase_s0_d'}(:);
am(2,3,2,:,:,:) = ncsd{'amp_e1_t'}(:);
ph(2,3,2,:,:,:) = ncsd{'phase_e1_t'}(:);
am(2,3,3,:,:,:) = ncsd{'amp_e1_u'}(:);
ph(2,3,3,:,:,:) = ncsd{'phase_e1_u'}(:);
am(2,3,4,:,:,:) = ncsd{'amp_e1_v'}(:);
ph(2,3,4,:,:,:) = ncsd{'phase_e1_v'}(:);
am(2,3,5,:,:,:) = ncsd{'amp_e1_w'}(:);
ph(2,3,5,:,:,:) = ncsd{'phase_e1_w'}(:);
am(2,3,1,:,:,:) = ncsd{'amp_e1_d'}(:);
ph(2,3,1,:,:,:) = ncsd{'phase_e1_d'}(:);
am(2,2,2,:,:,:) = ncsd{'amp_e2_t'}(:);
ph(2,2,2,:,:,:) = ncsd{'phase_e2_t'}(:);
am(2,2,3,:,:,:) = ncsd{'amp_e2_u'}(:);
ph(2,2,3,:,:,:) = ncsd{'phase_e2_u'}(:);
am(2,2,4,:,:,:) = ncsd{'amp_e2_v'}(:);
ph(2,2,4,:,:,:) = ncsd{'phase_e2_v'}(:);
am(2,2,5,:,:,:) = ncsd{'amp_e2_w'}(:);
ph(2,2,5,:,:,:) = ncsd{'phase_e2_w'}(:);
am(2,2,1,:,:,:) = ncsd{'amp_e2_d'}(:);
ph(2,2,1,:,:,:) = ncsd{'phase_e2_d'}(:);
am(2,1,2,:,:,:) = ncsd{'amp_e3_t'}(:);
ph(2,1,2,:,:,:) = ncsd{'phase_e3_t'}(:);
am(2,1,3,:,:,:) = ncsd{'amp_e3_u'}(:);
ph(2,1,3,:,:,:) = ncsd{'phase_e3_u'}(:);
am(2,1,4,:,:,:) = ncsd{'amp_e3_v'}(:);
ph(2,1,4,:,:,:) = ncsd{'phase_e3_v'}(:);
am(2,1,5,:,:,:) = ncsd{'amp_e3_w'}(:);
ph(2,1,5,:,:,:) = ncsd{'phase_e3_w'}(:);
am(2,1,1,:,:,:) = ncsd{'amp_e3_d'}(:);
ph(2,1,1,:,:,:) = ncsd{'phase_e3_d'}(:);

% 6 Diurnal Tides
lat = ncd{'lat'}(:);
lev = ncd{'lev'}(:);
month = ncd{'monrh'}(:);
am(1,6,2,:,:,:) = ncd{'amp_w2_t'}(:);
ph(1,6,2,:,:,:) = ncd{'phase_w2_t'}(:);
am(1,6,3,:,:,:) = ncd{'amp_w2_u'}(:);
ph(1,6,3,:,:,:) = ncd{'phase_w2_u'}(:);
am(1,6,4,:,:,:) = ncd{'amp_w2_v'}(:);
ph(1,6,4,:,:,:) = ncd{'phase_w2_v'}(:);
am(1,6,5,:,:,:) = ncd{'amp_w2_w'}(:);
ph(1,6,5,:,:,:) = ncd{'phase_w2_w'}(:);
am(1,6,1,:,:,:) = ncd{'amp_w2_d'}(:);
ph(1,6,1,:,:,:) = ncd{'phase_w2_d'}(:);
am(1,5,2,:,:,:) = ncd{'amp_w1_t'}(:);
ph(1,5,2,:,:,:) = ncd{'phase_w1_t'}(:);
am(1,5,3,:,:,:) = ncd{'amp_w1_u'}(:);
ph(1,5,3,:,:,:) = ncd{'phase_w1_u'}(:);
am(1,5,4,:,:,:) = ncd{'amp_w1_v'}(:);
ph(1,5,4,:,:,:) = ncd{'phase_w1_v'}(:);
am(1,5,5,:,:,:) = ncd{'amp_w1_w'}(:);
ph(1,5,5,:,:,:) = ncd{'phase_w1_w'}(:);
am(1,5,1,:,:,:) = ncd{'amp_w1_d'}(:);
ph(1,5,1,:,:,:) = ncd{'phase_w1_d'}(:);
am(1,4,2,:,:,:) = ncd{'amp_s0_t'}(:);
ph(1,4,2,:,:,:) = ncd{'phase_s0_t'}(:);
am(1,4,3,:,:,:) = ncd{'amp_s0_u'}(:);
ph(1,4,3,:,:,:) = ncd{'phase_s0_u'}(:);
am(1,4,4,:,:,:) = ncd{'amp_s0_v'}(:);
ph(1,4,4,:,:,:) = ncd{'phase_s0_v'}(:);
am(1,4,5,:,:,:) = ncd{'amp_s0_w'}(:);
ph(1,4,5,:,:,:) = ncd{'phase_s0_w'}(:);
am(1,4,1,:,:,:) = ncd{'amp_s0_d'}(:);
ph(1,4,1,:,:,:) = ncd{'phase_s0_d'}(:);
am(1,3,2,:,:,:) = ncd{'amp_e1_t'}(:);
ph(1,3,2,:,:,:) = ncd{'phase_e1_t'}(:);
am(1,3,3,:,:,:) = ncd{'amp_e1_u'}(:);
ph(1,3,3,:,:,:) = ncd{'phase_e1_u'}(:);
am(1,3,4,:,:,:) = ncd{'amp_e1_v'}(:);
ph(1,3,4,:,:,:) = ncd{'phase_e1_v'}(:);
am(1,3,5,:,:,:) = ncd{'amp_e1_w'}(:);
ph(1,3,5,:,:,:) = ncd{'phase_e1_w'}(:);
am(1,3,1,:,:,:) = ncd{'amp_e1_d'}(:);
ph(1,3,1,:,:,:) = ncd{'phase_e1_d'}(:);
am(1,2,2,:,:,:) = ncd{'amp_e2_t'}(:);
ph(1,2,2,:,:,:) = ncd{'phase_e2_t'}(:);
am(1,2,3,:,:,:) = ncd{'amp_e2_u'}(:);
ph(1,2,3,:,:,:) = ncd{'phase_e2_u'}(:);
am(1,2,4,:,:,:) = ncd{'amp_e2_v'}(:);
ph(1,2,4,:,:,:) = ncd{'phase_e2_v'}(:);
am(1,2,5,:,:,:) = ncd{'amp_e2_w'}(:);
ph(1,2,5,:,:,:) = ncd{'phase_e2_w'}(:);
am(1,2,1,:,:,:) = ncd{'amp_e2_d'}(:);
ph(1,2,1,:,:,:) = ncd{'phase_e2_d'}(:);
am(1,1,2,:,:,:) = ncd{'amp_e3_t'}(:);
ph(1,1,2,:,:,:) = ncd{'phase_e3_t'}(:);
am(1,1,3,:,:,:) = ncd{'amp_e3_u'}(:);
ph(1,1,3,:,:,:) = ncd{'phase_e3_u'}(:);
am(1,1,4,:,:,:) = ncd{'amp_e3_v'}(:);
ph(1,1,4,:,:,:) = ncd{'phase_e3_v'}(:);
am(1,1,5,:,:,:) = ncd{'amp_e3_w'}(:);
ph(1,1,5,:,:,:) = ncd{'phase_e3_w'}(:);
am(1,1,1,:,:,:) = ncd{'amp_e3_d'}(:);
ph(1,1,1,:,:,:) = ncd{'phase_e3_d'}(:);

% 04/13:  Print out amps and phases (UT of max for 0glon) for nh=40(97.5km),mo=3 and 12, il=19(0 glat)
OUTA = fopen('2_ctmt_97km_mar_dec_eq_amp_ph','w');
nh = 40;
il = 19;
for mo=3:9:12
 for nn=1:2
  nct = 6;
  if (nn==2)
   nct = 8;
  end 
  for nc=1:nct
  ss = nc - 4;
fprintf (OUTA,'%d %d %d %d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',mo,nn,nc,ss-nn,am(nn,nc,nvu,mo,nh,il),ph(nn,nc,nvu,mo,nh,il),am(nn,nc,nvv,mo,nh,il),ph(nn,nc,nvv,mo,nh,il),am(nn,nc,nvt,mo,nh,il),ph(nn,nc,nvt,mo,nh,il),am(nn,nc,nvd,mo,nh,il),ph(nn,nc,nvd,mo,nh,il));
  end	% nc=1:nct
 end	% nn=1:2  
end	% mo=3:9:12

% Print out for reading in plot_alt_planes.m March 105km (nh=43) Un,Vn for UT=0
mo = 3;	
nh=43;		% nh=43 for 105km
nh=39;		% nh=39 for 95km
ut0 = 0;
for jlon=1:73
 for il=1:37
  ut0_un(il,jlon) = 0.;
  ut0_vn(il,jlon) = 0.;
 for nn=1:2
  nct = 6;
  if (nn==2)
   nct = 8;
  end
  for nc=1:nct
  ss = nc - 4;
  ut0_un(il,jlon) =   ut0_un(il,jlon) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvu,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
  ut0_vn(il,jlon) = ut0_vn(il,jlon) + am(nn,nc,nvv,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvv,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
  end	% nc=1:nct
 end	% nn=1:2  
% Add ZM at 95km from HRDI (day 74) or TIDI (day 80)
 ut0_un_hrdi(il,jlon) = ut0_un(il,jlon) + un95h_mar(il);
 ut0_un_tidi(il,jlon) = ut0_un(il,jlon) + un95t_mar(il);
 ut0_vn_hrdi(il,jlon) = ut0_vn(il,jlon) + vn95h_mar(il);
 ut0_vn_tidi(il,jlon) = ut0_vn(il,jlon) + vn95t_mar(il);
 end	% for il=1:37 
 
fprintf (OUT,'%d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',glon(jlon),ut0_un(:,jlon));
fprintf (OUT,'%d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',glon(jlon),ut0_vn(:,jlon));

end	% for jlon=1:73

% Plot 95km CTMT plus zonal means
figure;
set (gcf,'position',[142 99 647 822])
H=gcf;
clf;
cmax = 60;
titalt = '0UT Mar CTMT 95km';
for nt=1:6
subplot(3,2,nt)
hold on;
if (nt==1)
 contourf(glon,glat,ut0_un);
 tit = sprintf('%s Un',titalt);
end
if (nt==3)
 contourf(glon,glat,ut0_un_tidi);
 tit = sprintf('Un+ZM TIDI');
end
if (nt==5)
 contourf(glon,glat,ut0_un_hrdi);
 tit = sprintf('Un+ZM HRDI');
end
if (nt==2)
 contourf(glon,glat,ut0_vn);
 tit = sprintf('%s Vn',titalt);
end
if (nt==4)
 contourf(glon,glat,ut0_vn_tidi);
 tit = sprintf('Vn+ZM TIDI');
end
if (nt==6)
 contourf(glon,glat,ut0_vn_hrdi);
 tit = sprintf('Vn+ZM HRDI');
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:4
hold off;

% 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) = am(2,4,4,9,nh,il);
  sep85_v_sd_phs_s0(il,nh-34) = ph(2,4,4,9,nh,il);
 end
 for nh=41:161
  alt100_400(nh-40) = alt(nh);
  sep100_v_sd_amp_s0(il,nh-40) = am(2,4,4,9,nh,il);
  sep100_v_sd_phs_s0(il,nh-40) = ph(2,4,4,9,nh,il);
 end
end
OUTAMP = fopen('2_mo390_amp_ph','w');
fprintf(OUTAMP,'den*100 E3 mo=1-10, nh=1,10, jl=1-10 \n');
for jl=1:10
for nh=1:10
fprintf(OUTAMP,'%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',am(1,1,1,:,nh,jl)*100.);
end
end
fprintf(OUTAMP,'v S0 amp 60S-60N 85-105km \n');
for il=7:31
fprintf(OUTAMP,'%6.1f  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',glat(il),sep85_v_sd_amp_s0(il,:));
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 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(1,2,1,mo,nh,il)*100.;
  mo390_phs_de2(mo,il) = ph(1,2,1,mo,nh,il);
  mo390_amp_de3(mo,il) = am(1,1,1,mo,nh,il)*100.;
  mo390_phs_de3(mo,il) = ph(1,1,1,mo,nh,il);
  mo390_amp_s0(mo,il) = am(2,4,1,mo,nh,il)*100.;
  mo390_phs_s0(mo,il) = ph(2,4,1,mo,nh,il);
  mo390_amp_se2(mo,il) = am(2,2,1,mo,nh,il)*100.;
  mo390_phs_se2(mo,il) = ph(2,2,1,mo,nh,il);
 end
end

fprintf(OUTAMP,'DE3 amp,ph \n');
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_amp_de3(:,il));
end
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_phs_de3(:,il));
end
fprintf(OUTAMP,'DE2 amp,ph \n');
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_amp_de2(:,il));
end
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_phs_de2(:,il));
end
fprintf(OUTAMP,'SE2 amp,ph \n');
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_amp_se2(:,il));
end
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_phs_se2(:,il));
end
fprintf(OUTAMP,'S0 amp,ph \n');
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_amp_s0(:,il));
end
for il=1:37
fprintf(OUTAMP,'0:  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',mo390_phs_s0(:,il));
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


mo = 9;				% September = mo #9
montit = 'Sep CTMT';
mo = 3;				% March = mo #3
montit = 'Mar CTMT';
nh = 40;
lt0 = 0;
lt12 = 12;
ut0 = 0;
ut12 = 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)= 0.;
tide100_00ut_d(jlon,il)= 0.;
tide100_00lt_u(jlon,il)= 0.;
tide100_00ut_u(jlon,il)= 0.;
tide100_12lt_d(jlon,il)= 0.;
tide100_12ut_d(jlon,il)= 0.;
tide100_12lt_u(jlon,il)= 0.;
tide100_12ut_u(jlon,il)= 0.;
 for nn=1:2
  nct = 6;
  if (nn==2)
   nct = 8;
  end
  for nc=1:nct
  ss = nc - 4;
tide100_00lt_d(jlon,il)= tide100_00lt_d(jlon,il) + 100.*am(nn,nc,nvd,mo,nh,il)*cos(nn*2.*pi*(lt0-ph(nn,nc,nvd,mo,nh,il))/24. + (ss-nn)*glon(jlon)*2*pi/360.);
tide100_00ut_d(jlon,il)= tide100_00ut_d(jlon,il) + 100.*am(nn,nc,nvd,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvd,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide100_00lt_u(jlon,il)= tide100_00lt_u(jlon,il) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(lt0-ph(nn,nc,nvu,mo,nh,il))/24. + (ss-nn)*glon(jlon)*2*pi/360.);
tide100_00ut_u(jlon,il)= tide100_00ut_u(jlon,il) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvu,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide100_12lt_d(jlon,il)= tide100_12lt_d(jlon,il)+ 100.*am(nn,nc,nvd,mo,nh,il)*cos(nn*2.*pi*(lt12-ph(nn,nc,nvd,mo,nh,il))/24. + (ss-nn)*glon(jlon)*2*pi/360.);
tide100_12ut_d(jlon,il)= tide100_12ut_d(jlon,il)+ 100.*am(nn,nc,nvd,mo,nh,il)*cos(nn*2.*pi*(ut12-ph(nn,nc,nvd,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide100_12lt_u(jlon,il)= tide100_12lt_u(jlon,il) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(lt12-ph(nn,nc,nvu,mo,nh,il))/24. + (ss-nn)*glon(jlon)*2*pi/360.);
tide100_12ut_u(jlon,il)= tide100_12ut_u(jlon,il) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(ut12-ph(nn,nc,nvu,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
  end	% nc=1:nct
 end	% nn=1:2   

 end
end


figure;
clf;
cmax = 60;
%cmax = 25;	% for Tn
for nt=1:4
subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glon,glat,tide100_00lt_u');
 tit = sprintf('%s 0LT Un 97.5km',montit);
end
if (nt==2)
 contourf(glon,glat,tide100_00ut_u');
 tit = sprintf('%s 0UT Un 97.5km',montit');
end
if (nt==3)
 contourf(glon,glat,tide100_12lt_u');
 tit = sprintf('%s 12LT Un 97.5km',montit);
end
if (nt==4)
 contourf(glon,glat,tide100_12ut_u');
 tit = sprintf('%s 12UT Un 97.5km',montit);
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:4
hold off;	


figure;
clf;
cmax = 20;
for nt=1:4
subplot(2,2,nt)
hold on;
if (nt==1)
 contourf(glon,glat,tide100_00lt_d');
 tit = sprintf('%s 0LT pcDen 97.5km',montit);
end
if (nt==2)
 contourf(glon,glat,tide100_00ut_d');
 tit = sprintf('%s 0UT pcDen 97.5km',montit);
end
if (nt==3)
 contourf(glon,glat,tide100_12lt_d');
 tit = sprintf('%s 12LT pcDen 97.5km',montit);
end
if (nt==4)
 contourf(glon,glat,tide100_12ut_d');
 tit = sprintf('%s 12UT pcDen 97.5km',montit);
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:4
hold off;	

% Plot March 15 for 97.5,105,112.5,120km (nh=40,43,46,49), then Sep
for nmo=1:2
mo = 3;	
motit = '0UT Mar CTMT';
if (nmo==2)
 mo = 9;
 motit = '0UT Sep CTMT';
end
ut0 = 0;
mh = 0;
for nh=40:3:49
mh = mh + 1;
for il=1:37
 for jlon=1:73
    tide_00ut_d(jlon,il,mh)= 0.;
    tide_00ut_u(jlon,il,mh)= 0.;
    tide_00ut_v(jlon,il,mh)= 0.;
    tide_00ut_t(jlon,il,mh)= 0.;

 for nn=1:2
  nct = 6;
  if (nn==2)
   nct = 8;
  end
  for nc=1:nct
  ss = nc - 4;
tide_00ut_d(jlon,il,mh) = tide_00ut_d(jlon,il,mh) + am(nn,nc,nvd,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvd,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide_00ut_u(jlon,il,mh) = tide_00ut_u(jlon,il,mh) + am(nn,nc,nvu,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvu,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide_00ut_v(jlon,il,mh) = tide_00ut_v(jlon,il,mh) + am(nn,nc,nvv,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvv,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
tide_00ut_t(jlon,il,mh) = tide_00ut_t(jlon,il,mh) + am(nn,nc,nvt,mo,nh,il)*cos(nn*2.*pi*(ut0-ph(nn,nc,nvt,mo,nh,il))/24. + ss*glon(jlon)*2*pi/360.);
  end	% nc=1:nct
 end	% nn=1:2   

 end
end
end	% for nh=40:3:49

for mh=1:4
 malt = 97.5 + (mh-1)*7.5;

figure;
clf;
for nt=1:4
subplot(2,2,nt)
hold on;
if (nt==1)
 cmax = 0.15;
 contourf(glon,glat,tide_00ut_d(:,:,mh)');
% tit = sprintf('0UT Mar CTMT pcDEN %5.1fkm',malt);
 tit = sprintf('%s pcDEN %5.1fkm',motit,malt);
end
if (nt==2)
 cmax = 35;
 contourf(glon,glat,tide_00ut_t(:,:,mh)');
% tit = sprintf('0UT Mar CTMT Tn %5.1fkm',malt);
 tit = sprintf('%s Tn %5.1fkm',motit,malt);
end
if (nt==3)
 cmax = 60;
 contourf(glon,glat,tide_00ut_u(:,:,mh)');
% tit = sprintf('0UT Mar CTMT Un %5.1fkm',malt);
 tit = sprintf('%s Un %5.1fkm',motit,malt);
end
if (nt==4)
 contourf(glon,glat,tide_00ut_v(:,:,mh)');
% tit = sprintf('0UT Mar CTMT Vn %5.1fkm',malt);
 tit = sprintf('%s Vn %5.1fkm',motit,malt);
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:4
hold off;

end	% for mh=1:4	
end	% for nmo=1:2
