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

 year1 = 2004;
 day1 = 314;
 day2 = 319;
 dayofmo1 = 9;
nf1 = 5;
it1nf1 = 61;
nf2 = 10;


% 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/tiegcm_dyncrit-linux/';

% 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_2004001-2004366.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


%%  Read JRO 5-min obs (some gaps)
jrofil = sprintf ('/hao/polar1/emery/aurora_d/jro_viz/jro041109b');
JRO=fopen(jrofil,'r');
nskip=2;
% Read over nskip lines
for n=1:nskip
 linjro = fgetl(JRO);
end

% JRO obs start at 0000 UT for the 20 min values
for n=1:nobs
 jro_viz_20min(n) = missing;
 jro_nan_20min(n) = NaN;
 time_jro_20min(n) = day1 + (n-1)/ntmpda;
% TEMP for something to plot in second plot
% timecon(n) = time_jro_20min(n);
% mlat1_con(n) = 75.;
% mlat2_con(n) = 60.;
end

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

% Specifics for Nov 9-13, 2004 (begin time 1150 UT and end time 1700 UT)
nmad1 = 12*3+1+midm1;
nmad2 = ntmpda*lenda - (24-17)*3 - midm1 + 1;
jrd = 1186;
iday0 = 314;

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) + 1;
	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,num_obs_20min
% 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

ifskp = 1

%%  Read TIEGCM .nc files starting at 0000 UT on 04309

itw5d10 = 0;
for num=nf1:nf2
 it1 = 1;
    if (num==nf1), it1=it1nf1; end
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 's_nov04_crb_5_10_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_5d10(itw5d10+it-it1+1) = jroav_wi(it);
 time5d10(itw5d10+it-it1+1) = dayfrac(it);
 mlat1_5d10(itw5d10+it-it1+1) = 90. - crit1(it);
 mlat2_5d10(itw5d10+it-it1+1) = 90. - crit2(it);
 crb(itw5d10+it-it1+1) = 90. - crit1(it) + 5.;
 mlat0pot(itw5d10+it-it1+1) = 90. - crit1(it) + 5. - crit1(it) + 5.;
end
itw5d10 = itw5d10 + nt - (it1-1);
end		
itw5d10

itwcon = 0;
for num=nf1:nf2
 it1 = 1;
    if (num==nf1), it1=it1nf1; end
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 's_nov04_wcon_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_con(itwcon+it-it1+1) = jroav_wi(it);
 timecon(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		
itwcon

itwd7p5 = 0;
for num=nf1:nf2
 it1 = 1;
    if (num==nf1), it1=it1nf1; end
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 's_nov04_crb_7_15_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=nf1:nf2 for 20 min saves 
itwd7p5

%if (ifskp==0)
itw5d15 = 0;
for num=nf1:nf2
 it1 = 1;
    if (num==nf1), it1=it1nf1; end
    str=num2str(num);
    if (num < 10), str=['0' str]; end
    Fname=[fpath 's_nov04_crb_5_15_0' str '.nc'];
    [nt,dayfrac,jroav_wi,crit1,crit2] = get_wi_nc(Fname,ii,jj,ht1,ht2);
for it=it1:nt
 jroav_5d15(itw5d15+it-it1+1) = jroav_wi(it);
 time5d15(itw5d15+it-it1+1) = dayfrac(it);
 mlat1_5d15(itw5d15+it-it1+1) = 90. - crit1(it);
 mlat2_5d15(itw5d15+it-it1+1) = 90. - crit2(it);
end
itw5d15 = itw5d15 + nt - (it1-1);
end		
itw5d15

%end	%ifskp


%% calculate rms and print from JRO aves every 20 min (from averaging 2 5-min values)
% [rmsw,nrms] = rms_jro(obs,model,n1,n2,missing)
filerms = sprintf('rms_nov04_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, 'Nov04  num   concrit  crb+5,10  crb+5,15  crb+7.5,15 \n');
iend = min(itw5d10,nobs);
%iend = min(iend,itwcon);
 [rmswcon,i20] = rms_jro(jro_viz_20min,jroav_con,1,iend,-999);
 [rmsw5d15,i20] = rms_jro(jro_viz_20min,jroav_5d15,1,iend,-999);
 [rmsw5d10,i20] = rms_jro(jro_viz_20min,jroav_5d10,1,iend,-999);
 [rmswd7p5,i20] = rms_jro(jro_viz_20min,jroav_7d15,1,iend,-999);
fprintf(RMS, '09-13 %4.0d %7.2f %7.2f %7.2f %7.2f \n \n', i20,rmswcon,rmsw5d10,rmsw5d15,rmswd7p5);
for j=1:lenda
 n1 = (j-1)*ntmpda + 1;
 n2 = min(iend,n1 + ntmpda - 1);
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_con,n1,n2,-999);
 dmswcon(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_5d10,n1,n2,-999);
 dmsw5d10(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_5d15,n1,n2,-999);
 dmsw5d15(j) = rmsw;
 [rmsw,nrms] = rms_jro(jro_viz_20min,jroav_7d15,n1,n2,-999);
  j20(j) = nrms;
 dmswd7p5(j) = rmsw;
fprintf(RMS,'%5.0d %4.0d %7.2f %7.2f %7.2f %7.2f \n', j+dayofmo1-1,j20(j),dmswcon(j),dmsw5d10(j),dmsw5d15(j),dmswd7p5(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
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_con,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);
end
iob,nmad
% Calculate daily and total mads:
mad_con = 0.;
mad_5d10 = 0.;
mad_5d15 = 0.;
mad_7d15 = 0.;
for i=1:nmad
 mad_con = mad_con + wimad_con(i);
 mad_5d10 = mad_5d10 + wimad_5d10(i);
 mad_5d15 = mad_5d15 + wimad_5d15(i);
 mad_7d15 = mad_7d15 + wimad_7d15(i);
end
mad_con = mad_con/nmad
mad_5d10 = mad_5d10/nmad
mad_5d15 = mad_5d15/nmad
mad_7d15 = mad_7d15/nmad
fprintf(RMS, '\n MAD (median absolute deviations) of 20-min JRO Viz and TIEGCM WI_VEL values \n');
fprintf(RMS, 'Nov04  num     obs   num   concrit  crb+5,10  crb+5,15  crb+7.5,15 \n');
fprintf(RMS, '09-13 %4.0d %7.2f %4.0d %7.2f %7.2f %7.2f %7.2f \n \n', iob,mad_obs,nmad,mad_con,mad_5d10,mad_5d15,mad_7d15);
j1 = 1;
for j=1:lenda
 dmad_con(j) = 0.;
 dmad_5d10(j) = 0.;
 dmad_5d15(j) = 0.;
 dmad_7d15(j) = 0.;
 j2 = j1 + ntmpda-1;
 if (j==1)
   j1 = j1 + midm1;
 end
 if (j==lenda)
   j2 = j2 - midm1 + 1;
 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);
 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);
fprintf(RMS,'%5.0d %4.0d %7.2f %4.0d %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_5d10(j),dmad_5d15(j),dmad_7d15(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=-190;
 y2=140;
y95 = y1 + 0.95*(y2-y1);
y90 = y1 + 0.90*(y2-y1);

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

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


%% Do lenda+1 plots:  each day separate, then all lenda days
for np=1:lenda+1

if (np == lenda+1)
 x1 = day1;
 x2 = day2;
 daytit = sprintf('TIEGCM W05 09-13 Nov 2004');
 plottit = sprintf('jro_09_13nov04');
 rmscon = sprintf('%4.2f',rmswcon);
 rms5d15 = sprintf('%4.2f',rmsw5d15);
 rms5d10 = sprintf('%4.2f',rmsw5d10);
 rms7d15 = sprintf('%4.2f',rmswd7p5);
else
x1 = day1 + (np-1);
x2 = x1 + 1;
daytit = sprintf('TIEGCM W05 %2.0d Nov 2004',np+dayofmo1-1);
plottit = sprintf('jro_%2.0dnov04',np+dayofmo1-1);
 rmscon = sprintf('%4.2f',dmswcon(np));
 rms5d10 = sprintf('%4.2f',dmsw5d10(np));
 rms5d15 = sprintf('%4.2f',dmsw5d15(np));
 rms7d15 = sprintf('%4.2f',dmswd7p5(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,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');

subplot('position',[0.08,0.480,0.87,0.240])	% 2/3
hold on;
plot (time5d10,crb,'Color','k','LineWidth',3.0);
text (x1+0.01*(x2-x1),ym95,'CRB','Color','k','FontSize',f9);
plot (time5d10,mlat0pot,'Color','k','LineWidth',2.0);
text (x1+wr1m*(x2-x1),ym95,'0Pot','Color','k','FontSize',f9);
plot (timecon,mlat1_con,'Color','b','LineStyle','--');
plot (timecon,mlat2_con,'b');
text (x1+(wr1m+pcwr)*(x2-x1),ym95,'concrit','Color','b','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 (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,0Pot', 'FontSize',fsize);
xlabel(' ','FontSize',fsize);
set (gca,'xlim',[x1,x2],'ylim',[ym1,ym2],'XTickLabel',{' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';' ';},'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 (timecon,jroav_con,'b');
text (x1+(wr1v+pcwr)*(x2-x1),y95,'concrit','Color','b','FontSize',f9);
text (x1+(wr1v+pcwr)*(x2-x1),y90,rmscon,'Color','b','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 (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_09_13nov04');
 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);
else
 plotmad = sprintf('jromad_%2.0dnov04',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));
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);
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(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);
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);
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

