% /hao/aim3/emery/cetic_climatology/calc_glob5_val.m
% 11/12 bae Calculate global value (den400, TEC, etc) from glat/glon grid UAM and TIME-AMIE

% ifuam = 0,4,5 for TIME-AMIE or UAM4,5
ifuam=4

% if2tec = 0,1 for TEC2000 and dens400 or for TEC500 and TECgps (21,000km) for UAM4
if2tec = 1


% nd = # of days 
nd=61;


if (ifuam>0)
% UAM in mlat and mlon 59x25 for dens400 and 87x25 for TEC2000
dglon = 15.;
lon=[0.:dglon:355.];
nlon=length(lon)
nlonp=nlon+1;
nhlatd=[90. 84.85 80.92 77.82 75.23 72.96 70.87 68.86 66.86 64.84 62.77 60.64 58.44 56.13 53.70 51.14 48.42 45.54 42.48 39.26 35.87 32.34 28.68 24.92 21.07 17.17 13.23 9.27 5.29 0.];
m = 60;
for n=1:29
 m = m - 1;
 latd(m) = -nhlatd(n);
 latd(n) = nhlatd(n);
 if (n>1)
  dlatd(n) = 0.5*(nhlatd(n-1)-nhlatd(n+1));
  dlatd(m) = dlatd(n);
 end
end
 latd(30) = 0.;
 nlatd=length(latd)
 dlatd(30) = latd(29);
 dlatd(1) = nhlatd(1)-nhlatd(2);
 dlatd(nlatd) = dlatd(1);

nhlatt=[90. 86.54 83.44 80.66 78.14 75.83 73.66 71.60 69.59 67.59 65.59 63.59 61.59 59.59 57.59 55.59 53.59 51.59 49.59 47.59 45.59 43.59 41.59 39.59 37.59 35.59 33.59 31.59 29.59 27.59 25.59 23.59 21.59 19.59 17.59 15.59 13.59 11.59 9.59 7.59 5.59 3.59 1.59 0.];
m = 88;
for n=1:43
 m = m - 1;
 latt(m) = -nhlatt(n);
 latt(n) = nhlatt(n);
 if (n>1)
  dlatt(n) = 0.5*(nhlatt(n-1)-nhlatt(n+1));
  dlatt(m) = dlatt(n);
 end
end
 latt(44) = 0.;
 nlatt=length(latt)
 dlatt(44) = latt(43);
 dlatt(1) = nhlatt(1)-nhlatt(2);
 dlatt(nlatt) = dlatt(1);

 if (ifuam==4)
  if (if2tec==0)
   filtec=sprintf('/hao/aim3/emery/ccmc_uam/v2/LongRun61_4_UAM/UAM4_07325_08020_TEC2000')
   filden=sprintf('/hao/aim3/emery/ccmc_uam/v2/LongRun61_4_UAM/UAM4_07325_08020_dens400')
   OUTG = fopen('daily_den_tec_UAM4','w');
  else
% Change nlatd to nlatt since are reading 2 TEC files
   nlatd = nlatt
   filden=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM/UAM4_07325_08020_TEC500')
   filtec=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_4_UAM/UAM4_07325_08020_TECgps')
   OUTG = fopen('daily_tec500_tecgps_UAM4','w');
  end
 else
  filtec=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_5_UAM/UAM5_07325_08020_TEC2000')
  filden=sprintf('/hao/aim3/emery/ccmc_uam/LongRun61_5_UAM/UAM5_07325_08020_dens400')
  OUTG = fopen('daily_den_tec_UAM5','w');
 end
  DTEC = fopen(filtec,'r');
  DDEN = fopen(filden,'r');

% Read 24 times to get daily values
%nd=1
 for md=1:nd
  tecd(1:nlon,1:nlatt,md) = 0.;
  dend(1:nlon,1:nlatd,md) = 0.;
  for nh=1:24
% Skip over 8 headers
   for ns=1:8
    rdhd = fgetl(DTEC);
    rdhd = fgetl(DDEN);
   end
   for nla=1:nlatt
%  Read lat and values for nlon (25) mlons
    lintec = fgetl(DTEC);
    intec = sscanf(lintec,'%f');
    for j=2:nlonp
     tech(j-1)=(10.^intec(j))*1.e-16;	% Convert log10(m-2) to TECU
    end		%     for j=2:nlonp
    for j=1:nlon
     tecd(j,nla,md) = tecd(j,nla,md) + tech(j);	
    end
   end	% for nla=1:nlatt
   for nla=1:nlatd
%  Read lat and values for nlon (25) mlons
    linden = fgetl(DDEN);
    inden = sscanf(linden,'%f');
    for j=2:nlonp
     if (if2tec==0)
      denh(j-1)=inden(j);
     else
      denh(j-1)=(10.^inden(j))*1.e-16;	% Convert log10(m-2) to TECU
     end
    end		%     for j=2:nlonp
    for j=1:nlon
     dend(j,nla,md) = dend(j,nla,md) + denh(j);
    end
   end	% for nla=1:nlatd
  end	% for nh=1:24
   for nla=1:nlatt
    for j=1:nlon
     tecd(j,nla,md) = tecd(j,nla,md)/24.;
    end		%     for j=2:nlonp
   end	% for nla=1:nlatt
   for nla=1:nlatd
    for j=1:nlon
     dend(j,nla,md) = dend(j,nla,md)/24.;
    end		%     for j=2:nlonp
   end	% for nla=1:nlatd
 end	% for md=1:nd

if (if2tec==0)
 fprintf(OUTG,'daynum(2007) 400km_log10(model_den(kg/m3)), tec global SP,SH,eq,NH,NP \n');
else
 fprintf(OUTG,'daynum(2007) tec500, tecgps(21,000km) global SP,SH,eq,NH,NP \n');
end

% Do DEN then TEC
for mpar=1:2
 if (mpar==1 & if2tec==0)
  lat=latd;
  dlat=dlatd;
 else
  lat=latt;
  dlat=dlatt;
 end
 nlat=length(lat)

% calculate global density (and tec) and compare with Emmert at 3 altitudes
re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
% Good for regular grids 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)*(dglon*dtr)*(re*dlat(ii)*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,area5,nlat5
%areat5=area5(1)+area5(2)+area5(3)+area5(4)+area5(5)

if (mpar==1) 
for it=1:nd
 gden(it) = 0.;
%  Assume no NaN for UAM densities
 area(it) = areatot;
 for kk=1:5
  area5k(it,kk) = area5(kk);
  gden5(it,kk) = 0.;
 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

% 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)
     gden(it) = gden(it) + da(ii)*dend(j,ii,it);
     gden5(it,k5) = gden5(it,k5) + da(ii)*dend(j,ii,it);
  end		% j=1:nlon
 end		% ii=1:nlat
end		% it=1:nd
%gden5
else
for it=1:nd
 gtec(it) = 0.;
 for kk=1:5
  gtec5(it,kk) = 0.;
 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

  for j=1:nlon
% Calculate global tec (assume always exists and is good)
  gtec(it) = gtec(it) + da(ii)*tecd(j,ii,it);
  gtec5(it,k5) = gtec5(it,k5) + da(ii)*tecd(j,ii,it);
  end		% j=1:nlon
 end		% ii=1:nlat
end		% it=1:nd
%gtec5
end	% if (mpar==1)  else

clear lat;
end 	% for mpar=1:2

else

% TIME-AMIE in glat and glon 36x72
 dglat = 5.;
 dglon = 5.;
 lat=[-87.5:dglat:87.5];
 lon=[-180.:dglon:175.]; 
 nlat=length(lat);
 nlon=length(lon)
 dattec = dlmread('astra_amie_time/61day_globalTEC.dat','',1,0);
 datden = dlmread('astra_amie_time/61day_globalRHO.dat','',2,0);
  OUTG = fopen('daily_den_tec_TIMA','w');
% Store 400 km den and TEC daily in glat/glon grid
 for md=1:nd
  for nla=1:nlat
   for nlo=1:nlon
    injl = nlo + (nla-1)*nlon + (md-1)*nlon*nlat;
    tecd(nlo,nla,md) = dattec(injl,5);
    dend(nlo,nla,md) = datden(injl,6); 
   end
  end
 end  

fprintf(OUTG,'daynum(2007) 400km_log10(model_den(kg/m3)), tec global SP,SH,eq,NH,NP \n');


% calculate global density (and tec) and compare with Emmert at 3 altitudes
re=6371.e5;   %  earth radius in cm
dtr = pi/180.;
% Good for regular grids 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)*(dglon*dtr)*(re*dglat*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




for it=1:nd
 gtec(it) = 0.;
 gden(it) = 0.;
 ngden(it) = 0;
 area(it) = 0.;
 for kk=1:5
  gtec5(it,kk) = 0.;
  gden5(it,kk) = 0.;
  ngden5(it,kk) = 0;
  area5k(it,kk) = 0.;
 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

% 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)*tecd(j,ii,it);
   gtec5(it,k5) = gtec5(it,k5) + da(ii)*tecd(j,ii,it);
% den400km not always present for TIME-AMIE!
   if (isnan(dend(j,ii,it))==0)
     ngden(it) = ngden(it) + 1;
     area(it) = area(it) + da(ii);
     area5k(it,k5) = area5k(it,k5) + da(ii);
     ngden5(it,k5) = ngden5(it,k5) + 1;
     gden(it) = gden(it) + da(ii)*dend(j,ii,it);   
     gden5(it,k5) = gden5(it,k5) + da(ii)*dend(j,ii,it);
   end
  end		% j=1:nlon
 end		% ii=1:nlat
end		% it=1:nd
 
end	% if ifuam>0 else

% Divide by total area
gtec = gtec/areatot;
avgtec=mean(gtec)
for it=1:nd
 gden(it) = gden(it)/area(it);
end
avgden = mean(gden)
if (if2tec==0)
 lg10avgden = log10(avgden)
 fprintf(OUTG,'%e %6.2f %7.2f mean 61-d 07325-08020 den400+lg10 kg/m3 and TEC2000 TECU \n',avgden,log10(avgden),avgtec);
else
 fprintf(OUTG,'%7.2f %7.2f mean 61-d 07325-08020 TEC 500km and TEC gps (21,000km) TECU \n',avgden,avgtec);
end
for it=1:nd
if (if2tec==0)
% get log10 of den400km
 gden(it) = log10(gden(it));
end
  for k5=1:5
   gtec5(it,k5) = gtec5(it,k5)/area5(k5);
   if (if2tec==0)
    gden5(it,k5) = log10(gden5(it,k5)/area5k(it,k5));
   else
    gden5(it,k5) = gden5(it,k5)/area5k(it,k5);
   end
  end	% for k5=1:5
ida = it + 324;
if (if2tec==0)
 fprintf(OUTG,'%5.1f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',ida+0.5,gden(it),gden5(it,1),gden5(it,2),gden5(it,3),gden5(it,4),gden5(it,5),gtec(it),gtec5(it,1),gtec5(it,2),gtec5(it,3),gtec5(it,4),gtec5(it,5));
else
 fprintf(OUTG,'%5.1f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n',ida+0.5,gden(it),gden5(it,1),gden5(it,2),gden5(it,3),gden5(it,4),gden5(it,5),gtec(it),gtec5(it,1),gtec5(it,2),gtec5(it,3),gtec5(it,4),gtec5(it,5));
end
end		% it=1:nd





