%%
%% This program is written for plotting the TIEGCM dres TEC data at 8 lonsra
%% BA Emery, October 2011 from /aim3/emery/tiegcm_dres/plots/tec8lon.m 
%% Revised for climatology output November 2011
%% Revised March 2012 to add 5 glat sections of global mean gden and gtec plus 140-170km Wi,Ui
clear all; close all

iftec = 1;
ifhmf2 = 0;
% iskpdata = 1 if skip the data and cf plots, and only read/plot model outputs
% iskp8lon = 1 if skip j=1:8 for 8lon plots and only do JRO Viz and global DEN
iskpdata = 1;
iskp8lon = 1;

% Data
 loncutsE=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
if iftec==0
dirname='/hao/aim3/emery/tiegcm_dres/data/cosmic_dec06/'; % directory with COSMIC NmF2,hmF2 data
 loncuts=[25,90,140,175,200,250,285,345]; % specify climatology longitudes to plot
 n15 = 4;
 n1 = 3;
 n2 = 95;
else
dirname='/hao/aim3/emery/tiegcm_dres/data/mit_climatology/'; % directory with TEC data
%dirname='/hao/aim3/emery/cetic_climatology/'; % directory with TEC data
 loncuts=[25,90,140,175,160,110,75,15]; % specify longitudes to plot 4 E and 4 W
 n15 = 1;
 n1 = 1;
 n2 = 96;
end
%latgps=-90:5:90;
latgps=-90:5:85;
min15=0:0.25:23.75;
min60=0.5:1:23.5;
length(min15),length(latgps)
% Dec 13-15, 2006 is daynumbers 347-349
%startday=datenum(2006,12,14,0,0,0); % specify first day to plot
%    [year1,month1,day1,hrs1,min1,sec1]=datevec(startday);

% Model
% longitude slices at begin lons:  25,90,140,175,-160,-110,-75,-15 for 2.5 deg dres
% TIE dres -180(#1) ... -157.5(#10)... -107.5(#30)... -72.5(#44) .. -12.5(#68) ..
%  27.5(#84) .. 92.5(#110) .. 142.5(#130) .. 177.5(#144 last)
jj8 = [84 110 130 144 10 30 44 68];
latdres=-90:2.5:87.5;
 year1 = 2007;
 nday1 = 325;
% nday1 = 366+23;
 nday2 = 366+20;
% nday2 = 366+26;
% nday2 = 369;
%nday1=366
%nday2=367
%month1=1
%day1=1
% TEMP
 nday2 = nday1+1

if iskp8lon~=1
% Loop over number of days
for ida=nday1:nday2-1
 idastr = sprintf('06%d',ida)
	day1 = day1 + ida-nday1;

if iskpdata~=1
 ratmd = NaN(length(min60),length(latgps),8);
 pcmd = NaN(length(min60),length(latgps),8);

% Data
        tecarrmean=NaN(length(min60),length(latgps),8);
        tecarrmed=NaN(length(min60),length(latgps),8);
	hmf2arrmed=NaN(length(min60),length(latgps),8);
	nmf2arrmed=NaN(length(min60),length(latgps),8);

figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
    clf
    for j=1:8
        clear year month day hrs mins lat lonmin lonmax tecmean tecmed dtec stderr npts hmf2med nmf2med stdhmf2 stdnmf2
        loncut=loncuts(j);
      if iftec==1
% TEC
        filename=[dirname,'long_',num2str(loncutsE(j)),...
            '_',num2str(year1),'_',num2str(month1),'_',num2str(day1),'.dat']
%            [year,month,day,hrs,mins,lat,lonmin,lonmax,tecmean,dtec,stderr,npts]=textread(filename,'%d %d %d %f %f %f %f %f %f %f %f %d');
            [year,month,day,hrs,mins,lat,lonmin,lonmax,tecmean,tecmed,dtec,stderr,npts]=textread(filename,'%d %d %d %f %f %f %f %f %f %f %f %f %d');
            %
            dates=datenum(year,month,day,hrs,mins,0);
            datesun=unique(dates);
            datestart=datenum(year(1),month(1),day(1),hrs(1),mins(1),0);
            datesun1=datesun-datestart;
            %
mt=length(tecmean)
            for i=1:length(tecmean)
% Emery:  Toss if dtec>4 or stderr>5 TECU - not much difference for 90E 22-Nov-2007 or 16-Dec-2007
	      if dtec(i)<4 & stderr(i)<5
                idates=find(datesun == dates(i));
                ilat=find(latgps == lat(i));
                tecarrmean(idates,ilat,j)=tecmean(i);
                tecarrmed(idates,ilat,j)=tecmed(i);
              end
            end		% of i=1:length(tecmean)
	else
% COSMIC
        filenamh=[dirname,'hmF2_long_',num2str(loncut),'_',idastr,'.dat'];
        filenamn=[dirname,'NmF2_long_',num2str(loncut),'_',idastr,'.dat'];
            [year,daynum,hrs,mins,secs, lat,lonmin,hmf2med,stdhmf2, npts]=textread(filenamh,'%d %d %d %d %d %f %f %f %f %d');
            [year,daynum,hrs,mins,secs, lat,lonmin,nmf2med,stdnmf2, npts]=textread(filenamn,'%d %d %d %d %d %f %f %f %f %d');
	    i = 0;
	    for n=1:length(min60)
%	     for k=1:length(latgps)-1
	     for k=1:length(latgps)
	      i = i + 1;
	      if (hmf2med(i)>50.)
	       hmf2arrmed(n,k,j) = hmf2med(i);
	      end
	      if (nmf2med(i)>1.)
% Convert from cm-3 to 1.e-12 m-3
	       nmf2arrmed(n,k,j) = nmf2med(i)*1.e-6;
	      end
	     end
	    end
	end
%% Make plots
            subplot(4,2,j)
       if iftec==1
            pcolor(datesun1,latgps,tecarrmed(:,:,j)')
            caxis([0 30]);		% for climatology
%            caxis([0 50]);		% for Dec06 storm
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
            title(['GPS TEC median, ',datestr(datestart,'dd-mmm-yyyy'),', LON=',num2str(loncutsE(j))]);
            set(gca,'TickDir','out');
            datetick('x','HH:MM');
       else
	 if ifhmf2==1
            pcolor(min60,latgps,hmf2arrmed(:,:,j)')
            caxis([0 500]);		% for Dec06 storm
            title(['COSMIC hmF2 (km) ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j))]);
	 else
            pcolor(min60,latgps,nmf2arrmed(:,:,j)')
            caxis([0 2.5]);		% for Dec06 storm
            title(['COSMIC NmF2*1.e-12 (m^3) ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j))]);
         end
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
	end
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);

hold off;

ida
stop
end		% iskpdata~=1

%%  Read TIEGCM .nc files starting at 0000 UT every 15 min and make lon slice plots
%Fname= sprintf('/hao/ftp/pub/emery/dres/dres_s_06%d.nc',ida)
% Read hourly TIEGCM_SABER_TIDI.nc files starting at 0000 UT every hour
%Fname = sprintf('/hao/aim3/emery/tiegcm_saber/tiegcm-linux/dres_saber_s_007325_366_%3.3d.nc',ida-324)
% Read 15min 08027 SSW day catted from 5 .nc files
Fname = sprintf('/hao/aim3/emery/tiegcm_dres/tiegcm-linux/saber_07027.nc')

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

    for j=1:8

     jj = jj8(j);
%  nmax in 1.e-12 m-3
%    [nt,dayfrac,avwi,nmax,hmax,tec,nz,nlat,ht,htne,lon,globden,globtec] = get_teclon_nc(Fname,jj);
   [nt,dayfrac,avwi,avwilo,avuilo,nmax,hmax,tec,nz,nlat,ht,htne,lon,globden,globtec,globden5,globtec5] = get_teclon_nc(Fname,jj);
if (j==1)
% Print out first time (should be 0 UT) in fraction of daynumber
 nt,dayfrac(1)
end
%  06347 0 UT is bad tec, nmax, hmax so delete
if (ida==347 | ida==350)
%hmax(1,1)
  for i=1:nlat
    avwi(1,i) = NaN;
    tec(1,i) = NaN;
    nmax(1,i) = NaN;
    hmax(1,i) = NaN;
  end
end
% Save 285 (j=7) wi at JRO (i=16 at ~12.5S)
if (j==7)
 for n=1:nt
  avwijro(n) = avwi(n,16);
 end
end

% have length 72 for some reason??? so do not print!
%maxnmf2 = max(nmax(:,:));
%minnmf2 = min(nmax(:,:));
%maxhmf2 = max(hmax(:,:));
%minhmf2 = min(hmax(:,:));

if iskpdata~=1
% Calculate ratio of model/data and percent: lat=1:2:nlat; time=0.5*(t+t+1)
     rmodat = 0;
     percentmd = 0;
     ndat = 0;
     for i=1:2:nlat
      i2 = (i+1)/2;
      m = 0;
      for n=n1:n15:n2
       m = m + 1;
       if iftec==1
	dat = tecarrmean(m,i2,j);
        modn = tec(n,i);
	if (n<nt)
	 modnp = tec(n+1,i);
	end
       else
	if ifhmf2==1
         dat = hmf2arrmed(m,i2,j);
	 modn = hmax(n,i);
	 modnp = hmax(n+1,i);
	else 
         dat = nmf2arrmed(m,i2,j);
	 modn = nmax(n,i);
	 modnp = nmax(n+1,i);
        end
       end       
        if (dat>=0.001)
	 ndat = ndat + 1;
         np1 = 0;
         if (ida==347 & n==1)
	   np1 = 1;
	 end
         if (ida==350 & n==1)
	   np1 = 1;
	 end
         if (np1==1)
	   tie = modnp;
	 else
	   if (n==nt)
	    tie = modn;
           else
	    tie = 0.5*(modn+modnp);
	   end
	 end
	ratmd(m,i2,j) = tie/dat;
        pcmd(m,i2,j) =100.*(tie/dat-1.);
        if (tie/dat < 1.)
          pcmd(m,i2,j) = -100.*(dat/tie-1.);
        end
	 rmodat = rmodat + tie/dat;
         percentmd = percentmd + pcmd(m,i2,j);
%if (tie/dat > 2. | tie/dat < 0.5)
%tie/dat,pcmd(m,i2,j)
%end
	end
      end
     end
	rmodat = rmodat/ndat
	percentmd = percentmd/ndat
rmd8(j) = rmodat;
pcmd8(j) = percentmd;
% 53-67 COSMIC npts vs ~1150-2650 TEC GPS npts
ndat
end		% iskpdata~=1

%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(dayfrac,latdres,tec')
            caxis([0 30]);		% for climatology
%            caxis([0 50]);		% for Dec06 storm
            title(['TIEGCM dres SABER/TIDI TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj))]);
%            title(['TIEGCM dres TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	else
	  if ifhmf2==1
            pcolor(dayfrac,latdres,hmax')
            caxis([0 500]);		% for Dec06 storm
            title(['TIEGCM dres SABER/TIDI hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj))]);
%            title(['TIEGCM dres hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	  else
            pcolor(dayfrac,latdres,nmax')
            caxis([0 2.5]);		% for Dec06 storm
            title(['TIEGCM dres SABER/TIDI NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj))]);
%            title(['TIEGCM dres NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(lon(jj)),' M/D=',num2str(rmodat)]);
	  end
	end
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[ida,ida+1],'Xtick',([ida:0.25:ida+1]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
if iftec==1
    plotname=['TIEdres_',idastr,'_TEC'];
else
  if ifhmf2==1
    plotname=['TIEdres_',idastr,'_hmf2'];
  else
    plotname=['TIEdres_',idastr,'_nmf2'];
  end
end
% This works, but matlab dies immediately so skip
%    print('-dpng',plotname);
hold off;

if (iskpdata~=1)
% Plot ratio of model/data
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
for j=1:8
%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(min15,latgps,ratmd(:,:,j)')		
            title(['model/data TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncutsE(j)),' M/D=',num2str(rmd8(j))]);
	else
            pcolor(min60,latgps,ratmd(:,:,j)')		
	  if ifhmf2==1
            title(['model/data hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' M/D=',num2str(rmd8(j))]);
	  else	
            title(['model/data NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' M/D=',num2str(rmd8(j))]);
	  end
	end
            caxis([0.5 2]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
hold off;

% Plot percent of model/data
figure;
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
for j=1:8
%% Make plots
            subplot(4,2,j)
	if iftec==1
            pcolor(min15,latgps,pcmd(:,:,j)')		
            title(['%model/data TEC ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncutsE(j)),' %M/D=',num2str(pcmd8(j))]);
	else
            pcolor(min60,latgps,pcmd(:,:,j)')		
	  if ifhmf2==1
            title(['%model/data hmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' %M/D=',num2str(pcmd8(j))]);
	  else	
            title(['%model/data NmF2 ',num2str(year1),' ',num2str(ida),' LON=',num2str(loncuts(j)),' %M/D=',num2str(pcmd8(j))]);
	  end
	end
            caxis([-100 100]);
            shading flat
            %xlabel('Date');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    end				% end of plotting 8 longitudes
    hc=colorbar('horiz');
    hc_pos=get(hc,'Position');
    hc_pos(1)=0.35;
    hc_pos(2)=0.065;
    set(hc,'Position', hc_pos);
    hold off;
end		% iskpdata~=1

% Plot avwijro
fsize = 12;
zero = [0.,0.];
xtot = [dayfrac(1),dayfrac(nt)];
 y1=-75;
 y2=85;
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
subplot('position',[0.08,0.050,0.87,0.850])	% 1/1
hold on;
title('TIEGCM dres SABER/TIDI 09027 SSW Jicamarca (11.9S, 284E)','FontSize',fsize);
plot (dayfrac,avwijro,'k');
plot(xtot,zero,'--k');
ylabel('Wi (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[dayfrac(1),dayfrac(nt)],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

end	% for ida=nday1:nday2-1

else 	% for iskp8lon==1 do following multi-day plots


% Read 15min 08027 SSW day catted from 5 .nc files
%Fname = sprintf('/hao/aim3/emery/tiegcm_dres/tiegcm-linux/saber_07027.nc')

% Open file with Scherliess quiet-time drift model for 285 glon 70 f10.7 flux
%   run from vdrift-model_jgr0499.f and located in file vdrift_07355_70fx_.. files
filwiq = sprintf('/home/emery/cetic/vdrift_07355_70fx_285')
datwiq = dlmread(filwiq,'',1,0);

% Open Emmert's MSIS and satellite drag files
f250km = '/home/emery/cetic/tledens_2007_2010_250km.txt'
f400km = '/home/emery/cetic/tledens_2007_2010_400km.txt'
f550km = '/home/emery/cetic/tledens_2007_2010_550km.txt'
%  Skip first line description for each file
%  Format: year daynum ut(0 for middday) ratio(MSIS/drag) uncert lg10MSIS lg10drag uncert
dat250km = dlmread(f250km,'',1,0);
dat400km = dlmread(f400km,'',1,0);
dat550km = dlmread(f550km,'',1,0);
% Open output files for daily global densities and medians of JRO Viz 
OUTDEN = fopen('daily_den_tec_test','w');
fprintf(OUTDEN,'daynum(2007) 250km_drag,msis,model 400km_drag,msis,model 550km_drag,msis,model tec tecSP,SH,eq,NH,NP \n');
OUTCLIM = fopen('climatology_viz_jro_test','w');
fprintf(OUTCLIM,'%d %d %d nday1,2 ndays \n UT JROWi(m/s)model hi,lo Uilo: median mean Schierliess\n',nday1,nday2-1,nday2-nday1);

% Make 8 lon files of TIEGCM SABER/TIDI output
% longitude slices at begin lons:  25,90,140,175,-160,-110,-75,-15 for 2.5 deg dres 
OUT025 = fopen('ccmc_output/lon_025_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT090 = fopen('ccmc_output/lon_090_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT140 = fopen('ccmc_output/lon_140_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT175 = fopen('ccmc_output/lon_175_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT200 = fopen('ccmc_output/lon_200_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT250 = fopen('ccmc_output/lon_250_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT285 = fopen('ccmc_output/lon_285_2_TIE-GCM_2007.325-2008.026e.txt','w');
OUT345 = fopen('ccmc_output/lon_345_2_TIE-GCM_2007.325-2008.026e.txt','w');


m = 0;
nd = 0;
nan3 = 0;
for ida=nday1:nday2-1

% Read TIEGCM SABER/TIDI climatology runs
if (ida<366)
Fname = sprintf('/hao/aim3/emery/tiegcm_saber/tiegcm-linux/dres_saber_s_007325_366_%3.3d.nc',ida-324)
else
Fname = sprintf('/hao/aim3/emery/tiegcm_saber/tiegcm-linux/dres_saber_s_08001_027_%3.3d.nc',ida-365)
end

yearwr = year1;
jda = ida;
if (ida > 365)
 yearwr = year1 + 1;
 jda = jda - 365;
end

 for j=1:8
     jj = jj8(j);
%  nmax in 1.e-12 m-3
%    [nt,dayfrac,avwi,nmax,hmax,tec,nz,nlat,ht,htne,lon,globden,globtec] = get_teclon_nc(Fname,jj);
[nt,dayfrac,avwi,avwilo,avuilo,nmax,hmax,tec,nz,nlat,ht,htne,lon,globden,globtec,globden5,globtec5] = get_teclon_nc(Fname,jj);

for n=1:24
for ii=1:nlat
ilat = -90 + (ii-1)*2.5;
% nmax in 1.e6 cm-3!
if j==1
fprintf(OUT025,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==2
fprintf(OUT090,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==3
fprintf(OUT140,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==4
fprintf(OUT175,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==5
fprintf(OUT200,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==6
fprintf(OUT250,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==7
fprintf(OUT285,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
if j==8
fprintf(OUT345,'%d %d %d 0 0 %5.1f %d %5.2f %e %7.2f %7.2f %7.2f %7.2f \n',yearwr,jda,n-1,ilat,loncutsE(j),tec(n,ii),nmax(n,ii),hmax(n,ii),avwi(n,ii),avwilo(n,ii),avuilo(n,ii));
end
end	% ii=1:nlat
end	% n=1:24

if j==7
% Print out first time (should be 0 UT) in fraction of daynumber
 dayfrac(1)

% Save 250,400,550 km global density and global daily TEC
  nd = nd + 1;
  midday(nd) = dayfrac(nt/2);
  if (ida > 365)
   midday(nd) =  dayfrac(nt/2) + 365.;
  end
  gtec(nd) = globtec;
  den3(nd,1) = globden(1);
  den3(nd,2) = globden(2);
  for k5=1:5
   gtec5(nd,k5) = globtec5(k5);
   den35(nd,1,k5) = globden5(1,k5);
   den35(nd,2,k5) = globden5(2,k5);
  end
% Put NaN in for missing model neutral densities at 550 km
  if (globden(3) < -98.)
   den3(nd,3) = NaN;
   for k5=1:5
    den35(nd,3,k5) = NaN;
   end
   nan3 = nan3 + 1;
   den3nonan = -99.;
  else
   den3(nd,3) = globden(3);
   den3nonan = globden(3);
   for k5=1:5
    den35(nd,3,k5) = globden5(3,k5);
   end
  end
% Save Emmert MSIS and drag value
dat250km(nday1+nd-1,2)
%  drag250(nd) = dat250km(nday1+nd-1,7);
  drag250(nd) = dat250km(nday1+nd,7);
  msis250(nd) = dat250km(nday1+nd-1,6);
%  drag400(nd) = dat400km(nday1+nd-1,7);
  drag400(nd) = dat400km(nday1+nd,7);
  msis400(nd) = dat400km(nday1+nd-1,6);
%  drag550(nd) = dat550km(nday1+nd-1,7);
  drag550(nd) = dat550km(nday1+nd,7);
  msis550(nd) = dat550km(nday1+nd-1,6);
%  print out Emmert daily values and model for later plots
fprintf(OUTDEN,'%5.1f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n',ida+0.5,drag250(nd),msis250(nd),den3(nd,1),drag400(nd),msis400(nd),den3(nd,2),drag550(nd),msis550(nd),den3nonan,gtec(nd),gtec5(nd,1),gtec5(nd,2),gtec5(nd,3),gtec5(nd,4),gtec5(nd,5),den35(nd,2,1),den35(nd,2,2),den35(nd,2,3),den35(nd,2,4),den35(nd,2,5));

% Save 285 (j=7) wi at JRO (i=16 at ~12.5S)
 for n=1:nt
  m = m + 1;
  dayfractot(m) = dayfrac(n);
  if (ida > 365)
   dayfractot(m) = dayfrac(n) + 365.;
  end
  avwijro(m) = avwi(n,16);
  utwi(nd,n) = avwi(n,16);
  avwijrolo(m) = avwilo(n,16);
  avuifrolo(m) = avuilo(n,16);
  utwilo(nd,n) = avwilo(n,16);
  utuilo(nd,n) = avuilo(n,16);
 end
%m,avwijro(m),dayfractot(m)
dayfractot(m)
end	%  for j=7
 end	%  for j=1:8

end	%  for ida=nday1:nday2-1

% Calculate median and mean values of density
med250d = median(drag250);
mean250d = mean(drag250);
med400d = median(drag400);
mean400d = mean(drag400);
med550d = median(drag550);
mean550d = mean(drag550);
med250m = median(msis250);
mean250m = mean(msis250);
med400m = median(msis400);
mean400m = mean(msis400);
med550m = median(msis550);
mean550m = mean(msis550);
med250t = median(den3(:,1));
mean250t = mean(den3(:,1));
med400t = median(den3(:,2));
mean400t = mean(den3(:,2));
if (nan3==nd)
med550t = -99;
mean550t = -99;
else
med550t = median(den3(:,3));
mean550t = mean(den3(:,3));
end
medgtec = median(gtec);
meangtec = mean(gtec);

%  print out median and mean values for lg10(den) at 3 alts
fprintf(OUTDEN,'%5.1f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %7.2f\n',0.5*(nday2-nday1)+nday1,med250d,med250m,med250t,med400d,med400m,med400t,med550d,med550m,med550t,medgtec);
fprintf(OUTDEN,'%5.1f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %7.2f\n',0.5*(nday2-nday1)+nday1,mean250d,mean250m,mean250t,mean400d,mean400m,mean400t,mean550d,mean550m,mean550t,meangtec);

% Calculate UT climatology of Viz at Jicamarca
for n=1:24
 nlt = n-5;
 if nlt<1
   nlt = nlt + 24;
 end
 ut24(n) = n-0.5;
 medwijro(n) = median(utwi(:,n));
 meanwijro(n) = mean(utwi(:,n));
 medwijrolo(n) = median(utwilo(:,n));
 meanwijrolo(n) = mean(utwilo(:,n));
 meduijrolo(n) = median(utuilo(:,n));
 meanuijrolo(n) = mean(utuilo(:,n));
 quietmod(n) = 0.5*(datwiq(nlt,2)+datwiq(nlt+1,2));
fprintf(OUTCLIM,' %6.2f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f\n',ut24(n),medwijro(n),meanwijro(n),medwijrolo(n),meanwijrolo(n),meduijrolo(n),meanuijrolo(n));
end



% Make plots
% Plot avwijro
fsize = 12;
zero = [0.,0.];
xtot = [dayfractot(1),dayfractot(m)];
 y1=-75;
 y2=85;
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
subplot('position',[0.08,0.050,0.87,0.850])	% 1/1
hold on;
title('TIEGCM dres SABER/TIDI Jicamarca (11.9S, 284E)','FontSize',fsize);
plot (dayfractot,avwijro,'k');
plot(xtot,zero,'--k');
ylabel('Wi (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[nday1,nday2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

% plot median and average Wi JRO - climatology
xt24 = [0,24];
xnoon = [17,17];
ytot = [-20,20];
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
subplot('position',[0.08,0.050,0.87,0.850])	% 1/1
hold on;
title('TIEGCM dres SABER/TIDI Jicamarca (11.9S, 284E)','FontSize',fsize);
plot (ut24,medwijro,'k');
plot (ut24,meanwijro,'--k');
plot (ut24,quietmod,'--r');
plot(xt24,zero,'--k');
plot(xnoon,ytot,':k');
text(16.5,19.,'noon','Color','k','FontSize',fsize,'FontWeight','bold');
text(9.5,19.,'quiet-time model','Color','r','FontSize',fsize,'FontWeight','bold');
text(0.5,19.,'__median --mean TIEGCM','Color','k','FontSize',fsize,'FontWeight','bold');
ylabel('Wi (m/s)','FontSize',fsize);
xlabel('UT(h)','FontSize',fsize);
set (gca,'xlim',[0,24],'ylim',[-20,20],'xtick',([0:2:24]),'FontSize',fsize,'box','on');
hold off;

% Plot global densities
% 09027 msis=[-10.49505 -12.15793 -13.44664] drag=[-10.62822 -12.38078 -13.71003]
% tiesaber g/cm3 = [3.88e-14 7.85e-16 NA], lg10(kg/m3) = [-10.4106 -12.1053 NA]  (like MSIS)
dark_green = [0 0.5 0];
den3
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;
clf
subplot('position',[0.08,0.500,0.87,0.450])	% 1/2
hold on;
title('TIEGCM dres SABER/TIDI lg10(DEN kg/m3)','FontSize',fsize);
plot (midday,den3(:,1),'s-r');
plot (midday,den3(:,2),'x-k');
%plot (midday,den3(:,2),'Color',dark_green,'x-');
plot (midday,den3(:,3),'d-c');
plot (midday,drag250,':r');
plot (midday,drag400,':k');
plot (midday,drag550,':c');
plot (midday,msis250,'--r');
plot (midday,msis400,'--k');
plot (midday,msis550,'--c');
text(nday1+0.01*(nday2-nday1),-10.1,'...drag ---MSIS ____TIEGCM  altitude:','FontSize',fsize,'FontWeight','bold');
text(nday1+0.5*(nday2-nday1),-10.1,'250','Color','r','FontSize',fsize,'FontWeight','bold');
text(nday1+0.6*(nday2-nday1),-10.1,'400','Color','k','FontSize',fsize,'FontWeight','bold');
text(nday1+0.7*(nday2-nday1),-10.1,'550','Color','c','FontSize',fsize,'FontWeight','bold');
ylabel('at 250,400,550km','FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[nday1,nday2],'ylim',[-14.5,-10],'FontSize',fsize,'box','on');

subplot('position',[0.08,0.050,0.87,0.400])	% 2/2
hold on;
plot (midday,gtec,'x-r');
ylabel('Global TEC (TECU)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[nday1,nday2],'FontSize',fsize,'box','on');
hold off;



end	% for iskp8lon==1




