   !&E---------------------------------------------------------------------

#ifdef key_palm_ww3
  SUBROUTINE mars2d3d()
#else
  PROGRAM mars2d3d
#endif

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE main  ***
   !&E
   !&E ** Purpose : main program
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : ionc4_init, init, bathy, grid_mask,grid_metric
   !&E                   : tool_decompdate, meteo, init_xe, init_uvz, init_uv, init_viscoef
   !&E                   : init_turbcoef, traceur, river, outflow, init_tssub, ts_stateeq_linear
   !&E                   : ts_stateeq_blumberg, sortie, diagnostic, track_points, step
   !&E                   : Agrif_Step
   !&E
   !&E ** Used ij-arrays : 
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History : 
   !&E        !  1987-02  (J.-C. Salomon, P. Lazure) original code
   !&E        !  2007-07  (M. Sourisseau) elimination var_read_obc
   !&E        !  2009-02  (V. Garnier) MPI2D
   !&E        !  2009-02  (P. Le Hir, F. Cayocca) introduction "new" sand/mud module
   !&E        !  2009-07  (N. Bruneau, F. Dumas) coupling with swan
   !&E        !  2009-09  (B. Thouvenin) coupling with MET&OR (contaminant)
   !&E        !  2009-10  (V. Garnier) reproductibility when using restart file
   !&E        !  2010-03  (A.-C. Bennis, T. Odaka) use of palm coupling
   !&E        !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
   !&E        !  2010-06  (F. Cayocca) coupling with the module of sedimentology
   !&E        !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
   !&E        !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E        !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
   !&E        !  2011-03  (M. Huret) Introduction of Individual Based Model
   !&E        !  2013-07 (A. Thevenin - CERFACS) coupling using OASIS3-MCT
   !&E        !  2014-07  (F. Ganthy) Introduction of obstruction module
   !&E        !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
   !&E---------------------------------------------------------------------

#include "toolcpp.h"
#ifdef key_sedim_MUSTANG
#include "../../SEDIM/coupleur_define_MUSTANG.h"
#endif
   !! * Modules used
   USE parameters
   USE ionc4, ONLY : ionc4_init,ionc4_closeall
   USE comvars2d
   USE comvars3d
   USE comvarp2d
   USE comvarp3d
   USE comcolig2d,    ONLY : igdv_old
   USE comcolig3d,    ONLY : intz_exp_vstar
   USE comdiag,       ONLY : l_diag,file_diag
   USE output,        ONLY : output_init,output_mng,l_out_demer,fileresu,nb_out
   USE comdem
   USE comfluxdo
   USE comturb,       ONLY : eps
   USE obcogcm,       ONLY : obc_tsref
#ifdef key_waveforcing
   USE waveforcing,   ONLY : init_wave_forcing_fric,compute_wave_forcing
#endif
#if defined key_wave && ! defined key_wave_monochromatic
!   USE wave_tools,    ONLY : wave_stokes_fill_tabkm
#endif
#ifdef key_substance
   USE subreaddat,    ONLY : sub_read_var,sub_read_vardiag
   USE subreadflxatm, ONLY : sub_sflx_read
   USE comsubstance,  ONLY : nv_dri, nv_state, cv_wat,l_cvrain_readfile,l_subflxatm_readfile,l_bilan
   USE usersubstance, ONLY : user_init_vdriv
   USE obcogcm,       ONLY : obc_subsref
   USE submassbalance,ONLY : sub_budget_readdomain
#ifdef key_biolo
   USE bioloinit,     ONLY : biolo_userinit
#endif
#ifdef key_eco3m
   USE DEF_TYPE
   USE VAR_USER
   USE VAR_GLOBAL
   USE COUPLEUR_PHY_BIO
   USE F_PROCESS
#endif
#ifdef key_contaminant                             
   USE contareaddat,  ONLY : conta_read_react  
#endif 
#endif
#ifdef key_sedim_MUSTANG
   USE comMUSTANG,        ONLY : MUSTANG_alloc,cv_sed,l_repsed,l_morphocoupl,coef_erolat
   USE initMUSTANG,       ONLY : MUSTANG_compatibility,MUSTANG_param,MUSTANG_sedinit
   USE sed_MUSTANG_MARS,  ONLY : sedinit_fromfile
   USE sedimorpho,        ONLY : sed_morphoinit,sed_alloc_morpho
   USE coupleur_MUSTANG,  ONLY : coupl_conv2MUSTANG
#endif
   USE sflxatm,     ONLY : drag,exch,imeteo_dragtype,imeteo_exchtype
   USE sflxturb,    ONLY : sflx_coef_init
   USE trackpoints, ONLY : l_points,track_points
#ifdef key_trajec3d
   USE trajinitsave,   ONLY : traj_init3d
#endif
#ifdef key_trajec2d
   USE traject2d,   ONLY : traj_init2d
#endif
#ifdef key_ibm
   USE ibm,         ONLY : ibm_init
#endif
#ifdef key_diagdyn
   USE diagdyn,     ONLY : diagdyn_init
#endif
#ifdef key_wave_qeuler_diag
   USE wavediag,    ONLY : wave_diag_qeuler_init
#endif
#ifdef key_nudging
   USE nudging,     ONLY : nud_init,l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean
#endif
   USE casinit,     ONLY : cas_init
   USE_MPI toolmpi, ONLY : MPI_START1,MPI_START2,MPI_START3,MPI_END1,NPROCS
   USE_MPI toolmpi, ONLY : ADD_TAG,ex_j_rsh,ex_i_rsh,MASTER
#ifdef key_oasis
   USE tooloasismct,ONLY : cpl_oasis_grid, cpl_oasis_define
   USE_MPI toolmpi, ONLY : MPI_COMM_MARS
#endif
   USE tooltiming,  ONLY : set_iluout_timing, time_stat_ll, time_header_ll, second_mars, &
                           timing_legend, timing_separator, xt_start, xt_diag, &
                           xt_store, xt_bound, xt_forcing, xt_nudging, xt_traj, &
                           xt_dynamics, xt_turb, xt_tsbuoy, xt_subs, xt_points, &
                           xt_misc, xt_oasis, xt_agrif
   USE_MPI toolmpi, ONLY : xt_mpi_com,xt_mpi_wait
#ifdef key_obstructions
   USE obstructions,    ONLY : obst_init,obst_dealloc
#endif

   !! * Arguments

   !! * Local declarations

   IMPLICIT NONE

#ifdef key_VT
   INCLUDE '/usr/local/vampir/vampirtraceold/include/VT.inc'
   INTEGER ierr
#endif

   INTEGER                            :: taille_germe,narg,iargc,nb
   INTEGER, ALLOCATABLE, DIMENSION(:) :: germe
   INTEGER                            :: iluout      ! Logical unit number for the output listing
#ifdef key_agrif
   EXTERNAL :: step, prepare_grids_connection, make_grids_connection
#endif
#ifdef key_MPI_2D
   CHARACTER(LEN=3)     :: tag
#endif
                                ! for computing time analysis
   REAL(KIND=rlg),DIMENSION(2)        :: ztime0,ztime1,ztime2,ztot,zall,xt_init,xt_step

!$   INTEGER  :: OMP_GET_NUM_THREADS

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

#ifdef key_VT
   CALL VTTRACEON(ierr)
#endif

   INST_AGRIF Agrif_Parallel_Sisters = .FALSE.

   CALL_MPI MPI_START1

   ! 0.0 Time logger initialization
   ! ------------------------------
   CALL second_mars(ztime0)
   xt_start     = 0.0_rlg ; xt_store     = 0.0_rlg ; xt_bound     = 0.0_rlg ;
   xt_forcing   = 0.0_rlg ; xt_nudging   = 0.0_rlg ; xt_dynamics  = 0.0_rlg ;
   xt_turb      = 0.0_rlg ; xt_tsbuoy    = 0.0_rlg ; xt_subs      = 0.0_rlg ;
   xt_diag      = 0.0_rlg ; xt_traj      = 0.0_rlg ; xt_points    = 0.0_rlg ;
   xt_misc      = 0.0_rlg ; xt_oasis     = 0.0_rlg ; xt_agrif     = 0.0_rlg ;
#ifdef key_MPI_2D
   xt_mpi_wait = 0.0_rlg
   xt_mpi_com = 0.0_rlg
#endif 


#ifdef key_agrif
   ! Initialisation de la grille mere (allocation des variables de la grille
   ! mere plus hierarchie de grille). Devrait s appeler Agrif_Init.
   ! Grilles filles n existent pas a ce stade
   CALL Agrif_Init_Grids()
#endif

   ! 1.0 Initialization
   ! ------------------

   !  Init des boucles
   !------------------------
   nbouc=0
   nbouc_macho=0

   ! initialisation de la procedure de tirage des nombres aleatoires
   !-----------------------------------------------------------------
   root='./'
   CALL RANDOM_SEED(size=taille_germe)
   ALLOCATE(germe(taille_germe))
   germe(:)=3
   CALL RANDOM_SEED(put=germe(:))

#ifdef key_oasis
   CALL ionc4_init(MPI_COMM_MARS)
#else
   CALL ionc4_init()
#endif
  
   narg=iargc()

   IF (l_testcase) THEN

   ! Initialize test case configuration
     CALL cas_init

   ! Initialize coefficients of turbulence
     IF (turb_nbeq==1 .OR. turb_nbeq==2) CALL init_turbcoef(ect,psi)

   ! Initialize xeadv
     xeadv(limin:limax,ljmin:ljmax)=ssh(limin:limax,ljmin:ljmax)

   ! Restart from a file ( restart file from MARS or climatological file )
     IF (l_initfromfile) THEN
       CALL init_fromfile(u,uz,v,vz,ssh,sal,temp,        &
                               intz_exp_vstar,igdv_old,  &
#ifdef key_substance
                                            cv_wat,      &
#endif
                                            ect,psi,eps)
     END IF

   ELSE

   ! Read namelists (paramain.txt, paracom.txt, paraspec.txt)
   ! Read head.conf
     CALL init(narg)

   ! Read parameters of XXXXX
#ifdef key_waveforcing
!     CALL lec_para_swan
#endif

   ! Initialize flux exchange coefficients (to prepare the bilinear interpolation)
     IF (imeteo_exchtype==1) CALL sflx_coef_init(drag,exch,imeteo_dragtype)

   ! For MARS_WEB
     IF (narg>=2) CALL init_marsweb_clim
     IF (narg==3) CALL init_marsweb_simu

   ! Read the bathymetry
     CALL bathy

   ! Initialize ig,id,jb,jh
     CALL grid_mask

   ! Divide the domain for MPI run
     CALL_MPI MPI_START2(1,1)
     CALL_MPI MPI_START3
     ! Toujours sur grille mere uniquement
     CALL_AGRIF Agrif_Declare_Variables


   ! Allocation of variables 
   ! -----------------------
     ! Dynamical variables
     CALL alloc
     ! Intermediate variables
     CALL alloc_inter
     ! Stationnary variables (grid metrics)
     CALL alloc_stat

   ! Initialize metrics (dx,dy,sig,lon,lat...)
     CALL grid_metric


#ifdef key_substance
   ! Read file variables.dat describing substances
   ! Allocate substance variables 
   ! Make variables in order
     CALL sub_read_var

   ! Initialize open boundary arrays for substance
     CALL obc_subsref

#ifdef key_sedim_MUSTANG

#if ! defined key_MARS
     ! coupleur definir des variables ou parametres dont a a besoin qui sont definis dans MARS ??
     CALL coupl_MUSTANG_hydro

     ! coupleur pour donner les nombres de variables de chaque type
     CALL coupl_MUSTANG_substance(0)

     ! initialisation des parametres pour la sedimento (appelle dans MARS par la subroutine init.F90)
     CALL MUSTANG_param('r')
#endif

     CALL MUSTANG_alloc

#if ! defined key_MARS
     ! LIRE les infos sur les variables a transporter
     ! pour l instant il faut que les variables soient bien rangees dans le bon ordre
     ! (a la lecture et numerotation des variables)
     ! ordre des variables iv : graviers, sable, vase, particulaires non constitutives, dissoutes,
     !                          fixees,forcantes,intermediaires

     ! coupleur pour definir et allouer certains tableaux apres qu'ils aient ete alloues 
     ! (remplace module substance)
      CALL coupl_MUSTANG_substance(1)
#endif
#endif

   ! Read file vardiag.dat describing diagnoses
   ! Allocate diagnoses variables
#ifdef key_biolo
     CALL sub_read_vardiag
#endif
#ifdef key_eco3m
! Initialization of ECO3M
     IF (nbouc == 0) THEN
       nbcallbio = -1
       CALL main_bio
       CALL alloc_var_eco3m
       !In order to compute the modulo of the call to ECO3M
       dt_phy_ini= dt
     END IF
#endif
#endif

   ! Initialize sea surface height
     CALL init_xe(ssh)
   ! Initialize height variables for vertical integration
     CALL init_height(ssh)
   ! Initialize 3D velocities
     CALL init_uvz(uz,vz)
   ! Initialize 2D velocities
     CALL init_uv(u,v,uz,vz)
   ! Initialize viscosity coefficients
     CALL init_viscoef
   ! Initialize coefficients of turbulence
     CALL init_turbcoef(ect,psi)
   ! Initialize sal, temp, cv_wat
     CALL init_tssub(sal,temp   &
#ifdef key_substance
                        ,cv_wat &
#endif
                             )
     CALL obc_tsref
   ! Initialize obstructions
#ifdef key_obstructions
     CALL obst_init
#endif

   ! Restart from a file ( restart file from MARS or climatological file )
     IF (l_initfromfile) THEN
       CALL init_fromfile(u,uz,v,vz,ssh,sal,temp,        &
                               intz_exp_vstar,igdv_old,  &
#ifdef key_substance
                                            cv_wat,      &
#endif
                                            ect,psi,eps)
     END IF


#ifdef key_sedim_MUSTANG
    ! coupleur pour recuperation des concentrations au fond 
    ! pour iappel=0 seuls WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD vont etre utilises
    CALL coupl_conv2MUSTANG(IMIN_BOUCL,IMAX_BOUCL,JMIN_BOUCL,JMAX_BOUCL,0,  &
                            BATHY_h0,WATER_ELEVATION,   &
#if ! defined key_MARS
                 THICK_LAYER_WAT,VERT_VISC,WAT_SETTL,SETTL_FLUXSUM_w2s, &
#endif
                       WATER_CONCENTRATION,SALINITY_MOD,TEMPERATURE_MOD )


     IF (l_repsed) THEN
#ifdef key_MARS
       CALL sedinit_fromfile(IMIN_BOUCL,IMAX_BOUCL,JMIN_BOUCL,JMAX_BOUCL,  &
                                 RHOREF,BATHY_h0)
#else
      !! a modifier selon modele hydro
      !! CALL sed_init_fromfile(RHOREF,BATHY_h0,TEMPERATURE,SALINITY,WATER_CONCENTRATION)
#endif
     ELSE
       CALL MUSTANG_sedinit(IMIN_BOUCL,IMAX_BOUCL,JMIN_BOUCL,JMAX_BOUCL,  &
                         RHOREF,BATHY_h0)
     END IF
     !CALL sed_settlveloc_ini
     IF(l_morphocoupl) THEN
        CALL sed_alloc_morpho
        IF(.NOT.l_initfromfile)CALL sed_morphoinit
     ENDIF
#endif
   ! Initialize wave fields
#ifdef key_waveforcing
     CALL init_wave_forcing_fric
!     CALL lect_wave
!cval ?     CALL bathy
     CALL compute_wave_forcing
#endif

   ! Initialize ecosystem variables
#ifdef key_substance
#ifdef key_biolo
   ! User special initialization of cv_wat.
     CALL biolo_userinit(cv_wat)
#endif
!!!! CAR INITIALISE DANS STEP POUR SEDIM REVOIR ?#if defined key_sedim_MUSTANG_anc || defined key_contaminant
#if defined key_contaminant
     CALL conta_read_react
!!!!!     CALL conta_userinit
#endif
#endif

#ifdef key_substance
   ! initialization of driving variables
     IF (nv_dri > 0) &
       CALL user_init_vdriv(cv_wat                       &
#ifdef key_sedim_MUSTANG                                
                           ,cv_sed                       &
#endif   
                                )

#endif
   ! Read output.dat and allocate variables with the correct number of output files
     CALL output_init

   END IF   ! ends test on l_testcase

#if defined key_chainop
   ! Read output.dat and allocate variables with the correct number of output files
     CALL output_init
#endif

#if defined key_oasis
   CALL cpl_oasis_grid(MASTER, MPI_COMM_MARS)
   CALL cpl_oasis_define
#endif

#ifdef key_nudging
   CALL nud_init(l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean)
#endif

#ifdef key_substance
   !!! introduction of submassbalance
     IF (l_bilan) CALL sub_budget_readdomain(ssh)
     PRINT_DBG*, 'call read_limits'
   !!! end submassbalance
#endif

   CALL_MPI ex_j_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,v(liminm2:limaxp2,ljminm2:ljmaxp2))   ! for dyn2d_dt val
   CALL_MPI ex_i_rsh(-2,2,1,liminm2,limaxp2,ljminm2,ljmaxp2,u(liminm2:limaxp2,ljminm2:ljmaxp2))   ! for dyn2d_dt val

   ! Save summary of simulation
   IF_MPI (MASTER) THEN
   CALL summary
#ifdef key_sedim_MUSTANG
   IF(.NOT. l_testcase)CALL MUSTANG_compatibility
#endif
   ENDIF_MPI

   ! Initialize listing file
   IF (narg == 0 .AND. .NOT. l_testcase) THEN
     filename= 'listing'
     CALL_MPI ADD_TAG(filename,LEN(filename))
     CALL_MPI ADD_TAG(file_diag,LEN(file_diag))
!     CALL edge_initialise_MPI
     IF_MPI (MASTER) THEN
     IF (iscreen/=6) OPEN(unit=iscreen,file=filename,status="replace")
     ENDIF_MPI
   END IF

   ! Initialize meteorological variables if l_meteo_stat or l_meteo_hom
   IF (.NOT. l_testcase) CALL meteo

   ! Read file river.dat, initialize river discharges inputs with time interpolation
   CALL river(ssh)

#ifdef key_substance
   ! Read file outflow.dat, file describing discharges
   ! initialize discharges inputs with time interpolation
   CALL outflow

   ! Read files sflx_cvrain and/or sflx_subflxatm 
   ! if temporal variations of cvrain or of atm. fluxes for at least one substance
   !  (l_cvrain_readfile=.true. or l_subflxatm_readfile=.true.)
   IF (l_cvrain_readfile .or. l_subflxatm_readfile) CALL sub_sflx_read
#endif


   ! Estimate buoyancy from sal and temp
   IF (l_modele2d) THEN
     IF (icon == 2) CALL ts_stateeq2d(sal)
   ELSE
     IF (icon >= 1 .AND. l_stateeq_lin) CALL ts_stateeq_linear (sal,temp,boz,pn2)
     IF (icon >= 1 .AND. .NOT.l_stateeq_lin) CALL ts_stateeq_blumberg(sal,temp,boz,pn2)
   END IF

   CALL second_mars(ztime2)
   xt_start = xt_start + ztime2 - ztime0
   ztime1 = ztime2

   ! Intialize variables for traking
   IF (l_points) CALL track_points(u,v,uz,vz,ssh,sal,temp, &
#ifdef key_substance
                                   cv_wat,                 &
#endif
                                   ect,boz)

   CALL second_mars(ztime2)
   xt_points = xt_points + ztime2 - ztime1
   ztime1 = ztime2

   ! Initialize diagnoses variables
   IF (l_diag) CALL diagnostic(u,v,uz,vz,ssh,sal,temp, &
#ifdef key_substance
                               cv_wat,                 &
#endif
                               ect,boz)

#ifdef key_wave_qeuler_diag
   CALL wave_diag_qeuler_init
#endif

#if defined key_wave && ! defined key_wave_monochromatic
!   CALL wave_stokes_fill_tabkm
#endif

   CALL second_mars(ztime2)
   xt_diag = xt_diag + ztime2 - ztime1
   ztime1 = ztime2

   ! manage outputs files
   CALL output_mng

   ! Initialize demerliac filter
   DO nb=1,nb_out

     IF (l_out_demer(nb)) &
        CALL out_demerliac(ssh(limin:limax,ljmin:ljmax),                                     &
                           u(limin:limax,ljmin:ljmax),uz(:,limin:limax,ljmin:ljmax),         &
                           v(limin:limax,ljmin:ljmax),vz(:,limin:limax,ljmin:ljmax),         &
                           sal(:,limin:limax,ljmin:ljmax),temp(:,limin:limax,ljmin:ljmax),   &
                           xed(limin:limax,ljmin:ljmax),                                     &
                           ud(limin:limax,ljmin:ljmax),uzd(:,limin:limax,ljmin:ljmax),       &
                           vd(limin:limax,ljmin:ljmax),vzd(:,limin:limax,ljmin:ljmax),       &
                           sald(:,limin:limax,ljmin:ljmax),tempd(:,limin:limax,ljmin:ljmax), &
#ifdef key_substance
                           cv_wat(:,:,limin:limax,ljmin:ljmax),cv_watd(:,:,limin:limax,ljmin:ljmax), &
#endif
                           nsdem)
   END DO

   CALL second_mars(ztime2)
   xt_store = xt_store + ztime2 - ztime1
   ztime1 = ztime2

   ! Initialize diagnoses variables relative to terms of the momentum equation
   ! after output_mng to be allowed to complete the main output file ( fileresu(1) )
#ifdef key_diagdyn
   CALL diagdyn_init
   xt_diag = xt_diag + ztime2 - ztime1
   ztime1 = ztime2
#endif

#ifdef key_trajec3d
   IF (.NOT.l_modele2d) CALL traj_init3d(ssh)
#endif
#ifdef key_trajec2d
   IF (l_modele2d) CALL traj_init2d
#endif

#ifdef key_ibm
   IF (.NOT.l_modele2d) CALL ibm_init(ssh,sal,temp, &
#ifdef key_substance
                                      cv_wat,       &
#endif
                                      boz)
#endif

#if defined key_trajec3d || defined key_trajec2d || defined key_ibm
   CALL second_mars(ztime2)
   xt_traj = xt_traj + ztime2 - ztime1
   ztime1 = ztime2
#endif

   xt_init = xt_start + xt_points + xt_diag + xt_store + xt_traj


   ! Save the repartition of the cpus
   ! --------------------------------

#ifdef key_MPI_2D
   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' NUMBER OF MPI CPUS' , nprocs
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   ENDIF_MPI
#endif
   iscreenlog=iscreenlog   ! because of bug in AGRIF conv
   IF_MPI (MASTER) THEN
!$OMP PARALLEL
!$OMP MASTER
!$   WRITE(iscreenlog,*) ' '
!$   WRITE(iscreenlog,*) '**************************************************'
!$   WRITE(iscreenlog,*) ' NUMBER OF OMP THREADS' , OMP_GET_NUM_THREADS()
!$   WRITE(iscreenlog,*) '**************************************************'
!$   WRITE(iscreenlog,*) ' '
!$OMP END MASTER
!$OMP END PARALLEL
   ENDIF_MPI

   ! 2.0 Time integration
   ! --------------------

   IF_MPI (MASTER) THEN
   WRITE(iscreen,*) ' '
   WRITE(iscreen,*) '**************************************************'
   WRITE(iscreen,*) '**************** Start time integration **********'
   WRITE(iscreen,*) '**************************************************'
   WRITE(iscreen,*) ' '
   ENDIF_MPI

   halfdt = dt*0.5_rlg

   ! Time integration
#ifdef key_agrif
   DO WHILE (t < tdebagrif .OR. MOD(nbouc,2) ==1 )  
     CALL step()
   END DO
   DO WHILE (t <= tfin .OR. MOD(nbouc,2)==0)
     ! call agrif_regrid au premier appel, donc Agrif_InitValues
     ! avant integration temporelle de la grille mere.
     CALL Agrif_Step(step)
   END DO
#else
   DO WHILE (t <= tfin .OR. MOD(nbouc,2)==0)
#ifdef key_waveforcing
     CALL compute_wave_forcing
#endif
#ifdef key_eco3m
! Call of ECO3M (in order to compute the dcdt term
! with a variable time step, modulo of dt ini of MARS
   IF (mod(nint(dt_bio/(2*dt_phy_ini)),2)==0) THEN
     IF (mod(nbouc,nint(dt_bio/(2*dt_phy_ini)))==0) THEN
       CALL transfert_marseco(ssh)
       CALL main_bio
       CALL transfert_ecomars
     END IF
   ELSE
     IF (mod(nbouc,1+nint(dt_bio/(2*dt_phy_ini)))==0) THEN
       CALL transfert_marseco(ssh)
       CALL main_bio
       CALL transfert_ecomars
     END IF
   END IF
#endif

     CALL step()
   END DO
#endif

   CALL second_mars(xt_step)
   !xt_step = xt_step - (xt_start + xt_diag + xt_store + xt_traj + xt_points)
   xt_step = xt_step - ztime1


   ! 3.0 Print statistics
   ! --------------------

   ! Set File Timing OUTPUT
   iluout = 1970
   CALL set_iluout_timing(iluout)

   ! Compute global time
   CALL time_stat_ll(xt_init+xt_step,ztot)  ! reference 100%

   CALL time_header_ll

   CALL time_stat_ll(xt_init,ztot,      'MAIN INITIALIZATION','=')
   CALL time_stat_ll(xt_step,ztot,       'MAIN TIME_STEP','=')
   CALL timing_separator('=')
   CALL timing_separator(' ')

   CALL timing_legend()
   CALL time_stat_ll(xt_start,ztot,      ' INITIALIZATION','=')
   CALL time_stat_ll(xt_bound,ztot,      ' BOUNDARIES','=')
   CALL time_stat_ll(xt_store,ztot,      ' STORE-FIELDS','=')
   CALL time_stat_ll(xt_dynamics,ztot,   ' DYNAMICS (mtm cty)','=')
   CALL time_stat_ll(xt_forcing,ztot,    ' FORCING','=')
   CALL time_stat_ll(xt_tsbuoy,ztot,     ' T_S_BUOY','=')
   CALL time_stat_ll(xt_turb,ztot,       ' TURBULENCE','=')
   CALL time_stat_ll(xt_nudging,ztot,    ' NUDGING','=')
   CALL time_stat_ll(xt_points,ztot,     ' tracking POINTS','=')
   CALL time_stat_ll(xt_diag,ztot,       ' DIAG','=')
   CALL time_stat_ll(xt_traj,ztot,       ' TRAJ','=')
   CALL time_stat_ll(xt_misc,ztot,       ' MISC','=')
   CALL time_stat_ll(xt_oasis,ztot,      ' OASIS','=')
   CALL time_stat_ll(xt_agrif,ztot,      ' AGRIF','=')
   CALL time_stat_ll(xt_subs,ztot,       ' SUBSTANCES','=')
#ifdef key_MPI_2D
   CALL time_stat_ll(xt_mpi_wait,ztot,   ' MPI_WAIT','=')
   CALL time_stat_ll(xt_mpi_com,ztot,    ' MPI_COM','=')
# endif
   CALL timing_separator('=')
   CALL timing_separator(' ')

   ! sum of call subroutine
   zall = xt_start + xt_bound  + xt_store  + xt_dynamics +  xt_forcing + &
          xt_tsbuoy+ xt_turb + xt_nudging + xt_points + xt_diag + xt_traj + &
          xt_misc + xt_oasis + xt_agrif + xt_subs
  
   CALL timing_legend()
   CALL time_stat_ll(zall,ztot,          ' SUM(CALL)','=')
   CALL second_mars(zall)
   CALL time_stat_ll(zall-ztime0,ztot,   ' TOTAL_RUN','=')
   CALL timing_separator('=')

   ! 4.0 Deallocations
   ! -----------------

   CALL dealloc_inter
   CALL dealloc_stat
   CALL ionc4_closeall()

   PRINT*, 'fin normale'

   CALL dealloc

#ifdef key_obstructions
   CALL obst_dealloc
#endif

#ifdef key_eco3m
   CALL dealloc_var_eco3m
#endif

   CALL_MPI MPI_END1

#ifdef key_VT
   CALL VTTRACEOFF(ierr)
#endif

#ifdef key_palm_ww3
  END subroutine mars2d3d
#else
  END PROGRAM mars2d3d 
#endif
