  SUBROUTINE cas_advspotmount

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_advspotmount  ***
   !&E
   !&E ** Purpose : advection of a spot of salinity by a 2D constant
   !&E              velocity field u=v=1
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : cas_init
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference : 
   !&E
   !&E ** History :
   !&E       !  2007-09  (T. Duhaut)  Original code
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : pi,rsh,rlg,riolg,imin,imax,jmin,jmax,kmax,   &
                          ljmin,ljminm1,ljminm2,limax,limaxp1,limaxp2, &
                          ljmax,ljmaxp1,ljmaxp2,limin,liminm1,liminm2, &
                          ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,suffix,l_maxu, &
                         l_obc_cycl_x,l_obc_cycl_y,l_modele2d,qmax,halfdt, &
                         ierrorlog,hex,hey,hx,hy,dhx,dhy,hx_g,hy_g,h0_g, &
                         h0fond,dx,t,tfin,h0,dx_u,ig,id,jb,jh
   USE comvarp2d, ONLY : ssh,u,v,temp
   USE comvarp3d, ONLY : uz,vz,boz
   USE comvars3d, ONLY : qmaxz,rhoref,tetaref,kzinit,nzinit,wz,dsigu
   USE output,    ONLY : pasor,l_out_sal,l_out_zaxis_w, &
                         l_out_ubt,l_out_ssh
   USE obccombine,ONLY : l_obc_cycl
   USE comfluxdo, ONLY : fxeau,fyeau
   USE toolgeom,  ONLY : f_zw,f_zu,f_dzu
#ifdef key_substance
   USE comsubstance, ONLY : nvpc,nvp,nv_cp,nv_dis,nv_fix,nv_dri,nv_int,nv_sorb,     &
                            name_var,long_name_var,standard_name_var,nv_ncp,        &
                            valid_min_var,valid_max_var,cv_wat,nb_var,date_startsub,&
                            name_filesubs,itypv,ws_free_opt,ws_free_min,ws_free_max
   USE casinit,   ONLY : cas_init_subs,cas_init_subs_disspart
#endif
   USE_MPI toolmpi, ONLY: ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh
   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
#ifdef key_tssub_casadvlense
   USE tssubeqtest, ONLY : uzc0,vzc0,uc0,vc0,wzc0,wzc0_phys
#endif
#ifdef key_diag_physwz
   USE comvars3d, ONLY : wz_phys
#endif

   IMPLICIT NONE

#ifdef key_MPI_2D
   INCLUDE 'mpif.h'
#endif

   !! * Arguments

   !! * Local declarations
   LOGICAL           :: l_nc4par
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   CHARACTER(LEN=19) :: tool_sectodat
   INTEGER           :: i,j,k,IERR_MPI
   REAL(KIND=rlg)    :: x,y,z,xini,yini,zini,radius,dist,dtemp,dz
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst
   REAL(KIND=rsh)    :: cff1,cff2,x0,cff,dh,drho,dhu,dhum,dhup,sigma
   REAL(KIND=rsh)    :: xx,xx0,yy,y0,rotation,limit
   REAL(KIND=rsh), PARAMETER    :: Lx = 300000.0_rsh, &  ! length of the domain
                                   cfl = 0.1_rsh,     &  ! CFL criterion
#ifdef key_agrif
                                   dxval = 1000.0_rsh, &
#else
                                   dxval = Lx / (imax-(ishift_obc_west+ishift_obc_east)-2), & ! horizontal resolution
#endif
                                   dyval = dxval,     &
                                   vel = 6.666666666666667E-002
#ifdef key_substance
   INTEGER           :: isubs
   LOGICAL           :: l_filesubs
#endif
   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT*, '!!!!!!!!!!!!!!!!!!!!'
   PRINT*, 'TEST CASE : AdvLense'
   PRINT*, '!!!!!!!!!!!!!!!!!!!!'

#if ! defined key_tssub_casadvlense
   WRITE(ierrorlog,*) 'For this test case, you must use key_tssub_casadvlense ccp key'
   WRITE(ierrorlog,*) 'Simulation stopped.'
   CALL_MPI MPI_FINALIZE(IERR_MPI)
   STOP
#else

   ! limit for the exponential
   IF (rsh==8) limit = -700.0_rsh
   IF (rsh==8) limit = -80.0_rsh
   IF (rsh==4) limit = -80.0_rsh


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

   ! dynamics
   CALL cas_init_param

   ! over-assignement
   l_maxu=.FALSE.

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

   ! over-assignement
   qmax=12.0
   qmaxz=10.0
   rhoref=1000.0_rsh
   tetaref=0.0_rsh


   ! Bathymetry
   ! ----------

   ! assignement : input variables for cas_init_bathy routine
   depthcst = 3000.0_rsh
   dh = 1100.0_rsh
   cff1 = (1.0_rsh / 5000.0_rsh)**2
   cff2 = (1.0_rsh / 8000.0_rsh)**2
   x0 = 150000.0_rsh
   y0 = (jmax-(jshift_obc_south+jshift_obc_north)-2-0.5_rsh)/2*dyval
   rotation = SIN(45.0_rsh*pi/180.0_rsh)   ! seamount rotated by 45 degrees

   DO j=jmin,jmax
     DO i=imin,imax
       xx = (i-0.5_rsh)*dxval
       yy = (j-0.5_rsh)*dyval
       xx0 = x0-rotation*(yy-y0)
       h0_g(i,j) = depthcst - dh*EXP(MAX( -cff1*(MIN(xx-xx0,0.0_rsh))**2 &
                                      -cff2*(MAX(xx-xx0,0.0_rsh))**2,limit) ) &
                            + h0fond
!IF (j==150) print*,h0_g(i,j),i,MIN((i-0.5_rsh)*dxval-x0,0.0_rsh),MAX((i-0.5_rsh)*dxval-x0,0.0_rsh)
       xx = (i        )*dxval
       yy = (j-0.5_rsh)*dyval
       xx0 = x0-rotation*(yy-y0)
       hx_g(i,j) = depthcst - dh*EXP(MAX( -cff1*(MIN(xx-xx0,0.0_rsh))**2 &
                                      -cff2*(MAX(xx-xx0,0.0_rsh))**2,limit) )
       xx = (i-0.5_rsh)*dxval
       yy = (j        )*dyval
       xx0 = x0-rotation*(yy-y0)
       hy_g(i,j) = depthcst - dh*EXP(MAX( -cff1*(MIN(xx-xx0,0.0_rsh))**2 &
                                      -cff2*(MAX(xx-xx0,0.0_rsh))**2,limit) )
       ! sans doute mieux de faire la moyenne
     END DO
   END DO

   depthcst = 0.0_rsh ! no change of the bathymetry

   ! no change of the bathymetry, no wall
   CALL cas_init_bathy(depthcst)


   ! over-assignement
   depthcst = 3000.0_rsh



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

   ! assignement : input variables for cas_init_grid routine
   fisac=1.0_riolg
   gwac=1.0_riolg

   CALL cas_init_grid(dxval,dyval)

   ! over-assignement
   fu(:,:)=0.0_rsh
   fv(:,:)=0.0_rsh


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

   ! assignement : input variables for cas_init_io routine
   dtmin=1500.0_rlg
   dtmax=1500.0_rlg
   dt=1500.0_rlg
   halfdt=0.5_rlg*dt

   date_start=   '01/01/1900 00:00:00'
!   tfin = 200.0_rlg *dxval**2*pi /(cfl*dt)
!   tfin = 200.0_rlg *dxval**2*pi /(cfl*dt) + dt
!   date_end= tool_sectodat( tfin )
   date_end= '18/02/1900 11:58:10'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'

   suffix='AdvLense'
   l_nc4par = .TRUE.

   CALL cas_init_io(date_start,date_end,date_start3d,date_startagrif,suffix,l_nc4par=l_nc4par)

   ! over-assignement
   pasor(1)=20.0_rlg*dt
   pasor(1)=200.0_rlg*1500.0_rlg
   l_out_sal(1)=.FALSE.
   l_out_ubt(1)=.FALSE.
   l_out_ssh(1)=.FALSE.
   l_out_zaxis_w(1)=.true.

   IF ((kmax==1) .AND. (l_modele2d.EQV..FALSE.) ) THEN
     l_modele2d=.TRUE.
     PRINT*,'variable l_modele2d should be set to .true. as kmax=1'
     PRINT*,'corrected directly in routine casAdvCst.F90'
   END IF


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

   ! no open boundary
   CALL cas_init_obc

   ! over-assignement
   ! used scheme to manage obc
   l_obc_cycl = .TRUE.
   l_obc_cycl_x=.TRUE.
   l_obc_cycl_y=.TRUE.

   ! viscosity
   !-----------
   ! none
   CALL cas_init_visc


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

   z0bottom   =0.0_rsh
   z0surface  =0.0_rsh
   kzinit=0.0_rsh
   nzinit=0.0_rsh

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

   ! over-assignement


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

   ! no forcing
   CALL cas_init_meteo


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

   ! no river
   CALL cas_init_river


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

   ! set to 0,tetaref,saliref
   CALL cas_init_progvar

   ! over-assignement

   ! maximum velocities
   ! ------------------
   ALLOCATE( uzc0(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2) )
   ALLOCATE( vzc0(1:kmax,liminm2:limaxp2,ljminm2:ljmaxp2) )
   ALLOCATE( uc0(liminm2:limaxp2,ljminm2:ljmaxp2) )
   ALLOCATE( vc0(liminm2:limaxp2,ljminm2:ljmaxp2) )
   ALLOCATE( wzc0(0:kmax,liminm1:limaxp1,ljminm1:ljmaxp1) )
   ALLOCATE( wzc0_phys(0:kmax,limin:limax,ljmin:ljmax) )

   zini = 930.0_rlg       ! center at 930 m depth
   dhu = 600.0_rlg        ! thickness woth horizontal velocity
!   vel = cfl * dxval / dt
   sigma = 80.0_rlg
   DO j=ljminm2,ljmaxp2
     y = (j-1.5d0)*dyval
     DO i=liminm2,limaxp2
       x = (i-1.5d0)*dxval
       DO k=1,kmax
         ! zonal velocity
         !z = f_zu(hx(i,j),0.5_rsh*(ssh(i,j)+ssh(i+1,j)),k,i,j)
         z = f_zu(hx(i,j),0.0_rsh,k,i,j)
         IF (ABS(z) > 200.0_rsh .AND. ABS(z) < 1600.0_rsh) THEN
           dhup = (dhu - z -zini)/sigma
           dhum = (-dhu - z -zini)/sigma
           uzc0(k,i,j) = vel * 0.5_rlg*(TANH(dhup) -TANH(dhum)) 
         END IF
         ! meridional velocity
         !z = f_zu(hy(i,j),0.5_rsh*(ssh(i,j)+ssh(i,j+1)),k,i,j)
         z = f_zu(hy(i,j),0.0_rsh,k,i,j)
         IF (ABS(z) > 200.0_rsh .AND. ABS(z) < 1600.0_rsh) THEN
           dhup = (dhu - z -zini)/sigma
           dhum = (-dhu - z -zini)/sigma
           vzc0(k,i,j) = vel * 0.5_rlg*(TANH(dhup) -TANH(dhum)) 
         END IF
       END DO
     END DO
   END DO

   ! initialization of the tracers
   ! -----------------------------
   temp(:,:,:) = tetaref
   xini = (75 + 2)*dxval  ! center at 75 km
   yini = (75 + 2)*dyval
   xini = (100 + 2 + 0.5_rlg)*dxval  ! center at 100 km
   yini = (100 + 2 + 0.5_rlg)*dyval
   xini = (100 + 0.5_rlg)*dxval  ! center at 100 km
   yini = (100 + 0.5_rlg)*dyval
   zini = 930.0_rlg       ! center at 930 m depth
   radius = 0.001_rlg
   DO j=ljmin,ljmax
     y = (j-1.5d0)*dyval
     DO i=limin,limax
       x = (i-1.5d0)*dxval
       DO k=1,kmax
         z = -f_zu(h0(i,j),ssh(i,j),k,i,j)
         dist = (x-xini)**2/(Lx**2) + (y-yini)**2/(Lx**2) + (z-zini)**2/(depthcst**2) 
         dist = ((x-xini)/Lx)**2 + ((y-yini)/Lx)**2 + ((z-zini)/depthcst)**2 
         dtemp = EXP(MAX( -(dist/radius ),limit) )
         IF (dtemp > 1.0e-12) THEN
         temp(k,i,j)=tetaref + EXP(MAX( -(dist/radius ),limit) )
         END IF
       END DO
     END DO
   END DO

   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))


#ifdef key_substance

   ! main initialization relative to the substances
   ! ----------------------------------------------

   ! assignement : input variables for cas_init_subs routine

   nvpc=0; nv_cp=0; nv_sorb=0; nv_fix=0; nv_dri=0; nv_int=0;
   nv_dis=1
   nv_ncp=1
   nvp=1
   name_var(2)='DISS'
   name_var(1)='PART'
   date_startsub=date_start

   CALL cas_init_subs(date_start)   ! nvp must be defined before

   ! over-assignement

   ! substances
   itypv(2)=0
   cv_wat(1,:,liminm1:limaxp2,ljminm1:ljmaxp2) = temp(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   cv_wat(2,:,liminm1:limaxp2,ljminm1:ljmaxp2) = temp(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   long_name_var(2)='nitrate'
   standard_name_var(2)='mole_concentration_of_nitrate_in_sea_water'
   valid_min_var(:)=-3.4e38
   valid_max_var(:)=3.4e38

 
   IF (nvp>0) THEN
      name_filesubs = 'bidon'
      l_filesubs=.FALSE.
      isubs=1
      name_var(isubs)='part_w0'
      long_name_var(isubs)='vase'
      standard_name_var(isubs)='vase'

      itypv(isubs)=4
      ws_free_opt(isubs)=0
      ws_free_min(isubs)=1.e-19_rsh
      ws_free_max(isubs)=1.e-19_rsh


      CALL cas_init_subs_disspart(l_filesubs)

   ! over-assignement
   ! particulate variable with settling velocity = 0 in order to compare results with dissolved variable
      ws_free_min(1)=0._rsh
      ws_free_max(1)=0._rsh


   ENDIF

  CALL_MPI ex_i_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
  CALL_MPI ex_j_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))

#endif


   ! density
   ! -------
   drho = 5.0_rsh
   DO j=ljminm1,ljmaxp2
     DO i=liminm1,limaxp2
       DO k=1,kmax
         z = -f_zu(hx(i,j),ssh(i,j),k,i,j)
         boz(k,i,j) = (drho*z/depthcst) / rhoref 
       END DO
     END DO
   END DO
!print*,'boz',boz(:,10,10)


   ! Miscellaneous
   ! ------------

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

   ! over-assignement

   ! shift of obc for non-linear advection terms
   id(:)=id(:)-ishift_obc_east
   ig(:)=ig(:)+ishift_obc_west
   jh(:)=jh(:)-jshift_obc_north
   jb(:)=jb(:)+jshift_obc_south


   ! Periodic conditions use ig,id,jb,jh which are defined in cas_init_misc

#ifdef key_MPI_2D
   CALL ex_peri_x_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_peri_y_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
   CALL ex_peri_x_rsh(2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_peri_y_rsh(2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#ifdef key_substance
#endif
#else
   DO j=ljminm1,ljmaxp2
     temp(:,ig(j)-1,j) = temp(:,id(j)-2,j)
     temp(:,ig(j),j) = temp(:,id(j)-1,j)
     temp(:,id(j),j)   = temp(:,ig(j)+1,j)
   END DO
   DO i=liminm1,limaxp2
     temp(:,i,jb(i)-1) = temp(:,i,jh(i)-2)
     temp(:,i,jb(i)) = temp(:,i,jh(i)-1)
     temp(:,i,jh(i))   = temp(:,i,jb(i)+1)
   END DO
#ifdef key_substance
   DO isubs=1,nb_var
    DO j=ljminm1,ljmaxp2
      cv_wat(isubs,:,ig(j)-1,j) = cv_wat(isubs,:,id(j)-2,j)
      cv_wat(isubs,:,ig(j),j) = cv_wat(isubs,:,id(j)-1,j)
      cv_wat(isubs,:,id(j),j)   = cv_wat(isubs,:,ig(j)+1,j)
    END DO
    DO i=liminm1,limaxp2
      cv_wat(isubs,:,i,jb(i)-1) = cv_wat(isubs,:,i,jh(i)-2)
      cv_wat(isubs,:,i,jb(i)) = cv_wat(isubs,:,i,jh(i)-1)
      cv_wat(isubs,:,i,jh(i))   = cv_wat(isubs,:,i,jb(i)+1)
    END DO
   END DO
#endif
#endif

!MANQUE PAS DES ECHANGES LA ???


   ! The call to subroutine but be done at the end of the test case
   ! because ig,id,jb,jh... must be defined
   ! --------------------------------------------------------------


   ! initialization of the barotropic velocities (need of hzex)
   ! -------------------------------------------
   DO j = ljmin,ljmax
     DO i = limin,limax
#if defined key_siggen || defined key_gencoord
       uc0(i,j)=SUM(uzc0(:,i,j)*hzex(:,i,j)*dsigu(:))/hex(i,j)
       vc0(i,j)=SUM(vzc0(:,i,j)*hzey(:,i,j)*dsigu(:))/hey(i,j)
#else
       uc0(i,j)=SUM(uzc0(:,i,j)*dsigu(:))
       vc0(i,j)=SUM(vzc0(:,i,j)*dsigu(:))
#endif
     END DO
   END DO
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,uc0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,uc0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,vc0(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,vc0(liminm2:limaxp2,ljminm2:ljmaxp2))



   ! initialization of the vertical velocity (need of ig,id,...,uc0,vc0)
   ! ---------------------------------------
#if defined key_siggen || defined key_gencoord
   fxeau(:,liminm1:limax,ljmin:ljmax)=hzex(:,liminm1:limax,ljmin:ljmax)*uzc0(:,liminm1:limax,ljmin:ljmax)*halfdt
   fyeau(:,limin:limax,ljminm1:ljmax)=hzey(:,limin:limax,ljminm1:ljmax)*vzc0(:,limin:limax,ljminm1:ljmax)*halfdt
#else
   DO k=1,kmax
     fxeau(k,liminm1:limax,ljmin:ljmax)=hex(liminm1:limax,ljmin:ljmax)*uzc0(k,liminm1:limax,ljmin:ljmax)*halfdt
     fyeau(k,limin:limax,ljminm1:ljmax)=hey(limin:limax,ljminm1:ljmax)*vzc0(k,limin:limax,ljminm1:ljmax)*halfdt
   END DO
#endif
!toseelater We can remove uz,vz,u,v in dyn3dwz
#if defined key_siggen || defined key_gencoord
   CALL dyn3d_wz(ssh,ssh,uzc0,vzc0,uc0,vc0)
#else
   CALL dyn3d_wz(ssh,ssh,uzc0,vzc0)
#endif
   wzc0(0:kmax,liminm1:limaxp1,ljminm1:ljmaxp1)=wz(0:kmax,liminm1:limaxp1,ljminm1:ljmaxp1)
#ifdef key_diag_physwz
   wzc0_phys(0:kmax,limin:limax,ljmin:ljmax)=wz_phys(0:kmax,limin:limax,ljmin:ljmax)
#endif


   ! Come back to correct initial conditions
   ! ---------------------------------------
#if defined key_siggen || defined key_gencoord
   fxeau(:,liminm1:limax,ljmin:ljmax)=hzex(:,liminm1:limax,ljmin:ljmax)*uz(:,liminm1:limax,ljmin:ljmax)*halfdt
   fyeau(:,limin:limax,ljminm1:ljmax)=hzey(:,limin:limax,ljminm1:ljmax)*vz(:,limin:limax,ljminm1:ljmax)*halfdt
#else
   DO k=1,kmax
     fxeau(k,liminm1:limax,ljmin:ljmax)=hex(liminm1:limax,ljmin:ljmax)*uz(k,liminm1:limax,ljmin:ljmax)*halfdt
     fyeau(k,limin:limax,ljminm1:ljmax)=hey(limin:limax,ljminm1:ljmax)*vz(k,limin:limax,ljminm1:ljmax)*halfdt
   END DO
#endif


#endif

  END SUBROUTINE cas_advspotmount
