function [nt,dayfrac,avwihi,avuihi,avwilo,avuilo,nmax,hmax,tec,nz,nlat,ht,htne,lon,globden,globtec,globden5,globtec5] = get_teclon_nc(Fname,jj,if2p5)

% open .nc file Fname
    nc=netcdf(Fname,'nowrite');

clear height
clear ne

    time=nc{'mtime'}(:);        % (day, ut_hour, ut_min)
    dayfrac=time(:,1)+time(:,2)/24.+time(:,3)/1440.;
    height=nc{'ZG'}(:);
    ilev=nc{'ilev'}(:);
    lat=nc{'lat'}(:);
    lon=nc{'lon'}(:);
% Units are cm/s for ExB components??
    den=nc{'DEN'}(:);
    ne=nc{'NE'}(:); 
   if (if2p5==1)
    dglatlon = 2.5;
    wi=nc{'WI_ExB'}(:);
    ui=nc{'UI_ExB'}(:);
    nmf2=nc{'NMF2'}(:);
    hmf2=nc{'HMF2'}(:);
   else
    dglatlon = 5.;
    wi=nc{'WI_VEL'}(:);
    ui=nc{'UI_VEL'}(:);
   end
    tecarr=nc{'TEC'}(:);
    close(nc)
    
    clear nc

nt = length(dayfrac);
nz = length(ilev);
nlat = length(lat);
nlon = length(lon);
% calculate global density (and tec) and compare with Emmert at 3 altitudes
re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
% Good for dres 2.5x2.5 deg only!
areatot = 0;
for kk=1:5
area5(kk) = 0;
nlat5(kk) = 0;
end
for ii=1:nlat
 da(ii)=re*cos(lat(ii)*dtr)*(dglatlon*dtr)*(re*dglatlon*dtr);
 areatot = areatot+da(ii)*nlon;
 if (lat(ii)<-50.) 
  area5(1) = area5(1) + da(ii)*nlon;
  nlat5(1) = nlat5(1) + 1;
 end
 if (lat(ii)>50.) 
  area5(5) = area5(5) + da(ii)*nlon;
  nlat5(5) = nlat5(5) + 1;
 end
 if (lat(ii)>=-50. & lat(ii)<-25.) 
  area5(2) = area5(2) + da(ii)*nlon;
  nlat5(2) = nlat5(2) + 1;
 end
 if (lat(ii)<=50. & lat(ii)>25.) 
  area5(4) = area5(4) + da(ii)*nlon;
  nlat5(4) = nlat5(4) + 1;
 end
 if (lat(ii)>=-25. & lat(ii)<=25.) 
  area5(3) = area5(3) + da(ii)*nlon;
  nlat5(3) = nlat5(3) + 1;
 end
end
%areatot


% TEC Jan 1, 2008 is missing so globtec,den will be bad - substitute 0UT for 1UT
if time(1,1)==1
 for ii=1:nlat
  for j=1:nlon
    tecarr(1,ii,j) = tecarr(2,ii,j);
   for kk=1:nz
    den(1,kk,ii,j) = den(2,kk,ii,j);
   end
  end
 end
end


htden(1) = 250.e+5;
htden(2) = 400.e+5;
htden(3) = 550.e+5;
for k=1:3
 noden(k) = 0;
 for k5=1:5
  noden5(k,k5)=0;
 end
end

%TEMP
%nt = 4

for it=1:nt
 gtec(it) = 0.;
 for kk=1:5
 gtec5(it,kk) = 0.;
  for k=1:3
  gden(it,k) = 0.;
  ngden(it,k) = 0;
  nlim(it,k) = 0;
  area(it,k) = 0.;
  gden5(it,k,kk) = 0.;
  ngden5(it,k,kk) = 0;
  nlim5k(it,k,kk) = 0;
  area5k(it,k,kk) = 0.;
  end
 end
 for ii=1:nlat
 if (lat(ii)<-50.) 
  k5 = 1;
 end
 if (lat(ii)>50.) 
  k5 = 5;
 end
 if (lat(ii)>=-50. & lat(ii)<-25.) 
  k5 = 2;
 end
 if (lat(ii)<=50. & lat(ii)>25.) 
  k5 = 4;
 end
 if (lat(ii)>=-25. & lat(ii)<=25.) 
  k5 = 3;
 end
temwhi=0.;
temuhi=0.;
nk = 0;
temwlo = 0.;
nklo = 0;
temulo = 0.;

if (if2p5==0)
% Went to nz-1 since Te ub e+35
   for kk=2:nz-1
% Calculate hmax (~290-340km) and nmax
     if (ne(it,kk,ii,jj)>=ne(it,kk-1,ii,jj) & ne(it,kk,ii,jj)>=ne(it,kk+1,ii,jj))
      h(1)=height(it,kk-1,ii,jj)*1.e-5;
      h(2)=height(it,kk,ii,jj)*1.e-5;
      h(3)=height(it,kk+1,ii,jj)*1.e-5;
      p(1)=ne(it,kk-1,ii,jj);
      p(2)=ne(it,kk,ii,jj);
      p(3)=ne(it,kk+1,ii,jj);
      deltax=h(1)*h(1)*(h(2)-h(3))+h(2)*h(2)*(h(3)-h(1))+h(3)*h(3)*(h(1)-h(2));
      atx=p(1)*(h(2)-h(3)) + p(2)*(h(3)-h(1)) + p(3)*(h(1)-h(2));
      ax=atx/deltax;
      btx=h(1)*h(1)*(p(2)-p(3))+h(2)*h(2)*(p(3)-p(1))+h(3)*h(3)*(p(1)-p(2));
      bx=btx/deltax;
      ctx=h(1)*h(1)*(h(2)*p(3)-h(3)*p(2))+h(2)*h(2)*(h(3)*p(1)-h(1)*p(3))+h(3)*h(3)*(h(1)*p(2)-h(2)*p(1));
      cx=ctx/deltax;
      hmf2(it,ii,jj) = -(bx/(2.*ax));
      nmf2(it,ii,jj) = -((bx*bx-4.*ax*cx)/(4.*ax));
     end
   end
end
% Convert from cm-3 to 1.e-12 m-3
%nmax(it,ii) = nmf2(it,ii,jj)*1.e-6;
% Leave nmf2 in cm-3
nmax(it,ii) = nmf2(it,ii,jj);
hmax(it,ii) = hmf2(it,ii,jj);
% Convert from cm-2 to m-2  PUT IN TECU!! or ~e+40!
tec(it,ii) = tecarr(it,ii,jj)*1.e+4;
% Convert from electrons/m2 to TEC units (1e+16 m-2 = 1 TECU)
tec(it,ii) = tec(it,ii)*1.e-16;
% Tn, Un, and composition are on midpts between ZG levels
   for kk=1:nz
% Convert from cm to km
    ht(it,kk,ii) = height(it,kk,ii,jj)*1.e-5;
% Convert from cm-3 to m-3
    htne(it,kk,ii) = ne(it,kk,ii,jj)*1.e+6;
    if (height(it,kk,ii,jj) >= 218.e+5 & height(it,kk,ii,jj) <= 577.e+5) 
     temwhi = temwhi + wi(it,kk,ii,jj);
     temuhi = temuhi + ui(it,kk,ii,jj);
     nk = nk + 1;
    end
    if (height(it,kk,ii,jj) >= 140.e+5 & height(it,kk,ii,jj) <= 170.e+5) 
     temwlo = temwlo + wi(it,kk,ii,jj);
     temulo = temulo + ui(it,kk,ii,jj);
     nklo = nklo + 1;
    end
   end
% Convert from cm/s to m/s
 avwihi(it,ii)=temwhi*0.01/nk;
 avuihi(it,ii)=temuhi*0.01/nk;
 avwilo(it,ii)=temwlo*0.01/nklo;
 avuilo(it,ii)=temulo*0.01/nklo;

% Calculate global density at 250, 400, 550 km
% DEN (g/cm3) is on ZG levels (ie, the interfaces, ilev pts)
  for j=1:nlon
% Calculate global tec (assume always exists and is good)
  gtec(it) = gtec(it) + da(ii)*tec(it,ii);
  gtec5(it,k5) = gtec5(it,k5) + da(ii)*tec(it,ii);

  for k=1:3
   if (height(it,nz,ii,j) < htden(k) | height(it,nz,ii,j) > 1.e+36)
    nlim(it,k) = nlim(it,k) + 1;
    nlim5k(it,k,k5) = nlim5k(it,k,k5) + 1;
   else
   for kk=1:nz-1
     if (height(it,kk,ii,j) <= htden(k) & height(it,kk+1,ii,j) >= htden(k)) 
     area(it,k) = area(it,k) + da(ii);
     lnden = ( (height(it,kk+1,ii,j)-htden(k))*log(den(it,kk,ii,j)) + (htden(k)-height(it,kk,ii,j))*log(den(it,kk+1,ii,j)) ) / (height(it,kk+1,ii,j)-height(it,kk,ii,j));
     gden(it,k) = gden(it,k) + da(ii)*exp(lnden);
     ngden(it,k) = ngden(it,k) + 1;
     area5k(it,k,k5) = area5k(it,k,k5) + da(ii);
     gden5(it,k,k5) = gden5(it,k,k5) + da(ii)*exp(lnden);
     ngden5(it,k,k5) = ngden5(it,k,k5) + 1;
% This confirmed the log extrapolation was correct
%if (it==nt & ngden(it,k) < 5)
%kk,j,ii,ngden(it,k),height(it,kk,ii,j),height(it,kk+1,ii,j),den(it,kk,ii,j),den(it,kk+1,ii,j),lnden,exp(lnden),da(ii)*exp(lnden)
%end
     end
    end		% kk=1:nz-1
    end		% if else
   end		% k=1:3
  end		% j=1:nlon
 end		% ii=1:nlat
% Check to see if area(it,k) is less than areatot
 for k=1:3
  if (nlim(it,k)<nlon*nlat)
   gden(it,k) = gden(it,k)/area(it,k);
  else
   noden(k) = noden(k) + 1;
  end
  for k5=1:5
  if (nlim5k(it,k,k5)<nlon*nlat5(k5))
   gden5(it,k,k5) = gden5(it,k,k5)/area5k(it,k,k5);
  else
   noden5(k,k5) = noden5(k,k5) + 1;
  end
  end	% for k5=1:5
%it,k,ngden(it,k),nlim(it,k),gden(it,k),area(it,k)/areatot
 end		% for k=1:3

% Divide global TEC by total area
 gtec(it) = gtec(it)/areatot;
  for k5=1:5
   gtec5(it,k5) = gtec5(it,k5)/area5(k5);
  end	% for k5=1:5
end		% it=1:nt

% Convert from g/cm3 to kg/m3 1.e-3kg/g*1.e6cm3/m3 = 1.e+3 and take log10 
% 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.4 -12.1 NA]  (like MSIS)
for k=1:3
 globden(k) = -99.;
 if (noden(k)<nt)
  globden(k) = log10(mean(gden(:,k))*1.e+3);
 end
 for k5=1:5
 globden5(k,k5) = -99.;
 if (noden5(k,k5)<nt)
  globden5(k,k5) = log10(mean(gden5(:,k,k5))*1.e+3);
 end
end
end

% Get mean of daily global TECs 
globtec = mean(gtec);
for k5=1:5
 globtec5(k5) = mean(gtec5(:,k5));
end

%noden(1),noden(2),noden(3),nlon*nlat,mean(gden(:,1)),mean(gden(:,2)),mean(gden(:,3))
%globden(1),globden(2),globden(3)

    clear time
%    stop

