clear all; close all;

   id=netcdf('/loceanfs/pulsation/cholod/peru12_nemo/inputs/grd/mesh_mask_peru12.L75_from_trop12.nc','r')
   lon=id{'nav_lon'}(:);
   lat=id{'nav_lat'}(:);
   bat=squeeze(id{'mbathy'}(1,:,:)); 


% Dans VertAvgSectionMmeanNemo.m 
%   vname='r2v';  prename='along_'; dist=3e5;
%   latpm=3;
%   latsec=-10;

%  lat1 = 10+3 = -7
%  lat2 = 10-3 = -13

%  A l'equateur, 110km ~= 1° => au 12eme, 1 pt ~= 9km
%  La section: dist = 3e5 = 300 km => 300/9km =>  33 points pour 300km

% sur la fenêtre matlab, en zommant, on cherche les 2 points de côte à -7 et -13° de lat
   lon1 = -79.67;
   lat1 = -7.065;
   lon2 = -76.33;
   lat2 = -13.05;

% puis on cherche les correspondants i et j:
% lon(1,245) = -79.67  i1 = 245
% lat(294,1) = -7.065  j1 = 294
% lon(1,285) = -76.33  i2 = 285
% lat(220,1) = -13.05  j2 = 220

%%%% on cherche la longitude pour les 2 points 28 points plus loin au large:
%%%% i1b =  i1-28 = 245-28 = 217
%%%% lon1b = lon(1,217) = -82
   %%%lon1b = -82;
%%%% i2b = i2-28 = 285-28 = 257
%%%% lon2b = lon(1,257) = -78.67 
   %%%lon2b = -78.67;

% on cherche la longitude pour les 2 points 33 points plus loin au large:
% i1b = i1-33 = 245-33 = 212
% lon1b = lon(1,212) = -82.42
   lon1b = -82.42;
% i2b = i2-33 = 285-33 = 252
% lon2b = lon(1,252) = -79.08 
   lon2b = -79.08;

   pcolor(lon,lat,bat); shading flat 
   line([lon1b lon1],[lat1 lat1],'color','k','LineWidth',2)
   line([lon2b lon2],[lat2 lat2],'color','k','LineWidth',2)
   line([lon1b lon2b],[lat1 lat2],'color','k','LineWidth',2)

   set(gca,'FontS',14); set(gca,'PlotBoxAspectRatio',[0.7 1 1]); box on;
   caxis ([0 75]); hcc=colorbar; set(hcc,'FontS',14);
   title('bathy (level) - Upwelling section area','FontSize',12);
   figname1='Upwelling_section_area';
   print ('-dpng',figname1)