  SUBROUTINE cas_bcvortex

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_bcvortex  ***
!&E
!&E ** Purpose : baroclinic vortex like Penven s one
!&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 : Penven P.,  Ocean modeling 12 (2006) 157-187
!&E
!&E ** History :
!&E       !  2007-03  (T. Duhaut)  Original code
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters,ONLY : rsh,rlg,riolg,kmax,jmax,imax,jmin,imin,             &
                         grav,limin,ljmin,ljmaxm1,limaxp1,ljmaxp1,ljminm1,   &
                         ljmaxp2,liminm1,liminp1,limaxp2,limaxm1,ljminp1,    &
                         liminm2,ljminm2,pi
   USE comvars2d, ONLY : fisac,gwac,dtmin,dtmax,dtmaxi,dt,halfdt,suffix, &
                         botstressmax,l_corio_cst,qmax,h0

   USE comvarp2d, ONLY : sal,temp,ssh,u,v
   USE comvars3d, ONLY : rhoref,saliref,kzinit,nzinit,qmaxz,sig,dsigu
   USE comvarp3d, ONLY : uz,vz
   USE output,    ONLY : pasor,l_out_sal

   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

   USE toolgeom,   ONLY : f_zu

   IMPLICIT NONE

!! * Arguments

!! * Local declarations
   LOGICAL           :: l_nc4par
   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
   INTEGER           :: ixl,iyl
   REAL(KIND=rsh)    :: dxval,dyval
   REAL(KIND=rlg)    :: Fzr,rhoz0,XYEXP,icentre,jcentre,Nfreq,zr,zru,zrv,disticentre,distjcentre
   REAL(KIND=rlg)    :: XYEXPx,XYEXPy,icentrex,icentrey,jcentrex,jcentrey,P0ref
   REAL(KIND=rlg)    :: disticentrex,disticentrey,distjcentrex,distjcentrey
   REAL(KIND=rlg)    :: lambda
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,limit,f0,xeu,umax,h1
   REAL(KIND=rlg)    :: rhox,rhoy

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

! 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_corio_cst=.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_rsh
   qmaxz=10.0_rsh
   rhoref=1024.4_rsh
   saliref=31.0_rsh


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

! assignement : input variables for cas_init_bathy routine
   depthcst = 5000.0_rsh
   ni_wall_e = 2  ! 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 = 2  ! 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
!dxval=20000.0_rsh
!dyval=20000.0_rsh
!dxval=30000.0_rsh
!dyval=30000.0_rsh
!dxval=5000.0_rsh
!dyval=5000.0_rsh
   fisac=38.5_riolg
   gwac=1.0_riolg





   CALL cas_init_grid(dxval,dyval)

! over-assignement
   

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

! assignement : input variables for cas_init_io routine
   date_start=   '01/01/1900 00:00:00'
   date_end=   '10/04/1900 00:00:00'

!   date_end=   '08/01/1900 00:00:00'    ! mh-20141125
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=960.0_rlg
   dtmax=960.0_rlg
   dt=960.0_rlg
!dtmin=480.0_rlg
!dtmax=480.0_rlg
!dt=480.0_rlg

   !IF (.NOT..TRUE.) THEN
     ! dtmin=Agrif_Parent(dtmin)/(Agrif_rhot())
     ! dtmax=Agrif_Parent(dtmax)/(Agrif_rhot())
     ! dt=Agrif_Parent(dt)/(Agrif_rhot())
   !ENDIF
   halfdt=dt*0.5_rlg

   suffix='Vortex'
   l_nc4par = .TRUE.

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

! over-assignement
   pasor(1) = 3.0_rlg*24.0_rlg*3600.0_rlg
!   pasor(1) = 24.0_rlg*3600.0_rlg    ! mh-20141125 : 1day
!   pasor(1) = 6.0_rlg*3600.0_rlg    ! mh-20141125 : 4 times a day
   l_out_sal(1)=.FALSE.
!tobs=24.0_rlg*3600.0_rlg


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

! no open boundary
   CALL cas_init_obc


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

! none
   CALL cas_init_visc


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

   z0bottom   =1.0e-13
   z0surface  =0.0_rsh
   kzinit=0.0000001_rsh
   nzinit=0.000013_rsh
   nzinit=0.0_rsh

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

! over-assignement
   botstressmax=1.0_rsh


! 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 (.TRUE.) THEN

   h1 = 2500.0_rsh
   umax = 1.0_rsh
   f0 = 2.0_rsh * 2.0_rsh*pi/(3600.0_rsh*24.0_rsh) * SIN(fisac*pi/180.0_rsh)
   lambda=60000.0_rlg
!P0ref=9127.4_rlg      ! Si Po>0 -> Vortex anticyclonic, et inversement!
   P0ref = rhoref * f0 * umax * lambda * SQRT(EXP(1.0_rlg))   ! rho0 * f0 * umax * lambda  *sqrt(e)

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

! domain centre
! (imax-imin+1) total number of points
! -1 because we start at 0 instead of 1
   icentre = REAL(imax-imin+1,rsh)/2.0_rsh*dxval
   jcentre = REAL(jmax-jmin+1,rsh)/2.0_rsh*dyval
   icentre = (REAL(imax-imin+1,rsh)/2.0_rsh - 1.0_rsh)*dxval
   jcentre = (REAL(jmax-jmin+1,rsh)/2.0_rsh - 1.0_rsh)*dyval
! shift of 0.5 to get the f - location
   icentre = icentre+0.5_rsh*dxval
   jcentre = jcentre+0.5_rsh*dyval
   DO j=ljmin,ljmaxp1
     distjcentre = REAL(j,rsh)*dyval-jcentre
     DO i=limin,limaxp1
       disticentre = REAL(i,rsh)*dxval-icentre
       XYEXP=EXP(MAX(-(disticentre**2+distjcentre**2)/(2.0_rsh*lambda*lambda),REAL(limit,rlg)))
       rhoz0=rhoref -(P0ref/grav)*XYEXP/(h1-1.0_rsh)
       ssh(i,j)=REAL((P0ref/grav)*XYEXP/(rhoz0),rsh)
     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=ljminp1,ljmaxm1
     distjcentre  = REAL(j,rsh)*dyval-jcentre
     distjcentrex = distjcentre
     distjcentrey = (REAL(j,rlg)+0.5_rlg)*dyval-jcentre
     DO i=liminp1,limaxm1
       disticentre  = REAL(i,rsh)*dxval-icentre
       disticentrex = (REAL(i,rlg)+0.5_rlg)*dxval-icentre
       disticentrey = disticentre
       XYEXP  = EXP(MAX(-((disticentre)**2+(distjcentre)**2)/(2.0_rsh*lambda*lambda),REAL(limit,rlg)))
       XYEXPx = EXP(MAX(-((disticentrex)**2+(distjcentrex)**2)/(2.0_rlg*lambda*lambda),REAL(limit,rlg)))
       XYEXPy = EXP(MAX(-((disticentrey)**2+(distjcentrey)**2)/(2.0_rlg*lambda*lambda),REAL(limit,rlg)))
       DO k=kmax,1,-1
         xeu = (ssh(i,j)+ssh(i+1,j))/2.0_rsh
         zru = f_zu(h0(i,j),xeu,k,i,j)
         xeu = (ssh(i,j)+ssh(i,j+1))/2.0_rsh
         zrv = f_zu(h0(i,j),xeu,k,i,j)
         IF (zru <= -h1) THEN
           uz(k,i,j) = 0.0_rsh
         ELSE
           Fzr=((h1-1.0_rsh)+zru+EXP(MAX((-zru-h1),REAL(limit,rlg))))/(h1-1.0_rsh)/(lambda*lambda)  ! F(z)/lambda**2
!uz(k,i,j) = REAL((P0ref/rhoref/fu(i,j))*Fzr*distjcentrex*XYEXPx,rsh)
           uz(k,i,j) = REAL((P0ref/rhoref/f0)*Fzr*distjcentrex*XYEXPx,rsh)
         END IF
         IF (zrv <= -h1) THEN
           vz(k,i,j) = 0.0_rsh
         ELSE
           Fzr=((h1-1.0_rsh)+zrv+EXP(MAX((-zrv-h1),REAL(limit,rlg))))/(h1-1.0_rsh)/(lambda*lambda)  ! F(z)/lambda**2
!vz(k,i,j) =-REAL((P0ref/rhoref/fv(i,j))*Fzr*disticentrey*XYEXPy,rsh)
           vz(k,i,j) =-REAL((P0ref/rhoref/f0)*Fzr*disticentrey*XYEXPy,rsh)
         END IF
       END DO

! barotropic velocities estimated after the call of cas_init_misc

     END DO
   END DO

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


! Tracers
!--------

   Nfreq=0.003_rsh
   sal(:,:,:)=saliref
   DO j=ljminm1,ljmaxp2
     distjcentre = REAL(j,rsh)*dyval-jcentre
     DO i=liminm1,limaxp2
       disticentre = REAL(i,rsh)*dxval-icentre
       DO k=1,kmax
         zr = f_zu(h0(i,j),ssh(i,j),k,i,j)
         IF (zr <= -h1) THEN
           temp(k,i,j) = (1030.0_rsh -rhoref*(1.0_rsh-zr*(Nfreq**2)/grav))/0.28_rsh
         ELSE
           XYEXP=EXP(MAX(-((disticentre)**2+(distjcentre)**2)/(2.0_rsh*lambda*lambda),REAL(limit,rlg)))
           temp(k,i,j)= (1030.0_rsh -rhoref*(1.0_rsh-zr*(Nfreq**2)/grav))/0.28_rsh &
             + ((P0ref/grav)*(1.0_rsh-EXP(MAX((-h1-zr),REAL(limit,rlg))))*XYEXP/(h1-1.0_rsh))/0.28_rsh 
         END IF
       END DO
     END DO
   END DO

   !ENDIF


! Trajectories
! ------------



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

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

! barotropic velocities require the estimate of hzex for generalized sigma

   DO j=ljminp1,ljmaxm1
     DO i=liminp1,limaxm1

       u(i,j)=SUM(uz(:,i,j)*dsigu(:))
       v(i,j)=SUM(vz(:,i,j)*dsigu(:))


     END DO
   END DO

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

  END SUBROUTINE cas_bcvortex
