  SUBROUTINE cas_btvortex

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_btvortex  ***
!&E
!&E ** Purpose : barotropic vortex in cyclo-geostrophic balance
!&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-03  (T. Duhaut)  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---------------------------------------------------------------------





























































!! * Modules used
   USE parameters,ONLY : pi,rsh,rlg,riolg,kmax,jmax,imax,jmin,imin, &
                         limin,ljmin,limaxp1,ljmaxp1,ljmaxp2,    &
                         liminp1,limaxp2,liminm2,ljminm2
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,halfdt,suffix,tobs, &
                         l_corio_cst,icon,l_modele2d,vish_xe,vish_phi
   USE comvarp2d, ONLY : ssh,u,v
   USE comvars3d, ONLY : rhoref,kzinit,nzinit,l_stateeq_lin
   USE comvarp3d, ONLY : uz,vz
   USE output,    ONLY : pasor,l_out_sal,l_out_temp
   USE toolmpi,ONLY : ex_i_rsh,ex_j_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

   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k
   INTEGER           :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   REAL(KIND=rsh)    :: dxval,dyval,z0bottom,z0surface,depthcst,fuval
   REAL(KIND=rlg)    :: r,eta0,rto,ito,jto,alpha1,vteta,test

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


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

! dynamics
   CALL cas_init_param

! over-assignement
!l_corio_cst=.FALSE.  ! beta-plane
   l_modele2d=.TRUE.

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

! over-assignement
   icon = 0
   rhoref=1000.0_rsh
   l_stateeq_lin=.FALSE.


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

! assignement : input variables for cas_init_bathy routine
   depthcst = 5000.0_rsh
   ni_wall_e = 1  ! land over ni_wall_e cells at the eastern part of the domain
   ni_wall_w = 2  ! land over ni_wall_w cells (from imin column) at the western 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

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

! over-assignement


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

! assignement : input variables for cas_init_grid routine
   dxval=10000.0_rsh
   dyval=10000.0_rsh
   fisac=38.5_riolg
   gwac=1.0_riolg

   CALL cas_init_grid(dxval,dyval)

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


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

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

   suffix='Vbar'

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

! over-assignement
   pasor(1)=24.0_rlg*3600.0_rlg
   l_out_sal(1)=.FALSE.
   l_out_temp(1)=.FALSE.
   tobs=24.0_rlg*3600.0_rlg


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

! no open boundary
   CALL cas_init_obc


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

! none
   CALL cas_init_visc

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


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

   z0bottom   =1.0e-13
   z0surface  =0.0_rsh
   kzinit=0.0000001_rsh
   nzinit=0.000013_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

! Sea surface height
!-------------------

   eta0      = -0.45_rlg ! Attention si tourbillon anti-cyclonique
! (valeur positive)il faut que eta0 soit
! plus petit ou egale a (rto*f)**2./(8*g)
   rto       = 60000.0_rlg
   ito = REAL(imax+imin,rlg)/2.0_rlg
   jto = REAL(jmax+jmin,rlg)/2.0_rlg
   DO j=ljmin,ljmaxp1
     DO i=limin,limaxp1
       r= SQRT((REAL(i,rlg)-ito)**2+(REAL(j,rlg)-jto)**2)*10000.0_rlg
       ssh(i,j) = eta0*EXP(-(r/rto)**2)
     END DO
   END DO
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,ssh(liminm2:limaxp2,ljminm2:ljmaxp2))

! Velocities
!-----------

   DO j=ljminm2,ljmaxp2
     DO i=liminm2,limaxp2

!   Calcul de U(i,j) Attention Vteta est different selon U ou V
!   r=r(u(i,j))

       r= SQRT((REAL(i,rlg)-ito+0.5_rlg)**2+(REAL(j,rlg)-jto)**2)*10000.0_rlg
       IF ((REAL(i,rlg)-ito+0.5_rlg).NE.0.0_rlg) THEN
         alpha1 = ATAN2(REAL(j,rlg)-jto,REAL(i,rlg)-ito+0.5_rlg)
       ELSEIF ((REAL(j,rlg)-jto).GE.0.0_rlg) THEN
         alpha1 = pi/2.0_rlg
       ELSE
         alpha1 = 3.0_rlg*pi/2.0_rlg
       END IF
       test=(fuval*r)**2/4.0_rlg-2.0_rlg*9.81_rlg*(r/rto)**2.0_rlg*eta0*EXP(-(r/rto)**2)
       IF (test >= 0.0_rlg ) THEN
         vteta=-fuval*r/2.0_rlg                                                                   &
               +SQRT((fuval*r)**2/4.0_rlg-2.0_rlg*9.81_rlg*(r/rto)**2.0_rlg*eta0*EXP(-(r/rto)**2))
       ELSE
         vteta=0.0_rlg
       END IF
       u(i,j)=-SIN(alpha1)*vteta
       IF (ABS(89-i)<=3.AND.ABS(86-j)<=3) THEN
!     print *,i,j,alpha1,vteta
       END IF

!   Calcul de V(i,j) Attention Vteta est different selon U ou V
!   r=r(v(i,j))
       r= SQRT((REAL(i,rlg)-ito)**2+(REAL(j,rlg)-jto+0.5_rlg)**2)*10000.0_rlg
       IF ((REAL(i,rlg)-ito).NE.0.0_rlg) THEN
         alpha1 = ATAN2(REAL(j,rlg)-jto+0.5_rlg,REAL(i,rlg)-ito)
       ELSEIF ((REAL(j,rlg)-jto+0.5_rlg).GE.0.0_rlg) THEN
         alpha1 = pi/2.0_rlg
       ELSE
         alpha1 = 3.0_rlg*pi/2.0_rlg
       END IF
       test=(fuval*r)**2/4.0_rlg-2*9.81_rlg*(r/rto)**2*eta0*EXP(-(r/rto)**2)
       IF (test >= 0.0_rlg) THEN
         vteta=-fuval*r/2.0_rlg                                                             &
               +SQRT((fuval*r)**2/4.0_rlg-2.0_rlg*9.81_rlg*(r/rto)**2*eta0*EXP(-(r/rto)**2))
       ELSE
         vteta=0.0_rlg
       END IF
       v(i,j)=COS(alpha1)*vteta
     END DO
   END DO

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


! initialization of vertical scale factor (part relative to dc*sds)
!------------------------------------------------------------------

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

  END SUBROUTINE cas_btvortex
