 MODULE casinit

!!======================================================================
!!                   ***  MODULE  casinit  ***
!! Casinit :  Initialization of the test cases
!!            cas_init : call the correct test case
!!            cas_init_* : global initialisation of each component
!!
!!   History :
!!        !  2012-07 (V. Garnier) Reorganisation of the implementation of the test cases
!!        !  2014-07 (F. Ganthy) Introduction of obstructions module
!!        !  2015-04 (B.Thouvenin) Introduction of general 1DV test case.
!!        !                        cas_init_sedim with options for dzs or cv_sed variable
!!
!!======================================================================





























































   IMPLICIT NONE
   PRIVATE

!! * Accessibility
   PUBLIC cas_init                    ! routine called by main.F90
   PUBLIC cas_init_param              ! routine called by each test case
   PUBLIC cas_init_grid               ! routine called by each test case
   PUBLIC cas_init_bathy              ! routine called by each test case
   PUBLIC cas_init_io                 ! routine called by each test case
   PUBLIC cas_init_obc                ! routine called by each test case
   PUBLIC cas_init_meteo              ! routine called by each test case
   PUBLIC cas_init_visc               ! routine called by each test case
   PUBLIC cas_init_turb               ! routine called by each test case
   PUBLIC cas_init_tspar              ! routine called by each test case
   PUBLIC cas_init_progvar            ! routine called by each test case
   PUBLIC cas_init_misc               ! routine called by each test case

   PUBLIC cas_init_subs               ! routine called by each test case
   PUBLIC cas_init_subs_budget        ! routine called by each test case and by cas_init_subs
   PUBLIC cas_init_subs_disspart      ! routine called by each test case
   PUBLIC cas_init_outflow            ! routine called by each test case

   PUBLIC cas_init_river              ! routine called by each test case
   PUBLIC cas_init_traj               ! routine called by each test case

   PUBLIC cas_init_wave               ! routine called by each test case


   PUBLIC cas_init_sedim



 CONTAINS
!!======================================================================

  SUBROUTINE cas_init
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (F. Dumas) Original code
!&E       !  2007-03  (F. Vandermeirch) Introduction of test cases 11 to 13
!&E       !  2007-03  (V. Garnier, T. Duhaut) Introduction of test cases 11 to 13
!&E       !  2010-03  (A.-C. Bennis) Introduction of test case (wave : Mellor approach)
!&E       !  2013-06  (F. Grasso) Introduction of test case "cas_conso" (1DV sediment consolidation)
!&E       !  2014-07  (F. Ganthy) Introduction of test case "casobsthydroflume2dv"
!&E       !  2014-08  (F. Ganthy) Introduction of test case "casobstsedflumecycl2dv"
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : num_testcase
   USE comvars2d,  ONLY : l_modele2d
   USE toolmpi,ONLY : MPI_START2,MPI_START3

!! * Arguments

!! * Local declarations

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

   CALL MPI_START2(0,0)
   !CALL Agrif_Declare_Variables

! Allocation of variables
! -----------------------

   l_modele2d = .FALSE.
!IF (num_testcase==10 .OR. num_testcase==11 .OR. num_testcase==24 .OR. &
   IF (num_testcase==10 .OR. &
       num_testcase==18 .OR. num_testcase==19 .OR. num_testcase==20 .OR. &
       num_testcase==16 .OR. num_testcase==17 .OR. num_testcase==25) THEN
     l_modele2d = .TRUE.
   END IF

! Dynamical variables
   CALL alloc
! Intermediate variables
   CALL alloc_inter
! Stationnary variable (grid metrics)
   CALL alloc_stat

! Initialize test case configuration

   IF (num_testcase==0)  CALL cas_gravadj
 
   IF (num_testcase==1)  CALL cas_cylinder

   IF (num_testcase==2)  CALL cas_lock

   IF (num_testcase==3)  CALL cas_kelvin

   IF (num_testcase==4)  CALL cas_upwelling

   IF (num_testcase==5)  CALL cas_evap

   IF (num_testcase==7)  CALL cas_deepening

   IF (num_testcase==8)  CALL cas_seamount

   IF (num_testcase==9)  CALL cas_profile

   IF (num_testcase==10) CALL cas_wetdrying

! cas test courant barotrope constant ...
!   IF (num_testcase==11) CALL cas_bt_current_anc

! cas test courant barotrope constant
   IF (num_testcase==11) CALL cas_bt_current

! cas test courant barocline constant ...
   IF (num_testcase==12) CALL cas_bc_current

! cas test jet ACC current (COMODO)
   IF (num_testcase==13) CALL cas_acc

! cas test vortex barocline  ...
   IF (num_testcase==14) CALL cas_bcvortex

! cas test maree interne
   IF (num_testcase==15) CALL cas_intwave

! barotropic vortex
   IF (num_testcase==16) CALL cas_btvortex

! barotropic current
   IF (num_testcase==17) CALL cas_instbtjet

! advection of a spot
   IF (num_testcase==18) CALL cas_advspot

! does advection converve constants ? (smolarkiewicz test case)
   IF (num_testcase==19) CALL cas_advcst

! smolarkiewicz test case (time dependent velocity field)
   IF (num_testcase==20) CALL cas_advsmolt

! advection of a lense over a seamount
   IF (num_testcase==24) CALL cas_advspotmount

! sedimetology
   IF (num_testcase==21) CALL cas_sedim

! turbulence
   IF (num_testcase==22) CALL cas_turbkatophil

   IF (num_testcase==25)  CALL cas_thacker

   IF (num_testcase==26)  CALL cas_innershelf

   IF (num_testcase==32) CALL cas_1DVfloc

   IF (num_testcase==33) CALL cas_tidalflat2dv

   IF (num_testcase==34) CALL cas_1DVconso

   IF (num_testcase==35) CALL cas_rip

   IF (num_testcase==36) CALL cas_1DVbio

   IF (num_testcase==37) CALL cas_bioplat

! Obstruction (vegetation, others...)
   IF (num_testcase==38) CALL cas_obsthydroflume

   IF (num_testcase==39) CALL cas_obstsedflumecycl

   IF (num_testcase==40) CALL cas_1DV

   IF (num_testcase==41) CALL cas_estuar


   CALL MPI_START3

   !PRINT*, 'exit CAS_INIT routine'

  END SUBROUTINE cas_init

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

  SUBROUTINE cas_init_param
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_param  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh
   USE comvars2d,  ONLY : l_corio_cst,hm,hminim,h0fond,l_tssubadv_over_smalldepth,dhj,fwetmin,l_maxu,l_modele2d
   USE comvars3d,  ONLY : l_stability
   USE cometemp,   ONLY : chp,coext
   USE cometdyn,   ONLY : rhoair
   USE comaree,    ONLY : l_genpot

!! * Arguments

!! * Local declarations

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

   l_corio_cst=.TRUE.    ! true : f-plane, false : beta-plane

   l_maxu=.TRUE.    ! to list the maximum of the velocity in the log file

   l_modele2d=.FALSE.  ! if TRUE barotropic test case

! bathymetrie
   h0fond=0.0_rsh
   hm=0.0_rsh
   hminim=0.0_rsh
   fwetmin=1.0_rsh
   dhj=0.0_rsh   ! > 0 if wet-drying
   l_tssubadv_over_smalldepth=.true.

! miscellaneous
! -------------

   l_genpot=.FALSE.
   l_stability=.FALSE.
   rhoair=1.25_rsh
   chp=3986.0_rsh
   coext=0.1_rsh

  END SUBROUTINE cas_init_param

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

  SUBROUTINE cas_init_grid(dxcst,dycst,betaplan)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_grid  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : pi,rsh,rlg,imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax, &
                          liminm1,ljmaxp1,ljminm1,limaxp1,liminm2,ljmaxp2,ljminm2,limaxp2
   USE comvars2d,  ONLY : fisac,gwac,dfiac,dgac,dx,dy,dx_u,dx_v,dy_u,dy_v,dx_f,dy_f,surf, &
                          fu,fv,tfiusr,tfivsr,h0,hx,hy, &
                          jb,jh,ig,id,goulet,govlet,radearth,l_corio_cst, &
                          lon2d_g,lat2d_g,lon2du_g,lon2dv_g,lon2df_g,lat2du_g,lat2dv_g,lat2df_g, &
                          lon2d,lat2d,lon2du,lon2dv,lon2df,lat2du,lat2dv,lat2df
   USE comvars3d,  ONLY : sig,sigw,dsigu,dsigw


!! * Arguments
   REAL(KIND=rsh), INTENT(IN)    :: dxcst,dycst
   REAL(KIND=rlg), INTENT(IN), OPTIONAL :: betaplan

!! * Local declarations
   INTEGER :: i,j,k
   REAL(KIND=rsh)    :: fi,f0,jf0,jj
   REAL(KIND=rlg)    :: beta

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

! Horizontal discretization
! -------------------------

   dx(:,:)=dxcst
   dx_u(:,:)=dxcst
   dx_v(:,:)=dxcst
   dx_f(:,:)=dxcst
   dy(:,:)=dycst
   dy_u(:,:)=dycst
   dy_v(:,:)=dycst
   dy_f(:,:)=dycst
   surf(:,:)=dxcst*dycst
   goulet(:,:)=1.0_rsh
   govlet(:,:)=1.0_rsh

! Vertical discretization
! -----------------------

   DO k=1,kmax
     sig(k)=(REAL(k-kmax,rsh)-0.5_rsh)/REAL(kmax,rsh)
   END DO
   sig(kmax+1)=-sig(kmax)
   sig(0)=ABS(sig(1))-2.0_rsh

   DO k=1,kmax
     dsigw(k)=sig(k+1)-sig(k)
     sigw(k)=0.5_rsh*(sig(k)+sig(k+1))
   END DO
   sigw(kmax+1)=-sigw(kmax)
   sigw(0)=-1.0_rsh

   DO k=1,kmax
     dsigu(k)=sigw(k)-sigw(k-1)
   END DO




! longitude latitude
! ------------------

   dfiac=dycst*360.0_rsh/(2*pi*radearth)
   dgac=dxcst*360.0_rsh/(2*pi*radearth)

   lon2d_g(:,0)=(/ (gwac + dgac *(i-1), i=0,imax) /)
   DO j=1,jmax
     lon2d_g(:,j)=lon2d_g(:,0)
   END DO
   lon2du_g(:,:)=lon2d_g(:,:)+dgac/2.0_rsh
   lon2dv_g(:,:)=lon2d_g(:,:)
   lon2df_g(:,:)=lon2du_g(:,:)
   lat2d_g(0,:)=(/ (fisac + dfiac *(j-1), j=0,jmax) /)
   DO i=1,imax
     lat2d_g(i,:)=lat2d_g(0,:)
   END DO
   lat2du_g(:,:)=lat2d_g(:,:)
   lat2dv_g(:,:)=lat2d_g(:,:)+dfiac/2.0_rsh
   lat2df_g(:,:)=lat2dv_g(:,:)

   lon2du(liminm2:limaxp2,ljminm2:ljmaxp2)=lon2du_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lon2d(liminm2:limaxp2,ljminm2:ljmaxp2)=lon2d_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lon2dv(liminm2:limaxp2,ljminm2:ljmaxp2)=lon2dv_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lon2df(liminm2:limaxp2,ljminm2:ljmaxp2)=lon2du_g(liminm2:limaxp2,ljminm2:ljmaxp2)

   lat2dv(liminm2:limaxp2,ljminm2:ljmaxp2)=lat2dv_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lat2d(liminm2:limaxp2,ljminm2:ljmaxp2)=lat2d_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lat2du(liminm2:limaxp2,ljminm2:ljmaxp2)=lat2du_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   lat2df(liminm2:limaxp2,ljminm2:ljmaxp2)=lat2dv_g(liminm2:limaxp2,ljminm2:ljmaxp2)

! coriolis parameter
! ------------------

   f0 = 2.0_rsh * 2.0_rsh*pi/(86164.0_rsh)   ! sidereal day = 86164 s
   fi = REAL(fisac,rsh)
   IF (l_corio_cst) THEN
! f-plane
     DO j=ljmin,ljmax
       fu(:,j)=f0*SIN(fi*pi/180.0_rsh)
       fv(:,j)=f0*SIN(fi*pi/180.0_rsh)
     END DO
   ELSE
! beta-plane
     IF ( PRESENT(betaplan) ) THEN
       beta=betaplan
     ELSE
       beta=f0*COS(fi*pi/180.0_rlg)/radearth
     END IF
     DO j=ljmin,ljmax
       jj = REAL(j,rsh)-REAL(jmin+jmax,rsh)/2.0_rsh
       !IF (.NOT..TRUE.) THEN
         ! fi = REAL(Agrif_Parent(fisac),rsh)
! j of the middle of the domain in the parent grid
         ! jf0 = REAL(Agrif_Parent(jmin)+Agrif_Parent(jmax),rsh)/2.0_rsh
! Remove the location of the middle of the parent grid in the child grid
         ! jj = REAL(j,rsh) - ( (jf0-Agrif_IY())*Agrif_irhoy()-INT(Agrif_irhoy()/2.0_rsh) )
         ! jj = jj - 1.0_rsh
       !ENDIF
       fu(:,j)=f0*SIN(fi*pi/180.0_rsh)+beta*dycst*jj
       fv(:,j)=f0*SIN(fi*pi/180.0_rsh)+beta*dycst*(jj+0.5_rsh)
     END DO
   END IF

! Oh la attention a revoir
   tfiusr(:,:)=0.0_rsh
   tfivsr(:,:)=0.0_rsh

  END SUBROUTINE cas_init_grid

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

  SUBROUTINE cas_init_bathy(depthcst,ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_bathy  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : valmanq,epsilon,rsh,imin,imax,jmin,jmax, &
                          limin,limax,ljmin,ljmax,liminm2,limaxp2,ljminm2,ljmaxp2, &
                          limin,limaxp1,ljmin,ljmaxp1,liminm1,ljminm1
   USE comvars2d,  ONLY : h0_g,hx_g,hy_g,h0,hx,hy,dhxy,dhx,dhy,dhx_g,dhy_g,dhj,h0fond,mslshift, &
                          dhxy_g,dhxy,hrdef
   USE comcolig2d, ONLY : me,mn,ms,mo,whx,why


   USE comMUSTANG,  ONLY : l_morphocoupl


!! * Arguments
   INTEGER, OPTIONAL, INTENT(IN) :: ni_wall_e,ni_wall_w,nj_wall_n,nj_wall_s
   REAL(KIND=rsh), INTENT(IN) :: depthcst

!! * Local declarations
   INTEGER :: i,j
   REAL(KIND=rsh),DIMENSION(imin:imax,jmin:jmax) :: hx_t,hy_t


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

   !PRINT*, 'enter CAS_INIT_BATHY routine'

   mslshift=0.0_rsh
   hrdef=0.0_rsh

! depth of the domain
! -------------------
   IF (depthcst/=0.0_rsh) THEN
     hx_g(:,:)=depthcst
     hy_g(:,:)=depthcst
     h0_g(:,:)=depthcst+h0fond
   END IF

! wall mask
! ---------
   me(:,:)=1
   mo(:,:)=1
   ms(:,:)=1
   mn(:,:)=1
   whx(:,:)=1
   why(:,:)=1

! walls of the domain
! -------------------

   !IF (.TRUE.) THEN

   IF (PRESENT(ni_wall_e)) THEN
     hx_g(imax-ni_wall_e  :imax,:)=-valmanq-h0fond
     hy_g(imax-ni_wall_e+1:imax,:)=-valmanq-h0fond
     h0_g(imax-ni_wall_e+1:imax,:)=-valmanq
     IF (limax==imax) THEN
     me(imax-ni_wall_e  :imax,:)=0
     mo(imax-ni_wall_e+1:imax,:)=0
     ms(imax-ni_wall_e  :imax,:)=0
     mn(imax-ni_wall_e  :imax,:)=0
     whx(imax-ni_wall_e  :imax,:)=0
     why(imax-ni_wall_e+1:imax,:)=0
     ENDIF
   END IF

   IF (PRESENT(ni_wall_w)) THEN
     hx_g(imin:imin+ni_wall_w-1,:)=-valmanq-h0fond
     hy_g(imin:imin+ni_wall_w-1,:)=-valmanq-h0fond
     h0_g(imin:imin+ni_wall_w-1,:)=-valmanq
     IF (limin==imin) THEN
     me(imin:imin+ni_wall_w-1,:)=0
     mo(imin:imin+ni_wall_w  ,:)=0
     ms(imin:imin+ni_wall_w-1,:)=0
     mn(imin:imin+ni_wall_w-1,:)=0
     whx(imin:imin+ni_wall_w-1,:)=0
     why(imin:imin+ni_wall_w-1,:)=0
     ENDIF
   END IF

   IF (PRESENT(nj_wall_n)) THEN
     hx_g(:,jmax-nj_wall_n+1:jmax)=-valmanq-h0fond
     hy_g(:,jmax-nj_wall_n  :jmax)=-valmanq-h0fond
     h0_g(:,jmax-nj_wall_n+1:jmax)=-valmanq
     IF (ljmax==jmax) THEN
     me(:,jmax-nj_wall_n  :jmax)=0
     mo(:,jmax-nj_wall_n  :jmax)=0
     ms(:,jmax-nj_wall_n+1:jmax)=0
     mn(:,jmax-nj_wall_n  :jmax)=0
     whx(:,jmax-nj_wall_n+1:jmax)=0
     why(:,jmax-nj_wall_n  :jmax)=0
     ENDIF
   END IF

   IF (PRESENT(nj_wall_s)) THEN
     hx_g(:,jmin:jmin+nj_wall_s-1)=-valmanq-h0fond
     hy_g(:,jmin:jmin+nj_wall_s-1)=-valmanq-h0fond
     h0_g(:,jmin:jmin+nj_wall_s-1)=-valmanq
     IF (ljmin==jmin) THEN
     me(:,jmin:jmin+nj_wall_s-1)=0
     mo(:,jmin:jmin+nj_wall_s-1)=0
     ms(:,jmin:jmin+nj_wall_s  )=0
     mn(:,jmin:jmin+nj_wall_s-1)=0
     whx(:,jmin:jmin+nj_wall_s-1)=0
     why(:,jmin:jmin+nj_wall_s-1)=0
     ENDIF
   END IF

   !ENDIF

!
! ----------------
   hx_t(:,:)=100000.0_rsh
   hy_t(:,:)=100000.0_rsh
   WHERE (hx_g(:,:) /= -valmanq -h0fond) hx_t(:,:)=hx_g(:,:)
   WHERE (hy_g(:,:) /= -valmanq -h0fond) hy_t(:,:)=hy_g(:,:)
   dhxy_g(:,:)=1.0_rsh
   DO j=jmin+1,jmax
      DO i=imin+1,imax
         dhxy_g(i,j)=MAX(hx_g(i,j),hx_g(i-1,j),hy_g(i,j),hy_g(i,j-1))- &
                     MIN(hx_t(i,j),hx_t(i-1,j),hy_t(i,j),hy_t(i,j-1))+epsilon
      END DO
   END DO

   dhx_g(:,:)=0.0_rsh
   dhy_g(:,:)=0.0_rsh
   WHERE (hx_g(:,:) < hrdef) dhx_g(:,:)=dhj
   WHERE (hy_g(:,:) < hrdef) dhy_g(:,:)=dhj




! local variables
! ---------------
   h0(:,:) = h0_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hx(:,:) = hx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   hy(:,:) = hy_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhx(:,:) = dhx_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhy(:,:) = dhy_g(liminm2:limaxp2,ljminm2:ljmaxp2)
   dhxy(:,:) = dhxy_g(limin:limax,ljmin:ljmax)


! initialization : no morphodynamic by default
! (l_morphocupl need to be initialize before init_subs because of MUSTANG_alloc)
   l_morphocoupl=.FALSE.


  END SUBROUTINE cas_init_bathy

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


  SUBROUTINE cas_init_io(date_start,date_end,date_start3d,date_startagrif, &
                         suffix,l_nc4par,nb_file,date_startshift_0)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_io  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg,riosh,lchain,imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax
   USE comvars2d, ONLY : date_ref,tdeb,t,tfin,tdeb3d,tdebagrif,date_startshift, &
                         cdate,ijour,imois,ian,iheure,iminu,isec, &
                         cflcrt,tobs,tdebobs,cflumaxmax,cflvmaxmax,l_obc2drank_save, &
                         l_saverestart_1file,l_initfromfile,l_saverestart_bydate,l_init_rtime, &
                         l_init_rdt,l_init_rssh,l_init_rbtvel,l_init_r3dvel,l_init_rwz, &
                         l_init_rturb,l_init_rsal,l_init_rtemp,out_dir, &
                         l_init_restart,file_saverestart,file_init,saverestart_step, &
                         name_bathy_hx,name_bathy_hy,name_bathy_h0,name_in_dt,name_in_xev,   &
                         name_in_u,name_in_v,name_in_uz,name_in_vz,      &
                         name_in_wz,name_in_temp,name_in_sal,name_in_psi,&
                         name_in_ect,name_in_kz,name_in_nz,name_in_vis2d,&
                         name_in_vis3d_xe,name_in_vis3d_phi, &
                         iscreen,iscreenlog,ierrorlog,iwarnlog
   USE comvars2d, ONLY : dx,dy,tsauv
   USE comvars2d, ONLY : cg_conv_title,cg_conv_comment,cg_conv_source,      &
                         cg_conv_area,cg_conv_institution,                  &
                         cg_conv_institution_references,cg_conv_data_centre,&
                         cg_conv_data_centre_references,                    &
                         cg_conv_forecast_range,cg_conv_forecast_type,      &
                         cg_conv_field_type,cg_conv_operational_status
   USE comvars2d, ONLY : riog_valid_min_h0,riog_valid_max_h0,         &
                         riog_valid_min_u,riog_valid_max_u,             &
                         riog_valid_min_temp,riog_valid_max_temp,       &
                         riog_valid_min_sal,riog_valid_max_sal,         &
                         riog_valid_min_xe,riog_valid_max_xe,           &
                         riog_valid_min_wstress,riog_valid_max_wstress, &
                         riog_valid_min_flux,riog_valid_max_flux,       &
                         riog_valid_min_visc,riog_valid_max_visc,       &
                         riog_valid_min_nz,riog_valid_max_nz,           &
                         riog_valid_min_kz,riog_valid_max_kz,           &
                         riog_valid_min_bz,riog_valid_max_bz, &
                         riog_valid_min_dx,riog_valid_min_dy,         &
                         riog_valid_max_dx,riog_valid_max_dy
   USE comvars2d, ONLY : name_out_h0,name_out_hx,name_out_hy,name_out_u, &
                         name_out_v,name_out_uz,name_out_vz,name_out_temp, &
                         name_out_sal,name_out_xe,name_out_tosx,name_out_tosy, &
                         name_out_vis2d,name_out_vishxe,name_out_vishphi, &
                         name_out_nz,name_out_kz,name_out_bz
   USE comvars2d, ONLY : name_out_mrad,name_out_mirlos,name_out_mlatlos,    &
                         name_out_msenslos,name_out_mevap,name_out_mrain,   &
                         name_out_mtosx,name_out_mtosy
   USE comvars2d, ONLY : nvalmax,perchrono,indchrono,l_chrono,tchrono

   USE output, ONLY : fileresu,tdebsor,tfinsor,nb_out,tinter,tsort,pasor,   &
                      l_out_precdate,l_out_demer,l_out_ssh,l_out_ubt,l_out_u3d, &
                      l_out_sal,l_out_temp,l_out_wind,l_out_buoy,l_out_wave,&
                      l_out_visc,l_out_turb,l_out_off,l_out_atmsflx,        &
                      iextmin,iextmax,l_out_nc4par,l_out_pack,l_out_zaxis_w,&
                      jextmin,jextmax,kextmin,kextmax,stepi,stepj,stepk,    &
                      igextmin,igextmax,jgextmin,jgextmax,l_filebydate,     &
                      l_out_timeave,nbave,file_output,fileresu_verif

   USE comdiag,     ONLY : l_diag,file_diag,tdebdiag,tfindiag
   USE toolmpi, ONLY : ADD_TAG,MASTER,min_all_mpi_rsh,max_all_mpi_rsh

   USE comMUSTANG,   ONLY : name_out_eptot,name_out_nbniv,name_out_dzs,name_out_tenfon,&
                           name_out_tenfonc,name_out_tenfonw,&
                           riog_valid_min_eptot,riog_valid_max_eptot,riog_valid_min_nbniv, &
                           riog_valid_max_nbniv,riog_valid_min_dzs,riog_valid_max_dzs,     &
                           riog_valid_min_tenfon,riog_valid_max_tenfon



!! * Arguments
   CHARACTER(LEN=19), INTENT(IN)     :: date_start,date_end,date_start3d,date_startagrif
   CHARACTER(len=lchain), INTENT(IN) :: suffix
   LOGICAL, OPTIONAL,     INTENT(IN) :: l_nc4par
   INTEGER, OPTIONAL,     INTENT(IN) :: nb_file
   REAL(KIND=rlg), OPTIONAL,INTENT(IN) :: date_startshift_0

!! * Local declarations
   CHARACTER(LEN=lchain) :: logfilename,chain
   INTEGER               :: i
   REAL(KIND=rsh)        :: value
   REAL(KIND=rlg)        :: tool_datosec,tdb,tfi
   CHARACTER(LEN=19)     :: tool_sectodat,charint


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

   !PRINT*, 'enter CAS_INIT_IO routine'

! dates
! -----

   date_ref = '01-Jan-1900 00:00:00'
   tdeb     = tool_datosec(date_start)
   t=tdeb
   tfin     = tool_datosec(date_end)
   tdeb3d   = tool_datosec(date_start3d)
   tdebagrif = tool_datosec(date_startagrif)
   IF (PRESENT(date_startshift_0) ) THEN
     date_startshift=date_startshift_0
   ELSE
     date_startshift= 0.0_rlg
   ENDIF
   cflcrt=0.6_rsh
   tobs=12.0_rlg*3600.0_rlg
   tdebobs=t
   cdate = tool_sectodat(t)
   CALL tool_decompdate(cdate,ijour,imois,ian,iheure,iminu,isec)

! restart file
! ------------

   l_initfromfile=.FALSE.
   file_init = 'save.nc'
   l_saverestart_1file=.FALSE.
   l_saverestart_bydate=.FALSE.
   file_saverestart='save.nc'
   saverestart_step=5.0_rlg*3600.0_rlg   ! in hours
   tsauv=tdeb + saverestart_step

   l_init_restart=.TRUE.
   l_init_rtime=.TRUE.
   l_init_rdt=.TRUE.
   l_init_rssh=.TRUE.
   l_init_rbtvel=.TRUE.
   l_init_r3dvel=.TRUE.
   l_init_rwz=.TRUE.
   l_init_rturb=.TRUE.
   l_init_rsal=.TRUE.
   l_init_rtemp=.TRUE.

! input parametrization
! ---------------------

   name_in_dt='dt'
   name_bathy_h0='h0'
   name_bathy_hx='hx'
   name_bathy_hy='hy'
   name_in_xev='xev'
   name_in_u='u'
   name_in_v='v'
   name_in_uz='uz'
   name_in_vz='vz'
   name_in_wz='wz'
   name_in_temp='temp'
   name_in_sal='sal'
   name_in_ect='ect'
   name_in_psi='psi'
   name_in_kz='kz'
   name_in_nz='nz'
   name_in_vis2d='vis2d'
   name_in_vis3d_xe='vish3d_xe'
   name_in_vis3d_phi='vish3d_phi'

! output parametrization
! ----------------------

   IF ( PRESENT(nb_file) ) THEN
     nb_out = nb_file
   ELSE
     nb_out = 1
   END IF
   ALLOCATE( fileresu(nb_out),tdebsor(nb_out),tfinsor(nb_out),tinter(nb_out) )
   ALLOCATE( tsort(nb_out),pasor(nb_out),l_filebydate(nb_out),l_out_precdate(nb_out) )
   ALLOCATE( l_out_off(nb_out),l_out_atmsflx(nb_out),l_out_visc(nb_out),l_out_turb(nb_out) )
   ALLOCATE( l_out_timeave(nb_out),l_out_buoy(nb_out), l_out_zaxis_w(nb_out) )
   ALLOCATE( l_out_demer(nb_out),l_out_ubt(nb_out),l_out_u3d(nb_out) )
   ALLOCATE( l_out_sal(nb_out),l_out_temp(nb_out),l_out_ssh(nb_out),l_out_wind(nb_out) )
   ALLOCATE( l_out_wave(nb_out))
   ALLOCATE( fileresu_verif(nb_out+2))   ! to manage diagnostic files (diagdyn.F90 and wavediagdyn.F90)
   ALLOCATE( iextmin(nb_out),iextmax(nb_out),jextmin(nb_out),jextmax(nb_out) )
   ALLOCATE( igextmin(nb_out),igextmax(nb_out),jgextmin(nb_out),jgextmax(nb_out))
   ALLOCATE( kextmin(nb_out),kextmax(nb_out),stepi(nb_out),stepj(nb_out),stepk(nb_out) )


   file_output = 'output.dat'
   fileresu(1)='results.nc'

   !IF (.NOT. .TRUE.) THEN
   ! fileresu(1)='results'//'_zoom'//TRIM(Agrif_Cfixed())//'.nc'
   !ENDIF
   IF (nb_out >= 2) THEN
     DO i=2,nb_out
       WRITE(chain,123) i
       fileresu(i)='results_'//TRIM(chain)//'.nc'
       !IF (.NOT. .TRUE.) THEN
       ! fileresu(i)='results_'//TRIM(chain)//'_agrif_'//TRIM(Agrif_Cfixed())//'.nc'
       !ENDIF
     END DO
   END IF
123 FORMAT(I1)
   tdebsor(1:nb_out)=tool_datosec(date_start)
   tfinsor(1:nb_out)=tool_datosec(date_end)
   pasor(1:nb_out)=3600.0_rlg
   l_filebydate(1:nb_out)=.FALSE.
   l_out_off(1:nb_out)=.FALSE.
   l_out_precdate(1:nb_out)=.FALSE.
   l_out_demer(1:nb_out)=.FALSE.
   igextmin(1:nb_out)=imin
   igextmax(1:nb_out)=imax
   jgextmin(1:nb_out)=jmin
   jgextmax(1:nb_out)=jmax
   iextmin(1:nb_out)=limin
   iextmax(1:nb_out)=limax
   jextmin(1:nb_out)=ljmin
   jextmax(1:nb_out)=ljmax
   kextmin(1:nb_out)=1
   kextmax(1:nb_out)=kmax
   stepi(1:nb_out)=1
   stepj(1:nb_out)=1
   stepk(1:nb_out)=1
   l_out_ssh(1:nb_out)=.TRUE.
   l_out_ubt(1:nb_out)=.TRUE.
   l_out_u3d(1:nb_out)=.TRUE.
   l_out_sal(1:nb_out)=.TRUE.
   l_out_temp(1:nb_out)=.TRUE.
   l_out_wind(1:nb_out)=.FALSE.
   l_out_atmsflx(1:nb_out)=.FALSE.
   l_out_timeave(1:nb_out)=.FALSE.
   l_out_visc(1:nb_out)=.FALSE.

   l_out_turb(1:nb_out)=.FALSE.

   l_out_buoy(1:nb_out)=.FALSE.
   l_out_wave(1:nb_out)=.FALSE.

   l_out_zaxis_w(1:nb_out)=.FALSE.


   out_dir="./"
   l_obc2drank_save=.FALSE.
   l_obc2drank_save=.FALSE.
   IF ( PRESENT(l_nc4par) ) THEN
     l_out_nc4par = l_nc4par
   ELSE
     l_out_nc4par = .FALSE.
   END IF
   IF ( .NOT. l_out_nc4par) THEN
     DO i=1,nb_out
      IF ( .NOT. l_filebydate(i) ) THEN
        CALL ADD_TAG(fileresu(i),LEN(fileresu(i)))
      ENDIF
      fileresu_verif(i)=fileresu(i)
     END DO
   END IF
   l_out_pack = .FALSE.
   nbave=0    ! no file for averaged variables

   name_out_h0='h0'
   name_out_hx='hx'
   name_out_hy='hy'
   name_out_u='u'
   name_out_v='v'
   name_out_uz='uz'
   name_out_vz='vz'
   name_out_temp='temp'
   name_out_sal='sal'
   name_out_xe='ssh'
   name_out_vis2d='vis2d'
   name_out_vishxe='vish_xe'
   name_out_vishphi='vish_phi'
   name_out_nz='nz'
   name_out_kz='kz'
   name_out_bz='bz'
   name_out_tosx='taux'
   name_out_tosy='tauy'

   name_out_eptot='EPTOT'
   name_out_nbniv='NBNIV'
   name_out_dzs='DZS'
   name_out_tenfon='TENFON'
   name_out_tenfonc='TENFONC'
   name_out_tenfonw='TENFONW'


   name_out_mrad='flsolar'
   name_out_mirlos='fltherm'
   name_out_mlatlos='fllat'
   name_out_msenslos='flsens'
   name_out_mevap='evap'
   name_out_mrain='rain'
   name_out_mtosx='ustress'
   name_out_mtosy='vstress'

! initialization of attributes (max and min valid values for each variable)
! -------------------------------------------------------------------------

   riog_valid_min_h0=-50.0_riosh
   riog_valid_max_h0=10000.0_riosh
   riog_valid_min_u=-20.0_riosh
   riog_valid_max_u=20.0_riosh
   riog_valid_min_temp=-5.0_riosh
   riog_valid_max_temp=50.0_riosh
   riog_valid_min_sal=-5.0_riosh
   riog_valid_max_sal=50.0_riosh
   riog_valid_min_xe=-50.0_riosh
   riog_valid_max_xe=50.0_riosh
   riog_valid_min_wstress=-2.0_riosh
   riog_valid_max_wstress=2.0_riosh
   riog_valid_min_flux=-1000.0_riosh
   riog_valid_max_flux=2000.0_riosh
   riog_valid_min_visc=0.0_riosh
   riog_valid_max_visc=1000.0_riosh
   riog_valid_min_nz=0.0_riosh
   riog_valid_max_nz=1000.0_riosh
   riog_valid_min_kz=0.0_riosh
   riog_valid_max_kz=1000.0_riosh
   riog_valid_min_bz=-1.0_riosh
   riog_valid_max_bz=1.0_riosh
   value=MINVAL(dx(limin:limax,ljmin:ljmax))
   CALL MIN_ALL_MPI_RSH(value)
   riog_valid_min_dx=value-1.0_riosh
   value=MAXVAL(dx(limin:limax,ljmin:ljmax))
   CALL MAX_ALL_MPI_RSH(value)
   riog_valid_max_dx=value+1.0_riosh
   value=MINVAL(dy(limin:limax,ljmin:ljmax))
   CALL MIN_ALL_MPI_RSH(value)
   riog_valid_min_dy=value-1.0_riosh
   value=MAXVAL(dy(limin:limax,ljmin:ljmax))
   CALL MAX_ALL_MPI_RSH(value)
   riog_valid_max_dy=value+1.0_riosh

   riog_valid_min_eptot=0.0_riosh
   riog_valid_max_eptot=100.0_riosh
   riog_valid_min_nbniv=0.0_riosh
   riog_valid_max_nbniv=200.0_riosh
   riog_valid_min_dzs=0.0_riosh
   riog_valid_max_dzs=100.0_riosh
   riog_valid_min_tenfon=0.0_riosh
   riog_valid_max_tenfon=50.0_riosh



! initialization of other attributes
! ----------------------------------

   cg_conv_title="name of area"
   cg_conv_comment=" "
   cg_conv_source="MARS CONF"
   cg_conv_area="area"
   cg_conv_institution="IFREMER"
   cg_conv_institution_references="http://www.ifremer.fr"
   cg_conv_data_centre="DYNECO/PHYSED"
   cg_conv_data_centre_references="http://wwz.ifremer.fr/mars3d"
   cg_conv_forecast_range="n/a"
   cg_conv_forecast_type="n/a"
   cg_conv_field_type="n/a"
   cg_conv_operational_status="experimental"

! log files
! ---------

   !IF (.TRUE.) THEN
     iscreen=12
     iscreenlog=18
     ierrorlog=17
     iwarnlog=16
   !ELSE
     ! charint = Agrif_Cfixed()
     ! READ(charint,'(I4)') i
     ! iscreen=(i+1)*10+2
     ! iscreenlog=(i+1)*10+8
     ! ierrorlog=(i+1)*10+7
     ! iwarnlog=(i+1)*10+6
   !ENDIF


   logfilename= 'listing'


   !IF (.NOT. .TRUE.) logfilename= 'listing_zoom'//TRIM(Agrif_Cfixed())
   CALL ADD_TAG(logfilename,LEN(logfilename))
   IF (MASTER) THEN
   IF (iscreen/=6) OPEN(unit=iscreen,file=logfilename,status="replace")
   ENDIF
   logfilename= 'simu_'//TRIM(suffix)//'.log'
   !IF (.NOT. .TRUE.) ! logfilename= 'simu_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
!CALL ADD_TAG(logfilename,LEN(logfilename))
   IF (MASTER) THEN
   IF (iscreenlog/=6) OPEN(unit=iscreenlog,file=logfilename,status="replace")
   ENDIF
   logfilename='error_'//TRIM(suffix)//'.log'
   !IF (.NOT. .TRUE.) logfilename= 'error_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
   IF (ierrorlog/=6) OPEN(unit=ierrorlog,file=logfilename,status="replace")
   logfilename='warning_'//TRIM(suffix)//'.log'
   !IF (.NOT. .TRUE.) logfilename= 'warning_'//TRIM(suffix)//'_'//TRIM(Agrif_Cfixed())//'.log'
   IF (iwarnlog/=6) OPEN(unit=iwarnlog,file=logfilename,status="replace")


! diagnoses
! ---------

   l_diag=.FALSE.
   file_diag='diagnostic.txt'
   tdebdiag=tdebsor(1)
   tfindiag=tfinsor(1)

! max of velocity
! ---------------

   cflumaxmax=0.0_rsh
   cflvmaxmax=0.0_rsh

! Initialize time vector tchrono relative to
! river dischage and outflow events
! (required for the restart file)
! --------------------------------

   perchrono=0.1_rlg
   perchrono=perchrono*3600.0_rlg
!   IF (l_meteo_hom .and. .not.l_meteo_stat) perchrono=3.0_rlg*3600.0_rlg
   IF (l_initfromfile .AND. l_init_rtime) THEN
     tdb=tdeb-perchrono
   ELSE
     tdb=tdeb+date_startshift*86400.0_rlg-perchrono
   END IF
   tfi=tfin+perchrono*10.0_rlg ! to make sure the run will end properly
   nvalmax=NINT((tfi-tdb)/perchrono+0.5_rlg)+1
   ALLOCATE(tchrono(nvalmax))
   DO i=1,nvalmax
     tchrono(i)=tdb+(i-1.0_rlg)*perchrono
   END DO
   indchrono=1
   l_chrono=.TRUE.
   t=tdeb+date_startshift*86400.0_rlg


   !PRINT*, 'exit CAS_INIT_IO routine'

  END SUBROUTINE cas_init_io

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

  SUBROUTINE cas_init_obc
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_obc  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh
   USE comvars2d,  ONLY : l_obc_cycl_x,l_obc_cycl_y
   USE comvars3d,  ONLY : saliref,tetaref
   USE obccombine, ONLY : l_obc_south, l_obc_north,l_obc_east,l_obc_west, &
                          l_obc_mars,l_obc_diri,l_obc_char,l_obc_cycl,    &
                          l_obc_tide,l_obc_ogcm, &
                          obc_tide_sshn,obc_tide_sshs,obc_tide_sshe,obc_tide_sshw, &
                          obc_tide_us,obc_tide_vs,obc_tide_uw,obc_tide_vw, &
                          obc_tide_un,obc_tide_vn, obc_tide_ue,obc_tide_ve
   USE obcogcm,    ONLY : obc_te,obc_tw,obc_ts,obc_tn,obc_se,obc_sw,obc_ss,obc_sn 
   USE obcogcm,    ONLY : obc_coefrel
   USE comaree,    ONLY : l_tide_harmcompo

!! * Arguments

!! * Local declarations

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

! geometry (where are the open boundaries ?)
   l_obc_south=.FALSE.
   l_obc_north=.FALSE.
   l_obc_west=.FALSE.
   l_obc_east=.FALSE.
   !IF (.NOT..TRUE.) THEN
     ! l_obc_south = .TRUE.
     ! l_obc_north = .TRUE.
     ! l_obc_east = .TRUE.
     ! l_obc_west = .TRUE.
   !ENDIF

! used scheme to manage obc

   l_obc_mars = .FALSE.
   l_obc_diri = .FALSE.
   l_obc_char = .FALSE.
   l_obc_cycl = .FALSE.
   l_obc_cycl_x=.FALSE.
   l_obc_cycl_y=.FALSE.

! physics introduced into the simulated area (tide or ogcm solution)

   l_obc_tide=.FALSE.
   l_tide_harmcompo=.FALSE.
   l_obc_ogcm=.FALSE.

! tide at the open boundaries

   obc_tide_sshs(:)=0.0_rsh
   obc_tide_sshn(:)=0.0_rsh
   obc_tide_sshe(:)=0.0_rsh
   obc_tide_sshw(:)=0.0_rsh

   obc_coefrel=0.0_rsh


  END SUBROUTINE cas_init_obc

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

  SUBROUTINE cas_init_meteo
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_meteo  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg
   USE cometdyn,   ONLY : l_meteo_stat,l_meteo_hom,wind_veloc,wind_dir,cds, &
                          windamp_duration,pa,paref
   USE cometemp,   ONLY : l_sflx_solarcst,sflx_solarcst,l_sflx_radlossbot
   USE output,     ONLY : nbatm
   USE sflxatm,    ONLY : imeteo_dragtype,imeteo_exchtype

!! * Arguments

!! * Local declarations

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

   l_meteo_stat=.TRUE.
   l_meteo_hom=.TRUE.
   wind_veloc=0.0_rsh
   wind_dir=0.0_rsh
   windamp_duration = 0.0_rlg
   l_sflx_solarcst=.TRUE.
   sflx_solarcst=0.0_rsh
   l_sflx_radlossbot = .FALSE.
   nbatm=0
   imeteo_dragtype=0
   imeteo_exchtype=0
   cds=0.0016_rsh

   paref=101500.0
   pa(:,:)=paref

  END SUBROUTINE cas_init_meteo

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

  SUBROUTINE cas_init_visc
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_visc  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh
   USE comvars2d,  ONLY : vis2d,vish_phi,vish_xe,vismin,vismax,cosmag,l_smagor, &
                          l_spongw,l_sponge,l_spongs,l_spongn,sponge_nbcell, &
                          fvisc,vismul

!! * Arguments

!! * Local declarations

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

   fvisc=0.0_rsh
   vis2d(:,:) = 0.0_rsh
   vish_xe(:,:,:) = 0.0_rsh
   vish_phi(:,:,:)= 0.0_rsh

   l_smagor=.FALSE.
   cosmag=0.1_rsh
   vismin=0.001_rsh
   vismax=4.0_rsh

   vismul=1.0_rsh
   sponge_nbcell=5
   l_spongs=.FALSE.
   l_spongn=.FALSE.
   l_sponge=.FALSE.
   l_spongw=.FALSE.

  END SUBROUTINE cas_init_visc

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

  SUBROUTINE cas_init_turb(kzcst,nzcst,z0scst,z0bcst)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_turb  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!!        !  2015-09 (F. Ganthy) Introduction of a spatialized bottom roughness length
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters,ONLY : rsh
   USE comvars2d, ONLY : z0b,z0s,l_incbotstress,botstressmax,frot_lin, &
                         hminfrot,hmaxfrot,l_z0b_hom,file_z0b
   USE comvars3d, ONLY : turb_nbeq,turb_0eq_option,turb_2eq_option,nz,kz, &
                         bgdiff
   USE comvarp3d, ONLY : ect,psi
   USE comturb,   ONLY : eps

!! * Arguments
   REAL(KIND=rsh), INTENT(IN) :: kzcst,nzcst,z0scst,z0bcst

!! * Local declarations

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

   turb_nbeq=0
   turb_0eq_option=1
   turb_2eq_option=0
   l_z0b_hom =.TRUE.
   file_z0b  ='./'
   z0b(:,:)  =z0bcst
   z0s(:,:)  =z0scst
   ect(:,:,:)=0.0_rsh
   eps(:,:,:)=0.0_rsh
   psi(:,:,:)=0.0_rsh
   nz(:,:,:) = nzcst
   kz(:,:,:) = kzcst
   l_incbotstress=.FALSE.
   frot_lin=0.0_rsh  ! used if key_linear_bottom_stress
   botstressmax=100.0_rsh
   hminfrot=0.5_rsh
   hmaxfrot=5.0_rsh
   bgdiff=0.0_rsh

  END SUBROUTINE cas_init_turb

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

  SUBROUTINE cas_init_tspar
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_tspar  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg
   USE comvars2d,  ONLY : icon,qmax,kx,ky,hminkxky
   USE comvars3d,  ONLY : saliref,tetaref,rhoref,saliref_lin,tetaref_lin,qmaxz, &
                          l_stateeq_lin
   USE cometemp,   ONLY : phitemp,phisal

!! * Arguments

!! * Local declarations

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

   icon=2
   kx=0.0_rsh
   ky=0.0_rsh
   hminkxky=0.0_rsh
   qmax=120000.0_rsh
   qmaxz=1000.0_rsh
   rhoref=1027.34_rsh
   saliref=35.5_rsh
   tetaref=10.0_rsh
   l_stateeq_lin=.TRUE.
   saliref_lin=saliref
   tetaref_lin=tetaref
   phisal(:,:,:) = 0.0_rlg
   phitemp(:,:,:) = 0.0_rlg

  END SUBROUTINE cas_init_tspar

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

  SUBROUTINE cas_init_progvar
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_progvar  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : valmanq,rsh,ljmin,ljmax,limin,limax,liminm2,limaxp2,ljminm2,ljmaxp2, &
                          liminm1,ljminm1,kmax
   USE comvars2d,  ONLY : h0
   USE comvarp2d,  ONLY : sal,temp,ssh,u,v
   USE comvarp3d,  ONLY : uz,vz
   USE comvars3d,  ONLY : wz,saliref,tetaref

!! * Arguments

!! * Local declarations
   INTEGER :: k

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

   ssh(:,:)  = MAX(0.0_rsh,-h0(liminm2:limaxp2,ljminm2:ljmaxp2))
   u(:,:)    = 0.0_rsh
   v(:,:)    = 0.0_rsh
   uz(:,:,:) = 0.0_rsh
   vz(:,:,:) = 0.0_rsh
   wz(:,:,:) = 0.0_rsh
   temp(:,:,:) = tetaref
   sal(:,:,:) = saliref

! over land
   WHERE( h0(liminm2:limaxp2,ljminm2:ljmaxp2)==-valmanq ) ssh(liminm2:limaxp2,ljminm2:ljmaxp2)=valmanq
   DO k=1,kmax
     WHERE( h0(liminm1:limaxp2,ljminm1:ljmaxp2)==-valmanq ) sal(k,liminm1:limaxp2,ljminm1:ljmaxp2)=valmanq
     WHERE( h0(liminm1:limaxp2,ljminm1:ljmaxp2)==-valmanq ) temp(k,liminm1:limaxp2,ljminm1:ljmaxp2)=valmanq
   END DO

  END SUBROUTINE cas_init_progvar

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

  SUBROUTINE cas_init_misc
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_misc  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg,kmax,ljmin,ljmax,limin,limax,imin,imax,jmin,jmax,valmanq, &
                          liminm1,limaxp1,ljminm1,ljmaxp1,liminm2,limaxp2,ljminm2,ljmaxp2, &
                          ljminp1,liminp2,limaxm1,ljminp2,ljmaxm1,epsilon
   USE comvars2d,  ONLY : fwet,fwetp,fwetmin,dhxy,h0,jb,jh,ig,id,h0,h0_g,hm, &
                          htx,hty,hex,hey,hx,hy,dhx,dhy,goulet,govlet,halfdt
   USE comvars2d,  ONLY : l_obc_cycl_x,l_obc_cycl_y
   USE comvarp2d,  ONLY : ssh
   USE comfluxdo, ONLY : fxeau,fyeau
   USE comvarp3d, ONLY : uz,vz

   USE comcolig2d,ONLY : igdu,igdv
   USE toolmpi, ONLY: ex_i_rsh,ex_j_rsh,ex_i_int,ex_j_int

!! * Arguments

!! * Local declarations
   INTEGER :: i,j,k

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

   !PRINT*, 'enter CAS_INIT_MISC routine'

! Domain of calculation
! ---------------------
   DO i=imin,imax
     jb(i)=jmin+1
     IF ( .NOT. l_obc_cycl_y) THEN
       DO j = jmin+1,jmax-2
         IF (h0_g(i,j) > -valmanq .OR. h0_g(i,j+1) > -valmanq .OR. h0_g(i,j+2) > -valmanq) THEN
           jb(i) = j
           GOTO 1050
         ELSE
           jb(i) = j
         END IF
       END DO
 1050 CONTINUE
     END IF
   END DO

   DO i=imin,imax
     jh(i)=jmax
     IF ( .NOT. l_obc_cycl_y) THEN
       DO j = jmax-2,jmin+1,-1
         IF (h0_g(i,j) > -valmanq .OR. h0_g(i,j+1) > -valmanq .OR. h0_g(i,j+2) > -valmanq) THEN
           jh(i) = j + 2
           GOTO 1100
         ELSE
           jh(i) = j
         END IF
       END DO
 1100 CONTINUE
     END IF
   END DO

   DO j=jmin,jmax
     ig(j)=imin+1
     IF ( .NOT. l_obc_cycl_x) THEN
       DO i = imin+1,imax-2
         IF (h0_g(i,j) > -valmanq .OR. h0_g(i+1,j) > -valmanq .OR. h0_g(i+2,j) > -valmanq) THEN
           ig(j) = i
           GOTO 1200
         ELSE
           ig(j) = i
         END IF
       END DO
 1200 CONTINUE
     END IF
   END DO

   DO j=jmin,jmax
     id(j)=imax
     IF ( .NOT. l_obc_cycl_x) THEN
       DO i = imax-2,imin+1,-1
         IF (h0_g(i,j) > -valmanq .OR. h0_g(i+1,j) > -valmanq .OR. h0_g(i+2,j) > -valmanq) THEN
           id(j) = i + 2
           GOTO 1300
         ELSE
           id(j) = i
         END IF
       END DO
 1300  CONTINUE
     END IF
   END DO


! Fraction of wet surface on each cell
! ------------------------------------

   DO j=ljmin,ljmax
     DO i=limin,limax
       IF (dhxy(i,j) == 0.0_rsh) THEN
          fwet(i,j)=1.0_rsh
         fwetp(i,j)=1.0_rsh
       ELSE
          fwet(i,j)=MIN(MAX(fwetmin,(ssh(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
         fwetp(i,j)=MIN(MAX(fwetmin,(ssh(i,j)+h0(i,j))/dhxy(i,j)),1.0_rsh)
       ENDIF
     END DO
   END DO

! Water depth at T and u locations
! --------------------------------

! Attention OK que si ssh=cst car on est a position u ou v
   hex(limin:limax,ljmin:ljmax)=(hx(limin:limax,ljmin:ljmax)+ssh(limin:limax,ljmin:ljmax)+ &
                                dhx(limin:limax,ljmin:ljmax))*goulet(limin:limax,ljmin:ljmax)
   hey(limin:limax,ljmin:ljmax)=(hy(limin:limax,ljmin:ljmax)+ssh(limin:limax,ljmin:ljmax)+ &
                                dhy(limin:limax,ljmin:ljmax))*govlet(limin:limax,ljmin:ljmax)
   htx(limin:limax,ljmin:ljmax)=hx(limin:limax,ljmin:ljmax)+ssh(limin:limax,ljmin:ljmax)
   hty(limin:limax,ljmin:ljmax)=hy(limin:limax,ljmin:ljmax)+ssh(limin:limax,ljmin:ljmax)
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljmin,ljmax,hex(liminm2:limaxp2,ljmin:ljmax))
   CALL ex_j_rsh(-2,2,1,limin,limax,ljminm2,ljmaxp2,hey(limin:limax,ljminm2:ljmaxp2))
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,htx(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,htx(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hty(liminm2:limaxp2,ljminm2:ljmaxp2))
   CALL ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,hty(liminm2:limaxp2,ljminm2:ljmaxp2))


   igdu(liminm1:limaxp1,ljminm1:ljmaxp1)=4
   DO j=ljminp2,ljmaxm1
     DO i=limin,limax
       IF ((i>=ig(j)) .AND. (i<=id(j)-1)) THEN
         IF (htx(i,j) > 0.0_rsh) THEN
           igdu(i,j)=1
           IF (htx(i,j) < hm) igdu(i,j)=5
           IF (htx(i,j) < MIN(epsilon,REAL(hm,rlg))) igdu(i,j)=4
         END IF
       END IF
     END DO
   END DO
   igdv(liminm1:limaxp1,ljminm1:ljmaxp1)=4
   DO j=ljminp1,ljmaxm1
     DO i=liminp2,limaxm1
       IF ( j>=jb(i) .AND. j<=jh(i)-1) THEN
         IF (hty(i,j) > 0.0_rsh) THEN
           igdv(i,j)=1
           IF (hty(i,j) < hm) igdv(i,j)=5
           IF (hty(i,j) < MIN(epsilon,REAL(hm,rlg))) igdv(i,j)=4
         END IF
       END IF
     END DO
   END DO
   CALL ex_i_int(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,igdu(liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL ex_j_int(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,igdu(liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL ex_i_int(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,igdv(liminm1:limaxp1,ljminm1:ljmaxp1))
   CALL ex_j_int(-1,1,1,liminm1,limaxp1,ljminm1,ljmaxp1,igdv(liminm1:limaxp1,ljminm1:ljmaxp1))



   DO k=1,kmax
     fxeau(k,liminm1:limax,ljmin:ljmax)=hex(liminm1:limax,ljmin:ljmax)*uz(k,liminm1:limax,ljmin:ljmax)*halfdt
     fyeau(k,limin:limax,ljminm1:ljmax)=hey(limin:limax,ljminm1:ljmax)*vz(k,limin:limax,ljminm1:ljmax)*halfdt
   END DO


   !PRINT*, 'exit CAS_INIT_MISC routine'

  END SUBROUTINE cas_init_misc

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

  SUBROUTINE cas_init_traj
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_traj  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules


!! * Arguments

!! * Local declarations

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




  END SUBROUTINE cas_init_traj

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

  SUBROUTINE cas_init_subs(date_start,date_start_sub)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_subs  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description : pour substances dissoutes seules
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E       !  2013-06  (F. Grasso) Modifications for using sediments
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax,obc_width, &
                          rsh,rlg,liminm1,limaxp2,ljminm1,ljmaxp2
   USE comvars2d, ONLY : hm,iscreenlog,ierrorlog
   USE comvars3d, ONLY : dsigu
   USE parameters, ONLY : nb_var
   USE parameters, ONLY : liminm2,ljminm2,limaxp3,ljmaxp3   ! VALERIE PROVISOIRE
   USE comsubstance





!! * Arguments
   CHARACTER(LEN=19), INTENT(IN) :: date_start
   CHARACTER(LEN=19), INTENT(IN),OPTIONAL :: date_start_sub

!! * Local declarations
   INTEGER :: n,k,isubs,IERR_MPI,iv,iiv,indgrp
   REAL(KIND=rlg)    :: tool_datosec

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

   nv_adv=nb_var
   nv_state=nb_var
   nv_dri=0
   nv_int=0
   nv_tot=nv_state+nv_dri+nv_int

   ALLOCATE(cv_wat(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   cv_wat(:,:,:,:)=0.0_rsh
   ALLOCATE(f90(nb_var))
   f90(1:nb_var) = 1.0_rsh
   ALLOCATE(t90(nv_state))
   t90(1:nv_state) = 0.0_rsh
   ALLOCATE(irk_fil(nv_tot))
   DO n=1,nv_tot
     irk_fil(n)=n
   END DO
   ALLOCATE(l_out_subs(nv_state))
   ALLOCATE(l_subs2D(nv_state))
   ALLOCATE(obc_cv_name(nv_state))
   obc_cv_name(:)='none'
   IF (nvp > 0) THEN
     ALLOCATE(ws_free_opt(nvp)) 
     ALLOCATE(ws_hind_opt(nvp)) 
     ALLOCATE(ws_free_para(4,nvp)) 
     ALLOCATE(ws_free_min(nvp)) 
     ALLOCATE(ws_free_max(nvp)) 
     ALLOCATE(ws_hind_para(2,nvp)) 
     ALLOCATE(tocd(nvp))
   ENDIF
   ALLOCATE(typart(nv_state))
   ALLOCATE(typdiss(nv_adv))
   ALLOCATE(itypv(nv_tot))

   unit_var(:)='unknown'
   long_name_var(:)=name_var(:)
   standard_name_var(:)=long_name_var(:)
   valid_min_var(:)=-3.4e12
   valid_max_var(:)=3.4e12
   IF ( PRESENT(date_start_sub) ) THEN
      tdebsubs=tool_datosec(date_start_sub)
   ELSE
      tdebsubs=tool_datosec(date_start)
   END IF
   l_out_subs(:)=.TRUE.

   ALLOCATE(sub_flx_atm(nv_adv))
   sub_flx_atm(:)=0.0_rsh
   ALLOCATE(cv_rain(nv_adv))
   cv_rain(:)=0.0_rsh
   ALLOCATE(cini_wat(nv_state))
   cini_wat(:)=0.0_rsh
   ALLOCATE(obc_cv_wat_e(nv_adv,kmax,imax-obc_width+1:imax,ljmin:ljmax))
   ALLOCATE(obc_cv_wat_w(nv_adv,kmax,imin:imin+obc_width,ljmin:ljmax))
   ALLOCATE(obc_cv_wat_n(nv_adv,kmax,limin:limax,jmax-obc_width+1:jmax))
   ALLOCATE(obc_cv_wat_s(nv_adv,kmax,limin:limax,jmin:jmin+obc_width))
   obc_cv_wat_e(:,:,:,:)=0.0_rsh
   obc_cv_wat_w(:,:,:,:)=0.0_rsh
   obc_cv_wat_n(:,:,:,:)=0.0_rsh
   obc_cv_wat_s(:,:,:,:)=0.0_rsh
   ALLOCATE(phicon(nv_tot,kmax,limin:limax,ljmin:ljmax))
   phicon(:,:,:,:) = 0.0_rsh
   ALLOCATE(phicon_drycell(nv_tot,kmax,limin:limax,ljmin:ljmax))
   phicon_drycell(:,:,:,:) = 0.0_rsh
   ALLOCATE(phicon_atm(nv_tot,limin:limax,ljmin:ljmax))
   phicon_atm(:,:,:)=0.0_rsh
   sflx_sub_atm_depth=1.0_rsh
! Allocate intermediate variable used in ultimate quickest macho scheme and macho3d organisation
   ALLOCATE(c2(nv_tot,kmax,liminm1:limaxp2,ljminm1:ljmaxp2))
   c2(:,:,:,:)=0.0_rsh



! initialialization as in cas_init_subs_disspart
     unit_modif_mudbio_N2dw(:)=1.0_rsh







   CALL alloc_fxfyfz

   l_bilan=.FALSE.
   call cas_init_subs_budget(date_start)

  END SUBROUTINE cas_init_subs


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


  SUBROUTINE cas_init_subs_budget(date_start)

   USE comsubstance

!! * Arguments
   CHARACTER(LEN=19), INTENT(IN) :: date_start

   REAL(KIND=rlg)    :: tool_datosec

! Initialization of budget parameters
! -----------------------------------
   nb_border=1
   dtout_budget=2.0_rlg  ! in hours
   tdeb_budget=tool_datosec(date_start)
   iout_budget=-1
   name_bilfil = ''

  END SUBROUTINE cas_init_subs_budget


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


  SUBROUTINE cas_init_subs_disspart(l_filesubs)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_subs_disspart  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E       !  2013-06  (F. Grasso) Modifications for using sediments
!&E       !  2014-01  (B. Thouvenin) Modifications for using sediments and biolo
!&E       !  2013-06  (F. Grasso) Modifications for using sediments
!&E       !  2014-01-05  (B. Thouvenin) Modifications for using sediments and biolo
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax,obc_width, &
                          rsh,rlg,liminm1,limaxp2,ljminm1,ljmaxp2
   USE comvars2d, ONLY : hm,iscreenlog,ierrorlog,dt
   USE comvars3d, ONLY : dsigu
   USE parameters, ONLY : nb_var
   USE comvarp2d,  ONLY : sal,temp
   USE comsubstance


   USE subreaddat, ONLY : sub_read_var,sub_read_vardiag

   USE comMUSTANG, ONLY : MUSTANG_alloc,ros,diam_sed


   USE comMUSTANG, ONLY : l_outsandrouse



   USE toolmpi, ONLY : MASTER

!! * Arguments
   LOGICAL, INTENT(IN) :: l_filesubs

!! * Local declarations
   INTEGER :: n,k,isubs,IERR_MPI,iv,iiv,indgrp,ivpc
   REAL(KIND=rlg)    :: tool_datosec
   LOGICAL           :: l_varassoc

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

   IF (l_filesubs) THEN   !we use "sub_read_var"
     write(*,*)'We are reading data file :',name_filesubs



     CALL sub_read_var
     CALL init_tssub(sal,temp,cv_wat)
     DO isubs=1,nv_adv
       obc_cv_wat_w(isubs,:,:,:)=cini_wat(isubs) 
       obc_cv_wat_e(isubs,:,:,:)=cini_wat(isubs) 
       obc_cv_wat_s(isubs,:,:,:)=cini_wat(isubs)
       obc_cv_wat_n(isubs,:,:,:)=cini_wat(isubs)
     ENDDO


     CALL MUSTANG_alloc(l_filesubs) 



     
   ELSE
     nvpc=nv_cp+nv_mud+nv_sand+nv_grav
     nvp=nvpc+nv_ncp+nv_sorb
     nv_adv=nvp+nv_dis
     nv_state=nv_adv+nv_fix
     nv_tot=nv_state+nv_dri+nv_int
! identification of igrav1,igrav2,isand1,isand2,imud1,imud2
     igrav1=1
     igrav2=nv_grav
     isand1=igrav2+1
     isand2=igrav2+nv_sand
     imud1=isand2+1
     imud2=isand2+nv_mud
     typdiss(nvp+1:nv_adv)=1.
     typdiss(1:nvp)=0.


     CALL MUSTANG_alloc(l_filesubs) 


! Mise en groupe des variables particulaires
! ------------------------------------------
! Allocation des tableaux se rapportant aux variables non dissoutes
     ALLOCATE(iv_ws_max(nvp,nvp))
     ALLOCATE(nv_ws_max(nvp))
     ALLOCATE(ndt_part(nvp))
     ALLOCATE(sv_part(nvp))
     ALLOCATE(ap_ws_gr(nvp))
     ALLOCATE(ws3(kmax,nv_adv,limin:limax,ljmin:ljmax))
     ALLOCATE(irkm_var_assoc(nvp))

     ALLOCATE(l_outsandrouse(nvp))



! Initialisation des différentes variables et tableaux
     nb_ws_max=0
     iv_ws_max(:,:)=0
     nv_ws_max(:)=0
     ap_ws_gr(:)=0
     ndt_part(:)=0
     sv_part(:)=0.0_rsh
     ws3(:,:,:,:) = 0.0_rsh
     irkm_var_assoc(:)=0
     l_subs2D(:)=.false.

     l_outsandrouse(:)=.false.

     unit_modif_mudbio_N2dw(:)=1.0_rsh
     

! If non constitutive particulate variables : Initialization of irkm_var_assoc
! ---------------------------------
! reperage de variables particulaires constitutives associees aux variables part. SORB
! ------------------------------------------------------------------------------------
   DO iv=1,nv_sorb
     isubs=nvpc+nv_ncp+iv
     irkm_var_assoc(isubs)=0
     l_varassoc=.FALSE.
     DO ivpc=1,nvpc
       IF ( TRIM(ADJUSTL(ADJUSTR(name_var(irk_fil(ivpc))))) == TRIM(ADJUSTL(ADJUSTR(name_varpc_assoc(irk_fil(isubs))))) )THEN
         irkm_var_assoc(isubs)=ivpc
         IF (MASTER) THEN
           WRITE(iscreenlog,*) ' '
           WRITE(iscreenlog,*)'constitutive part. variable where is sorbed'
           WRITE(iscreenlog,*)'the variable :',name_var(irk_fil(isubs))
           WRITE(iscreenlog,*)'is the variable ', name_var(irk_fil(ivpc))
         ENDIF
         l_varassoc=.TRUE.
         ws_free_para(:,irk_fil(isubs))=ws_free_para(:,irk_fil(ivpc))
         ws_hind_para(:,irk_fil(isubs))=ws_hind_para(:,irk_fil(ivpc))
         ws_free_opt(irk_fil(isubs))=ws_free_opt(irk_fil(ivpc))
         ws_free_opt(irk_fil(isubs))=ws_free_opt(irk_fil(ivpc))
         ws_free_max(irk_fil(isubs))=ws_free_max(irk_fil(ivpc))
         ws_free_min(irk_fil(isubs))=ws_free_min(irk_fil(ivpc))

       END IF
     END DO
     IF (.NOT.l_varassoc) THEN
       IF (MASTER) THEN
         WRITE(iscreenlog,*)' '
         WRITE(iscreenlog,*)'the SORB variable :',name_var(irk_fil(isubs))
         WRITE(iscreenlog,*)'does not have associated constitutive part. variable'
         WRITE(iscreenlog,*)'See variable.dat to give exactly the name of the constitutive associated variable'
         WRITE(iscreenlog,*)'otherwise, it is not a SORB variable, but a NoCP variable '
         WRITE(ierrorlog,*)' '
         WRITE(ierrorlog,*)'the SORB variable :',name_var(irk_fil(isubs))
         WRITE(ierrorlog,*)'does not have associated constitutive part. variable'
         WRITE(ierrorlog,*)'See variable.dat to give exactly the name of the constitutive associated variable'
         WRITE(ierrorlog,*)'otherwise, it is not a SORB variable, but a NoCP variable '
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       ENDIF
     END IF
   END DO

! Equation of substances with settling velocity
! is solved with a sub time step

! estimate of the minimum layer thickness
! ---------------------------------------

     dhmin=hm

     DO k=1,kmax
       dhmin=MIN(dhmin,hm*dsigu(k))
     ENDDO

     IF(dhmin<=0_rsh)dhmin=0.10_rsh



     nb_ws_max=0
     DO iv=1,nvp
       nb_ws_max=iv     ! number of variables with different settling velocities
! all of them for 1 because ws3 is a 4D variable
       nv_ws_max(iv)=1  ! each particulate variable is relative to a unique group
! because ws3 is a 4D variable
       sv_part(iv)=ws_free_max(iv)   ! important to determinate the number of iterations
! index of variable associated with a settling velocity
       iv_ws_max(iv,1)=iv
       IF(iv < imud1 ) THEN
           ndt_part(nb_ws_max)=1  ! will be calculated after estimation of settling velocity
       ELSE
         IF (ws_free_max(iv) == 0.0_rsh) THEN
           IF (MASTER) THEN
           WRITE(ierrorlog,*) ''
           WRITE(ierrorlog,*) ' ******* Error detected from SUB_READ_VAR from SUBREADDAT.F90 ******* '
           WRITE(ierrorlog,*) 'the variable :',name_var(irk_fil(iv))
           WRITE(ierrorlog,*) 'is a particulate one'
           WRITE(ierrorlog,*) 'but its settling velocity is set to 0.0 in ',TRIM(name_filesubs), ' file.'
           WRITE(ierrorlog,*) 'Consequently, the sub time step used to solve substance equation cannot be estimated'
           WRITE(ierrorlog,*) 'The simulation is stopped'
           CALL MPI_FINALIZE(IERR_MPI)
           STOP
           ENDIF
         END IF
         ndt_part(nb_ws_max)=NINT(dt/MIN((dhmin/ws_free_max(iv)),dt))
       END IF
       PRINT*, 'ndt_part dans castest',nb_ws_max,ndt_part(nb_ws_max)
     ENDDO




   nvpi=nvp
   IF (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'SUBSTANCES WITH A TIME CONSTANT SETTLING VELOCITY '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*)'NUMBER GROUP OF SUBSTANCES            : ',nb_ws_max
   ENDIF
   DO iiv=1,nb_ws_max
     IF (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*)'group number                   : ',iiv
     WRITE(iscreenlog,*)'number of substances           : ',nv_ws_max(iiv)
     WRITE(iscreenlog,*)'group settling velocity        : ',sv_part(iiv)
     WRITE(iscreenlog,*)'group number of iterations     : ',ndt_part(iiv)
     ENDIF
     DO iv=1,nv_ws_max(iiv)
       IF (MASTER) THEN
       WRITE(iscreenlog,*) 'substance number        name'
       WRITE(iscreenlog,'(4x,i4,6x,a18,2x)')  &
       iv,TRIM(name_var(irk_fil(iv_ws_max(iiv,iv))))
       ENDIF
     ENDDO
   ENDDO
  ENDIF

  END SUBROUTINE cas_init_subs_disspart


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


  SUBROUTINE cas_init_sedim(l_filesubs,l_csedinit_varz,l_dzsinit_varz,sinit_sed, &
                            tinit_sed,z0sedinit,dzsmaxinit,l_casinit_hsed,l_hsedinit_cste)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E       !  2014-05  (F. Grasso) Modifications for using sediments in test case consolidation
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : valmanq,imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax, &
                          rsh,rlg,ksdmax
   USE comvars2d, ONLY : h0,ig,id,iscreenlog,ierrorlog,dt
   USE comsubstance, ONLY : nvpc,nvp,nv_adv,nv_state,nv_tot,tdebsubs,irkm_var_assoc, &
                            irk_mod,isand1,isand2,igrav1,igrav2,imud1,imud2,         &
                            nb_ws_max,iv_ws_max,ndt_part,dhmin,irk_fil,unit_modif_mudbio_N2dw
   USE comaree
   USE obccombine
   USE toolmpi, ONLY : ADD_TAG,ex_i_rsh,ex_j_rsh,MASTER,min_all_mpi_rsh,max_all_mpi_rsh
   USE comMUSTANG
   USE initMUSTANG, ONLY : MUSTANG_init_ws_sand,MUSTANG_init_hsed
   USE comvars3d,    ONLY : rhoref



!! * Arguments
   LOGICAL, INTENT(IN)                :: l_filesubs
   LOGICAL, INTENT(IN), OPTIONAL      :: l_csedinit_varz,l_dzsinit_varz,l_casinit_hsed,l_hsedinit_cste
   REAL(KIND=rsh),INTENT(IN),OPTIONAL :: sinit_sed,tinit_sed,z0sedinit,dzsmaxinit

!! * Local declarations
   INTEGER :: i,j,k,iv,ivg,imud_nobio,IERR_MPI
   LOGICAL :: l_csedinit_Cste,l_dzsinit_Cste,l_hsed_cste
   REAL(KIND=rsh)    :: somalp,dzsf,cini_modif,tot_csed_init
   REAL(KIND=rlg)    :: tool_datosec

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

     CALL MUSTANG_init_ws_sand(rhoref)

! Initialization of sedimentary parameters
! ----------------------------------------
   l_fricwave=.FALSE.
   l_repsed=.FALSE.
   l_initsed_vardiss=.FALSE.

!pre-assignement    l_consolid=.FALSE.

   IF(PRESENT(l_csedinit_varz)) THEN
     l_csedinit_Cste=.NOT. l_csedinit_varz
   ELSE
     l_csedinit_Cste=.TRUE.
   ENDIF
   IF(PRESENT(l_dzsinit_varz)) THEN
     l_dzsinit_Cste=.NOT. l_dzsinit_varz
   ELSE
     l_dzsinit_Cste=.TRUE.
   ENDIF
  IF(PRESENT(sinit_sed)) THEN
    sini_sed=sinit_sed
  ELSE
    sini_sed=0.0_rsh
  ENDIF
  IF(PRESENT(tinit_sed)) THEN
    tini_sed=tinit_sed
  ELSE
    tini_sed=10.0_rsh
  ENDIF
  IF(PRESENT(z0sedinit)) THEN
    z0seduni=z0sedinit
  ELSE
    z0seduni = 0.0001_rsh
  ENDIF
  IF(PRESENT(dzsmaxinit)) THEN
    dzsmaxuni=dzsmaxinit
  ELSE
    dzsmaxuni = 1.0e-3
  ENDIF
  IF(PRESENT(l_casinit_hsed)) THEN
    l_init_hsed=l_casinit_hsed
  ELSE
    l_init_hsed = .TRUE.
  ENDIF
   IF(PRESENT(l_hsedinit_cste)) THEN
    l_hsed_cste=l_hsedinit_cste
  ELSE
    l_hsed_cste = .TRUE.
  ENDIF

   l_diffused=.FALSE.
   l_bioturb=.FALSE.
   l_dyn_insed=.FALSE.
   l_unised = .TRUE.
   fricwav=0.06_rsh
!pre-assignement    hseduni = 0.0_rsh
!pre-assignement    cseduni = 200.0_rsh
!pre-assignement    ksmiuni = 1
!pre-assignement    ksmauni = 0       ! no sediment layer
   dzsmin = 1.0e-5
   cfreshmud = 300.0_rsh
   xdifs1 = 1.e-8
   xdifs2 = 0.0_rsh
   xdifsi1 = 1.e-6
   xdifsi2 = 0.0_rsh
   epdifi = 0.01_rsh
   aref_sand=0.02_rsh
   

   csedmin=30.0_rsh
! csed_mud_ini : mud concentration in inital sediment (=cfreshmud if no consolidated sediment)
   csed_mud_ini=300.0_rsh
   xperm1=4.0e-12        ! best param for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
   xperm2=-6.0_rsh       ! best param for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
   xsigma1=6.0e5         ! best param for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
   xsigma2=6.0_rsh       ! best param for Co<120 kg/m3 (Grasso et al., Ocean Dynamics 2014)
   dt_sed=60.0_rsh
   dtcal_sed=30.0_rsh
   fexcs = 0.50_rsh
   cexcs=1._rsh-fexcs
   csegreg=250.0_rsh
   csandseg=1250.0_rsh
   activlayer=0.02_rsh
   frmudcr2=0.7_rsh
   coef_frmudcr1=1000._rsh
   x1toce_mud=0.3_rsh
   x2toce_mud=0._rsh
   E0_sand_option=0
   E0_sand_para=0.0002_rsh
   n_eros_sand=1.6_rsh
   E0_mud=0.0002_rsh
   n_eros_mud=1._rsh
   corfluer1=0.01_rsh
   corfluer2=350._rsh
   cmudcr=200._rsh
   ero_option = 1
   xexp_ero = 2.0_rsh
! florence
   cmudcr=600._rsh
   cvolmaxsort=0.58_rsh
   cvolmaxmel=0.67_rsh
   xsigma1sg=xsigma1/9.81_rsh
   MF=1.0_rsh
   htncrit_eros=0.0_rsh
   slopefac=0.01_rsh
   coef_erolat=0.002_rsh
   coef_tenfon_lat=5._rsh
   l_erolat_wet_cell=.FALSE.
   l_bathy_smoothing=.FALSE.
   choice_flxdiss_diffsed=3 ! formulation Boudreau/Shaw-Hanratty
   l_bioturb=.FALSE.
   xbioturbmax=1.157e-12_rsh
   xbioturbk=6.0_rsh
   dbiotu0=0.01_rsh
   dbiotum=0.002_rsh
   frmud_db_min=0.6_rsh
   frmud_db_max=0.8_rsh
   

   fws2=fricwav*0.5_rsh


! Initialization of the sediment layer
! ------------------------------------
   z0sed(:,:)=z0seduni
   ksmi(:,:)=1  !0
   WHERE (h0(limin:limax,ljmin:ljmax) /= -valmanq) ksmi(limin:limax,ljmin:ljmax)=ksmiuni
   IF(l_hsed_cste) THEN
     hsed(:,:)=0.0_rsh
     ksma(:,:)=0
     WHERE (h0(limin:limax,ljmin:ljmax) /= -valmanq) hsed(limin:limax,ljmin:ljmax)=hseduni
     WHERE (h0(limin:limax,ljmin:ljmax) /= -valmanq) ksma(limin:limax,ljmin:ljmax)=ksmauni
   ENDIF
   WHERE (h0(limin:limax,ljmin:ljmax) /= -valmanq) dzsmax(limin:limax,ljmin:ljmax)=dzsmaxuni

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! ATTTENTION: the total concentration is given as an input (cseduni), and the
! contribution of each variable (cini_sed(iv)) is given as a fraction [0-1] of cseduni
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



   IF(l_csedinit_Cste) THEN
     tot_csed_init=0.0_rsh
     DO iv=1,nvpc
! constitutive variables : cini_sed in fraction of cseduni
        tot_csed_init=tot_csed_init+cini_sed(iv)              
     ENDDO
     IF(tot_csed_init .NE. 1.0_rsh) THEN
        WRITE(ierrorlog,*)' '
        WRITE(ierrorlog,*)' inital concentrations of constitutive particulate matter'
        WRITE(ierrorlog,*)' have not been well defined'
        WRITE(ierrorlog,*)' See variable.dat to give fraction of total concentration '
        WRITE(ierrorlog,*)' for each constitutive variable (GRAV, SAND, MUD) (but not MUDB) '
        WRITE(ierrorlog,*)' : the sum of all must be =1 '
        CALL MPI_FINALIZE(IERR_MPI)
        STOP         
     ENDIF
   ENDIF
      
! ajustement du hsed et cv_sed
! ==============================

   IF (l_init_hsed .AND. l_csedinit_Cste ) THEN
    IF(hseduni > 0.0_rsh) THEN
        ALLOCATE(cv_sedini(1:nvpc))
!!! Concentration massique
        DO iv=1,nvpc
! constitutive variables : cini_sed in fraction of cseduni
           cv_sedini(iv)=cseduni*cini_sed(iv)
        ENDDO
        CALL MUSTANG_init_hsed(cv_sedini)
        cseduni=cseduni*hseduni/hsed_new         
        hseduni=hsed_new
        WRITE(*,*)'NEW TOTAL SEDIMENT THICKNESS AFTER ADJUSTEMENT = ',hsed_new
        WRITE(*,*)'NEW TOTAL CONCENTRATION SEDIMENT AFTER ADJUSTEMENT = ',cseduni
        DEALLOCATE(cv_sedini)
    ENDIF
   ENDIF
   WHERE (h0(limin:limax,ljmin:ljmax) /= -valmanq) hsed(limin:limaxp1,ljmin:ljmaxp1)=hseduni



   DO j=ljmin,ljmax
      DO i=MAX0(limin,ig(j)),MIN0(limax,id(j))
         IF(h0(i,j) /= -valmanq) THEN
            DO k=ksmi(i,j),ksma(i,j)
              somalp=0.
              c_sedtot(k,i,j)=0.0_rsh
              DO iv=1,nvpc
! constitutive variables : cini_sed in % of cseduni
                IF(l_csedinit_Cste)cv_sed(iv,k,i,j)=cseduni*cini_sed(iv)        
                c_sedtot(k,i,j)=c_sedtot(k,i,j)+cv_sed(iv,k,i,j)
                somalp=somalp+cv_sed(iv,k,i,j)/ros(iv)
              ENDDO
              poro(k,i,j)=1.0_rsh-somalp
              IF(l_csedinit_Cste) THEN
                DO iv=nvpc+1,nvp
! non constitutive variables : cini_sed in conc/kg of sediment
                  IF (irkm_var_assoc(iv) >0) THEN
                    cv_sed(iv,k,i,j)=cini_sed(iv)*cv_sed(irkm_var_assoc(iv),k,i,j)
                  ELSE
                    cv_sed(iv,k,i,j)=cini_sed(iv)*c_sedtot(k,i,j)
                  END IF
                ENDDO
                DO iv=nvp+1,nv_state
! dissolved variables : cini_sed in conc/m3 of pore water
! + fixed variables : cini_sed in conc/m3 of sediment
                  cv_sed(iv,k,i,j)=cini_sed(iv)
                ENDDO
                cv_sed(-1,k,i,j)=tini_sed
                cv_sed( 0,k,i,j)=sini_sed
              ENDIF
            ENDDO
         ENDIF
      ENDDO
   ENDDO

!   PRINT*, 'maxval(c_sedtot) ',maxval(c_sedtot)
!   PRINT*, 'maxval(cv_sed) ',maxval(cv_sed)


!   vertical discretization within the sediment
!   -------------------------------------------
   IF(l_dzsinit_Cste) THEN
     dzs(:,:,:)=0.0_rsh
     DO j=ljmin,ljmax
      DO i=MAX0(limin,ig(j)),MIN0(limax,id(j))
        IF (ksma(i,j).NE.0) THEN
          dzsf=hsed(i,j)/float(ksma(i,j)-ksmi(i,j)+1)
          DO k=ksmi(i,j),ksma(i,j)
            dzs(k,i,j)=dzsf
          END DO
        END IF
      END DO
     END DO
  ENDIF

!!flo IF(l_morphocoupl) CALL sed_morphoinit

!!CALL sed_settlveloc_ini  ! not used now


!  initialization of consolidation/diffusion/bioturbation
!  -------------------------------
   IF (l_consolid) THEN

     tstart_sed=tool_datosec(date_start_sed)
     tstart_sed=MAX(tstart_sed,tdebsubs)
     t_sed=tstart_sed+dt_sed

   END IF

! initalization of sedim output
!-------------------------------
  choice_nivsed_out=1
  nk_nivsed_out=ksdmax
  ep_nivsed_out(:)=0.0_rsh
  epmax_nivsed_out=0.0_rsh
  l_outsed_saltemp=.TRUE.

!ends key_flocmod

  END SUBROUTINE cas_init_sedim


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

  SUBROUTINE cas_init_river(nbriver)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : kmax,rsh
   USE comfluxdo,  ONLY : nbfl,file_river,tmfle,iriv,jriv,kriv,phiriver,amtfle,l_tempriv,l_flestat, &
                          debmoy,phiriver,river_thickness,nomriviere
   IMPLICIT NONE 

!! * Arguments
   INTEGER, INTENT(IN), OPTIONAL :: nbriver

!! * Local declarations

!!----------------------------------------------------------------------
!! * Executable part
   file_river='e'
   IF ( PRESENT(nbriver) ) THEN
     nbfl=nbriver
   ELSE
     nbfl=0
   END IF


  END SUBROUTINE cas_init_river

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


  SUBROUTINE cas_init_outflow

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg,nb_var
   USE comfluxdo,  ONLY : iriv,jriv,kriv,phiriver,debmoy
   USE comvars2d,  ONLY : tdeb,tfin,date_startshift,nvalmax
   USE comsubstance,ONLY : nbrej,xrej,yrej,krej,numriv,nucon,l_rejflux,valrej

   IMPLICIT NONE 

!! * Arguments

!! * Local declarations
   REAL(KIND=rlg)    :: tdb,tfi
   INTEGER           :: i,iv,irej
   REAL(KIND=rlg)    :: tool_datosec
!!----------------------------------------------------------------------
!! * Executable part

   IF(nbrej > 0 ) THEN                                         

     ALLOCATE(xrej(nbrej),yrej(nbrej),krej(nbrej))
     ALLOCATE(nucon(nbrej),numriv(nbrej))
     ALLOCATE(l_rejflux(nbrej),valrej(nvalmax,nbrej))
   

     xrej(:)=2
     yrej(:)=2
     krej(:)=1

! numero des variables rejetees
     numriv(:)=1
     l_rejflux(:)=.FALSE.
     nucon(1)=2              ! MES
! Flux apportes (en masse/s)
     valrej(:,1)=0.00_rsh
                                                   
     IF (nbrej > 1 .OR. nb_var > 1 ) THEN
       iv=0
! by default : each input corresponds to one variable
       DO irej=2,nbrej
         iv=min(iv+irej,nb_var)
         nucon(irej)=iv  
! Flux apportes (en masse/s)
         valrej(:,irej)=0.00_rsh
       ENDDO  
     ENDIF                                        



  ENDIF
  
  END SUBROUTINE cas_init_outflow


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

  SUBROUTINE cas_init_agrif
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules

!! * Arguments

!! * Local declarations

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

  END SUBROUTINE cas_init_agrif

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



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


  SUBROUTINE cas_init_wave
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE cas_init_wave  ***
!&E
!&E ** Purpose : initialize a test case
!&E
!&E ** Description :
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !           (V. Garnier) Original code : reorganization of the test cases
!&E
!&E---------------------------------------------------------------------
!! * Used modules
   USE parameters, ONLY : rsh,rlg
   USE wave,       ONLY : file_wave,longitude_name_wave,latitude_name_wave,l_wave_marsgrid, &
                          l_wave_rdir,l_wave_rfoc,l_wave_rubot,l_wave_ruubot,l_wave_rvubot, &
                          l_wave_rtp, l_wave_rhs,l_wave_rt0m1,l_wave_rt0m1,l_wave_rtauwx,l_wave_rtauwy, &
                          l_wave_rtauox,l_wave_rtauoy,l_wave_rwlm,l_wave_rbhd, &
                          l_wave_rfbb,l_wave_rtbb,l_wave_rtus,l_wave_russ
   USE wave,       ONLY : name_wave_hs,name_wave_ubot,name_wave_uubot,name_wave_vubot,name_wave_foc,name_wave_t0m1, &
                          name_wave_dir,name_wave_tp,name_wave_tauwx,name_wave_tauwy,name_wave_tauox,name_wave_tauoy, &
                          name_wave_fbb,name_wave_tbbx,name_wave_tbby, &
                          name_wave_wlm,name_wave_bhd,name_wave_ussx,name_wave_ussy, &
                          name_wave_tusx, name_wave_tusy
   USE wave,       ONLY : name_out_wave_dir,name_out_wave_hs,name_out_wave_ubot, &
                          name_out_wave_foc,name_out_wave_tp,name_out_wave_t0m1, &
                          name_out_wave_tauox,name_out_wave_tauoy,name_out_wave_wlm, &
                          name_out_wave_tauwx,name_out_wave_tauwy, &
                          name_out_wave_fbb,name_out_wave_tbbx,name_out_wave_tbby, &
                          name_out_wave_bhd,name_out_wave_ussx,name_out_wave_ussy, &
                          name_out_wave_usx,name_out_wave_usy,name_out_wave_usz, &
                          name_out_wave_tusx,name_out_wave_tusy, &
                          name_out_wave_ubotx,name_out_wave_uboty
   USE wave,       ONLY : riog_valid_min_wave_dir,riog_valid_max_wave_dir,   &
                          riog_valid_min_wave_ubot,riog_valid_max_wave_ubot, &
                          riog_valid_min_wave_hs,riog_valid_max_wave_hs,     &
                          riog_valid_min_wave_foc,riog_valid_max_wave_foc,   &
                          riog_valid_min_wave_tp,riog_valid_max_wave_tp,     &
                          riog_valid_min_wave_t0m,riog_valid_max_wave_t0m,   &
                          riog_valid_min_wave_two,riog_valid_max_wave_two,   &
                          riog_valid_min_wave_wlm,riog_valid_max_wave_wlm,   &
                          riog_valid_min_wave_bhd,riog_valid_max_wave_bhd,   &
                          riog_valid_min_wave_stokes,riog_valid_max_wave_stokes

!! * Arguments

!! * Local declarations

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

   file_wave = "wave.nc"
   longitude_name_wave="longitude";latitude_name_wave="latitude"

   l_wave_marsgrid=.TRUE.

   l_wave_rdir=.FALSE.;l_wave_rfoc=.FALSE.;l_wave_rfbb=.FALSE.;l_wave_rubot=.FALSE.;l_wave_ruubot=.FALSE.;l_wave_rvubot=.FALSE.
   l_wave_rtp=.FALSE.;l_wave_rhs=.FALSE.;l_wave_rt0m1=.FALSE.;l_wave_rt0m1=.FALSE.;l_wave_rtauwx=.FALSE.;l_wave_rtauwy=.FALSE.
   l_wave_rtauox=.FALSE.;l_wave_rtauoy=.FALSE.;l_wave_rtbb=.FALSE.;l_wave_rwlm=.FALSE.;l_wave_rbhd=.FALSE.
   l_wave_rtus=.FALSE.;l_wave_russ=.FALSE.

   name_wave_hs="hs";name_wave_ubot="ubr";name_wave_uubot="uubr";name_wave_vubot="vubr";name_wave_foc="foc";name_wave_t0m1="t0m1"
   name_wave_dir="dir";name_wave_tp="t0m1";name_wave_tauwx="utaw";name_wave_tauwy="vtaw"
   name_wave_tauox="utwo";name_wave_tauoy="vtwo";name_wave_tbbx="utbb";name_wave_tbby="vtbb"
   name_wave_wlm="lm";name_wave_bhd="bhd"
   name_wave_fbb="fbb";name_wave_tbbx="utbb";name_wave_tbby="vtbb"
   name_wave_ussx="uuss";name_wave_ussy="vuss"
   name_wave_tusx="utus";name_wave_tusy="vtus"

   name_out_wave_dir="dir"; name_out_wave_hs="hs"; name_out_wave_ubot="ubot"
   name_out_wave_ubotx="uubr";name_out_wave_uboty="vubr"
   name_out_wave_foc="foc"; name_out_wave_tp="t0m2"; name_out_wave_t0m1="t0m1"
   name_out_wave_fbb="fbb";name_out_wave_tbbx="utbb";name_out_wave_tbby="vtbb"
   name_out_wave_tauox="tauox";name_out_wave_tauoy="tauoy";name_out_wave_wlm="lm"
   name_out_wave_tauwx="tauwx";name_out_wave_tauwy="tauwy"
   name_out_wave_bhd="bhd";name_out_wave_ussx="uuss";name_out_wave_ussy="vuss"
   name_out_wave_usx="usz";name_out_wave_usy="vsz";name_out_wave_usz="wsz"
   name_out_wave_tusx="utus";name_out_wave_tusy="vtus"

   riog_valid_min_wave_dir=-0.5;riog_valid_max_wave_dir=360.5
   riog_valid_min_wave_ubot=-3.0;riog_valid_max_wave_ubot=3.0
   riog_valid_min_wave_hs=-0.1;riog_valid_max_wave_hs=20.0
   riog_valid_min_wave_foc=-10.0;riog_valid_max_wave_foc=+100.0
   riog_valid_min_wave_tp=-0.1;riog_valid_max_wave_tp=50.0
   riog_valid_min_wave_t0m=-1.0;riog_valid_max_wave_t0m=100.0
   riog_valid_min_wave_two=-1.0;riog_valid_max_wave_two=1.0
   riog_valid_min_wave_wlm=-1.0;riog_valid_max_wave_wlm=30001.0
   riog_valid_min_wave_bhd=-1.0;riog_valid_max_wave_bhd=10.0
   riog_valid_min_wave_stokes=-1.0;riog_valid_max_wave_stokes=1.0

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

  END SUBROUTINE cas_init_wave


  END MODULE casinit
