clear all;

% 11/11  Program to calculate the median of the median TEC values from CCMC model outputs
%        from 07325-08020 for 8 longitudes for the Emery et al., Fall AGU 2011 poster 

% 8 begin lons (for 5 deg glon): 25,90,140,175,200,250,285,345
glon=[25,90,140,175,200,250,285,345];
  
% IRI TEC 20,200km ONLY in LT not UT! correct it in median! LT=0:23, UT=LT-glon/15=LT+corrh;
corrh = [-2,-6,-9,-12,-13,-17,-19,-23];
%ich = 1;
ich = 0;
%  12 Dec 2011:  IRI TEC now in UT and up to 500 and 1500 km so have to read 2 files
ird2 = 1;
ird2 = 0;

%  The model longitudes are 25 deg wide like COSMIC - not like MIT GPS TEC
% lon_LLL_1_IRI(CTIPE,TIE-GCM)_2007.325-2008.020.txt where LLL=glon-10

% day 42 (Jan 1, 2008) is bad for 0 UT dres W05 TIMED lb
ndays = 61;
%ndays = 40
n36 = 36;
n36 = 72;


for lon=1:8
glon1 = glon(lon)-10;
glon2 = glon(lon)+15;
if (n36==72)
% NOT QUITE CORRECT (2.5deg), but close enough in integer form
 glon1 = glon(lon)+1;
 glon2 = glon(lon)+4;	
end

% Get name for output file of median TEC values for each glon
%outfile=sprintf('long_%3.3d_07325_08020_1_IRI_1500km.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_IRI.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_CTIPE.dat',glon(lon))
%outfile=sprintf('long_%3.3d_07325_08020_1_TIE-GCM.dat',glon(lon))
% Have double the number of lats (36*2=72)
%outfile=sprintf('long_%3.3d_07325_08020_2_TIE-GCMb.dat',glon(lon))
outfile=sprintf('long_%3.3d_07325_08020_2_TIE-GCMc.dat',glon(lon))
OUT = fopen(outfile,'w');

% Get name of 61-day output file (24h*61d*36glats=52704 lines each)
%infile=sprintf('ccmc_output/lon_%3.3d_1_IRI_2007.325-2008.020.txt',glon1)
%infilt=sprintf('ccmc_output/IRI_TEC_UT/TEC_lon_%3.3d_1_IRI_2007.325-2008.020.1500.txt',glon1)
%infile=sprintf('ccmc_output/lon_%3.3d_1_CTIPE_2007.325-2008.020.txt',glon1)
%infile=sprintf('ccmc_output/lon_%3.3d_1_TIE-GCM_2007.325-2008.020.txt',glon1)
%nfile=sprintf('ccmc_output/lon_%3.3d_2_TIE-GCM_2007.325-2008.026c.txt',glon(lon))
 infile=sprintf('ccmc_output/lon_%3.3d_07325_08026_2_TIE-GCMc.dat',glon(lon))

% Open file
INP = fopen(infile,'r');
if (ird2==1)
 INPT = fopen(infilt,'r');
end

for n=1:ndays
for mh=1:24
nh = mh;
% Go from LT to UT
%if ich==1
% nh = mh + corrh(lon);
% if nh<1
%  nh = nh + 24;
% end
%end
for nlat=1:n36
 line = fgetl(INP);
% year daynum beghr begmin begsec beg_glat beg_glon TEC(TECU) NmF2(cm-5?) hmF2(km)
% dres TIE adds wiav(m/s) on end
 input = sscanf(line,'%f');
if (ird2==1)
 lint = fgetl(INPT);
 inputt = sscanf(lint,'%f');
end

if (n36==36)
% Accumulate in arrays the values of TEC, NmF2, hmF2
tec(lon,n,nh,nlat) = input(8);
nmf2(lon,n,nh,nlat) = input(9);
hmf2(lon,n,nh,nlat) = input(10);
 if (ird2==1)
  tec(lon,n,nh,nlat) = inputt(8);
 end
else
% 2_TIE-GCM is dres so have to average 2 pts for single 5deg lat bin
tecb(nh,nlat) = input(8);
% Must correct nmf2 in 2_TIE-GCM by 1.e+6 before c.txt files
%nmf2b(nh,nlat) = input(9)*1.e+6;
nmf2b(nh,nlat) = input(9);
hmf2b(nh,nlat) = input(10);
wiavb(nh,nlat) = input(11);
end
% Model TEC was never negative
if (input(8) < 0.)
 input
 stop
end
% TIEGCM NmF2 was sometimes negative, but often hmF2 in low-lats <100km so skip
%if (input(9) < 0. | input(10)<100. | input(10)>800.)
%input
%end
end	% nlat=1:36
end	% nh=1:24
% Now combine 2 dres lat bins together (means have 55 to 60 from 53.75 to 58.75)
% Therefore, do NOT combine, but pick an exact latitude like 57.5 (ie, nl+1)!
if (n36==72)
for nh=1:24
for nl=1:2:72
nlat = (nl+1)/2;
%tec(lon,n,nh,nlat) = 0.5*(tecb(nh,nl)+tecb(nh,nl+1));
%wiav(lon,n,nh,nlat) = 0.5*(wiavb(nh,nl)+wiavb(nh,nl+1));
%if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800. & hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1)<800.)
% nmf2(lon,n,nh,nlat) = 0.5*(nmf2b(nh,nl)+nmf2b(nh,nl+1));
% hmf2(lon,n,nh,nlat) = 0.5*(hmf2b(nh,nl)+hmf2b(nh,nl+1));
%else
%  nmf2(lon,n,nh,nlat) = -1.;
%  hmf2(lon,n,nh,nlat) = 0.; 
% if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800.)
%  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl);
%  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl);
% end
% if (hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1) <800.)
%  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+1);
%  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+1);
% end
%end	% if,else
% Choose nl+1 latitude (center of 5deg bin in data)
tec(lon,n,nh,nlat) = tecb(nh,nl+1);
wiav(lon,n,nh,nlat) = wiavb(nh,nl+1);
if (hmf2b(nh,nl+1)>100. & hmf2b(nh,nl+1) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+1);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+1);
else
  nmf2(lon,n,nh,nlat) = -1.;
  hmf2(lon,n,nh,nlat) = 0.; 
 if (hmf2b(nh,nl)>100. & hmf2b(nh,nl) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl);
 end
 if (nl<71)
  if (hmf2b(nh,nl+2)>100. & hmf2b(nh,nl+2) <800.)
  nmf2(lon,n,nh,nlat) = nmf2b(nh,nl+2);
  hmf2(lon,n,nh,nlat) = hmf2b(nh,nl+2);
  end
 end
end		% if,else
end	% nl=1:2:72
end	% nh=1:24
end
end	% n=1:ndays

% Calculate the median TEC over n days for each longitude
for nh=1:24
 for nl=1:36
  clear tecnd;
  clear nmf2nd;
  clear hmf2nd;
  num = 0;
  for n=1:ndays
% revise 0 UT n=42 (1 Jan 2008) for dres - not necessary with c.dat files
%   if n36==72 & nh==1 & n==42
%    tec(lon,n,nh,nl) = tec(lon,n+1,nh,nl);
%    hmf2(lon,n,nh,nl) = 0.;
%    wiav(lon,n,nh,nl) = wiav(lon,n+1,nh,nl);
%   end
   tecnd(n) = tec(lon,n,nh,nl);
   if hmf2(lon,n,nh,nl)>100.
% 1_TIE-GCM bad hmnf2 at: lon=80, lat=+5; lon=140 lat=-5 to +10; lon=165 lat=-5 to +5 and 80-85N(13UT);
%   lon=190 lat=-10 to 0; lon=240 lat=-15 to -5
    num = num + 1;
    nmf2nd(num) = nmf2(lon,n,nh,nl);
    hmf2nd(num) = hmf2(lon,n,nh,nl);
   end    
  end
   medtec(lon,nh,nl) = median(tecnd);
   mednmf2(lon,nh,nl) = median(nmf2nd);
   medhmf2(lon,nh,nl) = median(hmf2nd);
   avetec(lon,nh,nl) = mean(tecnd);
   avenmf2(lon,nh,nl) = mean(nmf2nd);
   avehmf2(lon,nh,nl) = mean(hmf2nd);
   stdtec = 0.;
   stdnmf2 = 0.;
   stdhmf2 = 0.;
   for n=1:ndays
    stdtec = stdtec + (avetec(lon,nh,nl)-tec(lon,n,nh,nl))^2;
   end
   for n=1:num
    stdnmf2 = stdnmf2 + (avenmf2(lon,nh,nl)-nmf2nd(n))^2;
    stdhmf2 = stdhmf2 + (avehmf2(lon,nh,nl)-hmf2nd(n))^2;
   end
    stdtec = sqrt(stdtec/ndays);
    stdnmf2 = sqrt(stdnmf2/num);
    stdhmf2 = sqrt(stdhmf2/num);
if (n36==72)
  for n=1:ndays
   wind(n) = wiav(lon,n,nh,nl);
  end
   medwi(lon,nh,nl) = median(wind);
   avewi(lon,nh,nl) = mean(wind);   
   stdwi = 0.;
  for n=1:ndays
   stdwi = stdwi + (avewi(lon,nh,nl)-wiav(lon,n,nh,nl))^2;
  end
  stdwi = sqrt(stdwi/ndays); 
end
  if num<ndays
nl*5-95,nh-1,num
  end
if (n36==36)
% Print out median tec over ndays
  if nl<19	% allow for minus sign for glats -90 to -5
  fprintf(OUT,'%2.2d %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d\n',nh-1,nl*5-95,glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num);
  else
  fprintf(OUT,'%2.2d  %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d\n',nh-1,nl*5-95,glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num);
  end
else
  if nl<19	% allow for minus sign for glats -90 to -5
  fprintf(OUT,'%2.2d %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d %5.2f %5.2f %5.2f\n',nh-1,nl*5-95,glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num,medwi(lon,nh,nl),avewi(lon,nh,nl),stdwi);
  else
  fprintf(OUT,'%2.2d  %2.2d %3.3d %3.3d %5.2f %5.2f %5.2f %7.0f %7.0f %7.0f %5.1f %5.1f %5.1f %2.2d %2.2d %5.2f %5.2f %5.2f\n',nh-1,nl*5-95,glon1,glon2,medtec(lon,nh,nl),avetec(lon,nh,nl),stdtec,mednmf2(lon,nh,nl),avenmf2(lon,nh,nl),stdnmf2,medhmf2(lon,nh,nl),avehmf2(lon,nh,nl),stdhmf2,ndays,num,medwi(lon,nh,nl),avewi(lon,nh,nl),stdwi);
  end
end
 end	% nl=1:36
end	% nh=1:24   
   
end	% lon=1:8




