function [nt,dayfrac,ut,wiav,qjs,qjn] = get_wi_qj_nc(Fname,ht1,ht2,jj)

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

clear height
clear wi

    time=nc{'mtime'}(:);        % (day, ut_hour, ut_min)
    dayfrac=time(:,1)+time(:,2)/24.+time(:,3)/1440.;
    ut=time(:,2)+time(:,3)/60.;
    ilev=nc{'ilev'}(:);
    lat=nc{'lat'}(:);
    lon=nc{'lon'}(:);
    height=nc{'ZG'}(:);
    wi=nc{'WI_ExB'}(:);
    qj_integ=nc{'QJOULE_INTEG'}(:); 
    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.;
areatot = 0;
dlat = 180./nlat;
for j=1:nlat
 da(j)=re*cos(lat(j)*dtr)*(dlat*dtr)*(re*dlat*dtr);
 areatot = areatot+da(j)*nlon;
end
%areatot

% Calculate hemispheric Joule Heat in GW
for it=1:nt
 qjs(it) = 0.;
 qjn(it) = 0.;
 for j=1:nlat/3
  for ii=1:nlon
     qjs(it) = qjs(it) + da(j)*qj_integ(it,j,ii);
  end		% ii=1:nlon
 end		% j=1:nlat/3
 for j=2*nlat/3:nlat
  for ii=1:nlon
     qjn(it) = qjn(it) + da(j)*qj_integ(it,j,ii);
  end		% ii=1:nlon
 end		% j=2*nlat/3:nlat
% Convert from cm2*ergs/cm2-sec to GW (e-7J/s*ergs/s*e-9GW/W)
     qjs(it) = qjs(it)*1.e-16;
     qjn(it) = qjn(it)*1.e-16;
end		% it=1:nt
qjs,qjn

for it=1:nt
for ii=1:72
tem=0.;
nk = 0;
% Find mean between ht1 and ht2 (in cm)
% Tn, Un, and composition are on midpts between ZG levels, but for a mean, this is not too important
   for kk=1:nz
    if (height(it,kk,jj,ii) >= ht1 & height(it,kk,jj,ii) <= ht2) 
     tem = tem + wi(it,kk,jj,ii);
     nk = nk + 1;
    end
   end
% Convert from cm/s to m/s
 wiav(ii,it) = tem*0.01/nk;
 end
end

    clear time

