  SUBROUTINE cas_acc

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_acc  ***
   !&E
   !&E ** Purpose : definition of a test case to play with mixed layer instabilities
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : cas_init
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables : hx, hy, h0, sal, temp, dy, dhx, dhy, dsigu, dsigw
   !&E                         surf, fisac, gwac, dgac, dfiac, fu, fv, sig, sigw 
   !&E                         tfiusr, tfivsr, goulet, govlet, ig, id, jb, jh 
   !&E                         ssh,u,v,uz,vz
   !&E
   !&E
   !&E
   !&E Rq valerie : Revue du calcul de fu et fv non pas coherent avec la fiche COMODO
   !&E Rq valerie : Pb d ecriture pour le fichier nudging quand on utilise 160 proc
   !&E
   !&E
   !&E
   !&E
   !&E
   !&E
   !&E ** Reference : 
   !&E                Boccaletti et al, JPO 2007
   !&E
   !&E ** History :
   !&E       !  2013     (X. Couvelard) Original code
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : rsh,rlg,riolg,kmax,jmin,ljmin,ljminm2,jmax,ljmax,ljmaxp2, &
                          imin,limin,liminm2,imax,limax,limaxp2,obc_width,grav,pi,  &
                          liminm1,ljminm1,valmanq,liminp2,limaxm1,ljminp2,ljmaxm1,  &
                          liminp1,ljminp1,limaxp1,ljmaxp1,ishift_obc_east,ishift_obc_west, &
                          jshift_obc_north,jshift_obc_south
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dtmaxi,dt,suffix,qmax,l_modele2d, &
                         h0,h0fond,l_obc_cycl_x,frot_lin,hx,hy,dhx,dhy,tdeb,tsauv, &
                         halfdt,l_corio_cst,surf,ig,id,jb,jh,dx_u,dy_v,hex,hey, &
                         vish_xe,vish_phi,l_initfromfile,l_saverestart_bydate,saverestart_step
   USE comvarp2d, ONLY : sal,temp,ssh,u,v
#ifdef key_casacc
   USE modcasinit,ONLY : temp0,uz0,vz0,u0,v0,ssh0,cas_init_special
#endif
   USE comvarp3d, ONLY : uz,vz,boz
   USE comvars3d, ONLY : l_stateeq_lin,kzinit,nzinit,rhoref,dsigu, &
                         turb_nbeq,turb_2eq_option
   USE output,    ONLY : pasor,l_out_zaxis_w,l_out_turb,l_out_sal,stepi,l_out_zaxis_w, &
                         nbave,l_out_timeave,tdebsor,fileresu,l_filebydate
   USE obccombine,ONLY : l_obc_west,l_obc_east,l_obc_cycl
   USE comcolig2d,ONLY : igdu,igdv
#if defined key_siggen || defined key_gencoord
   USE comsiggen,  ONLY : b_sig,theta_sig,hc,hzex,hzey
   USE comsiggen,  ONLY : hcx_sig,hcy_sig,dcusds
#endif
   USE toolgeom,  ONLY : f_zw,f_zu
   USE pressure,  ONLY : pres_gradt
   USE obcapplyu, ONLY : obc_apply_u
   USE obcapplyuz, ONLY : obc_apply_uz
   USE obcapplyv, ONLY : obc_apply_v
   USE obcapplyvz, ONLY : obc_apply_vz
   USE traobc3dapply, ONLY : ts_obc3d_apply
   USE comcolig3d,ONLY : yz_exp
   USE casinit,   ONLY : cas_init_param,cas_init_grid,cas_init_bathy, &
                         cas_init_io,cas_init_obc,cas_init_meteo,     &
                         cas_init_visc,cas_init_turb,cas_init_tspar,  &
                         cas_init_progvar,cas_init_misc,cas_init_river
   USE_MPI toolmpi,ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ADD_ALL_MPI_RSH,add_tag
   USE output, ONLY : iextmin,iextmax,igextmin,igextmax
#ifdef key_nudging
   USE nudging, ONLY : l_out_nud
#endif
#ifdef key_chainop
   USE parameters, ONLY : lchain
   USE comvars2d, ONLY  : name_in_dt,l_initfromfile,file_init
#endif

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   LOGICAL           :: l_nc4par
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k,nj_wall_s,nj_wall_n,nb,nb_file
   INTEGER           :: istr,jstr,N,MMm,LLm,jet_perturb
   INTEGER           :: ixl,iyl
   REAL(KIND=rsh)    :: dxval,dyval,z0bottom,z0surface,depthcst
   REAL(KIND=rsh)    :: Ljet,rhomaxs,bgdrhos,zs1,dzs,drhos
   REAL(KIND=rsh)    :: rhomaxn,bgdrhon,zn1,dzn
   REAL(KIND=rsh)    :: drhosfs,drhosfn,z0,drhosf,z0p
   REAL(KIND=rsh)    :: cff_perturb
   REAL(KIND=rsh)    :: zk,dzs_a,dzn_a,drhon
   REAL(KIND=rsh)    :: x,y,Fyz,y0
   REAL(KIND=rsh)    :: mssh,area,cff1 
   REAL(KIND=rsh), DIMENSION(0:kmax,2)    :: rhoprof

   REAL(KIND=rsh), DIMENSION(1:kmax,liminm1:limaxp2,ljminm1:ljmaxp2) :: rho
   REAL(kind=rsh),ALLOCATABLE,DIMENSION(:,:,:) :: P,ipx,ipy
   REAL(KIND=rsh), DIMENSION(1:kmax)    :: hzex0,hzey0
   REAL(KIND=rsh)    :: f0
   REAL(KIND=rlg)    :: beta,el
   REAL(KIND=rlg)    :: rhox,rhoy
   REAL(KIND=rsh)    :: asymth
   REAL(KIND=rlg)    :: tool_datosec

#ifdef key_chainop
   CHARACTER(LEN=lchain) :: filenam
   NAMELIST/nmlrestart/ name_in_dt
   NAMELIST/namdate/ date_start,date_end
   NAMELIST/namchainop/ dtmin,dtmax,suffix,date_start3d,l_initfromfile,file_init,saverestart_step
#endif

   !!----------------------------------------------------------------------
   !! * Executable part

   !------------------------------------------------
   !------------------------------------------------
   !!     courant constant zonal et barotrope     !!
   !!     surface libre lineaire                  !!
   !!     densite constant                        !!
   !------------------------------------------------
   !------------------------------------------------

#ifdef key_casacc

   ! general parametrization (each parameter can be reassigned)
   ! -----------------------

   ! dynamics
   CALL cas_init_param

   ! over-assignement
   l_corio_cst=.false.

   ! tracers : among others, saliref=saliref_lin=35.5_rsh, tetaref=tetaref_lin=10.0_rsh
   CALL cas_init_tspar

   ! over-assignement
   rhoref=1024.4_rsh
   l_stateeq_lin=.TRUE.   ! oui et advection de densite avec -Dkey_equeta_temp


   ! bathymetry
   ! ----------

   ! assignement : input variables for cas_init_bathy routine
   depthcst = 4000.0_rsh
   nj_wall_s = 4  ! land over nj_wall_s cells (from jmin row) at the southern part of the domain
   nj_wall_n = 4  ! land over nj_wall_n cells at the northern part of the domain
   nj_wall_s = 2  ! land over nj_wall_s cells (from jmin row) at the southern part of the domain
   nj_wall_n = 1  ! land over nj_wall_n cells at the northern part of the domain
#if defined key_siggen || defined key_gencoord
   hc= 100.0_rsh
#endif



   ! flat bottom, plus walls
   CALL cas_init_bathy(depthcst,nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s)

   ! over-assignement


   ! grid parametrization
   ! --------------------

   ! assignement : input variables for cas_init_grid routine
#ifdef key_cas_20km
   dxval=20000.0_rsh
   dyval=20000.0_rsh
#elif key_cas_10km
   dxval=10000.0_rsh
   dyval=10000.0_rsh
#elif key_cas_5km
   dxval=5000.0_rsh
   dyval=5000.0_rsh
#elif key_cas_2km
   dxval=2000.0_rsh
   dyval=2000.0_rsh
#elif key_cas_1km
   dxval=1000.0_rsh
   dyval=1000.0_rsh
#endif
   fisac=0.0_riolg
   gwac=0.0_riolg
   f0 = 2.0_rsh * 2.0_rsh*pi/(86164.0_rsh)
   fisac=ASIN((1.1e-04-0.00001016_rsh)/f0)/(pi/180.0_rsh)  ! lat to get corio=1.1e-04 + shift (comparaison ROMS)
print*,'latitude',fisac,'pour corio = ',1.1e-04-0.00001016_rsh
   beta=1.6e-11
#if defined key_siggen || defined key_gencoord
   theta_sig= 5.0_rsh ! entre 0 et 20
   b_sig= 0.0_rsh     ! entre 0 et 1
#endif
#ifdef key_agrif
   IF (.NOT.agrif_root()) THEN
     rhox=Agrif_Rhox()
     rhoy=Agrif_Rhoy()

     ixl=Agrif_ix()
     iyl=Agrif_iy()

     dfiac = (Agrif_Parent(dfiac))/rhoy
     fisac = Agrif_Parent(fisac)+(iyl-0.5_rlg-0.5_rlg/rhoy) &
             *Agrif_Parent(dfiac)
     dgac = ( Agrif_Parent(dgac) )/rhox
     gwac = Agrif_Parent(gwac)+(ixl-0.5_rlg-0.5_rlg/rhox) &
             *Agrif_Parent(dgac)

     dxval=dxval * Agrif_Curgrid%agrif_dx(1)
     dyval=dyval * Agrif_Curgrid%agrif_dx(2)
   END IF
#endif


   CALL cas_init_grid(dxval,dyval,betaplan=beta)

   ! over-assignement

   i=10
   j=50
   IF_MPI ( i>=limin .AND. i>=limax .AND. j>=ljmin .AND. j>=ljmax) THEN
   DO k=0,kmax
!      PRINT*,'f_zw',f_zw(depthcst,0.0_rsh,k,i,j)
   END DO
   ENDIF_MPI


   ! open boundaries
   ! ---------------

   CALL cas_init_obc

   ! over-assignement
   ! geometry (where are the open boundaries ?)
   ! Do not specify for cyclic open boundaries
   ! used scheme to manage obc
   IF_AGRIF (agrif_root()) THEN
   l_obc_cycl = .TRUE.
   l_obc_cycl_x=.TRUE.
   ENDIF_AGRIF


   ! IO parametrization
   ! ------------------

   ! assignement : input variables for cas_init_io routine
   date_start=   '01/01/1900 00:00:00'
   date_end=   '06/01/1902 01:00:00'
   date_end=   '01/01/1902 01:00:00'
#ifdef key_chainvalid
   ! With AGRIF date_end must be a multiple of the mother time step
   date_end=   '06/01/1900 00:00:00'
#endif
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
#ifdef key_cas_20km
   dtmin=800.0_rlg
   dtmax=800.0_rlg
   dt=800.0_rlg
#elif key_cas_10km
   dtmin=400.0_rlg
   dtmax=400.0_rlg
   dt=400.0_rlg
#elif key_cas_5km
   dtmin=100.0_rlg
   dtmax=100.0_rlg
   dt=100.0_rlg
#elif key_cas_2km
   dtmin=50.0_rlg
   dtmax=50.0_rlg
   dt=50.0_rlg
#else
   dtmin=50.0_rlg
   dtmax=50.0_rlg
   dt=50.0_rlg
#endif
   suffix='acc'
#ifdef key_chainop
   filenam='./paracom.txt'
   OPEN(50,file=filenam,status='old',form='formatted',access='sequential')
   READ(50,namdate)
   CLOSE(50)
   filenam='./paraspec.txt'
   OPEN(50,file=filenam,status='old',form='formatted',access='sequential')
   READ(50,namchainop)
   CLOSE(50)
   dt = dtmin
#endif
#ifdef key_agrif
   dtmaxi = dtmax
#endif
   IF_AGRIF (.NOT.agrif_root()) THEN
     INST_AGRIF dtmin=Agrif_Parent(dtmin)/(Agrif_rhot())
     INST_AGRIF dtmax=Agrif_Parent(dtmax)/(Agrif_rhot())
     INST_AGRIF dt=Agrif_Parent(dt)/(Agrif_rhot())
   ENDIF_AGRIF
   halfdt=dt*0.5_rlg

   l_nc4par = .TRUE.

#if defined key_agrif
   nb_file=1
#else
   nb_file=2
#endif
   CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix,l_nc4par=l_nc4par,nb_file=nb_file)

   ! over-assignement
#if ! defined key_agrif
#if key_cas_20km
   fileresu(1)="champs_20km.nc"
#elif key_cas_10km
   fileresu(1)="champs_10km.nc"
#elif key_cas_5km
   fileresu(1)="champs_5km.nc"
#elif key_cas_2km
   fileresu(1)="champs_2km.nc"
#elif key_cas_1km
   fileresu(1)="champs_1km.nc"
#endif
#endif
   IF ( l_filebydate(1) ) fileresu(1)="results"
   IF ( .NOT. l_nc4par) THEN
     DO i=1,1   ! done for results_2.nc file
        CALL_MPI ADD_TAG(fileresu(i),LEN(fileresu(i)))
     END DO
   END IF

#if ! defined key_chainop
   pasor(:)= 5.0_rlg*24.0_rlg*3600.0_rlg ! 24*3600
   tdebsor(1)=tool_datosec('02/01/1900 00:00:00')
!   l_out_timeave(1)=.TRUE.
   IF (l_out_timeave(1)) THEN
     nbave=1
     tdebsor(1)=tdebsor(1)+pasor(1)
   END IF
#endif
#if defined key_agrif && defined key_chainvalid
   pasor(:)= 21600.0_rlg ! 6*3600
   tdebsor(1)=tool_datosec('01/01/1900 00:00:00')
#endif
   l_out_sal(:)=.FALSE.
   l_out_zaxis_w(:)=.TRUE.
   !l_initfromfile=.TRUE.
   l_saverestart_bydate=.TRUE.
#ifdef key_cas_20km
   saverestart_step=30.0_rlg
#elif key_cas_10km
   saverestart_step=30.0_rlg
#elif key_cas_5km
   saverestart_step=15.0_rlg
#elif key_cas_2km
   saverestart_step=10.0_rlg
#elif key_cas_1km
   saverestart_step=5.0_rlg
#endif
#ifdef key_chainop
   filenam='./paraspec.txt'
   OPEN(50,file=filenam,status='old',form='formatted',access='sequential')
   READ(50,namchainop)
   CLOSE(50)
#endif
   saverestart_step=saverestart_step*86400.0_rlg  ! days to seconds
   tsauv=tdeb + saverestart_step
#ifdef key_nudging
   l_out_nud(2)=.TRUE.
   !l_out_zaxis_w(1)=.TRUE.
   igextmin(2)=imax/2
   igextmax(2)=imax/2
   iextmin(2)=imax/2
   iextmax(2)=imax/2
#endif
#if defined key_MPI_2D
   DO nb=1,nb_file
     ! correct i indices
     DO WHILE ( iextmin(nb) < limin )
       iextmin(nb)=iextmin(nb)+stepi(nb)
     END DO
     iextmax(nb)=MIN(limax,iextmax(nb))
   END DO
#endif


   ! viscosity
   !-----------

   ! none
   CALL cas_init_visc

   ! over-assignement
   vish_xe(:,:,:) = 5.0_rsh
   vish_phi(:,:,:)= 5.0_rsh


   ! turbulence
   ! ----------

   z0bottom   =1.0e-13  ! pour permettre division mais useless -Dkey_no_bottom_stress or -Dkey_linear_bottom_stress
   z0surface  =0.0_rsh
   z0surface  =0.6_rsh
   kzinit=1.2e-5
   nzinit=2e-3
   !kzinit=0.0_rsh
   !nzinit=0.0_rsh

   ! Gaspard parametrization
   CALL cas_init_turb(kzinit,nzinit,z0surface,z0bottom)

   ! over-assignement
#ifdef key_linear_bottom_stress
   frot_lin=0.005_rsh ! Linear bottom drag coefficient (m/si) revoir unit
#endif
   turb_nbeq=2
   turb_2eq_option=2


   ! meteorology
   ! -----------

   ! no forcing
   CALL cas_init_meteo


   ! River and outflows
   ! -------------------

   ! no river
   CALL cas_init_river


   ! initialization of variables used by the nudging (special to this test case)
   ! -----------------------------------------
   CALL cas_init_special


   ! initialization of the pronostic variables
   ! -----------------------------------------

   ! set to 0,tetaref,saliref
   CALL cas_init_progvar

   IF_AGRIF (agrif_root()) THEN

   ! over-assignement
   jet_perturb = 0 ! no perturbation
   jet_perturb = 1

!
! --- Jet and density profile parameters ---
!
      Ljet        = 1600.e3;  ! Jet width in km
!
! southern profile
      rhomaxs=27.75     ! max density
      bgdrhos=9.8e-6    ! background stratification
      zs1=-1000
      dzs=700
      drhos=1.4         ! gives strength of Phillips mode (0 to 1.4)
! northern profile 
      rhomaxn=27.7573   ! max density
      bgdrhon=9.8e-6    ! background stratification
      zn1=-400
      dzn=300
! surface Charney mode 
      drhosfs=1.5       ! northern anomaly  (0 for Phillips mode only)
      drhosfn=0.0       ! southern anomaly  (0 for Phillips mode only)
      z0=-300.          ! vertical penetration (-600 for strong Charney mode)
      drhosf=0.00       ! weakens Charney mode (canceled for 1.05)
      z0p=-110.         ! only used for drhosf non zero
!
! zonal perturbation
!
      if (jet_perturb .eq. 0) then
        cff_perturb = 0.0    ! zonal perturbation coefficient
      else
        cff_perturb = 0.02
      endif
!
! global variables
!
      istr=(limax+limin)/2
      jstr=(ljmax+ljmin)/2
      N=kmax
      el=2000000.0_rsh   ! Length of the domain of interest (m)
      ! LLm,MMm  Number of the internal points of the PHYSICAL grid.
      !          in the XI- and ETA-directions [physical side boundary
      !          points and peroodic ghost points (if any) are excluded].
      MMm=el/dyval
      LLm=500000.0_rsh/dxval
print*,'LLm=',LLm,'MMm=',MMm
!
! --- Build northern and southern density profiles ---
!
! First compute background density profiles associated with
! a gentle stratification that does not depend 
! on jet side. It is there to ensure static stability. 
      do k=1,N
        zk = f_zu(depthcst,0.0_rsh,k,istr,jstr)
        rhoprof(k,1)=rhomaxs-bgdrhos*(zk+depthcst)
        rhoprof(k,2)=rhomaxn-bgdrhon*(zk+depthcst)
!print*,'k,rhoprof(k,1),rhoprof(k,2)',k,rhoprof(k,1),rhoprof(k,2)
      enddo
!
! Second, get main north/south contrast with a distorded 
! tanh shape of variable amplitude.Distorsion of the tanh 
! is done with the asym function that increases 
! stratification in the upper ocean
      dzs_a=1.3*dzs
      do k=1,N
        zk=asymth(f_zu(depthcst,0.0_rsh,k,istr,jstr),zs1,dzs_a)
        rhoprof(k,1)=rhoprof(k,1)-drhos*(0.5+0.5*tanh((zk-zs1)/dzs))
      enddo
      dzn_a=1.3*dzn
      drhon=-(rhoprof(N,1)-rhoprof(N,2))/(0.5+0.5*tanh((zk-zn1)/dzn))
      do k=1,N
        zk=asymth(f_zu(depthcst,0.0_rsh,k,istr,jstr),zn1,dzn_a)
        rhoprof(k,2)=rhoprof(k,2)-drhon*(0.5+0.5*tanh((zk-zn1)/dzn))
      enddo
      do k=1,N
        zk = f_zu(depthcst,0.0_rsh,k,istr,jstr)
        rhoprof(k,1)=rhoprof(k,1) &
     &             -drhosf*(exp((zk-z0p)/abs(z0p)))/(exp(1.))
        rhoprof(k,1)=rhoprof(k,1) &
     &             -drhosfs*0.5*(1+tanh((zk-z0)/abs(z0)))/tanh(1.)
        rhoprof(k,2)=rhoprof(k,2) &
     &             -drhosf*(exp((zk-z0p)/abs(z0p)))/(exp(1.))
        rhoprof(k,2)=rhoprof(k,2) &
     &             -drhosfn*0.5*(1+tanh((zk-z0)/abs(z0)))/tanh(1.)
      enddo
!
! --- Fit southern and northern profiles ---
!
       do k=1,N
         do j=ljminm1,ljmaxp2
           do i=liminm1,limaxp2
             y=j-jshift_obc_south
             x=i-1-ishift_obc_west
             !y=y/real(MMm)-0.5
             y=y/real(MMm)-0.51
             x=x/real(LLm-2)
!
! Initial state without the small perturbation
!
             y0 = y*pi*el/Ljet + pi/2

             if(y0.lt.0)then
               Fyz=1.0
             elseif(y0.gt.pi)then
               Fyz=0.0
             else
               Fyz=1.0-(y0-sin(y0)*cos(y0))/pi
             endif
             temp0(k,i,j)=Fyz*rhoprof(k,1)+(1.-Fyz)*rhoprof(k,2)
!
! Add vertical attenuation and XI perturbation
!
             y = y + cff_perturb*exp(f_zu(depthcst,0.0_rsh,k,i,j)/1000.) &
     &                          *(exp(-(x-0.5)**2/0.05)) &
     &                          *( 0.5*sin(2.*pi*x) &
     &                            +0.5*sin(6.*pi*x) )
!
! Set Jet structure and place its center at y=pi/2
!
             y = y*pi*el/Ljet + pi/2

             if(y.lt.0)then
               Fyz=1.0
             elseif(y.gt.pi)then
               Fyz=0.0
             else
               Fyz=1.0-(y-sin(y)*cos(y))/pi
             endif
             rho(k,i,j)=Fyz*rhoprof(k,1)+(1.-Fyz)*rhoprof(k,2)
             boz(k,i,j) = -rho(k,i,j)*grav/rhoref
             !boz(k,i,j) = rho(k,i,j)
           enddo
         enddo
       enddo

!print*,'init temp0',temp0(:,12,50)
!
! --- Initialization of the temperature ---
!
   temp(:,liminm1:limaxp2,ljminm1:ljmaxp2) = rho(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   !CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp0(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   !CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp0(:,liminm1:limaxp2,ljminm1:ljmaxp2))

!
! --- Initialization of the sea surface height ---
!

!
! --- Get surface pressure (and zeta) using pressure gradient routine ---
!
      ALLOCATE( P(1:kmax,limin:limaxp1,ljmin:ljmaxp1) )
      call pres_gradt(-boz,ssh0,P)
!ATTENTION ne marche que si ssh=0 initialement car la pression due a l elevation
!de surface n est pas pris en compte dans P
      do j=ljminp1,ljmax
        do i=liminp1,limax
          ssh(i,j)=-P(1,i,j)/grav ! cancels bottom pressure
        enddo
      enddo
!
! Remove mean sea level (mssh)
!
      mssh=0.0_rsh
      area=0.0_rsh
      do j=ljminp2,ljmaxm1
        do i=MAX0(limin,2+ishift_obc_west), MIN0(limax,imax-1-ishift_obc_east)
          IF (h0(i,j) > -valmanq) THEN
           cff1=surf(i,j)
           mssh=mssh+cff1*ssh(i,j)
           area=area+cff1
          END IF
        enddo
      enddo
#ifdef key_MPI_2D
      CALL_MPI ADD_ALL_MPI_RSH(mssh) 
      CALL_MPI ADD_ALL_MPI_RSH(area) 
! La parallelisation MPI change necessairement les resultats par rapport au sequentiel (car sommes effectuees dans ordres differents)
! Rq : mettre mssh=-100.485212797575_rsh (affichage de la valeur obtenue par le run sequentiel change tout autant (peu bien sur) les resultats
#endif
      mssh=mssh/area
#ifdef key_chainvalid
      mssh=-100.543328205658_rsh
#endif
print*,'mssh=mssh/area',mssh
      do j=ljminp1,ljmax    ! ATTENTION peux changer les resultats
        do i=MAX0(limin,2+ishift_obc_east), MIN0(limax,imax-1-ishift_obc_east)
          IF (h0(i,j) > -valmanq) THEN
           ssh(i,j)=ssh(i,j)-mssh
          END IF
        enddo
      enddo

   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))

   ENDIF_AGRIF

   ! Domain of calculation
   ! Percent of wet surface on each cell
   CALL cas_init_misc


   IF_AGRIF (agrif_root()) THEN

   ! over-assignement
   id(:)=id(:)-ishift_obc_east
   ig(:)=ig(:)+ishift_obc_west
   jh(:)=jh(:)-jshift_obc_north
   jb(:)=jb(:)+jshift_obc_south
!
! Apply ssh boundary conditions (must be done after the initialization of ig and id due to cyclic obc)
!
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     DO j=ljminm2,ljmaxp2
       ssh(ig(j)-1,j)=ssh(id(j)-2,j)
       ssh(ig(j),j)=ssh(id(j)-1,j)
       ssh(id(j),j)=ssh(ig(j)+1,j)
     END DO
#endif
   ENDIF
   ! Required for restart file
   DO j=ljminm2,ljmaxp2
     IF (limin==imin) ssh(imin:ig(j)-2,j)=0.0_rsh
     IF (limax==imax) ssh(id(j)+1:imax,j)=0.0_rsh
   END DO
   WHERE (h0(liminm2:limaxp2,ljminm2:ljmaxp2)==-valmanq) ssh(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh

!   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
!   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))

!
! Apply temp boundary conditions (must be done after the initialization of ig and id due to cyclic obc)
!
   CALL ts_obc3d_apply(uz0,vz0,sal,temp)

   !
   ! Estimate of the velocities from the pressure gradients
   !
!
! --- Get gressure gradient (ru,rv) using updated ssh field ---
!
      ALLOCATE( ipx(1:kmax,limin:limax,ljmin:ljmax) )
      ALLOCATE( ipy(1:kmax,limin:limax,ljmin:ljmax) )
      call pres_gradt(boz,ssh,P,ipx,ipy)
      ! warning P does not take into account the pressure introduced by the ssh ! hum not so sure...

!
! --- Compute geostrophic uz velocities ---
!
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (igdv(i,j) == 1) THEN
         DO k=1,kmax
           !attention u a la position v en fait
           uz(k,i,j)=-grav/dy_v(i,j)*(ssh(i,j+1)-ssh(i,j))/fv(i,j)+ipy(k,i,j)/fv(i,j)
         END DO
       END IF
     END DO
   END DO
!
! --- Compute geostrophic uz velocities ---
!
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF(igdu(i,j) == 1) THEN
         DO k=1,kmax
           !attention v a la position u en fait
           vz(k,i,j)=grav/dx_u(i,j)*(ssh(i+1,j)-ssh(i,j))/fu(i,j)-ipx(k,i,j)/fu(i,j)
         END DO
       END IF
     END DO
   END DO
!
! --- Compute geostrophic barotropic u velocity ---
!
   DO j=ljminp2,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igdu(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         u(i,j)=SUM(uz(:,i,j)*hzex(:,i,j)*dsigu(:))/hex(i,j)
#else
         u(i,j)=SUM(uz(:,i,j)*dsigu(:))
#endif
       END IF
     END DO
   END DO
!
! --- Compute geostrophic barotropic v velocity ---
!
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF (igdv(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         v(i,j)=SUM(vz(:,i,j)*hzey(:,i,j)*dsigu(:))/hey(i,j)
#else
         v(i,j)=SUM(vz(:,i,j)*dsigu(:))
#endif
       END IF
     END DO
   END DO

   ! For perfect restart
   DO j=ljminm2,ljmaxp2
   DO i=liminm2,limaxp2
     IF (hx(i,j)==-valmanq-h0fond) u(i,j)=0.0_rsh
     IF (hy(i,j)==-valmanq-h0fond) v(i,j)=0.0_rsh
     IF (hx(i,j)==-valmanq-h0fond) uz(:,i,j)=0.0_rsh
     IF (hy(i,j)==-valmanq-h0fond) vz(:,i,j)=0.0_rsh
   END DO
   END DO

   !
   ! Apply open boundary conditions (cyclic along x)
   !
   CALL obc_apply_u(u)
   CALL obc_apply_v(v)
   CALL obc_apply_uz(uz)
   CALL obc_apply_vz(vz)

   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz(:,liminm2:limaxp2,ljminm2:ljmaxp2))


   !
   ! Initialization without the perturbation (for the nudging)
   !
!
! --- Initialization of the sea surface height ---
!

!
! --- Get surface pressure (and zeta) using pressure gradient routine ---
!
      call pres_gradt(temp0*grav/rhoref,ssh0,P)
!ATTENTION ne marche que si ssh=0 initialement car la pression due a l elevation!de surface n est pas pris en compte dans P. Attention ssh utilise un peu quand meme
      do j=ljminp1,ljmax
        do i=liminp1,limax
          ssh0(i,j)=-P(1,i,j)/grav ! cancels bottom pressure
        enddo
      enddo
!
! Remove mean sea level (mssh)
!
      mssh=0.0_rsh
      area=0.0_rsh
      do j=ljminp2,ljmaxm1
        do i=MAX0(limin,2+ishift_obc_west), MIN0(limax,imax-1-ishift_obc_east)
          IF (h0(i,j) > -valmanq) THEN
           cff1=surf(i,j)
           mssh=mssh+cff1*ssh0(i,j)
           area=area+cff1
          END IF
        enddo
      enddo
#ifdef key_MPI_2D
      CALL_MPI ADD_ALL_MPI_RSH(mssh) 
      CALL_MPI ADD_ALL_MPI_RSH(area) 
! La parallelisation MPI change necessairement les resultats par rapport au sequentiel (car sommes effectuees dans ordres differents)
! Rq : mettre mssh=-100.485212797575_rsh (affichage de la valeur obtenue par le run sequentiel change tout autant (peu bien sur) les resultats
#endif
      mssh=mssh/area
print*,'mssh=mssh/area',mssh
      do j=ljminp1,ljmax    ! ATTENTION peux changer les resultats
        do i=MAX0(limin,2+ishift_obc_east), MIN0(limax,imax-1-ishift_obc_east)
          IF (h0(i,j) > -valmanq) THEN
           ssh0(i,j)=ssh0(i,j)-mssh
          END IF
        enddo
      enddo

   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh0(liminm2:limaxp2,ljminm2:ljmaxp2))

!
! Apply ssh boundary conditions (must be done after the initialization of ig and id due to cyclic obc)
!
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
     CALL_MPI ex_peri_x_rsh(2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh0(liminm2:limaxp2,ljminm2:ljmaxp2))
#else
     DO j=ljminm2,ljmaxp2
       ssh0(ig(j)-1,j)=ssh0(id(j)-2,j)
       ssh0(ig(j),j)=ssh0(id(j)-1,j)
       ssh0(id(j),j)=ssh0(ig(j)+1,j)
     END DO
#endif
   ENDIF
   ! Required for restart file (strange not so sure it is really required)
   DO j=ljminm2,ljmaxp2
     IF (limin==imin) ssh0(imin:ig(j)-2,j)=0.0_rsh
     IF (limax==imax) ssh0(id(j)+1:imax,j)=0.0_rsh
   END DO
   WHERE (h0(liminm2:limaxp2,ljminm2:ljmaxp2)==-valmanq) ssh0(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh


!   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh0(liminm2:limaxp2,ljminm2:ljmaxp2))
!   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh0(liminm2:limaxp2,ljminm2:ljmaxp2))

!
! Apply temp boundary conditions (must be done after the initialization of ig and id due to cyclic obc)
!
   CALL ts_obc3d_apply(uz0,vz0,sal,temp0)

   !
   ! Estimate of the velocities from the pressure gradients
   !
!
! --- Get gressure gradient (ru,rv) using updated ssh0 field ---
!
      call pres_gradt(-temp0*grav/rhoref,ssh0,P,ipx,ipy)
      ! warning P does not take into account the pressure introduced by the ssh0

!
! --- Compute geostrophic uz velocities ---
!
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
       IF (igdv(i,j) == 1) THEN
         DO k=1,kmax
           !attention u a la position v en fait
           uz0(k,i,j)=-grav/dy_v(i,j)*(ssh0(i,j+1)-ssh0(i,j))/fv(i,j)+ipy(k,i,j)/fv(i,j)
         END DO
       END IF
     END DO
   END DO
!
! --- Compute geostrophic uz velocities ---
!
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF(igdu(i,j) == 1) THEN
         DO k=1,kmax
           !attention v a la position u en fait
           vz0(k,i,j)=grav/dx_u(i,j)*(ssh0(i+1,j)-ssh0(i,j))/fu(i,j)-ipx(k,i,j)/fu(i,j)
         END DO
       END IF
     END DO
   END DO
!
! --- Compute geostrophic barotropic u velocity ---
!
   DO j=ljminp2,ljmaxm1
     DO i=liminp1,limaxm1
       IF (igdu(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         DO k=1,kmax
           hzex0(k)=(ssh0(i,j)+dhx(i,j)+hcx_sig(i,j)+(hx(i,j)-hcx_sig(i,j))*dcusds(k))
         END DO
         u0(i,j)=SUM(uz0(:,i,j)*hzex(:,i,j)*dsigu(:))/hex(i,j)
#else
         u0(i,j)=SUM(uz0(:,i,j)*dsigu(:))
#endif
       END IF
     END DO
   END DO
!
! --- Compute geostrophic barotropic v velocity ---
!

   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF (igdv(i,j)/=4) THEN
#if defined key_siggen || defined key_gencoord
         DO k=1,kmax
           hzey0(k)=(ssh0(i,j)+dhy(i,j)+hcy_sig(i,j)+(hy(i,j)-hcy_sig(i,j))*dcusds(k))
         END DO
         v0(i,j)=SUM(vz0(:,i,j)*hzey0(:)*dsigu(:))/hey(i,j)
#else
         v0(i,j)=SUM(vz0(:,i,j)*dsigu(:))
#endif
       END IF
     END DO
   END DO

   ! For perfect restart (strange not so sure it is really required)
   DO j=ljminm2,ljmaxp2
   DO i=liminm2,limaxp2
     IF (hx(i,j)==-valmanq-h0fond) u0(i,j)=0.0_rsh
     IF (hy(i,j)==-valmanq-h0fond) v0(i,j)=0.0_rsh
     IF (hx(i,j)==-valmanq-h0fond) uz0(:,i,j)=0.0_rsh
     IF (hy(i,j)==-valmanq-h0fond) vz0(:,i,j)=0.0_rsh
   END DO
   END DO

   !
   ! Apply open boundary conditions (cyclic along x) (for MPI, if limax-limin too small)
   !
   CALL obc_apply_u(u0)
   CALL obc_apply_v(v0)
   CALL obc_apply_uz(uz0)
   CALL obc_apply_vz(vz0)

   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz0(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,uz0(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz0(:,liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,kmax,liminm2,limaxp2,ljminm2,ljmaxp2,vz0(:,liminm2:limaxp2,ljminm2:ljmaxp2))


   ! deallocation
   DEALLOCATE(P,ipx,ipy)

   ENDIF_AGRIF
#endif

  END SUBROUTINE cas_acc

   !!======================================================================

  FUNCTION asymth(zz,zasym,dzasym)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION asymth  ***
   !&E
   !&E ** Purpose :
   !&E
   !&E ** Description : 
   !&E
   !&E ** Called by : casacc
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2013-10  (V. Garnier)  Extracted from ROMS
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : rsh

   !! * Arguments
   REAL(KIND=rsh), INTENT( in ) :: zz,zasym,dzasym

   !! * Local declarations
   REAL(KIND=rsh) :: asymth

   !!----------------------------------------------------------------------
   !! * Executable part

   asymth = (zz-zasym) &
     &       *sqrt(1.0_rsh+0.5_rsh*((zz-zasym)+abs(zz-zasym))**2/dzasym**2) &
     &       +zasym

  END FUNCTION asymth

   !!======================================================================

