  SUBROUTINE cas_thacker

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_thacker  ***
   !&E
   !&E ** Purpose : definition of a test case for a thacker
   !&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 ** Reference :
   !&E
   !&E ** History :
   !&E       !  2006-01  (F. Dumas, V. Garnier) Original code
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters,ONLY : rsh,rlg,riolg,imin,imax,jmin,jmax,kmax,  &
                         liminm1,limaxp2,ljminm1,ljmaxp2,liminm2,ljminm2,    &
                         valmanq,pi,limin,limax,ljmin,ljmax,limaxm1,ljmaxm1,grav
   USE comvars2d, ONLY : fisac,gwac,dtmin,dtmax,dt,halfdt,t,suffix,   &
                         qmax,h0,tfiusr,tfivsr,h0fond,radearth, &
                         hx_g,hy_g,h0_g,fu,h0,hx,hy,dhx,dhy, &
                         dhx_g,dhy_g,icon,l_modele2d,fwetmin
   USE comvarp2d, ONLY : ssh,u,v,temp
   USE comvars3d, ONLY : qmaxz,nzinit,kzinit,l_stateeq_lin,rhoref,saliref, &
                         tetaref,saliref_lin,tetaref_lin,sig
   USE output, ONLY : pasor,l_out_sal,l_out_temp
!   USE cometemp
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : theta_sig,b_sig,hc
#endif
   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

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k,icentre,jcentre
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst
   REAL(KIND=rsh)    :: L,hmax,dxp,dyp,dxval,dyval,eta,omega,fcorio

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


   PRINT_DBG*, 'enter CAS_THACKER routine'

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

   ! dynamics
   CALL cas_init_param

   ! over-assignement
   !h0fond=0.01_rsh
   !hm=0.3_rsh
   fwetmin=0.1


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

   ! over-assignement
   icon=0
   qmax=12.0_rsh
   qmaxz=10.0_rsh
   rhoref=1000.0_rsh
   saliref=35.0_rsh
   tetaref=5.0_rsh
   l_stateeq_lin=.FALSE.
   saliref_lin=0.0_rsh
   tetaref_lin=0.0_rsh


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

   ! assignement : input variables for cas_init_bathy routine
   depthcst = 0.0_rsh
   ni_wall_e = 5  ! land over ni_wall_e cells at the eastern part of the domain
   ni_wall_w = 4  ! land over ni_wall_w cells (from imin column) at the western part of the domain
   nj_wall_s = 4  ! land over nj_wall_s cells (from jmin row) at the southern part of the domain
   nj_wall_n = 5  ! land over nj_wall_n cells at the northern part of the domain
#if defined key_siggen || defined key_gencoord
   hc= 5.0_rsh
#endif


   ! As depthcst is null, specify your own bathymetry

   dxval = 1000.0_rsh
   dyval = 1000.0_rsh
   hmax=10.0_rsh
   L=80000.0_rsh
   icentre = 100  ! mars starts at imin=0
   jcentre = 100
   DO i=imin,imax
     DO j=jmin,jmax
       dxp=dxval*(REAL(i-icentre,rsh)+0.5_rsh)/L
       dyp=dyval*(REAL(j-jcentre,rsh))/L
       hx_g(i,j)=hmax*( 1.0_rsh - dxp**2 - dyp**2)
       dxp=dxval*(REAL(i-icentre,rsh))/L
       dyp=dyval*(REAL(j-jcentre,rsh)+0.5_rsh)/L
       hy_g(i,j)=hmax*( 1.0_rsh - dxp**2 - dyp**2)
     END DO
   END DO
   WHERE (hx_g(:,:) < -5.0_rsh) hx_g(:,:)=-valmanq-h0fond
   WHERE (hy_g(:,:) < -5.0_rsh) hy_g(:,:)=-valmanq-h0fond

   DO j=jmin+1,jmax
     DO i=imin+1,imax
       h0_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))+h0fond
     END DO
   END DO
   DO j=jmin+1,jmax
     h0_g(imin,j)=MAX(hx_g(imin,j),hy_g(imin,j),hy_g(imin,j-1))+h0fond
   END DO
   DO i=imin+1,imax
     h0_g(i,jmin)=MAX(hx_g(i,jmin),hx_g(i-1,jmin),hy_g(i,jmin))+h0fond
   END DO
   h0_g(imin,jmin)=MAX(hx_g(imin,jmin),hy_g(imin,jmin))+h0fond


   ! walls
   CALL cas_init_bathy(depthcst,ni_wall_w=ni_wall_w,ni_wall_e=ni_wall_e, &
                                nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s)

   ! over-assignement
   dhx(:,:)=0.002_rsh
   dhy(:,:)=0.002_rsh
   dhx_g(:,:)=0.002_rsh
   dhy_g(:,:)=0.002_rsh
   dhx(:,:)=0.2_rsh
   dhy(:,:)=0.2_rsh
   dhx_g(:,:)=0.2_rsh
   dhy_g(:,:)=0.2_rsh


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

   ! assignement : input variables for cas_init_grid routine
 
   fisac=ASIN(0.0001_rsh/0.0001458_rsh)/(pi/180.0_rsh)  ! lat to get corio=10-4
   gwac=1.0_riolg

#if defined key_siggen || defined key_gencoord
   theta_sig= 4.0_rsh ! entre 0 et 20
   theta_sig= 0.0000001_rsh ! entre 0 et 20
   b_sig= 1.0_rsh     ! entre 0 et 1
#endif

   CALL cas_init_grid(dxval,dyval)

   ! over-assignement
   !tfiusr(:,:)=TAN(fisac*pi/180.0_rsh)/radearth
   !tfivsr(:,:)=TAN(fisac*pi/180.0_rsh)/radearth


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

   ! assignement : input variables for cas_init_io routine
   date_start=   '01/01/1900 00:00:00'
   date_end=   '06/01/1900 00:06:00'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=20.0_rlg
   dtmax=20.0_rlg
   dt=20.0_rlg
   halfdt=dt*0.5_rlg

   suffix='thacker'

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

   ! over-assignement
   l_modele2d=.TRUE.
   pasor(1)=880.0_rlg
   l_out_sal(1)=.FALSE.
   l_out_temp(1)=.FALSE.


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

   ! no open boundary
   CALL cas_init_obc


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

   ! none
   CALL cas_init_visc

   ! over-assignement


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

   z0bottom   =0.0_rsh
   z0surface  =0.0_rsh
   nzinit=0.0_rsh
   kzinit=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

   fcorio = fu(limin,ljmax)
   omega = fcorio*.5 + SQRT((fcorio**2)*.25_rsh+2.0_rsh*grav*hmax/(L**2))

   ! sea surface height
   !-------------------

   eta = 8000.0_rsh ! 8 km, du meme ordre de grandeur que L sinon on ne voit rien
   DO j=ljminm2,ljmaxp2
     DO i=liminm2,limaxp2
       dxp=dxval*(REAL(i-icentre,rsh))/L
       dyp=dyval*(REAL(j-jcentre,rsh))/L
       ssh(i,j) = MAX( 2.0_rlg*eta*hmax/L * ( dxp*cos(omega*t)-dyp*sin(omega*t)-eta/(2.0_rsh*L) ), REAL(-h0(i,j),rlg))
       u(i,j) = - eta*omega*sin(omega*t)
       v(i,j) = - eta*omega*cos(omega*t)
     END DO
   END DO

   WHERE(h0(liminm2:limaxp2,ljminm2:ljmaxp2) == -valmanq) &
        ssh(liminm2:limaxp2,ljminm2:ljmaxp2)=valmanq
!   WHERE(h0(liminm2:limaxp2,ljminm2:ljmaxp2) + ssh(liminm2:limaxp2,ljminm2:ljmaxp2) <= 0.0_rsh) &
!        ssh(liminm2:limaxp2,ljminm2:ljmaxp2)=-h0(liminm2:limaxp2,ljminm2:ljmaxp2)
   WHERE(hx(liminm2:limaxp2,ljminm2:ljmaxp2) == -valmanq-h0fond) &
        u(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh
   WHERE(hy(liminm2:limaxp2,ljminm2:ljmaxp2) == -valmanq-h0fond) &
        v(liminm2:limaxp2,ljminm2:ljmaxp2)=0.0_rsh
! comment initialise t on les variables sur bancs decouvrants ?
   DO j=ljminm2,ljmaxm1
     DO i=liminm2,limaxm1
       IF (hx(i,j)+0.5_rsh*(ssh(i,j)+ssh(i+1,j)) <= 0.0_rsh) u(i,j)=0.0_rsh
       IF (hy(i,j)+0.5_rsh*(ssh(i,j)+ssh(i,j+1)) <= 0.0_rsh) v(i,j)=0.0_rsh
     END DO
   END DO

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


   ! tracers
   !--------

#ifdef key_virer
   DO j=ljminm1,ljmaxp2
     DO i=liminm1,limaxp2
       IF (h0(i,j) > -valmanq) THEN
         DO k=1,kmax
#if defined key_siggen || defined key_gencoord
           zr=hc_sig_g(i,j)*(sig(k)-1.0_rsh)+(h0(i,j)-hc_sig_g(i,j))*csu_sig(k)
           zr=zr*(-1)
#else
           zr=-sig(k)*h0(i,j)
#endif
           temp(k,i,j)= 5.0_rsh + 15.0_rsh*EXP(-zr/1000.0_rsh)
         END DO
       ELSE
         temp(:,i,j) = valmanq
       END IF
     END DO
   END DO
#endif


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

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


   PRINT_DBG*, 'exit CAS_THACKER routine'

  END SUBROUTINE cas_thacker
