function [my_latriv,my_lonriv,my_flow,my_rivername,number_river]=oct_runoff_glob_extract(grdfile,global_clim_river)
%
%  Read all the rivers from the netcdf global climatology file
%  Select the ones located into the model grid

disp(' ')
disp(['Reading the global monthly climatological run-off dataset... '])
disp(' ')

%  Read the grid
[latr,lonr,maskr]=read_latlonmask(grdfile,'r');
dl=0.0;
minlat=min(min(latr))-dl;
maxlat=max(max(latr))+dl;
minlon=min(min(lonr))-dl;
maxlon=max(max(lonr))+dl;

%  Read the global rivers using procedural netcdf
ncid = netcdf.open(global_clim_river,'NC_NOWRITE');
time = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'time'));
lonriv_mou = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lon_mou'));
latriv_mou = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'lat_mou'));
% name variables (char arrays) need transpose and char conversion
ct_name = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'ct_name'));
cn_name = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'cn_name'));
% Some name vars can be large char arrays; transpose to get [n x len]
try
  riv_name_tmp = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'riv_name'));
  riv_name = char(riv_name_tmp');
catch
  riv_name = '';
end
try
  ocn_name_tmp = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'ocn_name'));
  ocn_name = char(ocn_name_tmp');
catch
  ocn_name = '';
end
try
  stn_name_tmp = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'stn_name'));
  stn_name = char(stn_name_tmp');
catch
  stn_name = '';
end

FLOW_clm = netcdf.getVar(ncid,netcdf.inqVarID(ncid,'FLOW_clm'));
netcdf.close(ncid)

% Select the rivers in the domain grid
disp(' ')
disp('Selecting rivers inside the model domain...')

%River selection criteria
maxmaxflow=(max(max(FLOW_clm)));
maxflow=(max(FLOW_clm)');  % attention transpose

meanflowmax=mean(FLOW_clm(:,1));    % for the biggest
meanflow=nanmean(FLOW_clm(:,:),1)'; % for others

% by default : all rivers in the domain
rivdetectype='DEFAULT';
my_riv=find(latriv_mou>=minlat & latriv_mou<=maxlat & lonriv_mou >= minlon & lonriv_mou <= maxlon);

disp(['There are ',num2str(length(my_riv)),' rivers in the domain : '])
disp(['===='])
disp(' ')
disp([' => Rivers detection type : ',rivdetectype,' (see in runoff_glob_extract.m)'])
disp(' ')
disp(['===='])
disp(['There are ',num2str(length(my_riv)),' rivers in the domain : '])
disp(['Domain contains rivers :'])
for k=1:length(my_riv)
   if ischar(riv_name)
     nm = strtrim(riv_name(my_riv(k),:));
   else
     nm = '';
   end
   if ischar(ocn_name)
     oc = ocn_name(my_riv(k),1:min(4,size(ocn_name,2)));
   else
     oc = '';
   end
   disp([num2str(k),' - ',nm,' flowing in ocean ',oc])
end

my_flow=FLOW_clm(:,my_riv);
if ischar(riv_name)
  my_rivername=riv_name(my_riv,:);
else
  my_rivername = [];
end
number_river=length(my_riv);
my_latriv=latriv_mou(my_riv);
my_lonriv=lonriv_mou(my_riv);

% make a figure
if ~isempty(my_flow)
    figure(100)
    plot([1:12],my_flow)
    if ~isempty(my_rivername)
      legend(my_rivername,'location','northeastoutside')
    end
    box on, grid on
    title(['\bf Monthly clim of the domain run off'])
    xlabel(['\bf Month']);ylabel(['\bf Discharge in m3/s'])
    set(gca,'Xtick',[0.5:11.5],'XtickLabel',['J';'F';'M';'A';'M';'J';'J';'A';'S';'O';'N';'D']);
end

return
