function [Mvar,Mmask,Mh,distx]=get_section_meanlongshore_nemo(grd,ctlu,ctlv,ctlt,vname,rec,dist,latmin,latmax,NN);
%
%% input:
%  grd      Grid structure as loaded from grd=rnt_gridload_XA(model)
%  ctl      Output structure as loaded from ctl=ctlload(root,his,model,nstrt,nend,ninc);
%  vname    name of the variable (string) 'zeta', 'temp', 'u' ... and 'r1u'
%           where r stands for rotated, digit (1, 2 or 3 is for sharpness of rotation angle
%           from 1 sharp to 3 smooth and u is cross-shore positive toward it, v alongshore)
%  rec      record list in ctl eg as determined from get_recnumbers or get_reclist
%  dist     Number of points away from the mask (check what 0 means).
%  latmin   For longshore averaging. It is assumed that the coastline can be represented 
%           as a function of latitude. Regularize and send auxiliary grid if necessary. 
%  latmax   same
%  NN       number of vertical levels
%  output:
%   Mvar       Mean variable. 
%   Mmask      Mean mask for the section. Mask is determined by how many points are available
%   Mh         Mean bathymetry for the section. 
%   distx      x-coordinate (ready to use in pcolor(x,z,varextract))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%uu='vozocrtx';
 uu='uo';
%vv='vomecrty';
 vv='vo';

lon=grd.lonr;lat=grd.latr;pm=grd.pm;pn=grd.pn; h=grd.h;
gridinfo=rnt_gridinfo(grd.id);
grdfileH=gridinfo.grdfileH; % the grid file with the mask that we want to follow
                            % when doing alongshore averaging.
dL=0;dM=0;
if vname(1)=='u'
  mask=grd.masku;wmask='mask_u';lat=grd.latu;type='u';dM=-1;
elseif vname(1)=='v'
  mask=grd.maskv;wmask='mask_v';lat=grd.latv;type='v';dL=-1;
else
  mask=grd.maskr;wmask='mask_rho';lat=grd.latr;type='r';
end
  angle=grd.angle;
  for i=1:NN; angleg(i,:,:)=angle(:,:); end
  size(angleg)
M=grd.Mp+dM; L=grd.Lp+dL;

if vname(1:1)=='r'  % if rotated velocities.
  ncgrd=netcdf(grdfileH);
  if vname(2:2)=='1'
    ang=ncgrd{'ang_rho_sharp'}(:);
  elseif vname(2:2)=='2'
    ang=ncgrd{'ang_rho'}(:);
  elseif vname(2:2)=='3'
    ang=ncgrd{'ang_rho_smooth'}(:);
  end
  close(ncgrd)
end
%%%%angc=ones(size(ang,1),size(ang,2))*30.94; angc=angc*pi/180;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 LIM=follow_mask(grdfileH,wmask);
 for i=1:length(LIM); latlim(i)=lat(i,LIM(i)); end
 for i=1:length(LIM); lipm(i)=pm(i,LIM(i)); end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
jmax=min(find(latlim>latmax));
jmin=max(find(latlim<latmin));
%jmax
%jmin
mpm=mean(1./lipm(jmin:jmax)); nps=floor(dist/mpm);
%%%
[ll]=find(LIM(jmin:jmax)-nps < 0);
if isempty(ll)==0; display('ATTENTION ! dist trop long'); end
%%%
distx=[0:mpm:nps*mpm];distx=distx-distx(end);
%%%
for i=jmin:jmax; hsec(i-jmin+1,1:nps+1)=h(i,LIM(i)-nps-1:LIM(i)-1); end
%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% READ VARIABLE
if vname(1:1)=='r'  % then we want to compute a rotated velocity;
	u=rnt_loadvar_xa(ctlu,rec,uu); u=u(:,:,1:end-1); u=u2rho_3d(u);
	v=rnt_loadvar_xa(ctlv,rec,vv); v=v(:,1:end-1,:); v=v2rho_3d(v);
    size(u);
    size(v);
%%%    vns=v.*cos(angleg)+u.*sin(angleg);
%%%    uns=u.*cos(angleg)-v.*sin(angleg);
    var0=0.*u;
    if vname(end:end)=='v'  % alongshore velocity positive when the coast on the right
      for j=1:M
        var0(:,j,:)=cos(ang(j)).*v(:,j,:)-sin(ang(j)).*u(:,j,:);
      end
    elseif vname(end:end)=='u' % cross-shore velocity positive toward the coast
      for j=1:M
        var0(:,j,:)=sin(ang(j)).*v(:,j,:)+cos(ang(j)).*u(:,j,:);
      end
    end
else
  var0=rnt_loadvar_xa(ctlt,rec,vname);
end

%%%%hbl=rnt_loadvar_xa(ctl,rec,'hbl');
%%%%display('WARNING : here hbl is "hbl" in the file !!')

var=squeeze(var0(:,jmin:jmax,:)); %%%%mhbl=hbl(jmin:jmax,:);

for i=jmin:jmax; 
varsec(:,i-jmin+1,1:nps+1)=var(:,i-jmin+1,LIM(i)-nps-1:LIM(i)-1); 
%%%%mhblsec(i-jmin+1,1:nps+1)=mhbl(i-jmin+1,LIM(i)-nps-1:LIM(i)-1);
end
clear var
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% Estimation du "mask moyen" + Section Moyenne
%
Mmask=zeros(NN,length(distx));
Mvar=zeros(NN,length(distx));
for i=1:jmax-jmin+1
  vermask=squeeze(varsec(:,i,:));
  nn=zeros(size(vermask,1),size(vermask,2));
%%  nn(~isnan(vermask))=1;
  nn(vermask~=0)=1;
  Mmask=Mmask+nn;
%%  vermask(isnan(vermask))=0; Mvar=Mvar+vermask;
  vermask(vermask==0)=0; Mvar=Mvar+vermask;
end
Mvar=Mvar./Mmask;
Mmask(Mmask<(jmax-jmin+1)/2)=0;
Mmask(Mmask==0)=NaN;
Mmask(~isnan(Mmask))=1;
%% Mean section
Mvar=Mvar.*Mmask;
Mh=mean(hsec,1);
%%%%Mhbl=mean(mhblsec,1);
%%