function BP9_Field_RUNUP(filename1)
% -------------------------------------------------------------------------
%               NTHMP (National Tsunami Hazard Mitigation Program)
% -------------------------------------------------------------------------
%                          2011 NTHMP VALIDATION WORKSHOP
%              Venue: Galveston, TAMU TX
%                     Margh 30- Apr. 1
%                               ---- o ----
%               MATLAB SCRIPT:  Okushiri Island Runup
%               TYPE:           Field Benchmarking Problem (Field)
%                               NUMBER: 09
%               TSUNAMI SOURCE: Field fata Runup (RUNUP)                       
%               SCRIPT NAME   : BP9_Field_RUNUP.m
%                
%               AUTHOR(S)     :J. Horrillo, Gyeong-Bo Kim.  	        
%                              Copyright 2011 NTHMP/NOAA
%                              Revision: 0
%                              Date: Sep/2011
%               REFERENCE:
%               Synolakis, C.E. (1986): The Runup of Long Waves. Ph.D.
%                   Thesis, California Institute of Technology, Pasadena,
%                   California, 91125, 228 pp. 
% 
%                   NOAA/PMEL Center for Tsunami research
%               http://nctr.pmel.noaa.gov/benchmark/index.html
% 
%                   Synolakis, C.E. (1987): The runup of solitary waves.
%               J. Fluid Mech., 185, 523-545.  
%         
% 
%               CONTACT INFO:   Dr. Juan J. Horrillo, 
%                               Texas A&M University at Galveston 
%                               Maritime Systems Engineering (MASE) 
%                               200 Seawolf Park or P.O. Box 1675 
%                               Galveston, TX  77553-1675 
%                               Phone: (409) 740-4465
%                               horrillj@tamug.edu
%                               Website: http://www.tamug.edu/
% 
% -------------------------------------------------------------------------
%   > BP9_Field_RUNUP(filename1)
%
%   BP9_Field_RUNUP reads input files 'filename1' 
%   containing model results and compares with runup field data
%   (Prof Imamura and Shuto... before we used Prof. Tsuji data).
%    
%   Input files 'filename1' need to be named as follows:
%   XXXX_BP9_OKUSHIRI.dat, where XXXX is the name of the numerical model
%   to be tested,
%   eg. NEOWAVE_BP9_OKUSHIRI.dat
% 
%   The input file, 'filename1' must have the following
%   configuration:
%   'filename1' is a matrix ASCII files of size
%   19x5(row x column) as shown below (no heading).
%
%   filename1 
%   +-----------+-----------+-----------+--------------+  
%   |                       Model result Model Result
%   | long (Deg)  lat (Deg)    range(m)   Avg(m)
%   |                        (min - max)    
%   +-----------+-----------+-----------+--------------+ 
%    139.4292117 42.18818149   5.5  9.3   7.1
%    139.4111857 42.16276287   7.0 10.2   8.5
%    139.4182612 42.13740439   8.0 12.6  10.7
%    139.4280358 42.09301238  11.5 15.3  13.9
%    139.426245  42.11655479  11.0 15.4  12.9
%    139.4237147 42.10041415   9.6 17.1  11.8
%    139.4289018 42.07663658  10.6 15.8  13.2
%    139.4278534 42.06546152  10.4 17.2  13.6
%    139.4515399 42.04469655  11.6 16.0  12.9
%    139.4565284 42.05169226  5.2  12.0   6.7
%    139.4720138 42.05808988  8.1  11.8  10.3
%    139.5150461 42.21524909  5.1   5.9   5.4
%    139.5545494 42.22698164  4.1   8.6   6.8
%    139.4934307 42.06450128  7.4  14.7  10.8
%    139.5474599 42.18744879  2.2   5.3   3.8
%    139.5258982 42.17101221  6.4   9.0   7.4
%    139.5625242 42.21198369  2.6   3.8   3.1
%    139.5190997 42.11305805  6.2   8.8   7.7
%    139.5210766 42.15137635  2.9   4.7   3.6
%   +-----------+-----------+-----------+--------------+ 
%
%   NOTE:
%   Modelers need only to change columns 3 to 5 with their
%   model result data
%
%   Examples:
%   > BP9_Field_RUNUP
%     In this example the Matlab script just plots field data range.
%   > BP9_Field_RUNUP('NEOWAVE_BP9_OKUSHIRI.dat')
%     In this example the Matlab script plots model results
%     and compares against field data range.
% -------------------------------------------------------------------------
help BP9_Field_RUNUP
% Prof IMAMURA's RUNUP data from Yoshiki
%        Log         Lat    Min  Max  Avg
%----------------------------------------
Field_data=[
 139.4292117 42.18818149    4.9  6.5  5.8;
 139.4111857 42.16276287    5.8  6.3  6.1;
 139.4182612 42.13740439    6.4  6.4  6.4;
 139.4280358 42.09301238   16.2 16.2 16.2;
 139.426245  42.11655479   12.2 12.2 12.2;
 139.4237147 42.10041415   13.7 31.7 21.0;
 139.4289018 42.07663658   19.3 19.3 19.3;
 139.4278534 42.06546152   13.6 21.8 17.7;
 139.4515399 42.04469655    9.8 15.4 12.4;   %10.7 15.4;
 139.4565284 42.05169226    6.2 10.1  8.0;   % 4.8 11.9;
 139.4720138 42.05808988    7.1 11.6  9.9;
 139.5150461 42.21524909    5.7  6.8  6.2;   %5.4  6.8;
 139.5545494 42.22698164    6.8  9.8  8.7;   %6.8  8.0;
 139.4934307 42.06450128    3.8  6.9  5.2;   %3.8  8.4;
 139.5474599 42.18744879    3.0  3.0  3.0;
 139.5258982 42.17101221    5.9  7.6  6.7;
 139.5625242 42.21198369    3.9  3.9  3.9;
 139.5190997 42.11305805    3.9  5.4  4.8;
 139.5210766 42.15137635    3.1  3.3  3.2];
% Frank -> RUNUP Prof. Tsuji data
% Field_data=[
%  139.4292117	42.18818149  5.6  9.0 nan;
%  139.4111857	42.16276287  5.1  7.0 nan;
%  139.4182612	42.13740439 13.0 13.0 nan;
%  139.4280358	42.09301238 12.6 12.6 nan;
%  139.426245   42.11655479 16.3 16.3 nan;
%  139.4237147	42.10041415	30.6 30.6 nan;
%  139.4289018	42.07663658	18.7 18.7 nan;
%  139.4278534	42.06546152	22.8 22.8 nan;
%  139.4515399	42.04469655	12.4 12.4 nan;
%  139.4565284	42.05169226	 3.2 10.2 nan;
%  139.4720138	42.05808988  8.3 13.2 nan;
%  139.5150461	42.21524909	 5.7  9.6 nan;
%  139.5545494	42.22698164	 8.7 11.1 nan;
%  139.4934307	42.06450128	 4.8  4.8 nan;
%  139.5474599	42.18744879	 3.3  3.3 nan;
%  139.5258982	42.17101221	 3.4  3.4 nan;
%  139.5625242	42.21198369	 5.0  5.0 nan;
%  139.5190997	42.11305805	 3.4  3.4 nan;
%  139.5210766	42.15137635	 6.8  7.7];
%
str='OKUSHIRI ISLAND RUNUP';
%--------------------------------------------------------------------------
%   Time series of surface elevation at gauges around a conical island
%--------------------------------------------------------------------------
str1=[' '];
err=nan;
err1=nan;
if nargin < 1
   a1=ones(19,6)*nan;
   filename1=[];
else
   if isempty(filename1),
      a1=ones(19,6)*nan;
      filename1=[];
   else   
      load(filename1);
      len=strfind(filename1, '.');
      filename1=filename1(1:len-1);      
      eval(['a1=' filename1 ';']);
      eval(['clear ' filename1 ';']);
      len1=strfind(filename1, '_');
      str1=filename1(1:len1-1);    
   end
end
['RUNUP OKUSHIRI ISLAND']
fig1=figure;
set(fig1,'name',str,...
         'position',[100 300 1300 800]);
annotation(fig1,'textbox',...
    'String',[str],...
    'FontWeight','bold',...
    'FontSize',14,...
    'FitHeightToText','off',...
    'LineStyle','none',...
    'Position',[0.40 0.94 0.5 0.05]);
%--------------------------------------------------------------------------
% FIELD DATA
%--------------------------------------------------------------------------
   ax=[0 20 0 35];
   axis(ax);
   xloc=[1:1:19];
   a= Field_data(:,4)-Field_data(:,3);
   a(find(a==0))=0.2;
   Y=[Field_data(:,3) a];
   hndl=bar(Y,'stack');
   set(hndl(1), 'FaceColor', 'none', 'EdgeColor', 'none')
   set(hndl(2), 'FaceColor', [0.6 1.0 0.4], 'EdgeColor', 'k','linewidth',1)
%--------------------------------------------------------------------------
N=length(xloc);
err1=nan;
if ~isempty(filename1),
  % GET NRMSD error
   ina= find(~isnan(a1(:,5)));
   N=length(ina)
   % Field AVG
   % ----------------------------------------------------------------------
   Fmin=Field_data(ina,3);
   Mmin=a1(ina,3);
   Favg=Field_data(ina,5);
   Mavg=a1(ina,5);
   Fmax=Field_data(ina,4);
   Mmax=a1(ina,4);
   den_Field=[Fmin Favg Fmax];
   %[ina abs(Fmin-Mmin) abs(Favg-Mavg) abs(Fmax-Mmax)]
   % Juan Siggestion
   den_Field=[Fmin Favg Fmax];
   [num_Model, Ind]=min([abs(Fmin-Mmin) abs(Favg-Mavg) abs(Fmax-Mmax)],[],2);
   % Elena Suggestion
   %den_Field=[Fmin Favg Fmax Fmin Fmax];
   %[num_Model, Ind]=min([abs(Fmin-Mmin) abs(Favg-Mavg) abs(Fmax-Mmax) abs(Fmin-Mavg) abs(Fmax-Mavg)],[],2);  
   %[ina num_Model Ind]
   %return
   % GET MAX WAVE AMPLITUDE or RUNUP
   suma=0.0;
   TMP=0.0;
   for ii=1:N,
       if (Mavg(ii)>=Fmin(ii) & Mavg(ii)<=Fmax(ii)),
           TMP(ii)=0.0;
       else
           TMP(ii)=num_Model(ii,1)/den_Field(ii,Ind(ii))*100;
       end
       suma=suma +TMP(ii);
   end
   err1=suma/N;
end
%--------------------------------------------------------------------------
% MODEL RESULTS
%-----------------------------------------------------------------------
   hold
   xlocl1=num2str(Field_data(:,1),['%9.5f']);
   xlocl2=num2str(Field_data(:,2),['%9.5f']);
   a= a1(:,4)-a1(:,3);
   a(find(a==0))=0.2;
   Y1=[a1(:,3) a];
   hndl1=bar(Y1,'stack');
   set(hndl1(1), 'FaceColor', 'none', 'EdgeColor', 'none')
   set(hndl1(2), 'FaceColor', 'r', 'EdgeColor', 'r','BarWidth',0.2)
   hndl3=plot(Field_data(:,5),'*','Markersize',14,'linewidth',2,'color',[0.0 0.4 0.0]);
   hndl2=plot(a1(:,5),'k+','Markersize',14,'linewidth',2);
   grid on
   set(gca,'Layer','top',...
           'Xtick',xloc,...
           'XtickLabel',Field_data(:,1)') % display gridlines on top of graph
       
   set(gca,'XTickLabel','')
% rotate Xlabel  painful
% -----------------------------------------------
lx1 = text(xloc,ax(3)*ones(1,length(xloc))-0.1,xlocl1);
lx2 = text(xloc,ax(3)*ones(1,length(xloc))-0.1,xlocl2);
set(lx1,'HorizontalAlignment','right','VerticalAlignment','bottom', ...
'Rotation',90,'FontSize',10);
set(lx2,'HorizontalAlignment','right','VerticalAlignment','top', ...
'Rotation',90,'FontSize',10);
% Determine the lowest point. The X-label will be
% placed so that the top is aligned with this point.
% Place the axis label at this point
XMidPoint = xloc(1)+(xloc(end)-xloc(1))/2;
LowYPoint = -3.50;
text(XMidPoint,LowYPoint,'Runup Coordinate Around Okushiri (Degree)', ...
'VerticalAlignment','top', ...
'HorizontalAlignment','center','FontSize',14);
%---------------------------------
ylabel('Runup (m)','FontSize',14)
hasbehavior(hndl(1),'legend',false)
hasbehavior(hndl1(1),'legend',false)
if ~isempty(filename1),
   legend({'Field Data Range', [str1 ' Range'],'Field Data Avg.', [str1 ' Avg.']});
else
   legend({'Field Data Range'}); 
end
if ~isnan(err1)
    for ii=1:N,
      text(ina(ii), 1.01*ax(4),[num2str(TMP(ii),'%3.0f') '%'],'fontsize',10,...
          'HorizontalAlignment','center','VerticalAlignment','bottom')  
    end
text(0.77*ax(2),0.83*ax(4),[' - ERR Runup =' num2str(err1,'%3.0f') '%'],'fontsize',14)
end
return