%%
%% This program is written for comparing the TIEGCM and RFP FPI data
%% BA Emery, April 2011
clear all; close all

% notie=1 if have no TIE values to read (wd and wded have no TN,UN,VN)
notie = 1
% madskp = 1 means skip get_mad for obs in all directions
madskp = 0

 year1 = 2006;
 day1 = 348;
 day2 = 350;
nf1 = 1;
it1nf1 = 1;
nf2 = 1;

 iday0 = day1 - 1;
% 20 min saves
 ntmpda = 24*3;
lenda = day2 - day1;
nobs = ntmpda*lenda + 1
% Subtract 1 time so start at 0020 UT on first day
nobs = nobs - 1

% Number of 20 min points to look at the median absolute deviation (mad)
numhr = 4;
npmad = numhr*3+1;
% midm1=6 points in 4h of 13 20min obs
midm1 = (npmad-1)/2;
% Set the min number of nmad points small because is sparse data 
%   Therefore, for sparse data, may only find mad of combined data (and set madskp=0)
minnmad = 2;


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

% FPI red-line ranges from 210-300 km 
ht1 = 210.e+5;
ht2 = 300.e+5;

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];

% Read IMF file (TIEGCM run used hourly IMF) of 15-min aves every 5 min (t=0 from t=-5 to -20m)
  Fname=['/hao/aim1/tgcm/data/imf_OMNI_2006001-2006365.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:lenda
 ni=288*id
 dayfracimf(ni),vsw(ni)
end


% Have max of 68 obs per direction for 06347_350 Resolute Bay FPI (68 all dir but Z w 66), 45 for 348-349
nobsx = 68;
nobsx = 45;


% Resolute bay coordinates (74.73N, 94.89W)  36 lats from -87.5 to +87.5, 72 lons from -180 to +175(?)

% Z  74.73N, 94.89W,  j,i=33,18
% W  74.73N, 103.42W, j,i=33,16
% E  74.73N, 86.36W,  j,i=33,20
% N  76.98N, 94.89W,  j,i=34,18
% S  72.48N, 94.89W,  j,i=33,18 (same as Z)
ndir = 5;

dirs={'Z';'W';'E';'N';'S'};
jin = [33 33  33  34  33];
iin = [18 16  20  18  18];
color9={'k';yel_grn;'g';'b';'r'};

minmad1 = nobs;
maxmad2 = 0;

% Have RFP obs start at 0020 UT for the 20 min values to be the same as the TIEGCM
n9=ndir+1;
for n=1:nobs
 tn_obs_nan_20min(n,n9) = 0.;
% Start at 0000 UT
%time_obs_20min(n) = day1 + (n-1)/ntmpda;
% Start at 0020 UT (20 min is .0139 frac_days)
 time_obs_20min(n) = day1 + n/ntmpda;
 n9_obs_20min(n) = 0;
 octas_20min(n) = NaN;
end

for n9=1:ndir

obsfil = ['/hao/aim2/emery/cetic12/rbay_fpi_' dirs{n9}];

%%  Read RFP obs and assign each time to closest 20min model time
OBS=fopen(obsfil,'r');

% Have RFP obs start at 0020 UT for the 20 min values to be the same as the TIEGCM 
for n=1:nobs
 tn_obs_nan_20min(n,n9) = NaN;
 un_obs_nan_20min(n,n9) = NaN;
 vn_obs_nan_20min(n,n9) = NaN;
end


n = 0;
while (feof(OBS)==0)
 linobs = fgetl(OBS);
 obsdat = sscanf(linobs,'%f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %s');
 yrobs = obsdat(1);
 daynumobs = obsdat(2);
% Do not save if daynumobs < day1 or > day2-1
if (daynumobs >= day1 & daynumobs <= day2-1) 
 utobs = obsdat(4);
 n = n + 1;
 timeobs(n,n9) = daynumobs + utobs/24.;
 tnobs(n,n9) = NaN;
 unobs(n,n9) = NaN;
 vnobs(n,n9) = NaN;
% Put in NaN if gap over 0.2 day
   if (n > 1 & timeobs(n,n9)-timeobs(n-1,n9) > 0.2)
    timeobs(n,n9) = 0.5*(timeobs(n-1,n9)+timeobs(n,n9));
% Set next time
    n = n + 1;
    timeobs(n,n9) = daynumobs + utobs/24.;
    tnobs(n,n9) = NaN;
    unobs(n,n9) = NaN;
    vnobs(n,n9) = NaN;
   end
% Do not use Tn data if error bar > 100 K
 if (obsdat(5) > missing & obsdat(6) < 100.)
  tnobs(n,n9) = obsdat(5);
 end
 if (obsdat(7) > missing)
  unobs(n,n9) = obsdat(7);
 end
 if (obsdat(9) > missing)
  vnobs(n,n9) = obsdat(9);
 end
% Find nearest 20 min time 
 ntfrac = fix( ((timeobs(n,n9)-time_obs_20min(1))*1440. + 9.99)/20. ) + 1;
if (ntfrac >= 1)
 octas_20min(ntfrac) = obsdat(13);
 tn_obs_nan_20min(ntfrac,n9) = tnobs(n,n9);
 un_obs_nan_20min(ntfrac,n9) = unobs(n,n9);
 vn_obs_nan_20min(ntfrac,n9) = vnobs(n,n9); 
 if ~isnan(tnobs(n,n9))
  n9_obs_20min(ntfrac) = n9_obs_20min(ntfrac) + 1;
  tn_obs_nan_20min(ntfrac,ndir+1) = tn_obs_nan_20min(ntfrac,ndir+1)+tnobs(n,n9);
 end
end
end	% end of if daynumobs between day1 and day2-1
end	% end of while
n


nmad1 = fix((timeobs(1,n9)-time_obs_20min(1))*1440./20.) + 1;
nmad2 = min(nobs-npmad+1,fix((timeobs(n,n9)-time_obs_20min(1))*1440./20.) + 1);
minmad1 = min(minmad1,nmad1);
maxmad2 = max(maxmad2,nmad2);


if (n < nobsx)
% Add extra times
for n=n+1:nobsx
timeobs(n,n9) = timeobs(n-1,n9);
tnobs(n,n9) = NaN;
unobs(n,n9) = NaN;
vnobs(n,n9) = NaN;
end
end

if (madskp == 0)
% 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
iob = 0;
nmad = 0;
for i=nmad1:nmad2
 icen = i+midm1;
 [obsmed,obsresid,obsmad] = get_mad(tn_obs_nan_20min(:,n9),npmad,i,minnmad);
% if (obsmed ~= NaN)	% This fails
 if ~isnan(obsmed)
% Put in NaN if gap over 0.2 day
   if (iob > 0 & time_obs_20min(icen)-timemad_obs(iob,n9) > 0.2)
    iob = iob + 1;
    timemad_obs(iob,n9) = 0.5*(time_obs_20min(icen)+timemad_obs(iob-1,n9));
    tnmed_obs(iob,n9) = NaN;
    tnresid_obs(iob,n9) = NaN;
    tnmad_obs(iob,n9) = NaN;
   end
   iob = iob + 1;
   nmad = nmad + 1;
   timemad_obs(iob,n9) = time_obs_20min(icen);
   tnmed_obs(iob,n9) = obsmed;
   tnresid_obs(iob,n9) = obsresid;
   tnmad_obs(iob,n9) = obsmad;
 end
end
nmad
% Calculate daily and total mads:
mad_obs(n9) = 0.;
if (iob > 0)
 for i=1:iob
  if ~isnan(tnmad_obs(i,n9))
   mad_obs(n9) = mad_obs(n9) + tnmad_obs(i,n9);
  end
 end
 mad_obs(n9) = mad_obs(n9)/nmad;
end
nmad_obs(n9) = nmad;

if (n9==ndir)
% Calculate 4h median average of Tn over ndir positions
for it=1,nobs
tn_med_obs_20min(it) = NaN;
nmed = 0;
 for n9=1:ndir
  for n=1:nmad_obs(n9)
    if (abs(timemad_obs(n,n9)-time_obs_20min(it)) < 0.01) 
      nmed = nmed + 1
      tn_med_obs_20min(it) = tn_med_obs_20min(it) + tnmed(n,n9);
    end
  end
 end
if (nmed>0)
 tn_med_obs_20min(it) = tn_med_obs_20min(it)/nmed; 
end
end
end	% n9==ndir

end	% madskp

end	% end n9=1:ndir

% Calculate Tn obs to nearest 20 min from all ndir locations if present, or set to NaN if absent
for it=1:nobs
 if (n9_obs_20min(it) > 0)
  tn_obs_nan_20min(it,ndir+1) = tn_obs_nan_20min(it,ndir+1)/n9_obs_20min(it);
 else
  tn_obs_nan_20min(it,ndir+1) = NaN;
 end
end

% Calculate mad from tn_obs_nan_20min(:,ndir+1) from all obs
n9 = ndir+1;
iob = 0;
nmad = 0;
for i=minmad1:maxmad2
 icen = i+midm1;
 [obsmed,obsresid,obsmad] = get_mad(tn_obs_nan_20min(:,n9),npmad,i,midm1);
 if ~isnan(obsmed)
   if (iob > 0 & time_obs_20min(icen)-timemad_obs(iob,n9) > 0.2)
    iob = iob + 1;
    timemad_obs(iob,n9) = 0.5*(time_obs_20min(icen)+timemad_obs(iob-1,n9));
    tnmed_obs(iob,n9) = NaN;
    tnresid_obs(iob,n9) = NaN;
    tnmad_obs(iob,n9) = NaN;
   end
   iob = iob + 1;
   nmad = nmad + 1;
   timemad_obs(iob,n9) = time_obs_20min(icen);
   tnmed_obs(iob,n9) = obsmed;
   tnresid_obs(iob,n9) = obsresid;
   tnmad_obs(iob,n9) = obsmad;
 end
end
iob
% Calculate total mads:
mad_obs(n9) = 0.;
if (iob > 0)
 for i=1:iob
  if ~isnan(tnmad_obs(i,n9))
   mad_obs(n9) = mad_obs(n9) + tnmad_obs(i,n9);
  end
 end
 mad_obs(n9) = mad_obs(n9)/nmad;
end
nmad_obs(n9) = nmad;
% Extend n9=1:ndir timemad_obs to iob values for (:,ndir+1)
for n9=1:ndir
  for n=nmad_obs(n9)+1:iob
   timemad_obs(n,n9) = NaN;
   tnmed_obs(iob,n9) = NaN;
   tnresid_obs(iob,n9) = NaN;
   tnmad_obs(iob,n9) = NaN;
  end
end
   

if (notie == 0)
fpath = '/hao/ftp/pub/emery/';

%%  Read TIEGCM .nc files starting at 0020 UT every 20 min (so itwd = nobs-1)

itwd = 0;
it1 = 1;
    Fname=[fpath '/wd194/EMERY.tiegcm_dyncrit.s06348_350_001.nc'];
    [nt,dayfrac,tnav,unav,vnav] = get_tuvw_nc(Fname,ht1,ht2,iin,jin,ndir);
for it=it1:nt
 tn_wd(itwd+it-it1+1,ndir+1) = 0.;
for n9=1:ndir
 tn_wd(itwd+it-it1+1,n9) = tnav(it,n9);
 timewd(itwd+it-it1+1,n9) = dayfrac(it);
 un_wd(itwd+it-it1+1,n9) = unav(it,n9);
 vn_wd(itwd+it-it1+1,n9) = vnav(it,n9);
 tn_wd(itwd+it-it1+1,ndir+1) = tn_wd(itwd+it-it1+1,ndir+1) + tnav(it,n9);
end
 tn_wd(itwd+it-it1+1,ndir+1) = tn_wd(itwd+it-it1+1,ndir+1)/ndir;
end
itwd = itwd + nt - (it1-1);
    Fname=[fpath '/wd194/EMERY.tiegcm_dyncrit.s06348_350_002.nc'];
    [nt,dayfrac,tnav,unav,vnav] = get_tuvw_nc(Fname,ht1,ht2,iin,jin,ndir);
for it=it1:nt
 tn_wd(itwd+it-it1+1,ndir+1) = 0.;
for n9=1:ndir
 tn_wd(itwd+it-it1+1,n9) = tnav(it,n9);
 timewd(itwd+it-it1+1,n9) = dayfrac(it);
 un_wd(itwd+it-it1+1,n9) = unav(it,n9);
 vn_wd(itwd+it-it1+1,n9) = vnav(it,n9);
 tn_wd(itwd+it-it1+1,ndir+1) = tn_wd(itwd+it-it1+1,ndir+1) + tnav(it,n9);
end
 tn_wd(itwd+it-it1+1,ndir+1) = tn_wd(itwd+it-it1+1,ndir+1)/ndir;
end
itwd = itwd + nt - (it1-1);


%% calculate rms and print from obs every 20 min 

 filerms = sprintf('rms_dec06_rfp_tie');
RMS = fopen(filerms,'w');
fprintf(RMS, 'RMS of 20-min TIEGCM values of Tn from avs every 20 min of Resolute Bay FPI obs from 210-300 km \n');
fprintf(RMS, '%d    tot     Z       W     E       N     S \n',year1);
iend = min(itwd,nobs);
for n9=1:ndir
 rmswd(n9) = 0.;
 [rmswd(n9),n20(n9)] = get_rms(tn_obs_nan_20min(:,n9),tn_wd(:,n9),1,iend);
end
fprintf(RMS, '%d-%d %4.0d %7.2f %3.0d %7.2f %3.0d %7.2f %3.0d %7.2f %3.0d %7.2f \n \n', day1,day2,n20(1),rmswd(1),n20(2),rmswd(2),n20(3),rmswd(3),n20(4),rmswd(4),n20(5),rmswd(5));

% 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
nmad = itwd-npmad+1
for n9=1:ndir
for i=1:nmad
 icen = i+midm1;
 timemad(icen-midm1) = time_obs_20min(icen);
 [tnmed_wd(icen-midm1,n9),tnresid_wd(icen-midm1,n9),tnmad_wd(icen-midm1,n9)] = get_mad(tn_wd(:,n9),npmad,i,midm1);
end
% Calculate total mads:
mad_wd(n9) = 0.;
for i=1:nmad
 mad_wd(n9) = mad_wd(n9) + tnmad_wd(i,n9);
end
mad_wd(n9) = mad_wd(n9)/nmad;
end	% n9=1:ndir
fprintf(RMS, '\n MAD (median absolute deviations) of 20-min Resolute Bay Relative Tn obs and TIEGCM Tn values \n');
fprintf(RMS, '%     obs   wd1.94 1-5 \n',year1);
fprintf(RMS, '%d-%d %4.0d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f \n \n', day1,day2,nmad,mad_obs(1),nmad,mad_wd(1),mad_obs(2),nmad,mad_wd(2),mad_obs(3),nmad,mad_wd(3),mad_obs(4),nmad,mad_wd(4),mad_obs(5),nmad,mad_wd(5));
% Calculate 4h median average of Tn over ndir positions
for it=1,nobs
tn_med_tie_20min(it) = NaN;
nmed = 0;
 for n9=1:ndir
  for n=1:nmad
    if (abs(timemad(n)-time_obs_20min(it)) < 0.001) 
      nmed = nmed + 1;
      tn_med_tie_20min(it) = tn_tie_obs_20min(it) + tnmed_wd(n,n9);
    end
   end
 end
if (nmed>0)
 tn_med_tie_20min(it) = tn_med_tie_20min(it)/nmed; 
end
end

end	% Skip reading TIEGCM if notie=1

%% Make plots

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

yi1 = -20;
yi2 = 20;
yi95 = yi1 + 0.92*(yi2-yi1);

 y1=500;
 y2 = 1499;
y95 = y1 + 0.92*(y2-y1);
y90 = y1 + 0.85*(y2-y1);
y85 = y1 + 0.78*(y2-y1);

yu1 = -500;
yu2 = 499;
yu95 = yu1 + 0.92*(yu2-yu1);
yu90 = yu1 + 0.85*(yu2-yu1);
yu85 = yu1 + 0.78*(yu2-yu1);


 x1 = day1;
 x2 = day2;
 daytit = sprintf('Resolute Bay (74.73N, 94.89W) 14-15 Dec 2006');
 plottit = sprintf('rfp_06348_9');
 if (notie==1)
  plottit = sprintf('rfp_data_06348_9');
 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,bximf,'k');
plot (dayfracimf,byimf,'g');
plot (dayfracimf,bzimf,'r');
plot (dayfracimf,vsw/100.,'c');
plot(xtot,zero,'--k');
text (x1+0.01*(x2-x1),yi95,'Bx','Color','k','FontSize',f9);
text (x1+0.05*(x2-x1),yi95,'By','Color','g','FontSize',f9);
text (x1+0.09*(x2-x1),yi95,'Bz','Color','r','FontSize',f9);
text (x1+0.13*(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],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

% Obs in order N NW W SW S Z SE E NE for colors/lines: r m-- yel_grn c-- b k c g m
subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (time_obs_20min,tn_obs_nan_20min(:,ndir+1),'k','LineWidth',2.0);
text (x1+0.82*(x2-x1),y95,'All obs (err<100K)','Color','k','FontSize',f9,'LineWidth',2.0);
for n9=1:ndir
 plot (timeobs(:,n9),tnobs(:,n9),'Color',color9{n9});
 text (x1+(0.01+(n9-1)*0.09)*(x2-x1),y95,dirs{n9},'Color',color9{n9},'FontSize',f9);
end
if (notie ~=1)
text (x1+0.82*(x2-x1),y90,'RMS','Color','k','FontSize',f9);
text (x1+0.82*(x2-x1),y85,'--av TIEGCM ','Color',dark_green,'FontSize',f9);
plot (timewd,tn_wd(:,ndir+1),'Color',dark_green,'LineStyle','--','LineWidth',2.0);
for n9=1:ndir
 rms = sprintf('%4.2f',rmswd(n9));
  plot (timewd,tn_wd(:,n9),'Color',color9{n9});
 text (x1+(0.01+(n9-1)*0.09)*(x2-x1),y90,rms,'Color',color9{n9},'FontSize',f9);
end
end
ylabel('Relative Tn (K)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[y1,y2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

% dirs={'Z';'W';'E';'N';'S'} for colors: k yel_grn g b r
subplot('position',[0.08,0.265,0.87,0.215])	% 3/4
hold on;
plot (timeobs(:,2),unobs(:,2),'-o','Color',yel_grn);
plot (timeobs(:,3),unobs(:,3),'-og');
text (x1+0.01*(x2-x1),yu95,'W','Color',yel_grn,'FontSize',f9);
text (x1+0.10*(x2-x1),yu95,'E','Color','g','FontSize',f9);
text (x1+0.55*(x2-x1),yu95,'Contact Qian Wu for FPI data usage','Color','k','FontSize',f9);
if (notie ~= 1)
 plot (timewd,un_wd(:,2),'Color',yel_grn,'LineStyle','--');
 plot (timewd,un_wd(:,3),'g','LineStyle','--');
 text (x1+0.19*(x2-x1),yu95,'-- TIEGCM','Color','k','FontSize',f9);
end
plot(xtot,zero,':k');
ylabel('Un (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yu1,yu2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'FontSize',fsize,'box','on');

% Markers: +,o,.,x,*,s,d,p,h,^,v,<,>
subplot('position',[0.08,0.050,0.87,0.215])	% 4/4
hold on;
plot (timeobs(:,4),vnobs(:,4),'-xb');
plot (timeobs(:,5),vnobs(:,5),'-xr');
plot (time_obs_20min,octas_20min*50,'-k');
text (x1+0.01*(x2-x1),yu95,'N','Color','b','FontSize',f9);
text (x1+0.10*(x2-x1),yu95,'S','Color','r','FontSize',f9);
text (x1+0.60*(x2-x1),yu95,'cloud cover in tenths*100/2','Color','k','FontSize',f9);
if (notie ~= 1)
text (x1+0.19*(x2-x1),yu95,'-- TIEGCM','Color','k','FontSize',f9);
 plot (timewd,vn_wd(:,4),'b','LineStyle','--');
 plot (timewd,vn_wd(:,5),'r','LineStyle','--');
end
plot(xtot,zero,':k');
ylabel('Vn (m/s)','FontSize',fsize);
xlabel('Daynumber','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[yu1,yu2],'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')

if (notie ~= 1)
%% Plot mad, residuals, and median over npmad pts
ymadx = 100.;
yr1 = -79.;
yr2 = 79.;

 plotmad = sprintf('rfpmad_06348_9');

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);
n9 = ndir+1;
 plot (timemad_obs(:,n9),tnmad_obs(:,n9),'Color','k','LineWidth',2.0);
for n9=1:ndir
if (madskp==0)
 plot (timemad_obs(:,n9),tnmad_obs(:,n9),'Color',color9{n9},'LineStyle',':');
end
plot (timemad,tnmad_wd(:,n9),'Color',color9{n9});
end
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;
n9 = ndir+1;
 plot (timemad_obs(:,n9),tnresid_obs(:,n9),'Color','k','LineWidth',2.0);
for n9=1:ndir
if (madskp==0)
 plot (timemad_obs(:,n9),tnresid_obs(:,n9),'Color',color9{n9},'LineStyle',':');
end
plot (timemad,tnresid_wd(:,n9),'Color',color9{n9});
end
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;
n9=ndir+1;
%plot (time_obs_20min,tn_obs_nan_20min(:,n9),'-k','LineWidth',2.0);
 madobs = sprintf('%4.2f',mad_obs(n9));
 plot (timemad_obs(:,n9),tnmed_obs(:,n9),'Color','k','LineWidth',2.0);
text (x1+0.80*(x2-x1),y90,'20min 4h med Obs','Color','k','FontSize',f9);
text (x1+0.80*(x2-x1),y95,'All Obs mad','Color','k','FontSize',f9,'LineWidth',3.0);
text (x1+0.92*(x2-x1),y95,madobs,'Color','k','FontSize',f9);
plot (timewd,tn_wd(:,n9),'Color',dark_green,'LineStyle','--','LineWidth',2.0);
for n9=1:ndir
 text (x1+(0.01+(n9-1)*0.09)*(x2-x1),y95,dirs{n9},'Color',color9{n9},'FontSize',f9);
if (madskp==0)
 madobs = sprintf('%4.2f',mad_obs(n9));
 plot (timemad_obs(:,n9),tnmed_obs(:,n9),'Color',color9{n9},'LineStyle',':');
 text (x1+(0.01+(n9-1)*0.09)*(x2-x1),y90,madobs,'Color',color9{n9},'FontSize',f9);
end
 madwd = sprintf('%4.2f',mad_wd(n9));
 plot (timemad,tnmed_wd(:,n9),'Color',color9{n9});
text (x1+(0.01+(n9-1)*0.09)*(x2-x1),y85,madwd,'Color',color9{n9},'FontSize',f9);
end
text (x1+0.80*(x2-x1),y85,'--av TIEGCM ','Color',dark_green,'FontSize',f9);
plot(xtot,zero,'--k');
ytit = sprintf ('%d hour median Tn (K)',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


