%%
%% This program is written for comparing the TIEGCM Wi and JRO VperpN (Viz) data 
%% BA Emery, March 2011
clear all; close all

% Can choose fast option of reading (ifrdtie=1) the 20-min files for JRO and the 7 TIE runs for Sep 2005
%  If ifrdtie = 0 (~20min), then will print out values for 20-min JRO, and TIE, including crb, and different crit1,2
ifrdtie = 0;

 year1 = 2005;
 day1 = 256;
 day2 = 266;
dayofmo1 = 13;
nf1 = 1;
nf2 = 6;

% 20 min saves
 ntmpda = 24*3;
lenda = day2 - day1;
nobs = ntmpda*lenda + 1

% Number of 20 min points to look at the median absolute deviation (mad)
numhr = 4;
npmad = numhr*3+1;
midm1 = (npmad-1)/2;

zero = [0.,0.];
missing = -999.;

fpath = '/hao/aim2/emery/tiegcm_crit/';

% Jicamarca coordinates (11.9S, 76W)  36 lats from -87.5 to +87.5, 72 lons from -180 to +175(?)
jj=15+1;
ii=21+1;

% F reg aves are for 24 ranges from 218-577 km in the JRO Viz data (02Apr ilev=13-24)
% JULIA 150-km drifts are ~140-170 km (ilev=11-13? 200~15?) with kindat=1950 (<2006 and avg >=2006)
ht1 = 218.e+5;
ht2 = 577.e+5;


% Read IMF file used for TIEGCM run of 15-min aves every 5 min (t=0 from t=-5 to -20m)
  Fname=['/hao/aim1/tgcm/data/imf_OMNI_2005001-2005365.nc'];
  nc=netcdf(Fname,'nowrite');
  yrdayfracimf=nc{'date'}(:);
  yrbximf=nc{'bx'}(:);
  yrbyimf=nc{'by'}(:);
  yrbzimf=nc{'bz'}(:);
  yrvswimf=nc{'swvel'}(:);
  close(nc)
  clear nc 
  ntimf=length(yrdayfracimf);
  nt1 = (day1-1)*288 + 1
  nt2 = (day2-1)*288 + 1
  ntot = nt2-nt1 + 1
  for n=nt1:nt2
    dayfracimf(n-nt1+1) = yrdayfracimf(n)-year1*1000;
    bximf(n-nt1+1) = yrbximf(n);
    byimf(n-nt1+1) = yrbyimf(n);
    bzimf(n-nt1+1) = yrbzimf(n);
    btmag(n-nt1+1) = sqrt(yrbximf(n)^2+yrbyimf(n)^2+yrbzimf(n)^2);
    vsw(n-nt1+1) = yrvswimf(n);
  end
for id=1:10
 ni=288*id
 dayfracimf(ni),vsw(ni)
end


%%  Read JRO 5-min obs (some gaps - deleted 13 times on day 256 from jro_viz/jro050913)
jrofil = sprintf ('jro_av_viz_2005');
JRO=fopen(jrofil,'r');
nskip=2271;
% Read over nskip lines
for n=1:nskip
 linjro = fgetl(JRO);
end

% JRO obs start at 0000 UT
for n=1:nobs
 jro_viz_20min(n) = missing;
 jro_nan_20min(n) = NaN;
 time_jro_20min(n) = day1 + (n-1)/ntmpda;
end

for j=1:lenda
 num_obs_20min(j) = 0;
end

% Specifics for Sep 13-22, 2005 (begin time 0200 UT and end time 1550 UT)
nmad1 = 2*3+1+midm1;
nmad2 = ntmpda*lenda - fix((24.-15.6)*3) - midm1 + 1;
jrd = 2674;
iday0 = 256;

minendp = 55;
for i=1:jrd
 linjro = fgetl(JRO);
 jrodat = sscanf(linjro,'%f');
 hmbeg = jrodat(4);
 hbeg = fix(hmbeg/100.);
 minbeg = hmbeg - hbeg*100;
 hmend = jrodat(6);
 hend = fix(hmend/100.);
 minend = hmend - hend*100;
 yday(i) = jrodat(9);
 jro_viz(i) = jrodat(7);
 nummin = fix( (yday(i) - iday0)*24.*60. - 2.49);
% This works if minbeg=minendp or if minbeg function of 20 min 
%  - fails for minendp function of 20 min followed by a gap
 m20 = fix(minbeg/20);
 if (m20*20 == minbeg) 
	im1 = i;
	m20e = fix(minendp/20);
	if (m20e*20 == minendp)
          im1 = i-1;
        end
	i20 = fix((nummin+1)/20);
	jro_viz_20min(i20) = 0.5*(jro_viz(i)+jro_viz(im1));
	jro_nan_20min(i20) = 0.5*(jro_viz(i)+jro_viz(im1));
 end
 minendp = minend;
 i=i+1; 
end
time_jro_20min
jro_viz_20min

% Calculate the median absolute deviation (mad) sliding each 20 min for intervals of npmad pts
%  Best if npmad is odd number.  Usually is over an integer number of hours numhr.  
%  The median, the residual, and the mad are saved for the midpt of npmad pts or at (npmad+1)/2
% function [wimed,wiresid,wimad] = mad_jro(jro_wi,npmad,n1,missing);
iob = 0;
for i=nmad1:nmad2
 icen = i+midm1;
 [wimed,wiresid,wimad] = mad_jro(jro_viz_20min,npmad,i,missing);
 if (wimed > missing)
   iob = iob + 1;
   timemad_obs(iob) = time_jro_20min(icen);
   wimed_obs(iob) = wimed;
   wiresid_obs(iob) = wiresid;
   wimad_obs(iob) = wimad;
   jd = fix(timemad_obs(iob)-iday0+1);
   num_obs_20min(jd) = num_obs_20min(jd) + 1;
 end
end
iob
% Calculate daily and total mads:
mad_obs = 0.;
for i=1:iob
 mad_obs = mad_obs + wimad_obs(i);
end
mad_obs = mad_obs/iob

j1 = 1;
for j=1:lenda
 dmad_obs(j) = 0.;
 j2 = j1 + num_obs_20min(j)-1;
 for jd=j1:j2
   dmad_obs(j) = dmad_obs(j)+wimad_obs(jd);
 end
 dmad_obs(j) = dmad_obs(j)/(j2-j1+1);
 dnum_obs(j) = j2-j1+1;
 j1 = j2 + 1;
end
dmad_obs

%%  Read TIEGCM .nc files

if (ifrdtie == 1)

%% Read JRO 20 min values (from averaging 2 5-min values) from 1 of 2 sources 
% from rdwdjro.pro run where the 20-min aves were printed out
 jro20fil = sprintf ('wi_7_tie_sep05');
JRO20=fopen(jro20fil,'r');
iend = 719;
% Read over 4 lines
for n=1:4
 linjro20 = fgetl(JRO20);
end
for i=1:iend
 linjro20 = fgetl(JRO20);
 jro20dat = sscanf(linjro20,'%f');
 time5d15(i) = jro20dat(3)/24. + day1;
 time5d10(i) = jro20dat(3)/24. + day1;
 time5d05(i) = jro20dat(3)/24. + day1;
 time2d15(i) = jro20dat(3)/24. + day1;
 time7d15(i) = jro20dat(3)/24. + day1;
 time17d5(i) = jro20dat(3)/24. + day1;
 timewcon(i) = jro20dat(3)/24. + day1;
 timejro20(i) = jro20dat(3)/24. + day1;
 jro_av20(i) = jro20dat(5);
 jroav_wconcrit(i) = jro20dat(6);
 jroav_17d5(i) = jro20dat(7);
 jroav_7d15(i) = jro20dat(8);
 jroav_5d15(i) = jro20dat(9);
 jroav_2d15(i) = jro20dat(10);
 jroav_5d10(i) = jro20dat(11);
 jroav_5d05(i) = jro20dat(12);
% TEMP
 crb(i) = 78.;
 mlat1_con(i) = 75.;
 mlat2_con(i) = 60.;
 mlat1_5d05(i) = 73.;
 mlat2_5d05(i) = 68.;
 mlat1_5d10(i) = 73.;
 mlat2_5d10(i) = 63.;
 mlat1_2d15(i) = 75.5;
 mlat2_2d15(i) = 60.5;
 mlat1_5d15(i) = 73.;
 mlat2_5d15(i) = 58.;
 mlat1_7d15(i) = 70.5;
 mlat2_7d15(i) = 55.5;
 mlat1_17d5(i) = 60.5;
 mlat2_17d5(i) = 55.5;
end
itwd7p5 = iend;

else

%% Read TIEGCM runs.  All begin on 0000 with NaN in WI_VEL (which is OK, so read it)
itwd05 = 0;
for num=nf1:15
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dcrit05.swi10_sep05_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_10m_5d05(itwd05+it-it1+1) = jroav_wi(it);
 time5d05_10m(itwd05+it-it1+1) = dayfrac(it);
 mlat1_5d05_10m(itwd05+it-it1+1) = 90. - crit1(it);
 mlat2_5d05_10m(itwd05+it-it1+1) = 90. - crit2(it);
end
itwd05 = itwd05 + nt - (it1-1);
end		% end of for num=1,15 for 10-min saves of crb+5,+5
itwd05
% Now save only every 20 min or every second one where #1=0010 UT
it20 = 0;
for it=2:2:itwd05
 it20 = it20 + 1;
 jroav_5d05(it20) = jroav_10m_5d05(it);
 time5d05(it20) = time5d05_10m(it);
 mlat1_5d05(it20) = mlat1_5d05_10m(it);
 mlat2_5d05(it20) = mlat2_5d05_10m(it);
end
it20
if (it20 < 719)
 for it=it20+1:719
  jroav_5d05(it) = -999.;
  time5d05(it) = time5d05(it20) + (it-it20)*1./(3.*24.);
  mlat1_5d05(it) = mlat1_5d05(it20);
  mlat2_5d05(it) = mlat2_5d05(it20);
 end
end

itwd10 = 0;
for num=nf1:12
it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dcrit10.swi2_sep05_0' str '.nc'];
% Use 1 file made from ncrcat -v; Problem because matlab out of memory!
%it1 = 2;
%Fname=[fpath 'wi_crb5_10.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_2m_5d10(itwd10+it-it1+1) = jroav_wi(it);
 time5d10_2m(itwd10+it-it1+1) = dayfrac(it);
 mlat1_5d10_2m(itwd10+it-it1+1) = 90. - crit1(it);
 mlat2_5d10_2m(itwd10+it-it1+1) = 90. - crit2(it);
end
itwd10 = itwd10 + nt - (it1-1);
end		% end of for num=1,12 for 2-min saves of crb+5,+10 for 7199 times (takes ~10min!)
itwd10
% Now save only every 20 min or every 10th where #1=0002 UT
it20 = 0;
for it=10:10:itwd10
 it20 = it20 + 1;
 jroav_5d10(it20) = jroav_2m_5d10(it);
 time5d10(it20) = time5d10_2m(it);
 mlat1_5d10(it20) = mlat1_5d10_2m(it);
 mlat2_5d10(it20) = mlat2_5d10_2m(it);
end
it20

itwcon = 0;
for num=nf1:36
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_concrit.s_sep05_0' str '.nc'];

    %variable=var(nc); for i=1:length(variable); disp(name(variable{i}));end
% What does index do?    
%    [dayfrac, index]=unique(dayfrac);
%    height=height(index,:,:,:);
%    wi=wi(index,:,:,:);

[nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_wconcrit(itwcon+it-it1+1) = jroav_wi(it);
 timewcon(itwcon+it-it1+1) = dayfrac(it);
 mlat1_con(itwcon+it-it1+1) = 90.-crit1(it);
 mlat2_con(itwcon+it-it1+1) = 90.-crit2(it);
end
itwcon = itwcon + nt - (it1-1);
end		% end of for num=1,37
itwcon

itwd15 = 0;
for num=nf1:36
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dcrit15.s_sep05_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_5d15(itwd15+it-it1+1) = jroav_wi(it);
 time5d15(itwd15+it-it1+1) = dayfrac(it);
 mlat1_5d15(itwd15+it-it1+1) = 90. - crit1(it);
 mlat2_5d15(itwd15+it-it1+1) = 90. - crit2(it);
end
itwd15 = itwd15 + nt - (it1-1);
end		% end of for num=1,37
itwd15

itwd2p5 = 0;
for num=nf1:nf2
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dc1_2p5_mn10_diff15.s_sep05_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_2d15(itwd2p5+it-it1+1) = jroav_wi(it);
 time2d15(itwd2p5+it-it1+1) = dayfrac(it);
 mlat1_2d15(itwd2p5+it-it1+1) = 90. - crit1(it);
 mlat2_2d15(itwd2p5+it-it1+1) = 90. - crit2(it);
end
itwd2p5 = itwd2p5 + nt - (it1-1);
end		% end of for num=1,6 for 20 min saves of crb+2.5,15
itwd2p5

itwd7p5 = 0;
for num=nf1:nf2
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dc1_7p5_mn10_diff15.s_sep05_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_7d15(itwd7p5+it-it1+1) = jroav_wi(it);
 time7d15(itwd7p5+it-it1+1) = dayfrac(it);
 mlat1_7d15(itwd7p5+it-it1+1) = 90. - crit1(it);
 mlat2_7d15(itwd7p5+it-it1+1) = 90. - crit2(it);
 crb(itwd7p5+it-it1+1) = 90. - crit1(it) + 7.5;
 mlat0pot(itwd7p5+it-it1+1) = 90. - crit1(it) + 7.5 - crit1(it) + 7.5;
end
itwd7p5 = itwd7p5 + nt - (it1-1);
end		% end of for num=1,6 for 20 min saves of crb+7.5,15
itwd7p5

itwd17p5 = 0;
for num=nf1:nf2
 it1 = 1;
    %% read tiegcm output 
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 'EMERY.tiegcm_dyncrit.ssep05_crb_17p5_05_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_17d5(itwd17p5+it-it1+1) = jroav_wi(it);
 time17d5(itwd17p5+it-it1+1) = dayfrac(it);
 mlat1_17d5(itwd17p5+it-it1+1) = 90. - crit1(it);
 mlat2_17d5(itwd17p5+it-it1+1) = 90. - crit2(it);
end
itwd17p5 = itwd17p5 + nt - (it1-1);
end		% end of for num=1,6 for 20 min saves of crb+17.5,5 (same crit2 as crb+7.5,15)
itwd17p5

iend = itwd7p5
end		% end of if (ifrdtie==1) else (read TIEGCM .nc files) end

%% Print out values of Wi, crit1,2 from TIEGCM to a file (plus JRO, crb and IMF)
if (ifrdtie == 0)
filewi = sprintf('wi_sep05_jro_tie_b');
WIWR = fopen(filewi,'w');
fprintf(WIWR, 'Wi, crit1,2 values from 20-min TIEGCM \n Plus crb, JRO avs of 2 5-min bins every 20 min, \n  and IMF used (from 15-min aves lagged 5 min) \n');
fprintf(WIWR, 'year daynum uthincr uth   Vsw    Bx   By   Bz    Bt  VizJRO mlatCRB     wcon    mlatCrit1,2      wd05     mlatCrit1,2     wd10    mlatCrit1,2  w2.5,15    mlatCrit1,2     wd15  mlatCrit1,2  w7.5,15    mlatCrit1,2  w17.5,5    mlatCrit1,2 \n');
ida0 = fix(time7d15(1));
% i=1 at 0 UT 256 has -NaN WI_VEL at start
for i=1:itwd7p5 
 iday = fix(time7d15(i));
 utincr = (time7d15(i) - ida0)*24.;
 ndays = fix(utincr/24.);
 utda = utincr - ndays*24;
 i4=(i-1)*4+1;
% IMF start at 0000 UT every 5 min and JRO obs jro_viz_20min(i) start at 0000 UT every 20 min, and jro_20av(i) start at 0000 (first value NaN) every 20 min like TIE 
 fprintf(WIWR, '%4.0d %4.0d %7.2f %7.2f %6.1f %6.2f %6.2f %6.2f %6.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n', year1,iday,utincr,utda,vsw(i4),bximf(i4),byimf(i4),bzimf(i4),btmag(i4),jro_viz_20min(i),crb(i),jroav_wconcrit(i),mlat1_con(i),mlat2_con(i),jroav_5d05(i),mlat1_5d05(i),mlat2_5d05(i),jroav_5d10(i),mlat1_5d10(i),mlat2_5d10(i),jroav_2d15(i),mlat1_2d15(i),mlat2_2d15(i),jroav_5d15(i),mlat1_5d15(i),mlat2_5d15(i),jroav_7d15(i),mlat1_7d15(i),mlat2_7d15(i),jroav_17d5(i),mlat1_17d5(i),mlat2_17d5(i));
end
end


%% calculate rms and print from JRO aves every 20 min (from averaging 2 5-min values)

filerms = sprintf('rms_sep05_jro_tie');
RMS = fopen(filerms,'w');
fprintf(RMS, 'RMS of 20-min TIEGCM values of WI_VEL from avs every 20 min of JRO Viz obs averaged using 2 5-min JRO average F-region values from 218-577 km \n');
fprintf(RMS, 'Sep05  num   wcon    wd15    wd10    wd05 crb+2.5,15 +7.5,15* +17.5,5 \n');

 [rmswcon,i20] = rms_jro(jro_viz_20min,jroav_wconcrit,1,iend,-999);
 [rmswd15,i20] = rms_jro(jro_viz_20min,jroav_5d15,1,iend,-999);
 [rmswd10,i20] = rms_jro(jro_viz_20min,jroav_5d10,1,iend,-999);
 [rmswd05,i20] = rms_jro(jro_viz_20min,jroav_5d05,1,iend,-999);
 [rmswd2p5,i20] = rms_jro(jro_viz_20min,jroav_2d15,1,iend,-999);
 [rmswd7p5,i20] = rms_jro(jro_viz_20min,jroav_7d15,1,iend,-999);
 [rmswd17p5,i20] = rms_jro(jro_viz_20min,jroav_17d5,1,iend,-999);
fprintf(RMS, '13-22 %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n \n', i20,rmswcon,rmswd15,rmswd10,rmswd05,rmswd2p5,rmswd7p5,rmswd17p5);
for j=1:lenda
 n1 = (j-1)*ntmpda + 1;
 n2 = min(iend,n1 + ntmpda - 1);
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_wconcrit,n1,n2,-999);
 dmswcon(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_5d05,n1,n2,-999);
 dmswd05(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_5d10,n1,n2,-999);
 dmswd10(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_5d15,n1,n2,-999);
 dmswd15(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_7d15,n1,n2,-999);
 dmswd7p5(j) = rmsw;
 j20(j) = nrms;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_2d15,n1,n2,-999);
 dmswd2p5(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_17d5,n1,n2,-999);
 dmswd17p5(j) = rmsw;
%fprintf(RMS,'%5.0d %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f  \n', j+dayofmo1-1,j20(j),dmswcon(j),dmswd05(j),dmswd10(j),dmswd2p5(j),dmswd15(j),dmswd7p5(j),dmswd17p5(j));
fprintf(RMS,'%5.0d %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', j+dayofmo1-1,j20(j),dmswcon(j),dmswd15(j),dmswd10(j),dmswd05(j),dmswd2p5(j),dmswd7p5(j),dmswd17p5(j));
end


% Calculate the median absolute deviation (mad) sliding each 20 min for intervals of npmad pts
%  Best if npmad is odd number.  Usually is over an integer number of hours numhr.  
%  The median, the residual, and the mad are saved for the midpt of npmad pts or at (npmad+1)/2
% function [wimed,wiresid,wimad] = mad_jro(jro_wi,npmad,n1,missing);
nmad = nobs-npmad+1
nmad = iend-npmad+1
for i=1:nmad
 icen = i+midm1;
 timemad(icen-midm1) = time_jro_20min(icen);
 [wimed_con(icen-midm1),wiresid_con(icen-midm1),wimad_con(icen-midm1)] = mad_jro(jroav_wconcrit,npmad,i,missing);
 [wimed_5d05(icen-midm1),wiresid_5d05(icen-midm1),wimad_5d05(icen-midm1)] = mad_jro(jroav_5d05,npmad,i,missing);
 [wimed_2d15(icen-midm1),wiresid_2d15(icen-midm1),wimad_2d15(icen-midm1)] = mad_jro(jroav_2d15,npmad,i,missing);
 [wimed_5d10(icen-midm1),wiresid_5d10(icen-midm1),wimad_5d10(icen-midm1)] = mad_jro(jroav_5d10,npmad,i,missing);
 [wimed_5d15(icen-midm1),wiresid_5d15(icen-midm1),wimad_5d15(icen-midm1)] = mad_jro(jroav_5d15,npmad,i,missing);
 [wimed_7d15(icen-midm1),wiresid_7d15(icen-midm1),wimad_7d15(icen-midm1)] = mad_jro(jroav_7d15,npmad,i,missing);
 [wimed_17d5(icen-midm1),wiresid_17d5(icen-midm1),wimad_17d5(icen-midm1)] = mad_jro(jroav_17d5,npmad,i,missing);
end
iob,nmad
% Calculate daily and total mads:
mad_5d05 = 0.;
mad_2d15 = 0.;
mad_17d5 = 0.;
mad_con = 0.;
mad_5d10 = 0.;
mad_5d15 = 0.;
mad_7d15 = 0.;
for i=1:nmad
 mad_con = mad_con + wimad_con(i);
 mad_5d05 = mad_5d05 + wimad_5d05(i);
 mad_5d10 = mad_5d10 + wimad_5d10(i);
 mad_5d15 = mad_5d15 + wimad_5d15(i);
 mad_7d15 = mad_7d15 + wimad_7d15(i);
 mad_2d15 = mad_2d15 + wimad_2d15(i);
 mad_17d5 = mad_17d5 + wimad_17d5(i);
end
mad_con = mad_con/nmad
mad_5d10 = mad_5d10/nmad
mad_5d15 = mad_5d15/nmad
mad_7d15 = mad_7d15/nmad
mad_5d05 = mad_5d05/nmad
mad_2d15 = mad_2d15/nmad
mad_17d5 = mad_17d5/nmad
fprintf(RMS, '\n MAD (median absolute deviations) of 20-min JRO Viz and TIEGCM WI_VEL values \n');
fprintf(RMS, 'Sep05  num     obs   num   concrit  crb+5,5 +5,10 +2.5,15 +5,15 +7.5,15 +17.5,15 \n');
fprintf(RMS, '15-18 %4.0d %7.2f %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n \n', iob,mad_obs,nmad,mad_con,mad_5d05,mad_5d10,mad_2d15,mad_5d15,mad_7d15,mad_17d5);
j1 = 1;
for j=1:lenda
 dmad_con(j) = 0.;
 dmad_5d10(j) = 0.;
 dmad_5d15(j) = 0.;
 dmad_7d15(j) = 0.;
 dmad_5d05(j) = 0.;
 dmad_2d15(j) = 0.;
 dmad_17d5(j) = 0.;
 j2 = j1 + ntmpda-1;
 if (j==1)
   j1 = j1 + midm1;
 end
 if (j==lenda)
   j2 = j2 - midm1 + 1
   j2 = iend - npmad
 end
 for jd=j1:j2
   dmad_con(j) = dmad_con(j)+wimad_con(jd-midm1);
   dmad_5d10(j) = dmad_5d10(j)+wimad_5d10(jd-midm1);
   dmad_5d15(j) = dmad_5d15(j)+wimad_5d15(jd-midm1);
   dmad_7d15(j) = dmad_7d15(j)+wimad_7d15(jd-midm1);
   dmad_5d05(j) = dmad_5d05(j)+wimad_5d05(jd-midm1);
   dmad_2d15(j) = dmad_2d15(j)+wimad_2d15(jd-midm1);
   dmad_17d5(j) = dmad_17d5(j)+wimad_17d5(jd-midm1);
 end
 dmad_con(j) = dmad_con(j)/(j2-j1+1);
 dmad_5d10(j) = dmad_5d10(j)/(j2-j1+1);
 dmad_5d15(j) = dmad_5d15(j)/(j2-j1+1);
 dmad_7d15(j) = dmad_7d15(j)/(j2-j1+1);
 dmad_5d05(j) = dmad_5d05(j)/(j2-j1+1);
 dmad_2d15(j) = dmad_2d15(j)/(j2-j1+1);
 dmad_17d5(j) = dmad_17d5(j)/(j2-j1+1);
fprintf(RMS,'%5.0d %4.0d %7.2f %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n', j+dayofmo1-1,dnum_obs(j),dmad_obs(j),j2-j1+1,dmad_con(j),dmad_5d05(j),dmad_5d10(j),dmad_2d15(j),dmad_5d15(j),dmad_7d15(j),dmad_17d5(j));
 j1 = j2 + 1;
end


%% Make plots

% position [left bottom width height]
fsize = 12;
f9 = 10;
wr1c = 0.33;
wr1m = 0.06;
wr1v = 0.01;
pcwr = 0.10;

purple = [0.48, 0.06, 0.89];
dark_green = [0 0.5 0];
red_brown = [0.6 0.2 0];
yel_grn = [0.75 0.75 0];

 y1=-75;
 y2=85;
y95 = y1 + 0.95*(y2-y1);
y90 = y1 + 0.90*(y2-y1);

yi1 = -10;
yi2 = 10;
yi95 = yi1 + 0.95*(yi2-yi1);

ym1 = 50.;
ym2 = 85.;
ym95 = ym1 + 0.95*(ym2-ym1);


%% Do 11 plots:  each day separate, all 10 days
for np=1:11

if (np == 11)
 x1 = day1;
 x2 = day2;
 daytit = sprintf('TIEGCM W05 13-22 Sep 2005');
 plottit = sprintf('jro_13_22sep05');
rmscon = sprintf('%4.2f',rmswcon);
rms5d05 = sprintf('%4.2f',rmswd05);
rms5d10 = sprintf('%4.2f',rmswd10);
rms5d15 = sprintf('%4.2f',rmswd15);
rms2d15 = sprintf('%4.2f',rmswd2p5);
rms7d15 = sprintf('%4.2f',rmswd7p5);
rms17d5 = sprintf('%4.2f',rmswd17p5);
else
x1 = day1 + (np-1);
x2 = x1 + 1;
daytit = sprintf('TIEGCM W05 %2.0d Sep 2005',np+dayofmo1-1);
plottit = sprintf('jro_%2.0dsep05',np+dayofmo1-1);
rmscon = sprintf('%4.2f',dmswcon(np));
rms5d05 = sprintf('%4.2f',dmswd05(np));
rms5d10 = sprintf('%4.2f',dmswd10(np));
rms5d15 = sprintf('%4.2f',dmswd15(np));
rms2d15 = sprintf('%4.2f',dmswd2p5(np));
rms7d15 = sprintf('%4.2f',dmswd7p5(np));
rms17d5 = sprintf('%4.2f',dmswd17p5(np));
end
xtot = [x1,x2];

figure;
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

subplot('position',[0.08,0.720,0.87,0.240])	% 1/3
hold on;
title(daytit,'FontSize',fsize);
plot (dayfracimf,byimf,'g');
plot (dayfracimf,bzimf,'k');
plot (dayfracimf,vsw/100.,'c');
plot(xtot,zero,'--k');
text (x1+0.01*(x2-x1),yi95,'Bz','Color','k','FontSize',f9);
text (x1+0.05*(x2-x1),yi95,'By','Color','g','FontSize',f9);
text (x1+0.09*(x2-x1),yi95,'Vsw/100','Color','c','FontSize',f9);
ylabel('IMF (nT), Vsw/100 (km/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yi1,yi2],'FontSize',fsize,'box','on');

subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (time7d15,crb,'Color','k','LineWidth',3.0);
text (x1+0.01*(x2-x1),ym95,'CRB','Color','k','FontSize',f9);
plot (time7d15,mlat0pot,'Color','k','LineWidth',2.0);
text (x1+wr1m*(x2-x1),ym95,'0Pot','Color','k','FontSize',f9);
plot (timewcon,mlat1_con,'Color','b','LineStyle','--');
plot (timewcon,mlat2_con,'b');
text (x1+(wr1m+pcwr)*(x2-x1),ym95,'concrit','Color','b','FontSize',f9);
plot (time5d05,mlat1_5d05,'Color',purple,'LineStyle','--');
plot (time5d05,mlat2_5d05,'Color',purple);
text (x1+wr1m*(x2-x1),ym95,'crb5,05','Color',purple,'FontSize',f9);
plot (time5d10,mlat1_5d10,'Color',dark_green,'LineStyle','--');
plot (time5d10,mlat2_5d10,'Color',dark_green);
text (x1+(wr1m+2*pcwr)*(x2-x1),ym95,'crb5,10','Color',dark_green,'FontSize',f9);
plot (time2d15,mlat1_2d15,'Color','g','LineStyle','--');
plot (time2d15,mlat2_2d15,'g');
text (x1+(wr1m+3*pcwr)*(x2-x1),ym95,'crb2.5,15','Color','g','FontSize',f9);
plot (time17d5,mlat1_17d5,'Color',red_brown,'LineStyle','--');
plot (time17d5,mlat2_17d5,'Color',red_brown);
text (x1+(wr1m+6*pcwr)*(x2-x1),ym95,'crb17.5,5','Color',red_brown,'FontSize',f9);
plot (time5d15,mlat1_5d15,'Color','r','LineStyle','--');
plot (time5d15,mlat2_5d15,'r');
text (x1+(wr1m+4*pcwr)*(x2-x1),ym95,'crb5,15','Color','r','FontSize',f9);
plot (time7d15,mlat1_7d15,'Color',yel_grn,'LineStyle','--');
plot (time7d15,mlat2_7d15,'Color',yel_grn);
text (x1+(wr1m+5*pcwr)*(x2-x1),ym95,'crb7.5,15','Color',yel_grn,'FontSize',f9);
ylabel('mlat of CRB,crit1,crit2', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[ym1,ym2],'FontSize',fsize,'box','on');

subplot('position',[0.08,0.050,0.87,0.430])	% 3/3
hold on;
plot (yday,jro_viz,'k');
text (x1+wr1v*(x2-x1),y95,'Obs','Color','k','FontSize',f9);
text (x1+wr1v*(x2-x1),y90,'RMS m/s','Color','k','FontSize',f9);
plot (timewcon,jroav_wconcrit,'b');
text (x1+(wr1v+2*pcwr)*(x2-x1),y95,'concrit','Color','b','FontSize',f9);
text (x1+(wr1v+2*pcwr)*(x2-x1),y90,rmscon,'Color','b','FontSize',f9);
plot (time5d05,jroav_5d05,'Color',purple);
text (x1+(wr1v+pcwr)*(x2-x1),y95,'crb5,05','Color',purple,'FontSize',f9);
text (x1+(wr1v+pcwr)*(x2-x1),y90,rms5d05,'Color',purple,'FontSize',f9);
plot (time5d10,jroav_5d10,'Color',dark_green);
text (x1+(wr1v+3*pcwr)*(x2-x1),y95,'crb5,10','Color',dark_green,'FontSize',f9);
text (x1+(wr1v+3*pcwr)*(x2-x1),y90,rms5d10,'Color',dark_green,'FontSize',f9);
plot (time2d15,jroav_2d15,'g');
text (x1+(wr1v+4*pcwr)*(x2-x1),y95,'crb2.5,15','Color','g','FontSize',f9);
text (x1+(wr1v+4*pcwr)*(x2-x1),y90,rms2d15,'Color','g','FontSize',f9);
plot (time17d5,jroav_17d5,'Color',red_brown);
text (x1+(wr1v+7*pcwr)*(x2-x1),y95,'crb17.5,5','Color',red_brown,'FontSize',f9);
text (x1+(wr1v+7*pcwr)*(x2-x1),y90,rms17d5,'Color',red_brown,'FontSize',f9);
plot (time5d15,jroav_5d15,'Color','r');
text (x1+(wr1v+5*pcwr)*(x2-x1),y95,'crb5,15','Color','r','FontSize',f9);
text (x1+(wr1v+5*pcwr)*(x2-x1),y90,rms5d15,'Color','r','FontSize',f9);
plot (time7d15,jroav_7d15,'Color',yel_grn);
text (x1+(wr1v+6*pcwr)*(x2-x1),y95,'crb7.5,15','Color',yel_grn,'FontSize',f9);
text (x1+(wr1v+6*pcwr)*(x2-x1),y90,rms7d15,'Color',yel_grn,'FontSize',f9);
plot(xtot,zero,'--k');
ylabel('Wi (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

%set(gcf,'color',[1 1 1],'inverthardcopy','off');
%set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
%saveas(gcf,plottit,'tif')
saveas(gcf,plottit,'epsc')
%saveas(gcf,plottit,'jpg')



%% Plot mad, residuals, and median over npmad pts
ymadx = 40.;
yr1 = -29.;
yr2 = 29.;

mad5d15 = 'missing';
if (np == lenda+1)
 plotmad = sprintf('jromad_13_22sep05');
 madobs = sprintf('%4.2f',mad_obs);
 madcon = sprintf('%4.2f',mad_con);
 mad5d10 = sprintf('%4.2f',mad_5d10);
 mad5d15 = sprintf('%4.2f',mad_5d15);
 mad7d15 = sprintf('%4.2f',mad_7d15);
 mad5d05 = sprintf('%4.2f',mad_5d05);
 mad2d15 = sprintf('%4.2f',mad_2d15);
 mad17d5 = sprintf('%4.2f',mad_17d5);
else
 plotmad = sprintf('jromad_%2.0dsep05',np+dayofmo1-1);
 madobs = sprintf('%4.2f',dmad_obs(np));
 madcon = sprintf('%4.2f',dmad_con(np));
 mad5d10 = sprintf('%4.2f',dmad_5d10(np));
 mad5d15 = sprintf('%4.2f',dmad_5d15(np));
 mad7d15 = sprintf('%4.2f',dmad_7d15(np));
 mad5d05 = sprintf('%4.2f',dmad_5d05(np));
 mad2d15 = sprintf('%4.2f',dmad_2d15(np));
 mad17d5 = sprintf('%4.2f',dmad_17d5(np));
end

figure;
% Set size to be 8.5x11 inches
set(gcf,'position',[142    99   647   822])
H = gcf;

subplot('position',[0.08,0.720,0.87,0.240])	% 1/3
hold on;
title(daytit,'FontSize',fsize);
plot (timemad_obs,wimad_obs,'k');
plot (timemad,wimad_con,'b');
plot (timemad,wimad_5d10,'Color',dark_green);
plot (timemad,wimad_5d15,'r');
plot (timemad,wimad_7d15,'Color',yel_grn);
plot (timemad,wimad_5d05,'Color',purple);
plot (timemad,wimad_2d15,'g');
plot (timemad,wimad_17d5,'Color',red_brown);
ylabel('median abs deviation (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[0.,ymadx],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (timemad_obs,wiresid_obs,'k');
plot (timemad,wiresid_con,'b');
plot (timemad,wiresid_5d10,'Color',dark_green);
plot (timemad,wiresid_5d15,'r');
plot (timemad,wiresid_7d15,'Color',yel_grn);
plot (timemad,wiresid_5d05,'Color',purple);
plot (timemad,wiresid_2d15,'g');
plot (timemad,wiresid_17d5,'Color',red_brown);
plot(xtot,zero,'--k');
ylabel('Residual (m/s)', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yr1,yr2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

subplot('position',[0.08,0.050,0.87,0.430])	% 3/3
hold on;
plot (time_jro_20min,jro_nan_20min,':k');
plot (timemad_obs,wimed_obs,'k');
plot (timemad,wimed_con,'b');
plot (timemad,wimed_5d10,'Color',dark_green);
plot (timemad,wimed_5d15,'r');
plot (timemad,wimed_7d15,'Color',yel_grn);
plot (timemad,wimed_5d05,'Color',purple);
plot (timemad,wimed_2d15,'g');
plot (timemad,wimed_17d5,'Color',red_brown);
text (x1+0.85*(x2-x1),y95,'....20-min Obs','Color','k','FontSize',f9);
text (x1+wr1v*(x2-x1),y95,'Obs','Color','k','FontSize',f9);
text (x1+wr1v*(x2-x1),y90,madobs,'Color','k','FontSize',f9);
text (x1+(wr1v+pcwr)*(x2-x1),y95,'concrit','Color','b','FontSize',f9);
text (x1+(wr1v+pcwr)*(x2-x1),y90,madcon,'Color','b','FontSize',f9);
text (x1+(wr1v+3*pcwr)*(x2-x1),y95,'crb5,10','Color',dark_green,'FontSize',f9);
text (x1+(wr1v+3*pcwr)*(x2-x1),y90,mad5d10,'Color',dark_green,'FontSize',f9);
text (x1+(wr1v+5*pcwr)*(x2-x1),y95,'crb5,15','Color','r','FontSize',f9);
text (x1+(wr1v+5*pcwr)*(x2-x1),y90,mad5d15,'Color','r','FontSize',f9);
text (x1+(wr1v+6*pcwr)*(x2-x1),y95,'crb7.5,15','Color',yel_grn,'FontSize',f9);
text (x1+(wr1v+6*pcwr)*(x2-x1),y90,mad7d15,'Color',yel_grn,'FontSize',f9);
text (x1+(wr1v+2*pcwr)*(x2-x1),y95,'crb5,05','Color',purple,'FontSize',f9);
text (x1+(wr1v+2*pcwr)*(x2-x1),y90,mad5d05,'Color',purple,'FontSize',f9);
text (x1+(wr1v+4*pcwr)*(x2-x1),y95,'crb2,15','Color','g','FontSize',f9);
text (x1+(wr1v+4*pcwr)*(x2-x1),y90,mad2d15,'Color','g','FontSize',f9);
text (x1+(wr1v+7*pcwr)*(x2-x1),y95,'crb17.5,5','Color',red_brown,'FontSize',f9);
text (x1+(wr1v+7*pcwr)*(x2-x1),y90,mad17d5,'Color',red_brown,'FontSize',f9);
plot(xtot,zero,'--k');
ytit = sprintf ('%d hour median Wi (m/s)',numhr);
ylabel (ytit,'FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'FontSize',fsize,'box','on');
hold off;

set(gcf,'PaperPosition',[0.5, 0.5, 7.5, 10.5]);
set(gcf,'PaperPositionmode', 'auto');
saveas(gcf,plotmad,'epsc')

end	% end np=1:lenda+1 of plots

