clear all;
close all;

% Dec 2012  Program to calculate the ave of CHAMP PLP Te and Ne 
%    into 3-deg lat bins and 1-h LT bins for various days
%  02050-09355 PLP available from 00208-10052 CHAMP lifetime

fsize = 12;

nc1 = 1
nc2 = 14
%iy1 = 2
%iy2 = 9
iy1 = 7
iy2 = 8
% Days with F107>=185: max (275.4 on 03302)
n185 = 0;

  outbad = sprintf('champ_plp_bad2')
  OUTB = fopen(outbad,'w');
  fprintf (OUTB,'iyr nda n slt glat Te \n');

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

titseas1='spr';
titseas2='sum';
titseas3='fal';
titseas4='win';
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-125';
titcond6='Kp_2-4_F107_105-125';
titcond7='Kp_0-1_F107_125-145';
titcond8='Kp_2-4_F107_125-145';
titcond9='Kp_0-1_F107_145-165';
titcond10='Kp_2-4_F107_145-165';
titcond11='Kp_0-1_F107_165-185';
titcond12='Kp_2-4_F107_165-185';
titcond13='Kp_0-1_F107_185-285';
titcond14='Kp_2-4_F107_185-285';



nday1 = -999;

% apskp1 number of lines in 1972-2011 file to skip over for years 2001-2009
apskp1=[10593 10958 11323 11688 12054 12419 12784 13149 13515];
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)

% 02050-08007 none missing (do not have 06151-365, 07151-289, 2000-2002)
miss08 =[8 10 11 12 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 75 78 79 80 81 82 83 84 85 86 87 88 89 90 92 93 94 95 96 97 98 100 101 102 103 104 105 106 107 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 199 200 205 206 207 208 209 214 221 227 229 230 231 232 233 239 240 242 243 249 250 251 253 254 255 263 264 266 267 272 273 274 276 277 278 284 285 286 287 288 289 290 291 292 293 294 296 297 298 299 300 301 302 309 323 325 326 327 328 329 330 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354];	% 217 missing
miss09 =[ 10 11 12 13 14 15 16 17 18 19 28 29 30 31 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 71 72 73 74 76 77 78 79 80 90 91 92 94 95 96 97 113 114 115 116 117 119 120 121 122 123 124 125 126 127 128 129 130 131 142 143 144 145 146 147 148 149 150 151 152 153 155 160 161 162 163 164 165 166 168 169 170 171 172 173 174 175 176 177 178 179 180 182 211 212 213 214 215 216 217 219 220 236 237 238 239 240 241 242 243 244 245 246 248 249 250 251 252 253 262 263 264 265 266 267 268 269 270 271 272 273 275 276 277 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317  320 321 322  324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341  343 344 345 346 347 348 349];	% 200 missing

miss07_08=[8 10 11 12 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51];	% 42 missing
% Added 2 more days (09056-7) to 132 from 08290-09055 as missing
miss08_09=[290 291 292 293 294 296 297 298 299 300 301 302 309 323 325 326 327 328 329 330 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 10 11 12 13 14 15 16 17 18 19 28 29 30 31 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57];	% 82=43+39 missing
% 09290-355=66 ndt since then CHPLP missing to 10052 (end)
miss09_10=[294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317  320 321 322  324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341  343 344 345 346 347 348 349]; 	% 52 missing

% 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.;
   te(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.;
     tel(i,j,jl,iseas,k) = 0.;
     htl(i,j,jl,iseas,k) = 0.;
    end
   end
  end
 end
end

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

neqlt = 24;
xeqlt = 0;

nht = 5.e+7;
xht = 0.;


for iy=iy1:iy2
 iyear0 = iy+2000
 ndpy = 365;

nd1jan = 1;
nd1feb = 32;
nd1mar = 60;
nd1apr = 91;
nd1may = 121;
nd1jun = 152;
nd1jul = 182;
nd1aug = 213;
nd1sep = 244;
nd1oct = 274;
nd1nov = 305;
nd1dec = 335;
if (iyear0==2000 | iyear0==2004 | iyear0==2008)
  ndpy = 366;
 nd1mar = nd1mar + 1;
 nd1apr = nd1apr + 1;
 nd1may = nd1may + 1;
 nd1jun = nd1jun + 1;
 nd1jul = nd1jul + 1;
 nd1aug = nd1aug + 1;
 nd1sep = nd1sep + 1;
 nd1oct = nd1oct + 1;
 nd1nov = nd1nov + 1;
 nd1dec = nd1dec + 1;
end

if (iy==iy1)
 nskp=apskp1(iy1)+shortd1(1)-1;
 if (iy==2)
  nskp = apskp1(iy1)+50-1;
 end
% 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;
is2 = 4;
 for iseas=is1:is2
  nda1 = shortd1(iseas)
  if (iyear0==2002 & iseas==1)
    nda1 = 50;
  end
  nda2 = shortd1(iseas+1)-1;
%  Add extra day for leap year to season 4 (NH winter)
  if (iseas==4)   
   if (iyear0==2004 | iyear0==2008)
    nda2 = nda2 + 1;
   end
  end
  if (iy==9 & iseas==4)
   nda2 = 355;
  end

nbad = 0;
for nda=nda1:nda2
iyr = iyear0;
jda = nda;
if (nda>ndpy)
 iyr = iyear0+1;
 jda = nda - ndpy;
end
% Find month and daynumber
 if (nda<nd1feb)
  mo = 1;
  moday = jda;
 end
 if (nda>=nd1feb & nda<nd1mar)
  mo = 2;
  moday = jda - nd1feb + 1;
 end
 if (nda>=nd1mar & nda<nd1apr)
  mo = 3;
  moday = jda - nd1mar + 1;
 end
 if (nda>=nd1apr & nda<nd1may)
  mo = 4;
  moday = jda - nd1apr + 1;
 end
 if (nda>=nd1may & nda<nd1jun)
  mo = 5;
  moday = jda - nd1may + 1;
 end 
 if (nda>=nd1jun & nda<nd1jul)
  mo = 6;
  moday = nda - nd1jun + 1;
 end
 if (nda>=nd1jul & nda<nd1aug)
  mo = 7;
  moday = nda - nd1jul + 1;
 end
 if (nda>=nd1aug & nda<nd1sep)
  mo = 8;
  moday = nda - nd1aug + 1;
 end
 if (nda>=nd1sep & nda<nd1oct)
  mo = 9;
  moday = nda - nd1sep + 1;
 end
 if (nda>=nd1oct & nda<nd1nov)
  mo = 10;
  moday = nda - nd1oct + 1;
 end
 if (nda<nd1dec & nda>=nd1nov)
  mo = 11;
  moday = nda - nd1nov + 1;
 end
 if (nda>=nd1dec)
  mo = 12;
  moday = nda - nd1dec + 1;
 end
 if (nda>ndpy & jda<nd1feb)
  mo = 1;
  moday = jda;
 end
 if (nda>ndpy & jda>=nd1feb)
  mo = 2;
  moday = jda - nd1feb + 1;
 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);
%  Check to see if 165 found for later years
    if (f107>=185)
n185 = n185 + 1;
nday,f107
    end
    if (nday~=jda)
jyr,nday,iyr,jda,iyear0,nda,ndpy
stop
    end
   fx = max(f107,fx);

ifskp = 0;
% Skip missing days for PLP
if (iyr==2008)
 for nm=1:217
  if (jda==miss08(nm))
    ifskp=1;
  end
 end
end	% if iyr==2008
if (iyr==2009)
 for nm=1:200		
  if (jda==miss09(nm))
    ifskp=1;
  end
 end
end	% if iyr==2009

% Set nt from Ap and F107 conditions
nt = 0;
% Kp1=0-1Ap=0-4(7), Kp2=2-4-Ap=7-22(27);
% F1=65-85,F2=85-105,F3=105-125,F4=125-145,F5=145-165
 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<125)
  nt = 5;
 end
 if (apav>=7 & apav<27 & f107>=105 & f107<125)
  nt = 6;
 end
 if (apav<7 & f107>=125 & f107<145)
  nt = 7;
 end
 if (apav>=7 & apav<27 & f107>=125 & f107<145)
  nt = 8;
 end
 if (apav<7 & f107>=145 & f107<165)
  nt = 9;
 end
 if (apav>=7 & apav<27 & f107>=145 & f107<165)
  nt = 10;
 end
 if (apav<7 & f107>=165 & f107<185)
  nt = 11;
 end
 if (apav>=7 & apav<27 & f107>=165 & f107<185)
  nt = 12;
 end
 if (apav<7 & f107>=185 & f107<285)
  nt = 13;
 end
 if (apav>=7 & apav<27 & f107>=185 & f107<285)
  nt = 14;
 end

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

filin=sprintf('/hao/aim3/emery/cetic_climatology/champ_density/plp_ne_te/CH-ME-2-PLP+%d-%2.2d-%2.2d_21.dat',iyr,mo,moday);

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

% 12 inputs: gpssec year mo day hh mm ss radius glat glon Ne(cm-3) Te(K)
% calculate local time (LT = UT + glon/15.)
 slt = datkh(n,5)+datkh(n,6)/60.+datkh(n,7)/3600.+datkh(n,10)/15.;
 if (slt<0.)
  slt = slt + 24.;
 end
 if (slt>=24.)
  slt = slt - 24.;
 end
 nlt = fix(slt) + 1;

% Find bad night pts 2007 1.5LT 7.5glat, 21.5LT -25.5glat
ifbad=0;
 if (abs(datkh(n,9))<30. & slt<4. & datkh(n,12)>6000.)
  fprintf (OUTB,'%d %d %d %f %f %f \n',iyear0,nda,n,slt,datkh(n,9),datkh(n,12));
  ifbad=1;
 end
 if (abs(datkh(n,9))<30. & slt>20. & datkh(n,12)>6000.)
  fprintf (OUTB,'%d %d %d %f %f %f \n',iyear0,nda,n,slt,datkh(n,9),datkh(n,12));
  ifbad=1;
 end

if (ifbad==0)

 nlat = fix((datkh(n,9)+90)/3.) + 1;

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

% Accumulate in arrays Ne,Te values 
 num(nlt,nlat,iseas,nt) = num(nlt,nlat,iseas,nt) + 1;
 te(nlt,nlat,iseas,nt) = te(nlt,nlat,iseas,nt) + datkh(n,12);
 den(nlt,nlat,iseas,nt) = den(nlt,nlat,iseas,nt) + datkh(n,11);
 ht(nlt,nlat,iseas,nt) = ht(nlt,nlat,iseas,nt) + datkh(n,8) - 6371.;


else
 nbad = nbad + 1;
end	% if (ifbad==0)

end	% for n=1,nll
 
end	% if (ifskp==0)
end	% for nda=nda1:nda2
end	%  for iseas=is1:is2
end	% for iy=iy1:iy2
htn,htx		% 320-388km for low-lats from 2004-8 and 364-461km from 2001-4
fx,n185

% 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;
    te12(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;
     te12l(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
 n0 = 0;
% Calculate if there are any values per iseas,nt:
  for nh=1:24
   for jl=1:60
    if (num(nh,jl,iseas,nt)>0)
     n0 = n0 + 1;
    end
   end	% for jl=1:60
  end	% for nh=1:24
  num0(iseas,nt) = n0;

if (n0>0)
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
if (iseas==1 & nt==11)
 seasnt=sprintf('_%s_%s',titseas1,titcond11);
end
if (iseas==2 & nt==11)
 seasnt=sprintf('_%s_%s',titseas2,titcond11);
end
if (iseas==3 & nt==11)
 seasnt=sprintf('_%s_%s',titseas3,titcond11);
end
if (iseas==4 & nt==11)
 seasnt=sprintf('_%s_%s',titseas4,titcond11);
end
if (iseas==1 & nt==12)
 seasnt=sprintf('_%s_%s',titseas1,titcond12);
end
if (iseas==2 & nt==12)
 seasnt=sprintf('_%s_%s',titseas2,titcond12);
end
if (iseas==3 & nt==12)
 seasnt=sprintf('_%s_%s',titseas3,titcond12);
end
if (iseas==4 & nt==12)
 seasnt=sprintf('_%s_%s',titseas4,titcond12);
end
if (iseas==1 & nt==13)
 seasnt=sprintf('_%s_%s',titseas1,titcond13);
end
if (iseas==2 & nt==13)
 seasnt=sprintf('_%s_%s',titseas2,titcond13);
end
if (iseas==3 & nt==13)
 seasnt=sprintf('_%s_%s',titseas3,titcond13);
end
if (iseas==4 & nt==13)
 seasnt=sprintf('_%s_%s',titseas4,titcond13);
end
if (iseas==1 & nt==14)
 seasnt=sprintf('_%s_%s',titseas1,titcond14);
end
if (iseas==2 & nt==14)
 seasnt=sprintf('_%s_%s',titseas2,titcond14);
end
if (iseas==3 & nt==14)
 seasnt=sprintf('_%s_%s',titseas3,titcond14);
end
if (iseas==4 & nt==14)
 seasnt=sprintf('_%s_%s',titseas4,titcond14);
end
  outcnofs = sprintf('champ_plp_13glat_vs_lt%sb',seasnt)
  OUTC = fopen(outcnofs,'w');
  fprintf (OUTC,'SLT 13(# Te(K) NeDen(#/cm3) ht(km)) for glat=-90:15:90 %s \n',seasnt);
  outlon = sprintf('champ_plp_equ_vs_lt_glon%s',seasnt)
  OUTL = fopen(outlon,'w');
  fprintf (OUTL,'SLT equ(# Te(K) NeDen(#/cm3) 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);
    te12(nh,ilat,iseas,nt) = te12(nh,ilat,iseas,nt) + te(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);
     te12l(nh,ilat,ilon,iseas,nt) = te12l(nh,ilat,ilon,iseas,nt) + tel(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
    te(nh,nl,iseas,nt) = te(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)
     tel(nh,ilat,ilon,iseas,nt) = tel(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)
   te(nh,nl,iseas,nt) = NaN;
   den(nh,nl,iseas,nt) = NaN;
   ht(nh,nl,iseas,nt) = NaN;
  else
   den(nh,nl,iseas,nt) = log10(den(nh,nl,iseas,nt));
  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);
    te12(nh,ilat,iseas,nt) = te12(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)
     te12l(nh,ilat,ilon,iseas,nt) = te12l(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),te12(nh,1,iseas,nt),d400(nh,1,iseas,nt),ht12(nh,1,iseas,nt), ...
num12(nh,2,iseas,nt),te12(nh,2,iseas,nt),d400(nh,2,iseas,nt),ht12(nh,2,iseas,nt),num12(nh,3,iseas,nt),te12(nh,3,iseas,nt),d400(nh,3,iseas,nt),ht12(nh,3,iseas,nt), ...
num12(nh,4,iseas,nt),te12(nh,4,iseas,nt),d400(nh,4,iseas,nt),ht12(nh,4,iseas,nt),num12(nh,5,iseas,nt),te12(nh,5,iseas,nt),d400(nh,5,iseas,nt),ht12(nh,5,iseas,nt), ...
num12(nh,6,iseas,nt),te12(nh,6,iseas,nt),d400(nh,6,iseas,nt),ht12(nh,6,iseas,nt),num12(nh,7,iseas,nt),te12(nh,7,iseas,nt),d400(nh,7,iseas,nt),ht12(nh,7,iseas,nt), ...
num12(nh,8,iseas,nt),te12(nh,8,iseas,nt),d400(nh,8,iseas,nt),ht12(nh,8,iseas,nt),num12(nh,9,iseas,nt),te12(nh,9,iseas,nt),d400(nh,9,iseas,nt),ht12(nh,9,iseas,nt), ...
num12(nh,10,iseas,nt),te12(nh,10,iseas,nt),d400(nh,10,iseas,nt),ht12(nh,10,iseas,nt),num12(nh,11,iseas,nt),te12(nh,11,iseas,nt),d400(nh,11,iseas,nt),ht12(nh,11,iseas,nt), ...
num12(nh,12,iseas,nt),te12(nh,12,iseas,nt),d400(nh,12,iseas,nt),ht12(nh,12,iseas,nt),num12(nh,13,iseas,nt),te12(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),te12l(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	%  if (n0>0)
end 	%  for iseas=1:4
end	%  for nt=nc1:nc2

% Make plots
fsize = 12;
unx=300.;
unn=-300.;

for iseas=1:4
%  if (max(num(:,:,iseas,:))>0)

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.;

% 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;

tex=3000.;
ten=500.;

 alt = 350;
 nden = 3;
 xden1 = 6;
 xden2 = 6;

 for ifx=1:3

 if (ifx==1)
nc1 = 1;
nc2 = 2;
nc3 = 3;
nc4 = 4;
tit1ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond1);
tit2ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond2);
tit3ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond3);
tit4ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond4);
tit1un = sprintf('CHAMP Te (K) %s %s',seas,titcond1);
tit2un = sprintf('CHAMP Te (K) %s %s',seas,titcond2);
tit3un = sprintf('CHAMP Te (K) %s %s',seas,titcond3);
tit4un = sprintf('CHAMP Te (K) %s %s',seas,titcond4);
tit1ht = sprintf('CHAMP ht km %s %s',seas,titcond1);
tit2ht = sprintf('CHAMP ht km %s %s',seas,titcond2);
tit3ht = sprintf('CHAMP ht km %s %s',seas,titcond3);
tit4ht = sprintf('CHAMP ht km %s %s',seas,titcond4);
 end
 if (ifx==2)
nc1 = 5;
nc2 = 6;
nc3 = 7;
nc4 = 8;
tit1ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond5);
tit2ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond6);
tit3ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond7);
tit4ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond8);
tit1un = sprintf('CHAMP Te (K) %s %s',seas,titcond5);
tit2un = sprintf('CHAMP Te (K) %s %s',seas,titcond6);
tit3un = sprintf('CHAMP Te (K) %s %s',seas,titcond7);
tit4un = sprintf('CHAMP Te (K) %s %s',seas,titcond8);
tit1ht = sprintf('CHAMP ht km %s %s',seas,titcond5);
tit2ht = sprintf('CHAMP ht km %s %s',seas,titcond6);
tit3ht = sprintf('CHAMP ht km %s %s',seas,titcond7);
tit4ht = sprintf('CHAMP ht km %s %s',seas,titcond8);
 end
 if (ifx==3)
nc1 = 9;
nc2 = 10;
nc3 = 11;
nc4 = 12;
tit1ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond9);
tit2ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond10);
tit3ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond11);
tit4ave = sprintf('CHAMP lg10(Ne) in cm-3 %s %s',seas,titcond12);
tit1un = sprintf('CHAMP Te (K) %s %s',seas,titcond9);
tit2un = sprintf('CHAMP Te (K) %s %s',seas,titcond10);
tit3un = sprintf('CHAMP Te (K) %s %s',seas,titcond11);
tit4un = sprintf('CHAMP Te (K) %s %s',seas,titcond12);
tit1ht = sprintf('CHAMP ht km %s %s',seas,titcond9);
tit2ht = sprintf('CHAMP ht km %s %s',seas,titcond10);
tit3ht = sprintf('CHAMP ht km %s %s',seas,titcond11);
tit4ht = sprintf('CHAMP ht km %s %s',seas,titcond12);
 end

n0=num(iseas,nc1)+num(iseas,nc2)+num(iseas,nc3)+num(iseas,nc4);

if (n0>0)
figure
clf
hold on;
if (num0(iseas,nc1)>0)
subplot(2,2,1)
            pcolor(lt24,lat,te(:,:,iseas,nc1)')
            title([tit1un]);
            caxis([ten tex]);
            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');
end
if (num0(iseas,nc2)>0)
subplot(2,2,2)
            pcolor(lt24,lat,te(:,:,iseas,nc2)')
            title([tit2un]);
            caxis([ten tex]);
            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');
end
if (num0(iseas,nc3)>0)
subplot(2,2,3)
            pcolor(lt24,lat,te(:,:,iseas,nc3)')
            title([tit3un]);
            caxis([ten tex]);
            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');
end
if (num0(iseas,nc4)>0)
subplot(2,2,4)
            pcolor(lt24,lat,te(:,:,iseas,nc4)')
            title([tit4un]);
            caxis([ten tex]);
            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');
end
    hold off;	

figure
clf
hold on;
if (num0(iseas,nc1)>0)
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');
end
if (num0(iseas,nc2)>0)
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');
end
if (num0(iseas,nc3)>0)
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');
end
if (num0(iseas,nc4)>0)
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;
end	
	
figure
clf
if (num0(iseas,nc1)>0)
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');
end
if (num0(iseas,nc2)>0)
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');
end
if (num0(iseas,nc3)>0)
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');
end
if (num0(iseas,nc4)>0)
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');
end
hold off;

end	%  if (n0>0)

end	%  for ifx=1:3
end	% for iseas=1:4
