%%
%% 08/14 bae hao/aim3/emery/gswm/plot_mon_ut0.m 
%% Plot 95km UT0 GSWM09 perturbation gph(m) relrho(%) Tn(K) Un,Vn,Wn(m/s)
%  for dres TIEGCM glon glat for diurnal, semi-diurnal, or total values

clear all; close all

calc = 1
if (calc==1)
% pi defined or no?
 pi = 4.*atan(1.)

% Plot amps and phases in GSWM09 72 glons and 59 glats
glon09 = -180:5:175
glat09=-87:3:87
 for nn=1:2
% Read GSWM09 file used in TIEGCM
if (nn==1)
g9fil = sprintf('/hao/aim3/emery/gswm/gswm09_D_95km_pham_mar');
else
g9fil = sprintf('/hao/aim3/emery/gswm/gswm09_S_95km_pham_mar');
end
ING9 = fopen(g9fil,'r');
% Skip over headers and glon and glat
for i=1:16
 rdg9 = fgetl(ING9);
end
%  Read monthly file as made in /hao/aim3/emery/gswm/mkpham_gswm09.f
  for i=1:72
         for j=1:59
    rdg9 = fgetl(ING9);
    tmp = sscanf(rdg9,'%f %e %f %e %f %e %f %e');
            z_amp9(i,j,nn) = tmp(2);
            z_ph9(i,j,nn) = tmp(1);
            t_amp9(i,j,nn) = tmp(4);
            t_ph9(i,j,nn) = tmp(3);
            u_amp9(i,j,nn) = tmp(6);
            u_ph9(i,j,nn) = tmp(5);
            v_amp9(i,j,nn) = tmp(8);
            v_ph9(i,j,nn) = tmp(7);
         end
   end
end


figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;
 cmin=0;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glon09,glat09,z_amp9(:,:,1)',30);
 tit=sprintf('95km real GSWM09 gph Damp (km)');
 cmax=1;
end
if (nf==2)
 contourf(glon09,glat09,z_amp9(:,:,2)',30);
 tit=sprintf('95km rGSWM09 gph Samp (km)');
 cmax=cmax/2.;
end
if (nf==3)
 contourf(glon09,glat09,t_amp9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Tn Damp (K)');
 cmax=30;
end
if (nf==4)
 contourf(glon09,glat09,t_amp9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Tn Samp (K)');
 cmax=cmax/2.;
end
if (nf==5)
 contourf(glon09,glat09,u_amp9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Un Damp (m/s)');
 cmax=80;
end
if (nf==6)
 contourf(glon09,glat09,u_amp9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Un Samp (m/s)');
 cmax=cmax/2.;
end
if (nf==7)
 contourf(glon09,glat09,v_amp9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Vn Damp (m/s)');
 cmax=80;
end
if (nf==8)
 contourf(glon09,glat09,v_amp9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Vn Samp (m/s)');
 cmax=cmax/2.;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;
	
figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;
 cmin=0;
 cmax = 24.;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glon09,glat09,z_ph9(:,:,1)',30);
 tit=sprintf('95km real GSWM09 gph Dph (h)');
end
if (nf==2)
 contourf(glon09,glat09,z_ph9(:,:,2)',30);
 tit=sprintf('95km GSWM09 gph Sph (h)');
end
if (nf==3)
 contourf(glon09,glat09,t_ph9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Tn Dph (h)');
end
if (nf==4)
 contourf(glon09,glat09,t_ph9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Tn Sph (h)');
end
if (nf==5)
 contourf(glon09,glat09,u_ph9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Un Dph (h)');
end
if (nf==6)
 contourf(glon09,glat09,u_ph9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Un Sph (h)');
end
if (nf==7)
 contourf(glon09,glat09,v_ph9(:,:,1)',30);
 tit=sprintf('95km GSWM09 Vn Dph (h)');
end
if (nf==8)
 contourf(glon09,glat09,v_ph9(:,:,2)',30);
 tit=sprintf('95km GSWM09 Vn Sph (h)');
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	


% have 144 glons  and 72 glats
 glondres=-180.:2.5:177.5;
 glatdres=-88.75:2.5:88.75;


 for nn=1:2
% Read GSWM09 file used in TIEGCM
if (nn==1)
tie9fil = sprintf('/hao/aim3/emery/gswm/gswm09_dresD_90-120km_pham_mar');
else
tie9fil = sprintf('/hao/aim3/emery/gswm/gswm09_dresS_90-120km_pham_mar');
end
INTIE9 = fopen(tie9fil,'r');
%  Read monthly file as made in /hao/aim3/emery/gswm/mkpham_gswm09.f
 facamp = 1.;
 phsft = 0.;
 nlat = 72;
  for i=1:144
        for k=1:13	% 90:2.5:120km
         j = 0;
         for js=1:36
           j = j + 1;
    rdgtie9 = fgetl(INTIE9);
    tmp1 = sscanf(rdgtie9,'%f %e %f %e %f %e %f %e');
            z_amp(i,j,k,nn) = tmp1(2)*facamp;
            z_ph(i,j,k,nn) = tmp1(1)-phsft;
            t_amp(i,j,k,nn) = tmp1(4)*facamp;
            t_ph(i,j,k,nn) = tmp1(3)-phsft;
            u_amp(i,j,k,nn) = tmp1(6)*facamp;
            u_ph(i,j,k,nn) = tmp1(5)-phsft;
            v_amp(i,j,k,nn) = tmp1(8)*facamp;
            v_ph(i,j,k,nn) = tmp1(7)-phsft;
           j = j + 1;
    rdgtie9 = fgetl(INTIE9);
    tmp2 = sscanf(rdgtie9,'%f %e %f %e %f %e %f %e');
            z_amp(i,j,k,nn) = tmp2(2)*facamp;
            z_ph(i,j,k,nn) = tmp2(1)-phsft;
            t_amp(i,j,k,nn) = tmp2(4)*facamp;
            t_ph(i,j,k,nn) = tmp2(3)-phsft;
            u_amp(i,j,k,nn) = tmp2(6)*facamp;
            u_ph(i,j,k,nn) = tmp2(5)-phsft;
            v_amp(i,j,k,nn) = tmp2(8)*facamp;
            v_ph(i,j,k,nn) = tmp2(7)-phsft;
          end
         end
% TEMP - write out u_amp,ph for all i at 97.5km (k=4) at equator (nlat/2)
%nn,glondres(i),u_amp(i,nlat/2,4,nn),u_ph(i,nlat/2,4,nn)
       end
 end


figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;
 cmin=0;

iz = 4;
%iz = 3;
iz = 5;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,0.00001*z_amp(:,:,iz,1)',30);
 if (iz==3)
 tit=sprintf('95km GSWM09 gph Damp (km)');
 end
 if (iz==4)
 tit=sprintf('97.5km GSWM09 gph Damp (km)');
 end
 if (iz==5)
 tit=sprintf('100km GSWM09 gph Damp (km)');
 end
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,0.00001*z_amp(:,:,iz,2)',30);
 tit=sprintf('GSWM09 gph Samp (km)');
 cmax=cmax/2.;
end
if (nf==3)
 contourf(glondres,glatdres,t_amp(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Tn Damp (K)');
 cmax=30;
end
if (nf==4)
 contourf(glondres,glatdres,t_amp(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Tn Samp (K)');
 cmax=cmax/2.;
end
if (nf==5)
 contourf(glondres,glatdres,0.01*u_amp(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Un Damp (m/s)');
 cmax=80;
end
if (nf==6)
 contourf(glondres,glatdres,0.01*u_amp(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Un Samp (m/s)');
 cmax=cmax/2.;
end
if (nf==7)
 contourf(glondres,glatdres,0.01*v_amp(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Vn Damp (m/s)');
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,0.01*v_amp(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Vn Samp (m/s)');
 cmax=cmax/2.;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	


figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;
 cmin=0;
 cmax = 24.;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,z_ph(:,:,iz,1)',30);
 if (iz==3)
  tit=sprintf('95km GSWM09 gph Dph (h)');
 end
 if (iz==4)
  tit=sprintf('97.5km GSWM09 gph Dph (h)');
 end
 if (iz==5)
  tit=sprintf('100km GSWM09 gph Dph (h)');
 end
end
if (nf==2)
 contourf(glondres,glatdres,z_ph(:,:,iz,2)',30);
 tit=sprintf('GSWM09 gph Sph (h)');
end
if (nf==3)
 contourf(glondres,glatdres,t_ph(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Tn Dph (h)');
end
if (nf==4)
 contourf(glondres,glatdres,t_ph(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Tn Sph (h)');
end
if (nf==5)
 contourf(glondres,glatdres,u_ph(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Un Dph (h)');
end
if (nf==6)
 contourf(glondres,glatdres,u_ph(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Un Sph (h)');
end
if (nf==7)
 contourf(glondres,glatdres,v_ph(:,:,iz,1)',30);
 tit=sprintf('GSWM09 Vn Dph (h)');
end
if (nf==8)
 contourf(glondres,glatdres,v_ph(:,:,iz,2)',30);
 tit=sprintf('GSWM09 Vn Sph (h)');
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	

% Calculate nalt9 altitudes (90:2.5:120km)
secs = 60.;
      for k=1:13
% tide=amp*cos[n*omega*LT-phase] or amp*cos[n*omega*UT +s*glon -phase
      uth = secs/3600.;
      for il=1:144 
        slt1 = uth + glondres(il)/15.;
        if (slt1<0.) 
 slt1 = slt1 + 24.;
        end
        for jl=1:nlat
           gswm09_z(il,jl,k) = z_amp(il,jl,k,1)*cos(pi*(uth-z_ph(il,jl,k,1))/12.)+ z_amp(il,jl,k,2)*cos(2.*pi*(uth-z_ph(il,jl,k,2))/12.);
           gswm09_t(il,jl,k) = t_amp(il,jl,k,1)*cos(pi*(uth-t_ph(il,jl,k,1))/12.) + t_amp(il,jl,k,2)*cos(2.*pi*(uth-t_ph(il,jl,k,2))/12.);
           gswm09_u(il,jl,k) = u_amp(il,jl,k,1)*cos(pi*(uth-u_ph(il,jl,k,1))/12.)+ u_amp(il,jl,k,2)*cos(2.*pi*(uth-u_ph(il,jl,k,2))/12.);
           gswm09_v(il,jl,k) = v_amp(il,jl,k,1)*cos(pi*(uth-v_ph(il,jl,k,1))/12.) + v_amp(il,jl,k,2)*cos(2.*pi*(uth-v_ph(il,jl,k,2))/12.);
        end
      end
      end



figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,0.00001*gswm09_z(:,:,3)',30);
 tit=sprintf('95km 0001UT GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,gswm09_t(:,:,3)',30);
 tit=sprintf('GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==3)
 contourf(glondres,glatdres,0.01*gswm09_u(:,:,3)',30);
 tit=sprintf('GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==4)
 contourf(glondres,glatdres,0.01*gswm09_v(:,:,3)',30);
 tit=sprintf('GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==5)
 contourf(glondres,glatdres,0.00001*gswm09_z(:,:,4)',30);
 tit=sprintf('97.5km 0001UT GSWM09 gph (km)');
 cmin=95.5;
 cmax=97.5;
 cmin=-1;
 cmax=1;
end
if (nf==6)
 contourf(glondres,glatdres,gswm09_t(:,:,4)',30);
 tit=sprintf('GSWM09 Tn (K)');
 cmin=150;
 cmax=200;
 cmin=-30;
 cmax=30;
end
if (nf==7)
 contourf(glondres,glatdres,0.01*gswm09_u(:,:,4)',30);
 tit=sprintf('GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,0.01*gswm09_v(:,:,4)',30);
 tit=sprintf('97.5km GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	

%stop
end

% Read TIEGCM dres 148 glons -175:2.5:182.5, 72 glats -88.75:2.5:88.75, 
% 4 params at lbc ~97.5km and 0005 UT: gph Tn Un Vn w ZM, then same for GSWM09 perts
%gswmfil = sprintf('/hao/aim3/emery/gswm/gswm09_lbc_0000_5ut');
%gswmfil = sprintf('/hao/aim3/emery/gswm/gswm09_lbc_0000_zk');
%gswmfil = sprintf('/hao/aim3/emery/gswm/gswm09_0001_zk');
%gswmfil = sprintf('/hao/aim3/emery/gswm/gswm09_pert_0001ut');
%gswmfil = sprintf('/hao/aim3/emery/gswm/tie_0001ut_gswm09_mar_lbc_pert');
%if40gswm09=0
gswmfil = sprintf('/hao/aim3/emery/gswm/tie_0001ut_gswm09_mar_alt1');
if40gswm09=1
INGSWM = fopen(gswmfil,'r');
% Read 148 glons  and 72 glats, but skip first 2 glons and don't read last 2 glons
 glondres=-180.:2.5:177.5;
 glatdres=-88.75:2.5:88.75;
if (if40gswm09==1)
% Skip top line of hts, glon and glat lines
 rdgswm = fgetl(INGSWM);	% alt for gswm09_pert_0001ut from gswm09.F etc
% Skip glon and glat lines
  for n=1:23
    rdgswm = fgetl(INGSWM);
  end
% See glitch of more positive gph at 41.25N (#52), esp ~140E 
  for i=1:144
   for j=1:72
    rdgswm = fgetl(INGSWM);
    datgswm = sscanf(rdgswm,'%e');
    gph_09_95(i,j) = datgswm(1);
    tn_09_95(i,j) = datgswm(2);
    un_09_95(i,j) = datgswm(3);
    vn_09_95(i,j) = datgswm(4);
    gph_09_97(i,j) = datgswm(5);
    tn_09_97(i,j) = datgswm(6);
    un_09_97(i,j) = datgswm(7);
    vn_09_97(i,j) = datgswm(8);
   end
  end
figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,gph_09_95',30);
 tit=sprintf('95km 0001UT TIE GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,tn_09_95',30);
 tit=sprintf('TIE GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==3)
 contourf(glondres,glatdres,un_09_95',30);
 tit=sprintf('TIE GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==4)
 contourf(glondres,glatdres,vn_09_95',30);
 tit=sprintf('TIE GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==5)
 contourf(glondres,glatdres,gph_09_97',30);
 tit=sprintf('97.5km noZM TIE GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==6)
 contourf(glondres,glatdres,tn_09_97',30);
 tit=sprintf('TIE GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==7)
 contourf(glondres,glatdres,un_09_97',30);
 tit=sprintf('TIE GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,vn_09_97',30);
 tit=sprintf('v GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;
end	% if (if40gswm09==1)	

% Skip glon and glat lines, then first 2 glons
  for n=1:23
    rdgswm = fgetl(INGSWM);
  end
  for i=1:4	 % 2 fails?	% for tie_0001ut_gswm09_mar_lbc_pert etc from lbc.F
   for j=1:72
    rdgswm = fgetl(INGSWM);
   end
  end
  for i=1:144
   for j=1:72
    rdgswm = fgetl(INGSWM);
    datgswm = sscanf(rdgswm,'%e');
    gphlbc(i,j) = datgswm(1);
%    gphlbc(i,j) = datgswm(1)*0.1;	% 0.0001 instead of 0.00001
    tnlbc(i,j) = datgswm(2);
    unlbc(i,j) = datgswm(3);
    vnlbc(i,j) = datgswm(4);
    gph_09(i,j) = datgswm(5);
%    gph_09(i,j) = datgswm(5)*0.1;	
    tn_09(i,j) = datgswm(6);
    un_09(i,j) = datgswm(7);
    vn_09(i,j) = datgswm(8);
   end
  end

figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,gph_09',30);
 tit=sprintf('0001UT TIE GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,tn_09',30);
 tit=sprintf('TIE GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==3)
 contourf(glondres,glatdres,un_09',30);
 tit=sprintf('TIE GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==4)
 contourf(glondres,glatdres,vn_09',30);
 tit=sprintf('TIE GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==5)
 contourf(glondres,glatdres,gphlbc',30);
 tit=sprintf('TIE GSWM09+96.3723 gph (km)');
 cmin=95.5;
 cmax=97.5;
end
if (nf==6)
 contourf(glondres,glatdres,tnlbc',30);
 tit=sprintf('TIE GSWM09+ZM Tn (K)');
 cmin=150;
 cmax=200;
end
if (nf==7)
 contourf(glondres,glatdres,unlbc',30);
 tit=sprintf('TIE GSWM09+ZM Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,vnlbc',30);
 tit=sprintf('v GSWM09+ZM Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	

stop

% Read TIEGCM dres 144 glons -180:2.5:177.5, 72 glats -88.75:2.5:88.75, 
% D (then SD) 6 params at 95km and 0 UT: gph relrho Tn Un,Vn,Wn
ut0fil = sprintf('/hao/aim3/emery/gswm/gswm09_gph_relrho_tn_un_vn_wn_95km_0ut_mar_ut');
INUT0 = fopen(ut0fil,'r');
% Skip first line
 rdhead = fgetl(INUT0);
% Read 144 glons at 10/line, 72 glats at 10/line
 for n=1:15
  rdloc = fgetl(INUT0);
  datloc = sscanf(rdloc,'%f');
  m0 = (n-1)*10;
  m10 = 10;
  if (n==15)
   m10 = 4;
  end
  for m=1:m10
   glondres(m0+m) = datloc(m);
  end
 end
 for n=1:8
  rdloc = fgetl(INUT0);
  datloc = sscanf(rdloc,'%f');
  m0 = (n-1)*10;
  m10 = 10;
  if (n==8)
   m10 = 2;
  end
  for m=1:m10
   glatdres(m0+m) = datloc(m);
  end
 end
% Skip over params as characters
 rdhead = fgetl(INUT0);
% Read diurnal values then sd
 for nn=1:2
  for i=1:144
   for j=1:72
    rdut0 = fgetl(INUT0);
    datut0 = sscanf(rdut0,'%e');
    gph(i,j,nn) = datut0(1);
    relrho(i,j,nn) = datut0(2);
    tn(i,j,nn) = datut0(3);
    un(i,j,nn) = datut0(4);
    vn(i,j,nn) = datut0(5);
    wn(i,j,nn) = datut0(6);
   end
  end
 end	% for nn=1;2
  for i=1:144
   for j=1:72
      gph(i,j,3) = gph(i,j,1) + gph(i,j,2);
      relrho(i,j,3) = relrho(i,j,1) + relrho(i,j,2);
      tn(i,j,3) = tn(i,j,1) + tn(i,j,2);
      un(i,j,3) = un(i,j,1) + un(i,j,2);
      vn(i,j,3) = vn(i,j,1) + vn(i,j,2);
      wn(i,j,3) = wn(i,j,1) + wn(i,j,2);
   end
  end

% Plot D and SD GSWM perturbations
figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,0.001*gph(:,:,1)',30);
 tit=sprintf('D 95km GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,0.001*gph(:,:,2)',30);
 tit=sprintf('SD 95km GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==3)
 contourf(glondres,glatdres,tn(:,:,1)',30);
 tit=sprintf('D 95km GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==4)
 contourf(glondres,glatdres,tn(:,:,2)',30);
 tit=sprintf('SD 95km GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==5)
 contourf(glondres,glatdres,un(:,:,1)',30);
 tit=sprintf('D 95km GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==6)
 contourf(glondres,glatdres,un(:,:,2)',30);
 tit=sprintf('SD 95km GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==7)
 contourf(glondres,glatdres,vn(:,:,1)',30);
 tit=sprintf('95km GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,vn(:,:,2)',30);
 tit=sprintf('95km GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:8
hold off;	

%stop


% Read Mack Jones ZM file for 12 months, 72 glats -88.75:2.5:88.75, 
% 8 params: Un WINDII-HRDI,HWM07; Vn W-H,HWM07; Tn SABER,MSIS00; geop_ht 20-95km S,M00
mackfil = sprintf('/home/emery/cnofs/mackjr/zm_mon_vs_lat.dat');
INMACK = fopen(mackfil,'r');
% Read 72 glats at 12/line
 for n=1:6
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  m0 = (n-1)*12;
  for m=1:12
   glat72(m0+m) = datmack(m);
  end
 end
 nj = 6*12;
% Skip 12 months in character form
 rdmon = fgetl(INMACK);
 month=1:1:12;
% Now read 72 glats, 12 months per glat, 8 params
 for j=1:nj
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_un_wh(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_un_h07(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_vn_wh(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_vn_h07(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_tn_s(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_tn_m00(:,j) = datmack(:);
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_gph_s(:,j) = datmack(:)*0.001;
  rdmack = fgetl(INMACK);
  datmack = sscanf(rdmack,'%f');
  zm_gph_m00(:,j) = datmack(:)*0.001;
 end

% Plot month vs glat
figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
%set (gcf,'position',[482 32 668 912]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(month,glat72,zm_un_wh',30);
 tit=sprintf('95km WINDII-HRDI Zonal Mean Un (m/s)');
 cmin=-25;
 cmax=55;
end
if (nf==2)
 contourf(month,glat72,zm_un_h07',30);
 tit=sprintf('95km HWM07 Zonal Mean Un (m/s)');
 cmin=-25;
 cmax=55;
end
if (nf==3)
 contourf(month,glat72,zm_vn_wh',30);
 tit=sprintf('95km WINDII-HRDI Zonal Mean Vn (m/s)');
 cmin=-15;
 cmax=25;
end
if (nf==4)
 contourf(month,glat72,zm_vn_h07',30);
 tit=sprintf('95km HWM07 Zonal Mean Vn (m/s)');
 cmin=-15;
 cmax=25;
end
if (nf==5)
 contourf(month,glat72,zm_tn_s',30);
 tit=sprintf('95km SABER Zonal Mean Tn (K)');
 cmin=160;
 cmax=210;
end
if (nf==6)
 contourf(month,glat72,zm_tn_m00',30);
 tit=sprintf('95km MSIS00 Zonal Mean Tn (K)');
 cmin=160;
 cmax=210;
end
if (nf==7)
 contourf(month,glat72,zm_gph_s',30);
 tit=sprintf('95km SABER Zonal GeoP Ht (km)');
 cmin=94;
 cmax=96;
end
if (nf==8)
 contourf(month,glat72,zm_gph_m00',30);
 tit=sprintf('95km MSIS00 Zonal GeoP Ht (km)');
 cmin=87;
 cmax=92;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('month','fontsize',12,'fontweight','bold');	

end	% for nf=1:4
hold off;	

%stop


 f107='120F107';
 daygip = '21Mar'; 

altlb=92.5:2.5:120;


% Plot zonal and UT means (and cf w TIDI)

% Read zm file
 zmfil = sprintf('/home/emery/tidi/outzav_04080_120f107_90to120km');
 INZM = fopen(zmfil,'r');
%  read 13 alts 90:2.5:120km
 for nalt=1:13
  rdzm1 = fgetl(INZM);
  datzm1 = sscanf(rdzm1,'%f %d %d %f %f %f %f %d %s');
% datzm1(1) or altmsis is altlb=95:2.5:105 or 1 more alt than used above
  altmsis(nalt) = datzm1(1);
  ndmsis = datzm1(2);
  ndtidi = datzm1(3);
  utsecmsis = datzm1(4);
  f107msis = datzm1(5);
  yrtidi = datzm1(8);
  rdzm2 = fgetl(INZM);
  datzm2 = sscanf(rdzm2,'%f %e %f %f');
  globtn(nalt) = datzm2(1);
  globden(nalt) = datzm2(2);
  globun(nalt) = datzm2(3);
  globvn(nalt) = datzm2(4);
  rdzmhd = fgetl(INZM);
% Read 72 glat -88.75:1.25:88.75 (same as dres TIEGCM)
  for jl=1:72
   rdzm = fgetl(INZM);
   datzm = sscanf(rdzm,'%f %f %e %f %f');
   glat(jl) = datzm(1);
   t_zm(jl,nalt) = datzm(2);
   d_zm(jl,nalt) = datzm(3);
   u_zm(jl,nalt) = datzm(4);
   v_zm(jl,nalt) = datzm(5);
% Only store 2-5
   if (nalt>1)
   t_zm4(jl,nalt-1) = datzm(2);
   d_zm4(jl,nalt-1) = datzm(3);
   u_zm4(jl,nalt-1) = datzm(4);
   v_zm4(jl,nalt-1) = datzm(5);
   end
  end
 end	% for nalt=1:13
n95 = 3;

% Plot glon (SLT=0-24) vs glat at 0UT

% Add ZM to D+SD
for i=1:144
 for j=1:72
  gphzm(i,j) = 0.001*(gph(i,j,1) + gph(i,j,2)) + zm_gph_s(3,j);
  tnzm(i,j) = tn(i,j,1) + tn(i,j,2) + t_zm(j,n95);
  unzm(i,j) = un(i,j,1) + un(i,j,2) + u_zm(j,n95);
  vnzm(i,j) = vn(i,j,1) + vn(i,j,2) + v_zm(j,n95);
  tnzmm(i,j) = tn(i,j,1) + tn(i,j,2) + zm_tn_s(3,j);
  unzmm(i,j) = un(i,j,1) + un(i,j,2) + zm_un_wh(3,j);
  vnzmm(i,j) = vn(i,j,1) + vn(i,j,2) + zm_vn_wh(3,j);
 end
end

figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,0.001*gph(:,:,3)',30);
 tit=sprintf('95km GSWM09 gph (km)');
 cmin=-1;
 cmax=1;
end
if (nf==2)
 contourf(glondres,glatdres,tn(:,:,3)',30);
 tit=sprintf('95km GSWM09 Tn (K)');
 cmin=-30;
 cmax=30;
end
if (nf==3)
 contourf(glondres,glatdres,un(:,:,3)',30);
 tit=sprintf('95km GSWM09 Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==4)
 contourf(glondres,glatdres,vn(:,:,3)',30);
 tit=sprintf('95km GSWM09 Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==5)
 contourf(glondres,glatdres,gphzm',30);
 tit=sprintf('95km GSWM09+ZMM gph (km)');
 cmin=94;
 cmax=96;
end
if (nf==6)
 contourf(glondres,glatdres,tnzm',30);
 tit=sprintf('95km GSWM09+ZM Tn (K)');
 cmin=150;
 cmax=200;
end
if (nf==7)
 contourf(glondres,glatdres,unzm',30);
 tit=sprintf('95km GSWM09+ZM Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,vnzm',30);
 tit=sprintf('95km GSWM09+ZM Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:4
hold off;	


figure;
% Set size to be 8.5x11 inches
set (gcf,'position',[142 99 647 822]);
H = gcf;
clf;

for nf=1:8
subplot(4,2,nf)
hold on;
if (nf==1)
 contourf(glondres,glatdres,gphzm',30);
 tit=sprintf('95km GSWM09+ZMM gph (km)');
 cmin=94;
 cmax=96;
end
if (nf==2)
 contourf(glondres,glatdres,tnzmm',30);
 tit=sprintf('95km GSWM09+ZMM Tn (K)');
 cmin=150;
 cmax=200;
end
if (nf==3)
 contourf(glondres,glatdres,unzmm',30);
 tit=sprintf('95km GSWM09+ZMM Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==4)
 contourf(glondres,glatdres,vnzmm',30);
 tit=sprintf('95km GSWM09+ZMM Vn (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==5)
 contourf(glondres,glatdres,gphzm',30);
 tit=sprintf('95km GSWM09+ZMM gph (km)');
 cmin=94;
 cmax=96;
end
if (nf==6)
 contourf(glondres,glatdres,tnzm',30);
 tit=sprintf('95km GSWM09+ZM Tn (K)');
 cmin=150;
 cmax=200;
end
if (nf==7)
 contourf(glondres,glatdres,unzm',30);
 tit=sprintf('95km GSWM09+ZM Un (m/s)');
 cmin=-80;
 cmax=80;
end
if (nf==8)
 contourf(glondres,glatdres,vnzm',30);
 tit=sprintf('95km GSWM09+ZM Vn (m/s)');
 cmin=-80;
 cmax=80;
end
h=title([tit]);
shading flat
caxis([cmin cmax]);	
colorbar
set(gca,'layer','top')
set(h,'fontsize',12,'fontweight','bold')
ylabel('latitude','fontsize',12,'fontweight','bold');
xlabel('longitude','fontsize',12,'fontweight','bold');	

end	% for nf=1:4
hold off;	
stop

% Plot zonal means of MSIS00 (at 0UT) and of TIDI u,v

figure;
clf;

for nf=1:4
subplot(2,2,nf)
hold on;
if (nf==1)
 contourf(glat,altlb,log10(d_zm4)',30);
 tit=sprintf('LgDEN gm/cm^3 %s %s MSIS00',daygip,f107);
 cmin=-11;
 cmax=-9;
 cmin=-10;
end
if (nf==2)
 contourf(glat,altlb,t_zm4',30);
 tit=sprintf('Tn K %s %s MSIS00',daygip,f107);
 cmin = 165.;
 cmax = 415;
 cmin = 150;
 cmax = 220;
end
if (nf==3)
 contourf(glat,altlb,u_zm4',30);
 tit=sprintf('Un m/s %s %s TIDI',daygip,f107);
 cmin = -60;
 cmax = 20;
 cmin = -30;
end
if (nf==4)
 contourf(glat,altlb,v_zm4',30);
 tit=sprintf('Vn m/s %s %s TIDI',daygip,f107);
 cmin = -45;
 cmax = 35;
end
h=title([tit]);
shading flat
caxis([cmin 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 nf=1:4
hold off;	

%stop

figure;
clf;

for nf=1:4
subplot(2,2,nf)
hold on;
if (nf==1)
 contourf(glat,altmsis,log10(d_zm)',30);
 tit=sprintf('LgDEN gm/cm^3 %s %s MSIS00',daygip,f107);
 cmin=-11;
 cmax=-9;
 cmin=-10;
end
if (nf==2)
 contourf(glat,altmsis,t_zm',30);
 tit=sprintf('Tn K %s %s MSIS00',daygip,f107);
 cmin = 165.;
 cmax = 415;
 cmin = 150;
 cmax = 220;
end
if (nf==3)
 contourf(glat,altmsis,u_zm',30);
 tit=sprintf('Un m/s %s %s TIDI',daygip,f107);
 cmin = -60;
 cmax = 20;
 cmin = -30;
end
if (nf==4)
 contourf(glat,altmsis,v_zm',30);
 tit=sprintf('Vn m/s %s %s TIDI',daygip,f107);
 cmin = -45;
 cmax = 35;
end
h=title([tit]);
shading flat
caxis([cmin 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 nf=1:4
hold off;	

