clear all;

% Feb 2012  Program to calculate the median of the average 3-deg interpolated values 
%    of 400km density from CHAMP into 3-deg lat bins and 1-h LT bins for various days
%    for the Emery et al., ISEA13 March 2012 talk in Peru

fsize = 12;

%ifall = 1,0 if get all 61 days or hi/lo days
ifall = 1
% Density in /345km/ directory has 345km in 13th place and 400 km in 14th place
n13 = 13;
n13 = 14;
n13 = 12;
alt = 345;
nden = 7.5;
xden = 40;
msis = 16;
if (n13==14)
 alt = 400;
 nden = 0;
 xden = 15;
end
if (n13==12)
 msis = 15;
 xden = 50;
end

if (ifall==1)
% Name 'global' den 400km output file for daily global means (from equal wt ascent/descent)
if (n13>12)
 outglob=sprintf('champ_density/global_400km_LTden_07325_08020d.dat2')
else
% outglob=sprintf('champ_density/global_sat_LTden_07300_08045_a.dat')
 outglob=sprintf('champ_density/global_sat_LTden_07325_08020_b.dat2')
end
OUTG = fopen(outglob,'w');
end	% ifall==1

if (n13>12)
 kmdantit = sprintf('at %dkm 07325-08020',alt);
else
% kmdantit = sprintf('~337-368km 07300-08045');
 kmdantit = sprintf('~337-368km 07325-08020');
end
hitit = sprintf('Kp>=2 Vsw>=500km/s');
lotit = sprintf('Kp<=1 Vsw<=450km/s');

nt1 = 1;
nt2 = 2;
if (ifall==1)
 nt1 = 3;
 nt2 = 3;
end

% CHAMP densities missing for days 08_006-8

lt24 = 0:1:23
lat = -87:3:87

% 25 daynumbers for Kp>=2 and Vsw>=500km/s: 07325-330, 345-6, 351-7, 08005-8, 14-19
ndahi = [325,326,327,328,329,330,345,346,351,352,353,354,355,356,357, 5,6,7,8,14,15,16,17,18,19]
% 25 daynumbers for Kp=<1 and Vsw=<450km/s: 07332-343, 349-350, 358-360, 362-4, 08001-4, 11
ndalo = [332,333,334,335,336,337,338,339,340,341,342,343,349,350,358,359,360,362,363,364, 1,2,3,4,11]


% Zero out LT,glat arrays for accumulations
for i=1:24
 for j=1:59
  for k=1:3
   num(i,j,k) = 0;
   modden(i,j,k) = 0.;
   denmdmod(i,j,k) = 0.;
   den(i,j,k) = 0.;
   denmd(i,j,k) = 0.;
   den2(i,j,k) = 0.;
   denmd2(i,j,k) = 0.;
   ht(i,j,k) = 0.;
  end
 end
end

neqlt = 24;
xeqlt = 0;

for nt=nt1:nt2
nht = 5000.;
xht = 0.;
if (nt<3)
 ndt = 25;
else
 ndt = 61;
%if (n13==12)
%   ndt = 61+50;
% end
end
for nd=1:ndt
 nda = 325+nd-1;
% if (n13==12)
%  nda = 300+nd-1;
% end
 ida = nda;
 iyr = 2007;
 if (nda>365)
  nda = nda - 365;
  iyr = 2008;
 end
 if nt==1
   nda = ndahi(nd);
 end
 if nt==2
   nda = ndalo(nd);
 end
iy1 = 8;
if (nda>250)
 iy1 = 7;
end
% TEMP for missing CHAMP days!
 if (nda>5 & nda<9) 
 else

% Get name for CHAMP file
if (n13==12)
 denfile=sprintf('champ_density/200%d/ascii/Density_3deg_0%d_%3.3d.ascii',iy1,iy1,nda);
else
 denfile=sprintf('champ_density/345km/345km/200%d/ascii/Density_3deg_0%d_%3.3d.ascii',iy1,iy1,nda);
end
% Skip over first 2 lines of denfile and load rest into datden(length,19)
datden = dlmread(denfile,'',2,0);
nl=length(datden);

% nad = 1 for ascent (8.3 to 2.8LT) and 2 for descent (20.4 to 14.8LT) 
 nad = 2;
if (datden(1,4)<datden(2,4))
 nad = 1;
end
if (datden(1,4)==datden(2,4))
  if (datden(2,4)>datden(3,4))
    nad = 1;
  end
end
nad,datden(1,4),datden(2,4),datden(3,4)

if (nt==3)
% Zero out glat arrays for ascent/descent accumulations
for i=1:59
 for k=1:2
 numad(i,k) = 0;
 satalt(i,k) = 0.;
 denad(i,k) = 0.;
 denmdad(i,k) = 0.;
 denadmod(i,k) = 0.;
 denmdadmod(i,k) = 0.;
 denad2(i,k) = 0.;
 denmdad2(i,k) = 0.;
 end
end
end	% nt==3

% Read hourly TIEGCM_SABER_TIDI.nc files starting at 0000 UT every hour
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
% Read entire neutral density values for TIEGCM runs each day (nt=24 for 0-23UT)
[ntm,dayfrac,nz,height,numlat,glat,nlon,glon,dentie,globden,globden5] = get_den_nc(Fname);

for n=1:nl
 uth = datden(n,3)/3600.;
 nuth = fix(datden(n,3)/3600.)+1;
 glat = datden(n,5);
 glon = datden(n,6);
 htcm = datden(n,7)*1.e+5;
% Find TIEGCM den at (glat,glon) htkm for hour 0-23UT (it=nuth+1) at 400km and CHAMP alt
%  TIMED lb TIEGCM dres (72 lats -88.75 to +88.75 and 144 lons -180 to 177.5) 
% Assume that log interp in ht most important and get nearest glat/glon
 ilat = fix( (glat+88.75+2.5)/2.5);
 ilon = fix( (glon+180.0+2.5)/2.5);
 cm400 = 400.e+5;
 if (n13==13)
  cm400 = 345.e+5;
 end
 for nnz=15:nz-1
  if (height(nuth,nnz,ilat,ilon) <= cm400 & height(nuth,nnz+1,ilat,ilon) > cm400)
   nz400 = nnz;
  end
  if (height(nuth,nnz,ilat,ilon) <= htcm & height(nuth,nnz+1,ilat,ilon) > htcm)
   nzsat = nnz;
  end
 end
 if (nuth<ntm)
 nuthp = nuth+1;
 for nnz=15:nz-1
  if (height(nuth+1,nnz,ilat,ilon) <= cm400 & height(nuth+1,nnz+1,ilat,ilon) > cm400)
   nz400p = nnz;
  end
  if (height(nuth+1,nnz,ilat,ilon) <= htcm & height(nuth+1,nnz+1,ilat,ilon) > htcm)
   nzsatp = nnz;
  end
 end
 else
  nuthp = nuth;
  nz400p = nz400;
  nzsatp = nzsat;
 end
   lnden400 = ( (height(nuth,nz400+1,ilat,ilon)-cm400)*log(dentie(nuth,nz400,ilat,ilon)) + (cm400-height(nuth,nz400,ilat,ilon))*log(dentie(nuth,nz400+1,ilat,ilon)) ) / (height(nuth,nz400+1,ilat,ilon)-height(nuth,nz400,ilat,ilon));
 lndensat = ( (height(nuth,nzsat+1,ilat,ilon)-htcm)*log(dentie(nuth,nzsat,ilat,ilon)) + (htcm-height(nuth,nzsat,ilat,ilon))*log(dentie(nuth,nzsat+1,ilat,ilon)) ) / (height(nuth,nzsat+1,ilat,ilon)-height(nuth,nzsat,ilat,ilon));
 if (nuthp>nuth)
  lnden400p = ( (height(nuthp,nz400p+1,ilat,ilon)-cm400)*log(dentie(nuthp,nz400p,ilat,ilon)) + (cm400-height(nuthp,nz400p,ilat,ilon))*log(dentie(nuthp,nz400p+1,ilat,ilon)) ) / (height(nuthp,nz400p+1,ilat,ilon)-height(nuthp,nz400p,ilat,ilon));
  lndensatp = ( (height(nuthp,nzsatp+1,ilat,ilon)-htcm)*log(dentie(nuthp,nzsatp,ilat,ilon)) + (htcm-height(nuthp,nzsatp,ilat,ilon))*log(dentie(nuthp,nzsatp+1,ilat,ilon)) ) / (height(nuthp,nzsatp+1,ilat,ilon)-height(nuthp,nzsatp,ilat,ilon));
  denmodsat = (uth-nuth+1)*exp(lndensatp) + (nuth-uth)*exp(lndensat);
  denmod400 = (uth-nuth+1)*exp(lnden400p) + (nuth-uth)*exp(lnden400);
 else
  denmodsat = exp(lndensat);
  denmod400 = exp(lnden400);
 end
% Convert from g/cm3 to kg/m3 (e-3/e-6 or e+3)
  denmod = denmod400*1.e+3;
  if (n13==12)
   denmod = denmodsat*1.e+3;
  end

 nlat = (datden(n,4)+90)/3;
 nlt = fix(datden(n,8)) + 1;
if (nt==3)
% Accumulate ascent/descent arrays each day
 numad(nlat,nad) = numad(nlat,nad) + 1;
 mad = numad(nlat,nad);
 denadmod(nlat,nad) = denadmod(nlat,nad) + denmod;
 alldenadmod(nlat,nad,mad) = denmod;
 satalt(nlat,nad) = satalt(nlat,nad) + datden(n,7);
 denad(nlat,nad) = denad(nlat,nad) + datden(n,n13);
 alldenad(nlat,nad,mad) = datden(n,n13);
 denad2(nlat,nad) = denad2(nlat,nad) + datden(n,msis);
 alldenad2(nlat,nad,mad) = datden(n,msis);
% Save LT ascent,descent at equator
   if (nlat==30)
     lteq(nad) = datden(n,8);
   end
% Find if next point is change from ascent/descent over the poles
 if (n<nl & datden(n,4)==datden(n+1,4))
   nad = nad + 1;
   if (nad==3)
    nad = 1;
   end
 end
end	% if nt==3 (ifall==1)

% TEMP TO FIND MIN/MAX LT at equator is 2.7547 and 20.3931 LT
 if (nlat==30)
   if (neqlt>datden(n,8))
    neqlt = datden(n,8);
   end
   if (xeqlt<datden(n,8))
    xeqlt = datden(n,8);
   end
 end
% Accumulate in arrays the density values at 400km or at 345km
 num(nlt,nlat,nt) = num(nlt,nlat,nt) + 1;
 modden(nlt,nlat,nt) = modden(nlt,nlat,nt) + denmod;
 den(nlt,nlat,nt) = den(nlt,nlat,nt) + datden(n,n13);
 den2(nlt,nlat,nt) = den2(nlt,nlat,nt) + datden(n,msis);
 ht(nlt,nlat,nt) = ht(nlt,nlat,nt) + datden(n,7);
 m = num(nlt,nlat,nt);
 allden(nlt,nlat,nt,m) = datden(n,n13);
 alldenmod(nlt,nlat,nt,m) = denmod;
 allden2(nlt,nlat,nt,m) = datden(n,msis);
 if (datden(n,7)<nht)
  nht = datden(n,7);
 end
 if (datden(n,7)>xht)
  xht = datden(n,7);
 end
end	% n=1:nl
if (nt==3)
% Calculate and print out 'global' ascent/descent aves where all points are equal in s/c track
 for k=1:2
  for i=1:59
% Assume ALWAYS have all points in ascent/descent - if not - stop!
  if (numad(i,k)==0)
nl,i,k,numad(i,k)
	stop
  end
    clear dennhnlmod
    clear dennhnl
    clear dennhnl2
    for j=1:numad(i,k)
      dennhnlmod(j) = alldenadmod(i,k,j);
      dennhnl(j) = alldenad(i,k,j);
      dennhnl2(j) = alldenad2(i,k,j);
    end
    denmdadmod(i,k) = median(dennhnlmod);
    denadmod(i,k) = denadmod(i,k)/numad(i,k);
    denmdad(i,k) = median(dennhnl);
    denad(i,k) = denad(i,k)/numad(i,k);
    satalt(i,k) = satalt(i,k)/numad(i,k);
    denmdad2(i,k) = median(dennhnl2);
    denad2(i,k) = denad2(i,k)/numad(i,k);
  end	% i=1:59
  gadmod(k) = 0.;
  salt(k) = 0.;
  gad(k) = 0.;
  gad2(k) = 0.;
  for i=1:59
        gadmod(k) = gadmod(k) + denmdadmod(i,k);
        salt(k) = salt(k) + satalt(i,k);
	gad(k) = gad(k) + denmdad(i,k);
	gad2(k) = gad2(k) + denmdad2(i,k);
  end	% i=1/59
   gadmod(k) = gadmod(k)/59.;
   salt(k) = salt(k)/59.;
   gad(k) = gad(k)/59.;
   gad2(k) = gad2(k)/59.;
 end	% k=1:2
nda
   globdenadmod = 0.5*(gadmod(1)+gadmod(2))
   globdenad = 0.5*(gad(1)+gad(2))
   globdenad2 = 0.5*(gad2(1)+gad2(2))
  fprintf(OUTG,'%4.4d %3.3d %5.2f %4.1f %12.4e %12.4e %12.4e %5.2f %4.1f %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e \n',iyr,nda,lteq(1),salt(1),gad(1),gad2(1),gadmod(1),lteq(2),salt(2),gad(2),gad2(2),gadmod(2),globdenad,globdenad2,globdenadmod);
end	% if nt==3 (ifall==1)
end	% TEMP to skip
end	% nd=1:ndt
neqlt,fix(neqlt)+1.5,xeqlt,fix(xeqlt)+1.5
end	% nt=nt1:nt2


% Print out 2 files of den over 25 days 
for nt=nt1:nt2
if (n13>12)
  outfile = sprintf('champ_density/all_den_%dkm_lat_vs_lt_d',alt)
else
%  outfile = sprintf('champ_density/all_den_sat_lat_vs_lt_111d')
  outfile = sprintf('champ_density/all_den_sat_lat_vs_lt_b')
end
 if nt==1
  outfile = sprintf('champ_density/hi_den_%dkm_lat_vs_lt_c',alt)
 end
 if nt==2
  outfile = sprintf('champ_density/lo_den_%dkm_lat_vs_lt_c',alt)
 end
OUT = fopen(outfile,'w');
nx = 0;
max = 0.;
maxmd = 0.;
nx2 = 0;
max2 = 0.;
maxmd2 = 0.;
nxmod = 0;
maxmod = 0.;
maxmdmod = 0.;
for nh=1:24
 for nl=1:59
% Divide by non-zero num
  if (num(nh,nl,nt)>0)
    clear dennhnlmod
    clear dennhnl
    clear dennhnl2
    for i=1:num(nh,nl,nt)
      dennhnlmod(i) = alldenmod(nh,nl,nt,i);
      dennhnl(i) = allden(nh,nl,nt,i);
      dennhnl2(i) = allden2(nh,nl,nt,i);
    end
    denmdmod(nh,nl,nt) = median(dennhnlmod);
    modden(nh,nl,nt) = modden(nh,nl,nt)/num(nh,nl,nt);
    denmd(nh,nl,nt) = median(dennhnl);
    den(nh,nl,nt) = den(nh,nl,nt)/num(nh,nl,nt);
    if (denmd(nh,nl,nt)>maxmd)
     maxmd = denmd(nh,nl,nt);
    end
    if (den(nh,nl,nt)>max)
     max = den(nh,nl,nt);
    end
    if (num(nh,nl,nt)>nx)
     nx = num(nh,nl,nt);
    end
    denmd2(nh,nl,nt) = median(dennhnl2);
    den2(nh,nl,nt) = den2(nh,nl,nt)/num(nh,nl,nt);
    if (denmd2(nh,nl,nt)>maxmd2)
     maxmd2 = denmd2(nh,nl,nt);
    end
    if (den2(nh,nl,nt)>max2)
     max2 = den2(nh,nl,nt);
    end
    if (num(nh,nl,nt)>nx2)
     nx2 = num(nh,nl,nt);
    end
    if (denmdmod(nh,nl,nt)>maxmdmod)
     maxmdmod = denmdmod(nh,nl,nt);
    end
    if (modden(nh,nl,nt)>maxmod)
     maxmod = modden(nh,nl,nt);
    end
    if (num(nh,nl,nt)>nxmod)
     nxmod = num(nh,nl,nt);
    end
  end
  fprintf(OUT,'%5.1f %2.2d %d %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n', ...
   lt24(nh)+0.5,lat(nl),num(nh,nl,nt),den(nh,nl,nt),denmd(nh,nl,nt), ...
   den2(nh,nl,nt),denmd2(nh,nl,nt),modden(nh,nl,nt),denmdmod(nh,nl,nt));
% Set 0's to NaN for plots
  if (num(nh,nl,nt)==0)
   den(nh,nl,nt) = NaN;
   denmd(nh,nl,nt) = NaN;
   den2(nh,nl,nt) = NaN;
   denmd2(nh,nl,nt) = NaN;
   modden(nh,nl,nt) = NaN;
   denmdmod(nh,nl,nt) = NaN;
  else
   den(nh,nl,nt) = den(nh,nl,nt)*1.e+13;
   denmd(nh,nl,nt) = denmd(nh,nl,nt)*1.e+13;
   den2(nh,nl,nt) = den2(nh,nl,nt)*1.e+13;
   denmd2(nh,nl,nt) = denmd2(nh,nl,nt)*1.e+13;
   modden(nh,nl,nt) = modden(nh,nl,nt)*1.e+13;
   denmdmod(nh,nl,nt) = denmdmod(nh,nl,nt)*1.e+13;
  end
 end	% nl=1:36
end	% nh=1:24   
nt,max,maxmd,nx,nht,xht
max2,maxmd2,nx2
end	% nt=1:2

if (nt2>nt1)
tit1num = [hitit,' CHAMP number in bins ',kmdantit];
tit1ave = [hitit,' CHAMP average density*1.e+13 in kg/m3 ',kmdantit];
tit1med = [hitit,' CHAMP median density*1.e+13 in kg/m3 ',kmdantit];
tit2num = [lotit,' CHAMP number in bins ',kmdantit];
tit2ave = [lotit,' CHAMP average density*1.e+13 in kg/m3 ',kmdantit];
tit2med = [lotit,' CHAMP median density*1.e+13 in kg/m3 ',kmdantit];


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

subplot(2,1,1)
            pcolor(lt24,lat,den(:,:,1)')
            title([tit1ave]);
            caxis([nden xden]);
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,den(:,:,2)')
            title([tit2ave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

    hold off;	


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

subplot(2,1,1)
            pcolor(lt24,lat,denmd(:,:,1)')
            title([tit1med]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,denmd(:,:,2)')
            title([tit2med]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

    hold off;	

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

subplot(2,1,1)
            pcolor(lt24,lat,num(:,:,1)')
            title([tit1num]);
            caxis([0 150]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(2,1,2)
            pcolor(lt24,lat,num(:,:,2)')
            title([tit2num]);
            caxis([0 150]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

    hold off;	

else

% Plot All CHAMP densities
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf
titnum = ['All CHAMP number in bins ',kmdantit];
titave = ['All CHAMP average density*1.e+13 in kg/m3 ',kmdantit];
titmed = ['All CHAMP median density*1.e+13 in kg/m3 ',kmdantit];
subplot(3,1,1)
            pcolor(lt24,lat,num(:,:,3)')
            title([titnum]);
            caxis([0 200]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(3,1,2)
            pcolor(lt24,lat,den(:,:,3)')
            title([titave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

subplot(3,1,3)
            pcolor(lt24,lat,denmd(:,:,3)')
            title([titmed]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

    hold off;	

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.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,denmd(:,:,3)')
            title([titmed],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','FontSize',fsize,'Fontweight','bold');
            ylabel('Latitude','FontSize',fsize,'Fontweight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	

% Plot All MSIS00 densities
figure
% position [left bottom width height]
% Set size to be 8.5x11 inches
%set(gcf,'position',[142    99   647   822])
%H = gcf;
clf
titnum = ['All MSIS00 number in bins ',kmdantit];
titave = ['All MSIS00 average density*1.e+13 in kg/m3 ',kmdantit];
titmed = ['All MSIS00 median density*1.e+13 in kg/m3 ',kmdantit];
subplot(3,1,1)
            pcolor(lt24,lat,num(:,:,3)')
            title([titnum]);
            caxis([0 200]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});

    hc=colorbar('horiz');

subplot(3,1,2)
            pcolor(lt24,lat,den2(:,:,3)')
            title([titave]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

subplot(3,1,3)
            pcolor(lt24,lat,denmd2(:,:,3)')
            title([titmed]);
            caxis([nden xden]);
            shading flat
            xlabel('LT');
            ylabel('Latitude');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'});
    hc=colorbar('horiz');

    hold off;	

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.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,denmd2(:,:,3)')
            title([titmed],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','FontSize',fsize,'Fontweight','bold');
            ylabel('Latitude','FontSize',fsize,'Fontweight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	

% Plot model densities
titmed = ['All TIEGCM-WT median density*1.e+13 in kg/m3 ',kmdantit];
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.10,0.10,0.87,0.85])
hold on;
            pcolor(lt24,lat,denmdmod(:,:,3)')
            title([titmed],'FontSize',fsize,'Fontweight','bold');
            caxis([nden xden]);
            shading flat
	    grid;
            xlabel('LT','FontSize',fsize,'Fontweight','bold');
            ylabel('Latitude','FontSize',fsize,'Fontweight','bold');
	    set (gca,'xlim',[0,24],'Xtick',([0:6:24]),'XTickLabel',{'00:00','06:00','12:00','18:00','24:00'},'ylim',[-90,90]);
    hc=colorbar('horiz');

    hold off;	
end	% for nt=3:3 or all 61-3=58 days
