clear all;
close all;

% Nov 2012  Program to calculate the ave of older legacy CHAMP winds from 
%  Kathrin Haeusler (using calc of spherical earth so alt ~5km lower)
%    into 3-deg lat bins and 1-h LT bins for various days
%    for the Emery et al., CCMC talk and Poster (w values +/-13.1 deg for C/NOFS)
%  Revise to separate days of various Ap and F107 levels
% Kp1=0-1Ap=0-4(7), Kp2=2-4-Ap=7-22(27); F1=<75,F2=75-95,F3=95-125,F4=125-175,F5>175
%  10 cases:  Kp1F1, Kp2F1, Kp1F2, Kp2F2, Kp1F3, Kp2F3, Kp1F4, Kp2F4, Kp1F5, Kp2F5
%  4 seasons of 131 days:  NHwin=290-055; NHspr=15-145; NHsum=107-237; NHfal=199-329
% Max F107 is 285.5


%nc1 = 1;
%nc2 = 4;
%iy1 = 4;
%iy2 = 8;

%nc1 = 5;
%nc2 = 10
%iy1 = 1
%iy2 = 4

nc1 = 1
nc2 = 10
iy1 = 1
iy2 = 8

%  outbad = sprintf('champ_winds_bad')
%  OUTB = fopen(outbad,'w');
%  fprintf (OUTB,'SLT 9.5 and 21.5LT bad Un (>+400 and <-400m/s) \n iyr nda slat un ht den \n');

% Midpoints 24h and 60 lats
lt24 = 0.5:1:23.5;
lat = -88.5:3:88.5;
lon=[0:15:355];

titseas1='spr';
titseas2='sum';
titseas3='fal';
titseas4='win';
%titcond1='Kp_0-1_F107_65-75';
%titcond2='Kp_2-4_F107_65-75';
%titcond3='Kp_0-1_F107_75-95';
%titcond4='Kp_2-4_F107_75-95';
%titcond5='Kp_0-1_F107_95-125';
%titcond6='Kp_2-4_F107_95-125';
%titcond7='Kp_0-1_F107_125-175';
%titcond8='Kp_2-4_F107_125-175';
%titcond9='Kp_0-1_F107_175-max';
%titcond10='Kp_2-4_F107_175-max';
titcond1='Kp_0-1_F107_65-85';
titcond2='Kp_2-4_F107_65-85';
titcond3='Kp_0-1_F107_85-105';
titcond4='Kp_2-4_F107_85-105';
titcond5='Kp_0-1_F107_105-145';
titcond6='Kp_2-4_F107_105-145';
titcond7='Kp_0-1_F107_145-195';
titcond8='Kp_2-4_F107_145-195';
titcond9='Kp_0-1_F107_195-285';
titcond10='Kp_2-4_F107_195-285';



% MJD = model Julian Day: First=0517=1 June 2001 (01152) and Last=3285=08364
% MJD 001: 2001=366 2002=731 2003=1096, 2004=1461, 2005=1827, 2006=2192, 2007=2557, 2008=2922
%                 sum01(136-179), fal01(136-236), win01(200-228)
% spr02(119-139), sum02(153-191), fal02(165-188), win02(121-165)
% spr03(118-139), sum03(119-133), fal03(113-150), win03(104-150)
% spr04(102-111), sum04(100-122), fal04(104-122), win04(92-105)
% spr05(87-102),  sum05(87-102),  fal05(75-88),   win05(75-88)
% spr06(75-90),   sum06(78-90),   fal06(74-84),   win06(74-84)
% spr07(71-81),   sum07(68-71),   fal07(67-74),   win07(67-72)
% spr08(69-72),   sum08(68-71),   fal08(67-68),   win08(67-68)
% spr09(68-72),   sum09(69-71),   fal09(69-72),   win09(72-83)
% spr10(75-83)

mjd01=[366 731 1096 1461 1827 2192 2557 2922 3287];
apskp1=[10593 10958 11323 11688 12054 12419 12784 13149];
longd1=[15 107 199 290];	% for ndt=131 days for full LT orbit
shortd1=[36 127 218 309 401];	% for ndt=91 days for shorter seasons (win=35+57,8=92,3d)

miss01=[0528 0529 0536 0544 0547 0558 0560 0593 0605 0606 0607 0608 0619 0620 0623 0651 0656 0657 0659 0668 0672 0673 0693 0696 0697 0699 0700 0701 0702 0712 0717 0720];
miss02=[0745 0816 0830 0889 0890 0891 0892 0924 1073 1074 1092 1093 1094 1095];
miss03=[1123 1148 1276 1280 1281 1444 1445 1446];
miss04=[1478 1491 1492 1493 1494 1495 1506 1540 1559 1626 1661 1680 1701 1707 1708 1778 1779 1780 1783 1784 1785 1786 1796 1810 1817 1826];
miss05=[1853 1854 1855 1856 1857 1858 1871 1915 1917 1932 1933 1954 1972 1986 2058 2096 2105 2106 2107 2111 2112 2113 2147 2148 2149 2154 2155 2156 2165 2191];
miss06=[2209 2217 2240 2253 2254 2258 2276 2277 2278 2321 2322 2329 2330 2374 2375 2376 2377 2378 2399 2403 2407 2408 2409 2444 2467 2468 2469 2475 2478 2479 2480 2484 2517 2536 2542 2543 2544 2545 2546 2552 2553];
miss07=[2561 2579 2585 2587 2588 2597 2599 2600 2601 2617 2645 2647 2651 2659 2756 2766 2767 2770 2782 2860 2874 2912];
miss08=[2927 2928 2929 2979 2980 2981 2982 3091 3092 3199 3200 3204];
miss_sum01 = [0528 0529 0536 0544 0547 0558 0560 0593 0605 0606 0607 0608 0619 0620 0623 0651 0656 0657 0659 0668 0672 0673 0693 0696 0697 0699 0700 0701 0702 0712 ...
0717  0720 0745 0816 0830 0889 0890 0891 0892 0924 1073 1074 1092 1093 1094 1095 1123 1148 1276 1280 1281 1444 1445 1446 1478 1491 1492 1493 1494 1495 1506 1540 1559 ...
1626 1661 1680 1701 1707 1708 1778 1779 1780 1783 1784 1785 1786 1796 1810 1817 1826 1853 1854 1855 1856 1857 1858 9999];	% 87 total
miss_spr04 = [1506 1540 1559 1626 1661 1680 1701 1707 1708 1778 1779 1780 1783 1784 1785 1786 1796 1810 1817 1826 1853 1854 1855 1856 1857 1858 1871 1915 1917 1932 ...
 1933 1954 1972 1986 2058 2096 2105 2106 2107 2111 2112 2113 2147 2148 2149 2154 2155 2156 2165 2191 2209 2217 2240 2253 2254 2258 2276 2277 2278 2321 2322 2329 ...
 2330 2374 2375 2376 2377 2378 2399 2403 2407 2408 2409 2444 2467 2468 2469 2475 2478 2479 2480 2484 2517 2536 2542 2543 2544 2545 2546 2552 2553 2561 2579 2585 ...
 2587 2588 2597 2599 2600 2601 2617 2645 2647 2651 2659 2756 2766 2767 2770 2782 2860 2874 2912 2927 2928 2929 2979 2980 2981 2982 3091 3092 3199 3200 3204 9999];	% 126 total
miss_all = [0528 0529 0536 0544 0547 0558 0560 0593 0605 0606 0607 0608 0619 0620 0623 0651 0656 0657 0659 0668 0672 0673 0693 0696 0697 0699 0700 0701 0702 0712 ...
0717  0720 0745 0816 0830 0889 0890 0891 0892 0924 1073 1074 1092 1093 1094 1095 1123 1148 1276 1280 1281 1444 1445 1446 1478 1491 1492 1493 1494 1495 1506 1540 1559 ...
1626 1661 1680 1701 1707 1708 1778 1779 1780 1783 1784 1785 1786 1796 1810 1817 1826 1853 1854 1855 1856 1857 1858 1871 1915 1917 1932 ...
 1933 1954 1972 1986 2058 2096 2105 2106 2107 2111 2112 2113 2147 2148 2149 2154 2155 2156 2165 2191 2209 2217 2240 2253 2254 2258 2276 2277 2278 2321 2322 2329 ...
 2330 2374 2375 2376 2377 2378 2399 2403 2407 2408 2409 2444 2467 2468 2469 2475 2478 2479 2480 2484 2517 2536 2542 2543 2544 2545 2546 2552 2553 2561 2579 2585 ...
 2587 2588 2597 2599 2600 2601 2617 2645 2647 2651 2659 2756 2766 2767 2770 2782 2860 2874 2912 2927 2928 2929 2979 2980 2981 2982 3091 3092 3199 3200 3204 9999];  % 186 total


% Zero out LT,glat arrays for accumulations
for i=1:24
 for j=1:60
  for iseas=1:4
   for k=nc1:nc2
   num(i,j,iseas,k) = 0;
   den(i,j,iseas,k) = 0.;
   un(i,j,iseas,k) = 0.;
   ht(i,j,iseas,k) = 0.;
    for jl=1:24
     numl(i,j,jl,iseas,k) = 0;
     denl(i,j,jl,iseas,k) = 0.;
     unl(i,j,jl,iseas,k) = 0.;
     htl(i,j,jl,iseas,k) = 0.;
    end
   end
  end
 end
end

htn = 999.;
htx = 0.;
fx = 0.;

nskp=apskp1(iy1)+shortd1(1)-1;
if (iy1==1)
 nskp = apskp1(iy1) + 517 - mjd01(1);
end
for iy=iy1:iy2
 iyear0 = iy+2000
ndpy = 365;
if (iyear0==2000 | iyear0==2004 | iyear0==2008)
  ndpy = 366;
end

if (iy==iy1)
% Get daily av Kp(5), Ap(8), Vsw(37) from eln_tot_72_11_combnx7_daily (daynum=2)
filap = sprintf('/hao/polar1/emery/aurora_i/xuw/combined_output_daily/eln_tot_72_11_combnx7_daily')
INAP = fopen(filap,'r');
% Skip over lines
 for n=1:nskp
     rdhd = fgetl(INAP);
 end
rdhd
end	% if (iy==iy1)

is1 = 1;
if (iy==1)
 is1 = 2;
end
 for iseas=is1:4
  nda1 = shortd1(iseas)
  nda2 = shortd1(iseas+1)-1;
  if (iseas==4)
   if (iyear0==2004 | iyear0==2008)
    nda2 = nda2 + 1;
   end
  end
  mjd1 = mjd01(iy) + nda1 - 1;
  if (mjd1 < 517)
   nda1 = 152
   mjd1 = 517;
  end
  mjd2 = mjd1 + nda2 - nda1;
  if (iseas==4)
   mjd2 = min(3285,mjd01(iy+1) + shortd1(1) - 2);
   nda2 = nda1 + mjd2 - mjd1;
  end

for mjd=mjd1:mjd2
 nda = nda1 + mjd - mjd1;
iyr = iyear0;
jda = nda;
if (nda>ndpy)
 iyr = iyear0+1;
 jda = nda - ndpy;
end
% Read indices file
    rdap = fgetl(INAP);
    datap = sscanf(rdap,'%f');
    kpav = datap(5);
    apav = datap(8);
    vswav = datap(37);
    f107 = datap(10);
    jyr = datap(1);
    nday = datap(2);
    if (nday~=jda)
jyr,nday,iyr,jda,iyear0,nda,ndpy
stop
    end
   fx = max(f107,fx);

ifskp = 0;
if (iy1==1 & iy2==8)
 for nm=1:186
  if (mjd==miss_all(nm))
    ifskp=1;
  end
 end
end
if (iy1==1 & iy2==4)
 for nm=1:87
  if (mjd==miss_sum01(nm))
    ifskp=1;
  end
 end
end
if (iy1==4)
 for nm=1:126
  if (mjd==miss_spr04(nm))
    ifskp=1;
  end
 end
end


% Set nt from Ap and F107 conditions
nt = 0;
% Kp1=0-1Ap=0-4(7), Kp2=2-4-Ap=7-22(27);
% F1=65-75,F2=75-95,F3=95-125,F4=125-175,F5>175
% F1=65-85,F2=85-105,F3=105-145,F4=145-195,F5>195
 if (apav<7 & f107<85)
  nt = 1;
 end
 if (apav>=7 & apav<27 & f107<85)
  nt = 2;
 end
 if (apav<7 & f107>=85 & f107<105)
  nt = 3;
 end
 if (apav>=7 & apav<27 & f107>=85 & f107<105)
  nt = 4;
 end
 if (apav<7 & f107>=105 & f107<145)
  nt = 5;
 end
 if (apav>=7 & apav<27 & f107>=105 & f107<145)
  nt = 6;
 end
 if (apav<7 & f107>=145 & f107<195)
  nt = 7;
 end
 if (apav>=7 & apav<27 & f107>=145 & f107<195)
  nt = 8;
 end
 if (apav<7 & f107>=195)
  nt = 9;
 end
 if (apav>=7 & apav<27 & f107>=195)
  nt = 10;
 end

% Skip 01279-1,289
if (iyr==2001 & nda>=279 & nda<=281)
 ifskp=1;
nda,apav,f107
end


if (ifskp==0 & nt>0)
clear datkh;

filin=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/khwinds/MJD/CHAMP_MJD_%4.4d.txt',mjd);

datkh = dlmread(filin,'',0,0);
nll = length(datkh);
for n=1:nll

% 13 inputs: mjd,slt,mlt,glat,glat360,glon,mlat,mlon,alt,den,Vn,Un,Wn
% vn and wn are components of the cross-track obs in the +N and +up directions (are small so ignore them)
 nlt = fix(datkh(n,2)) + 1;
 nlat = fix((datkh(n,4)+90)/3.) + 1;

% Accumulate in 24 15-deg arrays with glon
 ilon = fix((datkh(n,6)+180)/15.) + 1;
 numl(nlt,nlat,ilon,iseas,nt) = numl(nlt,nlat,ilon,iseas,nt) + 1;
 unl(nlt,nlat,ilon,iseas,nt) = unl(nlt,nlat,ilon,iseas,nt) + datkh(n,12);
 denl(nlt,nlat,ilon,iseas,nt) = denl(nlt,nlat,ilon,iseas,nt) + datkh(n,10);
 htl(nlt,nlat,ilon,iseas,nt) = htl(nlt,nlat,ilon,iseas,nt) + datkh(n,9) - 6371.;

% Accumulate in arrays the density values etc
 num(nlt,nlat,iseas,nt) = num(nlt,nlat,iseas,nt) + 1;
 un(nlt,nlat,iseas,nt) = un(nlt,nlat,iseas,nt) + datkh(n,12);
% Find problem fall nt=7-8 9.5LT Un>0, 21.5LT Un<0 
% 01280 n=?6167-7164
% 01281 n= 926-931 ,1899-7457-8415
%if (nlt==10 & datkh(n,12)>400. & abs(datkh(n,4))<50.)
%  fprintf (OUTB,'%d %d %d %5.1f %5.1f %7.2f %7.2f %e \n',iyr,nda,n,nlt+1.5,datkh(n,4),datkh(n,12),datkh(n,9),datkh(n,10));
%end
%if (nlt==22 & datkh(n,12)<-400. & abs(datkh(n,4))<50.)
%  fprintf (OUTB,'%d %d %d %5.1f %7.2f %7.2f %e \n',iyr,nda,n,nlt+1.5,datkh(n,4),datkh(n,12),datkh(n,9),datkh(n,10));
%end
 den(nlt,nlat,iseas,nt) = den(nlt,nlat,iseas,nt) + datkh(n,10);
 ht(nlt,nlat,iseas,nt) = ht(nlt,nlat,iseas,nt) + datkh(n,9) - 6371.;
   if (nlat>=27 & nlat<=34)
    htn = min(htn,datkh(n,9)-6378.);
    htx = max(htx,datkh(n,9)-6378.);
   end

end	% for n=1,nll

end	% if (ifskp==0)
end	% for mjd=mjd1:mjd2
end	%  for iseas=1:4
end	% for iy=iy1:iy2
htn,htx		% 320-388km for low-lats from 2004-8 and 364-461km from 2001-4
fx

% Print out LT files 
for nh=1:24
 for ilat=1:13
  for iseas=1:4
   for nt=nc1:nc2
    num12(nh,ilat,iseas,nt) = 0;
    un12(nh,ilat,iseas,nt) = 0.;
    d400(nh,ilat,iseas,nt) = 0.;
    ht12(nh,ilat,iseas,nt) = 0.;
    for jl=1:24
     num12l(nh,ilat,jl,iseas,nt) = 0;
     un12l(nh,ilat,jl,iseas,nt) = 0.;
     d400l(nh,ilat,jl,iseas,nt) = 0.;
     ht12l(nh,ilat,jl,iseas,nt) = 0.;
    end
   end
  end
 end
end

for nt=nc1:nc2
 for iseas=1:4

clear seasnt;
if (iseas==1 & nt==1)
 seasnt=sprintf('_%s_%s',titseas1,titcond1);
end
if (iseas==2 & nt==1)
 seasnt=sprintf('_%s_%s',titseas2,titcond1);
end
if (iseas==3 & nt==1)
 seasnt=sprintf('_%s_%s',titseas3,titcond1);
end
if (iseas==4 & nt==1)
 seasnt=sprintf('_%s_%s',titseas4,titcond1);
end
if (iseas==1 & nt==2)
 seasnt=sprintf('_%s_%s',titseas1,titcond2);
end
if (iseas==2 & nt==2)
 seasnt=sprintf('_%s_%s',titseas2,titcond2);
end
if (iseas==3 & nt==2)
 seasnt=sprintf('_%s_%s',titseas3,titcond2);
end
if (iseas==4 & nt==2)
 seasnt=sprintf('_%s_%s',titseas4,titcond2);
end
if (iseas==1 & nt==3)
 seasnt=sprintf('_%s_%s',titseas1,titcond3);
end
if (iseas==2 & nt==3)
 seasnt=sprintf('_%s_%s',titseas2,titcond3);
end
if (iseas==3 & nt==3)
 seasnt=sprintf('_%s_%s',titseas3,titcond3);
end
if (iseas==4 & nt==3)
 seasnt=sprintf('_%s_%s',titseas4,titcond3);
end
if (iseas==1 & nt==4)
 seasnt=sprintf('_%s_%s',titseas1,titcond4);
end
if (iseas==2 & nt==4)
 seasnt=sprintf('_%s_%s',titseas2,titcond4);
end
if (iseas==3 & nt==4)
 seasnt=sprintf('_%s_%s',titseas3,titcond4);
end
if (iseas==4 & nt==4)
 seasnt=sprintf('_%s_%s',titseas4,titcond4);
end
if (iseas==1 & nt==5)
 seasnt=sprintf('_%s_%s',titseas1,titcond5);
end
if (iseas==2 & nt==5)
 seasnt=sprintf('_%s_%s',titseas2,titcond5);
end
if (iseas==3 & nt==5)
 seasnt=sprintf('_%s_%s',titseas3,titcond5);
end
if (iseas==4 & nt==5)
 seasnt=sprintf('_%s_%s',titseas4,titcond5);
end
if (iseas==1 & nt==6)
 seasnt=sprintf('_%s_%s',titseas1,titcond6);
end
if (iseas==2 & nt==6)
 seasnt=sprintf('_%s_%s',titseas2,titcond6);
end
if (iseas==3 & nt==6)
 seasnt=sprintf('_%s_%s',titseas3,titcond6);
end
if (iseas==4 & nt==6)
 seasnt=sprintf('_%s_%s',titseas4,titcond6);
end
if (iseas==1 & nt==7)
 seasnt=sprintf('_%s_%s',titseas1,titcond7);
end
if (iseas==2 & nt==7)
 seasnt=sprintf('_%s_%s',titseas2,titcond7);
end
if (iseas==3 & nt==7)
 seasnt=sprintf('_%s_%s',titseas3,titcond7);
end
if (iseas==4 & nt==7)
 seasnt=sprintf('_%s_%s',titseas4,titcond7);
end
if (iseas==1 & nt==8)
 seasnt=sprintf('_%s_%s',titseas1,titcond8);
end
if (iseas==2 & nt==8)
 seasnt=sprintf('_%s_%s',titseas2,titcond8);
end
if (iseas==3 & nt==8)
 seasnt=sprintf('_%s_%s',titseas3,titcond8);
end
if (iseas==4 & nt==8)
 seasnt=sprintf('_%s_%s',titseas4,titcond8);
end
if (iseas==1 & nt==9)
 seasnt=sprintf('_%s_%s',titseas1,titcond9);
end
if (iseas==2 & nt==9)
 seasnt=sprintf('_%s_%s',titseas2,titcond9);
end
if (iseas==3 & nt==9)
 seasnt=sprintf('_%s_%s',titseas3,titcond9);
end
if (iseas==4 & nt==9)
 seasnt=sprintf('_%s_%s',titseas4,titcond9);
end
if (iseas==1 & nt==10)
 seasnt=sprintf('_%s_%s',titseas1,titcond10);
end
if (iseas==2 & nt==10)
 seasnt=sprintf('_%s_%s',titseas2,titcond10);
end
if (iseas==3 & nt==10)
 seasnt=sprintf('_%s_%s',titseas3,titcond10);
end
if (iseas==4 & nt==10)
 seasnt=sprintf('_%s_%s',titseas4,titcond10);
end
  outcnofs = sprintf('champ_winds_13glat_vs_lt%sb',seasnt)
  OUTC = fopen(outcnofs,'w');
  fprintf (OUTC,'SLT 13(# Un(m/s) NDen(kg/m3) ht(km)) for glat=-90:15:90 %s \n',seasnt);
  outlon = sprintf('champ_winds_equ_vs_lt_glon%s',seasnt)
  OUTL = fopen(outlon,'w');
  fprintf (OUTL,'SLT equ(# Un(m/s) NDen(kg/m3) ht(km)) for glon=7.5:15:352.5 %s \n',seasnt);

for nh=1:24

 for nl=1:60
  if (num(nh,nl,iseas,nt)>0)

% 60 3-deg glats from -88.5:3:88.5 (collect 15deg -60:15:60)
% Overlap 1-2 ilat values!
   nls = 0;
   if (nl>=27 & nl<=34) 
% Average over nl=27-34 (+/-12glat) for C/NOFS
    nls = nls + 1;
    ilatsv(nls) = 7;
   end
   if (nl>=23 & nl<=28) 
    nls = nls + 1;
    ilatsv(nls) = 6;
   end
   if (nl>=18 & nl<=23) 
    nls = nls + 1;
    ilatsv(nls) = 5;
   end
   if (nl>=13 & nl<=18) 
    nls = nls + 1;
    ilatsv(nls) = 4;
   end
   if (nl>=8 & nl<=13) 
    nls = nls + 1;
    ilatsv(nls) = 3;
   end
   if (nl>=3 & nl<=8) 
    nls = nls + 1;
    ilatsv(nls) = 2;
   end
   if (nl>=1 & nl<=3) 
    nls = nls + 1;
    ilatsv(nls) = 1;
   end
   if (nl>=33 & nl<=38) 
    nls = nls + 1;
    ilatsv(nls) = 8;
   end
   if (nl>=38 & nl<=43) 
    nls = nls + 1;
    ilatsv(nls) = 9;
   end
   if (nl>=43 & nl<=48) 
    nls = nls + 1;
    ilatsv(nls) = 10;
   end
   if (nl>=48 & nl<=53) 
    nls = nls + 1;
    ilatsv(nls) = 11;
   end
   if (nl>=53 & nl<=58) 
    nls = nls + 1;
    ilatsv(nls) = 12;
   end
   if (nl>=58 & nl<=60) 
    nls = nls + 1;
    ilatsv(nls) = 13;
   end
   for jlat=1:nls
    ilat = ilatsv(jlat);
    num12(nh,ilat,iseas,nt) = num12(nh,ilat,iseas,nt) + num(nh,nl,iseas,nt);
    un12(nh,ilat,iseas,nt) = un12(nh,ilat,iseas,nt) + un(nh,nl,iseas,nt);
    d400(nh,ilat,iseas,nt) = d400(nh,ilat,iseas,nt) + den(nh,nl,iseas,nt);
    ht12(nh,ilat,iseas,nt) = ht12(nh,ilat,iseas,nt) + ht(nh,nl,iseas,nt);
    for ilon=1:24
     if (numl(nh,nl,ilon,iseas,nt)>0)
     num12l(nh,ilat,ilon,iseas,nt) = num12l(nh,ilat,ilon,iseas,nt) + numl(nh,nl,ilon,iseas,nt);
     un12l(nh,ilat,ilon,iseas,nt) = un12l(nh,ilat,ilon,iseas,nt) + unl(nh,nl,ilon,iseas,nt);
     d400l(nh,ilat,ilon,iseas,nt) = d400l(nh,ilat,ilon,iseas,nt) + denl(nh,nl,ilon,iseas,nt);
     ht12l(nh,ilat,ilon,iseas,nt) = ht12l(nh,ilat,ilon,iseas,nt) + htl(nh,nl,ilon,iseas,nt);
     end
    end  
   end	%  for jlat=1:nls 

% Now divide by num to get averages
    un(nh,nl,iseas,nt) = un(nh,nl,iseas,nt)/num(nh,nl,iseas,nt);
    den(nh,nl,iseas,nt) = den(nh,nl,iseas,nt)/num(nh,nl,iseas,nt);
    ht(nh,nl,iseas,nt) = ht(nh,nl,iseas,nt)/num(nh,nl,iseas,nt);
    for ilon=1:24
     if (numl(nh,nl,ilon,iseas,nt)>0)
     unl(nh,ilat,ilon,iseas,nt) = unl(nh,ilat,ilon,iseas,nt)/numl(nh,nl,ilon,iseas,nt);
     denl(nh,ilat,ilon,iseas,nt) = denl(nh,ilat,ilon,iseas,nt)/numl(nh,nl,ilon,iseas,nt);
     htl(nh,ilat,ilon,iseas,nt) = htl(nh,ilat,ilon,iseas,nt)/numl(nh,nl,ilon,iseas,nt);
     end
    end   

  end	% if (num(nh,nl,iseas,nt)>0)

% Set 0's to NaN for plots
  if (num(nh,nl,iseas,nt)==0)
   un(nh,nl,iseas,nt) = NaN;
   den(nh,nl,iseas,nt) = NaN;
   ht(nh,nl,iseas,nt) = NaN;
  else
   den(nh,nl,iseas,nt) = den(nh,nl,iseas,nt)*1.e+13;
  end
 end	% nl=1:60

% Print out C/NOFS info as func of LT only
  for ilat=1:13
  if (num12(nh,ilat,iseas,nt)>0);
    un12(nh,ilat,iseas,nt) = un12(nh,ilat,iseas,nt)/num12(nh,ilat,iseas,nt);
    d400(nh,ilat,iseas,nt) = d400(nh,ilat,iseas,nt)/num12(nh,ilat,iseas,nt);
    ht12(nh,ilat,iseas,nt) = ht12(nh,ilat,iseas,nt)/num12(nh,ilat,iseas,nt);
    for ilon=1:24
     if (num12l(nh,ilat,ilon,iseas,nt)>0)
     un12l(nh,ilat,ilon,iseas,nt) = un12l(nh,ilat,ilon,iseas,nt)/num12l(nh,ilat,ilon,iseas,nt);
     d400l(nh,ilat,ilon,iseas,nt) = d400l(nh,ilat,ilon,iseas,nt)/num12l(nh,ilat,ilon,iseas,nt);
     ht12l(nh,ilat,ilon,iseas,nt) = ht12l(nh,ilat,ilon,iseas,nt)/num12l(nh,ilat,ilon,iseas,nt);
     end
    end  
  end
  end
   fprintf(OUTC,'%5.1f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f %d %8.2f %12.4e %5.0f\n',lt24(nh),num12(nh,1,iseas,nt),un12(nh,1,iseas,nt),d400(nh,1,iseas,nt),ht12(nh,1,iseas,nt), ...
num12(nh,2,iseas,nt),un12(nh,2,iseas,nt),d400(nh,2,iseas,nt),ht12(nh,2,iseas,nt),num12(nh,3,iseas,nt),un12(nh,3,iseas,nt),d400(nh,3,iseas,nt),ht12(nh,3,iseas,nt), ...
num12(nh,4,iseas,nt),un12(nh,4,iseas,nt),d400(nh,4,iseas,nt),ht12(nh,4,iseas,nt),num12(nh,5,iseas,nt),un12(nh,5,iseas,nt),d400(nh,5,iseas,nt),ht12(nh,5,iseas,nt), ...
num12(nh,6,iseas,nt),un12(nh,6,iseas,nt),d400(nh,6,iseas,nt),ht12(nh,6,iseas,nt),num12(nh,7,iseas,nt),un12(nh,7,iseas,nt),d400(nh,7,iseas,nt),ht12(nh,7,iseas,nt), ...
num12(nh,8,iseas,nt),un12(nh,8,iseas,nt),d400(nh,8,iseas,nt),ht12(nh,8,iseas,nt),num12(nh,9,iseas,nt),un12(nh,9,iseas,nt),d400(nh,9,iseas,nt),ht12(nh,9,iseas,nt), ...
num12(nh,10,iseas,nt),un12(nh,10,iseas,nt),d400(nh,10,iseas,nt),ht12(nh,10,iseas,nt),num12(nh,11,iseas,nt),un12(nh,11,iseas,nt),d400(nh,11,iseas,nt),ht12(nh,11,iseas,nt), ...
num12(nh,12,iseas,nt),un12(nh,12,iseas,nt),d400(nh,12,iseas,nt),ht12(nh,12,iseas,nt),num12(nh,13,iseas,nt),un12(nh,13,iseas,nt),d400(nh,13,iseas,nt),ht12(nh,13,iseas,nt));
ilat = 7;
fprintf(OUTL,'%5.1f',lt24(nh));
 for ilon=1:24
     fprintf(OUTL,' %d %8.2f %12.4e %5.0f',num12l(nh,ilat,ilon,iseas,nt),un12l(nh,ilat,ilon,iseas,nt),d400l(nh,ilat,ilon,iseas,nt),ht12l(nh,ilat,ilon,iseas,nt));
 end
 fprintf(OUTL,'\n');
end	%  for nh=1:24   

end 	%  for iseas=1:4
end	%  for nt=nc1:nc2

% Cannot go from mid-pt to one extra edge for plotting - bad dimensions!
%  But still only see 23 hours, not 24??? on plots?? is because is edging?
%lt24 = 0:1:24;
%lat = -90:3:90;

fsize = 12;
unx=300.;
unn=-300.;

for iseas=1:4

if (iseas==1)
 seas=sprintf('%s',titseas1);
end
if (iseas==2)
 seas=sprintf('%s',titseas2);
end
if (iseas==3)
 seas=sprintf('%s',titseas3);
end
if (iseas==4)
 seas=sprintf('%s',titseas4);
end

% Plot lo and hi Kp days together
hn = 320.;
hx = 460.;
nden = 0.;

for jy=1:2
iy1 = 4;
if (jy==2)
 iy1 = 1;
 iy2 = 4
end
 
if (iy1==4)
xden1 = 10;
xden2 = 15;
nc1 = 1;
nc2 = 2;
nc3 = 3;
nc4 = 4;
tit1ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond1);
tit2ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond2);
tit3ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond3);
tit4ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond4);
tit1un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond1);
tit2un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond2);
tit3un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond3);
tit4un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond4);
tit1ht = sprintf('KH CHAMP ht km %s %s',seas,titcond1);
tit2ht = sprintf('KH CHAMP ht km %s %s',seas,titcond2);
tit3ht = sprintf('KH CHAMP ht km %s %s',seas,titcond3);
tit4ht = sprintf('KH CHAMP ht km %s %s',seas,titcond4);
end
if (iy1==1 & iy2==4)
xden1 = 30;
xden2 = 60;
xden3 = 90;
nc1 = 5;
nc2 = 6;
nc3 = 7;
nc4 = 8;
nc5 = 9;
nc6 = 10;
tit1ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond5);
tit2ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond6);
tit3ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond7);
tit4ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond8);
tit1un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond5);
tit2un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond6);
tit3un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond7);
tit4un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond8);
tit1ht = sprintf('KH CHAMP ht km %s %s',seas,titcond5);
tit2ht = sprintf('KH CHAMP ht km %s %s',seas,titcond6);
tit3ht = sprintf('KH CHAMP ht km %s %s',seas,titcond7);
tit4ht = sprintf('KH CHAMP ht km %s %s',seas,titcond8);
tit5ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond9);
tit6ave = sprintf('KH CHAMP density*1.e+13 in kg/m3 %s %s',seas,titcond10);
tit5un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond9);
tit6un = sprintf('KH CHAMP Un(+E) m/s %s %s',seas,titcond10);
tit5ht = sprintf('KH CHAMP ht km %s %s',seas,titcond9);
tit6ht = sprintf('KH CHAMP ht km %s %s',seas,titcond10);
end

figure
clf
subplot(2,2,1)
            pcolor(lt24,lat,un(:,:,iseas,nc1)')
            title([tit1un]);
            caxis([unn unx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,2)
            pcolor(lt24,lat,un(:,:,iseas,nc2)')
            title([tit2un]);
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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');
subplot(2,2,3)
            pcolor(lt24,lat,un(:,:,iseas,nc3)')
            title([tit3un]);
            caxis([unn unx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,4)
            pcolor(lt24,lat,un(:,:,iseas,nc4)')
            title([tit4un]);
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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;	

figure
clf
subplot(2,2,1)
            pcolor(lt24,lat,den(:,:,iseas,nc1)')
            title([tit1ave]);
            caxis([nden xden1]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,2,2)
            pcolor(lt24,lat,den(:,:,iseas,nc2)')
            title([tit2ave]);
            caxis([nden xden1]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,2,3)
            pcolor(lt24,lat,den(:,:,iseas,nc3)')
            title([tit3ave]);
            caxis([nden xden2]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,2,4)
            pcolor(lt24,lat,den(:,:,iseas,nc4)')
            title([tit4ave]);
            caxis([nden xden2]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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;	
	
figure
clf
subplot(2,2,1)
            pcolor(lt24,lat,ht(:,:,iseas,nc1)')
            title([tit1ht]);
            caxis([hn hx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,2)
            pcolor(lt24,lat,ht(:,:,iseas,nc2)')
            title([tit2ht]);
            caxis([hn hx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,3)
            pcolor(lt24,lat,ht(:,:,iseas,nc3)')
            title([tit3ht]);
            caxis([hn hx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,4)
            pcolor(lt24,lat,ht(:,:,iseas,nc4)')
            title([tit4ht]);
            caxis([hn hx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');

if (iy1==1 & iy2==4)
figure
clf
subplot(2,2,1)
            pcolor(lt24,lat,un(:,:,iseas,nc5)')
            title([tit5un]);
            caxis([unn unx]);
            shading flat
	    grid;
            ylabel('Latitude');
            xlabel('LT');
	    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');
subplot(2,2,2)
            pcolor(lt24,lat,un(:,:,iseas,nc6)')
            title([tit6un]);
            caxis([unn unx]);
            shading flat
	    grid;
            xlabel('LT');
            ylabel('Latitude');
	    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');
subplot(2,2,3)
            pcolor(lt24,lat,den(:,:,iseas,nc5)')
            title([tit5ave]);
            caxis([nden xden3]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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');
subplot(2,2,4)
            pcolor(lt24,lat,den(:,:,iseas,nc6)')
            title([tit6ave]);
            caxis([nden xden3]);
            shading flat
	    grid;
            ylabel('Latitude');
	    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	% if (iy1==1 & iy2==4)

end	% for (jy=1:2

end	% for iseas=1:4
