  SUBROUTINE cas_bt_current

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_advcst  ***
   !&E
   !&E ** Purpose : definition of a test case with constant zonal and baroctropic current
   !&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 ** Reference : 
   !&E
   !&E ** History :
   !&E       !  2006     (F. Vandermeirsch) module -Original code
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&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,  &
                          ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north, &
                          valmanq
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,halfdt,suffix,qmax,l_modele2d, &
                         h0fond,hx_g,hy_g,h0_g,hx,hy,h0,icon,jh,jb,ig,id
   USE comvarp2d, ONLY : sal,temp,ssh,u,v
   USE comvarp3d, ONLY : uz,vz
   USE comvars3d, ONLY : l_stateeq_lin,tetaref_lin,saliref_lin,kzinit,nzinit, &
                         saliref,tetaref,rhoref,qmaxz
   USE output,    ONLY : pasor
   USE comcfo,    ONLY : xext,uext,vext
   USE obccombine,ONLY : l_obc_south,l_obc_north,l_obc_west,l_obc_east, &
                          obc_ue,obc_uw,obc_us,obc_un,obc_uze,obc_uzw,obc_uzs,obc_uzn, &
                          obc_ve,obc_vw,obc_vs,obc_vn,obc_vze,obc_vzw,obc_vzs,obc_vzn, &
                          obc_sshe,obc_sshw,obc_sshs,obc_sshn,l_obc_char
   USE obcogcm,   ONLY : obc_tw,obc_se,obc_te,obc_ss,obc_ts,obc_sn,obc_tn,obc_sw,obc_tsref
   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

   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k,test_bathy
   INTEGER,PARAMETER :: N1=10
   REAL(KIND=rsh)    :: uinit,vinit,hcst,ubar1,ubar,vbar,theta_C,fact_perturb
   REAL(KIND=rsh)    :: f0,s,A2,B2,B3,coef,L1,ha,hp,fuval
   REAL(KIND=rsh)    :: d,radi,x,y,tetaref_perturbation
   REAL(KIND=rsh)    :: dxval,dyval,z0bottom,z0surface,depthcst

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

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


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

   ! dynamics
   CALL cas_init_param

   ! over-assignement
   !l_corio_cst=.FALSE.  ! beta-plane
   l_modele2d=.FALSE.
   h0fond=0.01_rsh

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

   ! over-assignement
   qmax=12.0_rsh
   qmaxz=10.0_rsh
   rhoref=1024.0_rsh
   saliref=31.0_rsh
   IF (icon==2) l_stateeq_lin=.TRUE.
   saliref_lin=saliref
   tetaref_lin=tetaref


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

   !-----------------------------------------------------
   ! Choix du test
   ! Sans bathy : test_bathy = 0 
   ! Avec bathy : test_bathy = 1

   test_bathy = 0

   ! clef a mettre en 3D key_no_bottom_stress
   IF (test_bathy == 0) THEN
     !-----------------------------------------------------
     ! Parametres specifiques au cas test courant_barotrope
     !-----------------------------------------------------

     WRITE(*,*)'####################################'
     WRITE(*,*)'        TEST CASE n 11.1            '
     WRITE(*,*)'       barotropic current           '
     WRITE(*,*)'####################################'
     WRITE(*,*)'icon = ',icon 

     ! bathymetrie
     !-------------

     hcst = 4000.0_rsh
     hx_g(:,:)=hcst
     hy_g(:,:)=hcst
     h0_g(:,:)=hcst

   ELSE

     WRITE(*,*)'####################################'
     WRITE(*,*)'          TEST CASE n 11.2           '
     WRITE(*,*)' barotropic current with topography '
     WRITE(*,*)'####################################'
     WRITE(*,*)'icon = ',icon


     hp = -200.0_rsh;
     s =0.05_rsh
     L1 = N1*dyval
     ha = -s*L1+ hp

     !----------------------
     ! premier domaine : 0 a N1 -1 , 0 a L1-dy 
     !----------------------
     DO j=jmin,N1-1
       DO i=imin,imax
         hx_g(i,j) = -ha
         hy_g(i,j) = hx_g(i,j)
         h0_g(i,j) = hx_g(i,j)
       END DO
     END DO
     !----------------------      
     ! deuxieme domaine : N1 a 2N1-1 : L1 a 2L1-dy
     !----------------------
     DO j=N1,2*N1-1
       DO i=imin,imax
         y = (j-jmin)*dyval
         hx_g(i,j) =  -(s*y+2.0_rsh*ha-hp)
         hy_g(i,j) =  hx_g(i,j)
         h0_g(i,j) =  hx_g(i,j)
       END DO
     END DO
     !----------------------
     ! troisieme domaine : 2N1 a 2N1+ N3 -1 : 2L1 a L
     !----------------------
     DO j=2*N1,jmax
       DO i=imin,imax
         hx_g(i,j) = -hp
         hy_g(i,j) = hx_g(i,j)
         h0_g(i,j) = hx_g(i,j)
       END DO
     END DO

   END IF
   h0(:,:) = h0_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hx(:,:) = hx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hy(:,:) = hy_g(liminm2:limaxp2,ljminm2:ljmaxp2)


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

   ! assignement : input variables for cas_init_grid routine
   dxval=20000.0_rsh
   dyval=20000.0_rsh
   fisac=0.0_riolg
   gwac=0.0_riolg

   CALL cas_init_grid(dxval,dyval)

   ! over-assignement
#ifdef key_symetrical_uv 
   fuval=0_rsh
   fu(:,:)=0_rsh
   fv(:,:)=0_rsh
#else
   fuval=0.0001_rsh
   fu(:,:)=0.0001_rsh
   fv(:,:)=0.0001_rsh
#endif  


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

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

   suffix='courant'

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

   ! over-assignement
   pasor(1)= 24.0_rlg*3600.0_rlg ! 24*3600


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

   ! none
   CALL cas_init_visc


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

   z0bottom   =1.0e-13
   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

   IF (test_bathy == 0) THEN
     !-----------------------------------------------------
     ! Parametres specifiques au cas test courant_barotrope

     ! grandeurs dynamiques
     !----------------------
     ! L angle theta_C par rapport a l axe des x
     ! permet d obtenir un courant oblique.
     ! theta_C variant entre 0 et 2*pi
     ! Dans le cas d un courant oblique
     ! theta_C # 0 ou # pi/2 (zonal ou meridien)
     ! alors les 4 frontieres sont ouvertes
     ! possibilite de génraliser avec des murs 
     ! et une bathy non-constante
     
       theta_C = pi/18. ! + pi
     
     ! Solution exterieure 
       ubar = 1.0_rsh*cos(theta_c)
       vbar = 1.0_rsh*sin(theta_c)
       ! perturbation de la condition initiale pour voir si le champ final
       ! converge vers la solution (donnee par conditions aux limites non perturbees)
       fact_perturb = 0.9_rsh  
       
       WRITE(*,*) 'Courant barotrope faisant un angle de'
       WRITE(*,*) theta_C*180/pi,"degre avec l'axe des x"
       WRITE(*,*) 'le facteur de perturbation est',fact_perturb
          
     ! vitesse initiale 
       uinit = ubar *fact_perturb
       vinit = vbar *fact_perturb
          
     ! Solution exterieure
     !--------------------    
     DO j=jmin,jmax
       DO i=imin,imax
         xext(i,j) = -fuval/grav*(ubar*(j-jshift_obc_south-jmin)*dyval - vbar*(i-ishift_obc_west-imin)*dyval)
       END DO
     END DO

     !------------------------------------------------------------------------------  
     ! Condition aux limites  
     ! devrait etre obc_ogcm_ssh ou obc_tide_ssh 
     obc_sshe(:)=valmanq
     obc_sshw(:)=valmanq
     obc_sshn(:)=valmanq
     obc_sshs(:)=valmanq
     DO j=MAX(jshift_obc_south+1,ljmin),MIN(jmax-jshift_obc_north,ljmax)
       obc_sshe(j)=xext(imax-ishift_obc_east,j)
       obc_sshw(j)=xext(imin+1+ishift_obc_west,j)
     END DO
     DO i=MAX(ishift_obc_west+1,limin),MIN(imax-ishift_obc_east,limax)
       obc_sshn(i)=xext(i,jmax-jshift_obc_north)
       obc_sshs(i)=xext(i,jmin+1+jshift_obc_south)
     END DO

     ! meme remarque 
     obc_ue(:,:) = ubar
     obc_uw(:,:) = ubar
     obc_us(:,:) = ubar
     obc_un(:,:) = ubar
     IF (.NOT. l_modele2d) THEN
       obc_uze(:,:,:) = ubar
       obc_uzw(:,:,:) = ubar
       obc_uzs(:,:,:) = ubar
       obc_uzn(:,:,:) = ubar
     END IF

     obc_ve(:,:) = vbar
     obc_vw(:,:) = vbar
     obc_vs(:,:) = vbar
     obc_vn(:,:) = vbar
     IF (.NOT. l_modele2d) THEN
       obc_vze(:,:,:) = vbar
       obc_vzw(:,:,:) = vbar
       obc_vzs(:,:,:) = vbar
       obc_vzn(:,:,:) = vbar
     END IF

     ! initialisation 
     !----------------------   
     ssh(liminm2:limaxp2,ljminm2:ljmaxp2) = xext(liminm2:limaxp2,ljminm2:ljmaxp2)*fact_perturb

     u(:,:)    = uinit
     v(:,:)    = vinit 

     uz(:,:,:) = uinit
     vz(:,:,:) = vinit

   ELSE
     !-----------------------------------------------------
     ! Parametres specifiques au cas test courant_barotrope avec bathy selon y
     !-----------------------------------------------------

     WRITE(*,*)'####################################'
     WRITE(*,*)'          TEST CASE n 11.2           '
     WRITE(*,*)' barotropic current with topography '
     WRITE(*,*)'####################################'
     WRITE(*,*)'icon = ',icon 

     ! grandeurs dynamiques
     !----------------------

     ! Solution exterieure 

     s =0.05_rsh
     ubar1 = 0.1_rsh
     vbar = 0.0_rsh
     L1 = N1*dyval
     f0 = fuval
     A2 = f0*s*L1 - ubar1
     B2 = -f0*s*L1*L1/2
     B3 = 3.0_rsh/2.0_rsh*f0*s*L1*L1
     coef = f0/grav

     !----------------------
     ! premier domaine : 0 a N1 -1 , 0 a L1-dy 
     !----------------------
     DO j=jmin,N1-1
       DO i=imin,imax
         y=(j-jmin)*dyval
         uext(i,j) = ubar1
         xext(i,j) = -coef*ubar1*y
       END DO
     END DO
     !----------------------      
     ! deuxieme domaine : N1 a 2N1-1 : L1 a 2L1-dy
     !----------------------
     DO j=N1,2*N1-1
       DO i=imin,imax
         y = (j-jmin)*dyval
         uext(i,j) = f0*s*y - A2;
         xext(i,j) = coef*(-f0*s/2.0_rsh*y*y + A2*y  +B2)
       END DO
     END DO
     !----------------------
     ! troisieme domaine : 2N1 a 2N1+ N3 -1 : 2L1 a L
     !----------------------
     DO j=2*N1,jmax
       DO i=imin,imax
         y=(j-jmin)*dyval
         uext(i,j) = ubar1+f0*s*L1
         xext(i,j) = coef*(-uext(i,j)*y + B3)
       END DO
     END DO

     !------------------------------------------------------------------------------  
     ! Condition aux limites 
     obc_sshe(ljminm2:ljmaxp2)=xext(imax,ljminm2:ljmaxp2)
     obc_sshw(ljminm2:ljmaxp2)=xext(imin+1,ljminm2:ljmaxp2)
     obc_sshn(liminm2:limaxp2)=xext(liminm2:limaxp2,jmax)
     obc_sshs(liminm2:limaxp2)=xext(liminm2:limaxp2,jmin+1)

     obc_ue(imax-obc_width+1:imax,ljmin:ljmax) = uext(imax-obc_width+1:imax,ljmin:ljmax)
     obc_uw(imin:imin+obc_width,ljmin:ljmax) = uext(imin:imin+obc_width,ljmin:ljmax)
     obc_us(limin:limax,jmin:jmin+obc_width) = uext(limin:limax,jmin:jmin+obc_width)
     obc_un(limin:limax,jmax-obc_width+1:jmax) = uext(limin:limax,jmax-obc_width+1:jmax)

     IF (.NOT. l_modele2d) THEN
       DO k=1,kmax
         obc_uze(k,imax-obc_width+1:imax,ljmin:ljmax)    = uext(imax-obc_width+1:imax,ljmin:ljmax) 
         obc_uzw(k,imin:imin+obc_width,ljmin:ljmax)      = uext(imin:imin+obc_width,ljmin:ljmax)
         obc_uzs(k,limin:limax,jmin:jmin+obc_width)      = uext(limin:limax,jmin:jmin+obc_width)
         obc_uzn(k,limin:limax,jmax-obc_width+1:jmax)    = uext(limin:limax,jmax-obc_width+1:jmax) 
       END DO
     END IF

     obc_ve(:,:) = vbar
     obc_vw(:,:) = vbar
     obc_vs(:,:) = vbar
     obc_vn(:,:) = vbar

     IF (.NOT. l_modele2d) THEN
       obc_vze(:,:,:) = vbar
       obc_vzw(:,:,:) = vbar
       obc_vzs(:,:,:) = vbar
       obc_vzn(:,:,:) = vbar
     END IF

     ! initialisation 
     !----------------------   

     ssh(liminm2:limaxp2,ljminm2:ljmaxp2) = xext(liminm2:limaxp2,ljminm2:ljmaxp2) 
     ! sans perturbation : a t=0 uinit = solution analytique
     u(liminm2:limaxp2,ljminm2:ljmaxp2)    = uext(liminm2:limaxp2,ljminm2:ljmaxp2)
     ! avec perturbation : a t=0 uinit = 0
     !  u(:,:)    = 0 

     DO k=1,kmax
       uz(k,:,:) = u(:,:)
     END DO

   ENDIF

   ! initialisation 
   !----------------------  
   tetaref_perturbation = 0
   sal(:,:,:)  = saliref
   temp(:,:,:) = tetaref+tetaref_perturbation


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

   CALL cas_init_obc

   ! over-assignement
   ! geometry (where are the open boundaries ?)
   l_obc_south=.TRUE.
   l_obc_north=.TRUE.
   l_obc_west=.TRUE.
   l_obc_east=.TRUE.
   ! used scheme to manage obc
   l_obc_char = .TRUE.
   ! values at the open boundaries
   CALL obc_tsref

   obc_te(:,:,:) = tetaref+tetaref_perturbation
   obc_tw(:,:,:) = tetaref+tetaref_perturbation
   obc_ts(:,:,:) = tetaref+tetaref_perturbation
   obc_tn(:,:,:) = tetaref+tetaref_perturbation
   obc_se(:,:,:) = saliref
   obc_sw(:,:,:) = saliref
   obc_ss(:,:,:) = saliref
   obc_sn(:,:,:) = saliref

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

   ! over-assignement
   id(:)=id(:)-ishift_obc_east
   ig(:)=ig(:)+ishift_obc_west
   jh(:)=jh(:)-jshift_obc_north
   jb(:)=jb(:)+jshift_obc_south

  END SUBROUTINE cas_bt_current
