   !!======================================================================

  SUBROUTINE cas_upwelling

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_upwelling  ***
   !&E
   !&E ** Purpose : define the upwelling test case
   !&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       !  2005-10-13  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,grav,pi,imin,imax,jmin,jmax, &
                          limin,liminp1,limax,ljmin,ljminp1,ljminp2,ljmax, &
                          ljmaxm1,ljminm2,ljmaxp2,liminm2,limaxp2, &
                          limaxp1,ljmaxp1,liminm1,ljminm1,valmanq
   USE comvars2d, ONLY : fisac,gwac,fu,fv,dtmin,dtmax,dt,halfdt,suffix,l_obc_cycl_x, &
                         l_obc_cycl_y,icon,qmax,t,h0, &
                         h0_g,hx_g,hy_g,vish_phi,vish_xe,frot_lin, &
                         vismin,vismax,cosmag,l_smagor
   USE comvarp2d, ONLY : ssh,temp
   !USE comvarp3d, ONLY : uz,vz
   USE comvars3d, ONLY : kzinit,nzinit,qmaxz,saliref,rhoref,nz,l_stateeq_lin, &
                         turb_0eq_option,bgdiff
   USE obccombine, ONLY : l_obc_cycl
   USE_MPI toolmpi, ONLY : ex_i_rsh,ex_j_rsh,ex_peri_x_rsh,ex_peri_y_rsh
   USE output,    ONLY : pasor,l_out_sal,l_out_wind,l_out_turb,l_out_zaxis_w
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : theta_sig,b_sig,hc
#endif
   USE toolgeom,  ONLY : f_zw,f_zu
   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 cometdyn,  ONLY : wind_veloc,windamp_duration


   IMPLICIT NONE

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=19) :: date_start,date_end,date_start3d,date_startagrif
   INTEGER           :: i,j,k
   INTEGER           :: nj_wall_n,nj_wall_s
   REAL(KIND=rsh)    :: dxval,dyval
   REAL(KIND=rsh)    :: y,z,zsurf
   REAL(KIND=rsh)    :: z0bottom,z0surface,depthcst,depth
   REAL(KIND=rsh)    :: Length_ETA,Length_XI,cff
   REAL(KIND=rsh)    :: Ekdepth,zd0,zd1,strat,thknss
   REAL(KIND=rsh), DIMENSION(imin:imax) :: Eknum


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

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

   ! dynamics
   CALL cas_init_param

   ! over-assignement

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

   ! over-assignement
   saliref=35.0_rsh
   rhoref=1000.0_rsh
   icon=2
   qmax=8.0_rsh   ! dire a quoi et ou cela sert
   qmaxz=10.0_rsh
   l_stateeq_lin=.FALSE.


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

   ! assignement : input variables for cas_init_bathy routine
   depthcst = 0.0_rsh
   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
#if defined key_siggen || defined key_gencoord
   hc= 5.0_rsh
#endif

   ! As depthcst is null, specify your own bathymetry
   Length_ETA=80000.0_rsh
   Length_XI = 20000.0_rsh
   dyval = 1250.0_rsh
   depth = 150.0_rsh
   DO j=jmin,jmax
     y=dyval*REAL(j-2,rsh)  ! -2 pour se caler avec Roms
     IF (y > Length_ETA/2.0_rsh) y=Length_ETA-y+dyval
     cff=MIN(depth,84.5_rsh+66.526_rsh*TANH(0.00015_rsh*(y-0.125_rsh*Length_ETA)))
     h0_g(:,j)=cff
     hx_g(:,j)=cff
     y=dyval*(REAL(j-2,rsh)+0.5_rsh)  ! -2 pour se caler avec Roms
     IF (y > Length_ETA/2.0_rsh) y=Length_ETA-y+dyval
     cff=MIN(depth,84.5_rsh+66.526_rsh*TANH(0.00015_rsh*(y-0.125_rsh*Length_ETA)))
     hy_g(:,j)=cff
   END DO
#if defined key_siggen || defined key_gencoord
   hc= h0_g(10,2)
#endif

   ! walls
   CALL cas_init_bathy(depthcst,nj_wall_n=nj_wall_n,nj_wall_s=nj_wall_s)

   ! over-assignement


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

   ! assignement : input variables for cas_init_grid routine
 
   ! parametres de forme
   !---------------------
   dxval = 1250.0_rsh
   dyval = 1250.0_rsh
   fisac=ASIN(-8.26e-05/0.0001458_rsh)/(pi/180.0_rsh)  ! lat to get corio=-8.26e-05
   gwac=1.0_riolg

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

   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=   '05/01/1900 14:00:00'
   date_start3d=  '01/01/1900 00:00:00'
   date_startagrif='01/01/1900 00:00:00'
   dtmin=120.0_rlg
   dtmax=120.0_rlg
   dt=120.0_rlg
   dtmin=82.0_rlg
   dtmax=82.0_rlg
   dt=82.0_rlg
   dtmin=164.0_rlg
   dtmax=164.0_rlg
   dt=164.0_rlg
   halfdt=dt*0.5_rlg

   suffix='updownwel'

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

   ! over-assignement
   pasor(1)=3600.0_rlg
   !pasor(1)=78720.0_rlg
   pasor(1)=1312.0_rlg
   l_out_sal(1)=.FALSE.
   l_out_wind(1)=.TRUE.
   l_out_turb(1)=.TRUE.
   l_out_zaxis_w(1)=.TRUE.


   ! 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.
   ! physics introduced into the simulated area (tide or ogcm solution)

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

   ! none
   CALL cas_init_visc

   ! over-assignement
   !vish_xe(:,:,:) = 100.0_rsh
   !vish_phi(:,:,:)= 100.0_rsh
   l_smagor=.TRUE.
   cosmag=0.2_rsh
   vismin=10.0_rsh
   vismax=700.0_rsh


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

   z0bottom   =0.0035_rsh
   z0surface  =0.0_rsh
   nzinit=1.0e-4
   kzinit=1.0e-4

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

   ! over-assignement
   turb_0eq_option=4
   bgdiff=1.0e-6
#ifdef key_virer
!   Si ANA_VMIX dans Romms seulement
   cff = 1.0_rsh/150.0_rsh
   DO j = ljmin,ljmaxp1
     DO i = limin,limaxp1
       zsurf = f_zw(h0(i,j),ssh(i,j),kmax,i,j)
       DO k = 1,kmax
         z = f_zw(h0(i,j),ssh(i,j),k,i,j)
         nz(k,i,j) = 0.002_rsh+0.008_rsh*EXP(cff*(z-zsurf))
            IF (i==10 .AND. j==30) then
              print*,'z_w(i,j,k)',z
              print*,'z_w(i,j,N)',zsurf
              print*,'Akv(i,j,k)',nz(k,i,j)
            end if
       END DO
     END DO
   END DO
#endif
#ifdef key_linear_bottom_stress
   frot_lin=0.0003_rsh ! Linear bottom drag coefficient (m/si) revoir unit
#endif

   ! Ekman depth
   Ekdepth = pi* SQRT( 2*nzinit/ABS(fu(1,1)) )
   ! Ekman number
   Eknum = nzinit / ABS(fu(1,1)) / hx_g(:,1)**2


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

   ! no forcing
   CALL cas_init_meteo

   ! over-assignement
   wind_veloc = -SQRT(0.1_rsh)  ! to get a wind stress of 0.1N/m2, assuming rhoair=1
   windamp_duration = 2.0_rlg * 86400.0_rlg ! 2 days
   ! toseelater pour instant vent selon x seulement. Generaliser ben comme dans
   ! meteo.F90 en fait


   ! 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

   ! tracers
   !--------

   zd0 = -35.0_rsh
   zd1 = -75.0_rsh
   thknss = 6.5_rsh
   strat = 150.0_rsh
   DO j=ljminm1,ljmaxp2
     DO i=liminm1,limaxp2
       DO k=1,kmax
         temp(k,i,j)=14.0_rsh &
                   +4.0_rsh*THKNSS*log( &
                       cosh((f_zw(h0(i,j),ssh(i,j),k,i,j)-Zd0)/THKNSS) &
                      /cosh((f_zw(h0(i,j),ssh(i,j),k-1,i,j)-Zd0)/THKNSS) ) &
                     /(f_zw(h0(i,j),ssh(i,j),k,i,j)-f_zw(h0(i,j),ssh(i,j),k-1,i,j)) &
                   +((f_zw(h0(i,j),ssh(i,j),k,i,j)+ &
                      f_zw(h0(i,j),ssh(i,j),k-1,i,j))/2.-Zd1)/STRAT

       END DO
     END DO
   END DO


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

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

  END SUBROUTINE cas_upwelling

   !!======================================================================

  SUBROUTINE cas_upwel_structure(h,D)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE cas_upwel_structure  ***
   !&E
   !&E ** Purpose : define the horizontally-varying structure which extends
   !&E              1D Ekman s model to horizontal dimensions using Welander s method
   !&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 :  Welander (1975), Brink (1983), Mitchum and Clarke (1986)
   !&E
   !&E ** History :
   !&E       !  2012-07    (V. Garnier)  Original code
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters, ONLY : imin,imax,pi,rsh
   USE comvars2d,  ONLY : fu

   !! * Arguments
   REAL(rsh), DIMENSION(imin:imax), INTENT(in) :: h    ! topography
   REAL(rsh),                       INTENT(in) :: D    ! Ekman depth  

   !! * Local declarations
   REAL(KIND=rsh), DIMENSION(imin:imax) :: alpha,S1,S2,T1,T2   !  struture functions
                                       !   S1p,S2p,T1p,T2p        !  derivative of struture functions
                                           !Eknum                  !  Ekman number
   INTEGER                :: i
   REAL(rsh)              :: phd,hemis

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

   hemis = fu(1,1)/ABS(fu(1,1))

   DO i=imin,imax
     phd = pi*h(i)/D
     alpha(i) = ( cosh(phd)*cos(phd) )**2 + ( sinh(phd)*sin(phd) )**2
     S1(i) = cosh(phd)*cos(phd) / alpha(i)
     S2(i) = sinh(phd)*sin(hemis*phd) / alpha(i)
     T1(i) = cosh(phd)*sinh(phd) / alpha(i)
     T2(i) = cos(phd)*sin(hemis*phd) / alpha(i)
   END DO


  END SUBROUTINE cas_upwel_structure
