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

  SUBROUTINE step

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE step  ***
   !&E
   !&E ** Purpose : manage the estimate of all pronostic variables
   !&E
   !&E ** Description : - save restart file
   !&E                  - update forcings
   !&E                  - update the time step
   !&E                  - call the step?dxy or step?dyx subroutines, depending on if the time-step
   !&E                    is integer or half-integer (i.e. depending on the value of nbouc)
   !&E                  - call diagnoses
   !&E                  - estimate maximum velocityies at the sea surface
   !&E
   !&E ** Called by : main
   !&E
   !&E ** External calls : step3dxy, step3dyx, tool_sectodat, tide_genpot, meteo, sflx_read, sub_sflx_read
   !&E                   : out_saverestart, outflow, river, tool_decompdate
   !&E                   : tool_julien, sflx_surf, step2dxy, step2dyx, diagnostic, track_points, output_mng
   !&E                   : tide_astroarg
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables : cdate, ijour, imois, ian, iheure, iminu, isec, jjulien
   !&E                         l_boolsauv, l_boolsort
   !&E                         indchrono, nbouc, t
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E        !           (P. Lazure)  Original code
   !&E        !  2005-02-23
   !&E        !  2007-07-25  (M. Sourisseau) elimination var_read_obc
   !&E        !  2009-03  (V. Garnier) MPI2D
   !&E        !  2009-02  (P. Le Hir, F. Cayocca) introduction "new" sand/mud module
   !&E        !  2009-09  (B. Thouvenin) introduction of module MET&OR (contaminant)
   !&E        !  2009-10  (F. Cayocca) correction for sedimentology
   !&E        !  2009-10  (V. Garnier) reproductibility when using restart file
   !&E        !  2010-02  (G. Simon) introduction of module wave
   !&E        !  2010-03  (A.-C. Bennis) coupling between ww3 and mars through palm
   !&E        !  2010-06  (F. Cayocca) coupling with the module of sedimentology
   !&E        !  2010-07  (R. Verney) Introduction of orbital velocity variable
   !&E        !  2010-10   (B. Thouvenin) Distribution of the river discharge throughout a thickness
   !&E        !  2013-07 (A. Thevenin - CERFACS) MARS-WRF coupling with OASIS3-MCT
   !&E        !  2013-12 (A. Thevenin - CERFACS) MARS-WW3 coupling with OASIS3-MCT
   !&E        !  2014-03  (A.-C. Bennis, F. Dumas, V. Garnier) Review of the wave-current quasi-eulerian formalism
   !&E        !  2014    (B. Thouvenin) : Modification for key_biolo and key_sedim
   !&E        !  2014-05 (R. Verney) : integration of module flocmod (flocculation)
   !&E        !  2014-01  (F. Ganthy) Introduction of obstruction module
   !&E        !  2015-08  (M. Honnorat) Correct ending of the run with AGRIF
   !&E        !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"
#ifdef key_sedim_MUSTANG
#include "../../SEDIM/coupleur_define_MUSTANG.h"
#endif

   !! * Modules used
   USE parameters
   USE comvars2d
   USE comvars3d
   USE comvarp2d
   USE comvarp3d
   USE comcolig2d,   ONLY : igdv_old
   USE comcolig3d,   ONLY : intz_exp_vstar
   USE cometdyn
   USE cometemp,     ONLY : l_sflx_solarcst
   USE comfluxdo
   USE comdiag,      ONLY : l_diag
   USE output,       ONLY : output_mng
   USE comaree
   USE sflxread             ! read sea surface data or atmospherical fluxes
   USE trackpoints,  ONLY : l_points,track_points
   USE tssubeqtest,  ONLY : tssub_eq_test
   USE_MPI toolmpi,  ONLY : findmax_mpi_3,MASTER,MAX_ALL_MPI_RSH
   USE tooltiming,   ONLY : second_mars,xt_store,xt_forcing,xt_nudging,xt_misc,xt_diag,xt_subs,xt_points
#ifdef key_palm_ww3
   USE palmlib,      ONLY : PL_LNAME,PL_NO_TAG
#endif
#if defined key_substance
   USE comsubstance,    ONLY : cv_wat
   USE comsubstance,    ONLY : nv_dri,l_driv_cst,t_driv,nv_state,nv_int,nv_adv,sv_part, &
                               l_cvrain_readfile,l_subflxatm_readfile,tdebsubs,l_bilan
   USE usersubstance,   ONLY : user_vintermediate,user_ftime_vdriv
   USE subreadflxatm,   ONLY : sub_sflx_read
   USE submassbalance,  ONLY : t_out_budget,sub_budget_out
#if defined key_sedim_MUSTANG
   USE comMUSTANG,      ONLY : cv_sed,phieau_s2w,flx_s2w,flx_w2s,l_morphocoupl
   USE comMUSTANG,      ONLY : hsed,dzs,flx_w2s_sum,ksma
   USE comsubstance,    ONLY : isand2
   USE sed_MUSTANG,     ONLY : sed_MUSTANG_first,sed_MUSTANG_deposition
   USE sedimorpho,      ONLY : sed_morpho
   USE comturb,         ONLY : eps
#ifdef key_MARS
   USE sed_MUSTANG_MARS,    ONLY : sed_outsaverestart
#endif
#if defined key_flocmod
   USE flocmod
#endif
#if defined key_biolo && ! defined key_sedim_partial
   USE comsubstance,   ONLY : dtbio
   USE biolosinksource,ONLY : biolo_sksc_sed
#endif
#if defined key_1DV_force_towfile || defined key_1DV_force_ufile
   USE readforcingk1DV,  ONLY :  cas1DV_read_forcing
#endif   
#endif
#if defined key_parsub_newdt && defined key_substance
   USE comsubstance,    ONLY : dhmin,ndt_part,nb_ws_max,iv_ws_max,ws3max,itypv,l_subs2D
#if defined key_biolo
   USE bioloinit,       ONLY : ws3max_bio
#endif
#endif

#if defined key_contaminant 
   USE contareaddat,    ONLY : nreacmax,l_treat_re_after
   USE contasinksource, ONLY : conta_reac_after
   USE comsubstance,    ONLY : dcdt,tdebsubs
#endif
#endif
#ifdef key_wave
   USE wave,            ONLY : wave_read
#endif
#ifdef key_wave_qeuler_dynamics
#ifdef key_wave_qeuler_diag
   USE wavediag,        ONLY : wave_diag_qeuler_output
#endif
#ifdef key_wave_readww3txt
   USE wave, ONLY : wave_read_ww3txt
#endif
#endif
#ifdef key_nudging
   USE nudging,         ONLY : nud_correct_fields,l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean
#endif
#ifdef key_oasis_mars_atm
   USE toolmarsatmcpl,  ONLY : snd_fields_to_atm, rcv_fields_from_atm
#endif
#ifdef key_oasis_mars_ww3
   USE toolmarsww3cpl,  ONLY : snd_fields_to_ww3, rcv_fields_from_ww3
#endif
#ifdef key_diagdyn
   USE diagdyn,    ONLY : diagdyn_output
#endif
#ifdef key_obstructions
   USE obstructions, ONLY : obst_comp_botstress,obst_update
#endif

   IMPLICIT NONE
   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=19)  :: tool_sectodat
   INTEGER            :: i,j,k,tool_julien,iumx,jumx,ivmx,jvmx,iv
   INTEGER            :: liumx,ljumx,livmx,ljvmx
   REAL(KIND=rsh)     :: umx,vmx,lumx,lvmx
#if defined key_substance && defined key_parsub_newdt
   INTEGER            :: num_ws,ndt_part_anc
#ifdef key_biolo
  REAL(KIND=rsh)      :: ww
#endif
#endif
#ifdef key_palm_ww3
   INTEGER :: il_err,ifile,timea=0,timeb=0
   INTEGER :: timea_nbouc,timeb_nbouc,timea_start,timeb_start
   INTEGER :: min_palm_time_a=1,min_palm_time_b=1,max_palm_time_a=30000,max_palm_time_b=30000
   CHARACTER(LEN=PL_LNAME) :: cl_object, cl_space
   REAL(KIND=rsh),DIMENSION(limin:limax,ljmin:ljmax) :: tauoxym
   REAL(KIND=rsh),DIMENSION(limin:limax,ljmin:ljmax) :: usend,usendxe
   REAL(KIND=rsh),DIMENSION(1:kmax) :: usendsig
#endif
#ifdef key_oasis
   INTEGER :: il_oasis_time 
#endif
   REAL(KIND=rlg),DIMENSION(2)         :: ztime1,ztime2


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

   PRINT_DBG*, 'enter STEP routine'
   !mh because t might be equal to tfin + epsilon
   IF_AGRIF ( t >= tfin .AND. .NOT.Agrif_Root() ) return

!MD 27/06/18
!   cv_sed=max(0.0_rsh,cv_sed)
!   cv_wat=max(0.0_rsh,cv_wat)
!END MD 27/06/2017

   ! increment the time
   ! ------------------

   t=t+0.5*dt
   cdate = tool_sectodat(t)
   CALL tool_decompdate(cdate,ijour,imois,ian,iheure,iminu,isec)
   jjulien=tool_julien(ijour,imois,ian)-tool_julien(1,1,ian)+1
!   IF (MOD(nbouc,10)==0 .AND. iscreen/=6) WRITE(*,*) cdate, ' dt=', dt
!   IF (iscreen==6)        WRITE(iscreen,*) cdate, ' dt=', dt

   ! increment nbouc
   nbouc=nbouc+1

   CALL second_mars(ztime1)

   ! update forcings
   ! ---------------

   ! for time interpolation of river discharges...
   IF (l_chrono) THEN
     DO WHILE((t-tchrono(indchrono))*(t-tchrono(indchrono+1)) > 0.0_rlg)
       indchrono=indchrono+1
     END DO
   END IF

   ! Potential generator at t
   IF (l_genpot) CALL tide_astroarg
   IF (l_genpot) CALL tide_genpot


#ifdef key_oasis_mars_atm
   il_oasis_time = int(t-tdeb-dt)
   IF (il_oasis_time .GE. 0) CALL rcv_fields_from_atm(il_oasis_time)
#else
   ! Fluxes at t
   IF (.NOT. l_testcase) CALL meteo        ! only if l_meteo_stat or l_meteo_hom
   IF (.NOT.l_meteo_stat .AND. .NOT.l_meteo_hom) CALL sflx_read
   IF (.NOT.l_sflx_rflux) CALL sflx_flux(temp)
#endif


#ifdef key_wave
   ! Waves at t

#ifdef key_wave_readww3txt
   IF (nbouc == 1) THEN
     ! lecture des champs de forcage venant de WW3
     CALL wave_read_ww3txt
   END IF
#elif key_oasis_mars_ww3
   il_oasis_time = int(t-tdeb-dt)
   IF (il_oasis_time .GE. 0) CALL rcv_fields_from_ww3(il_oasis_time)
#else
   CALL wave_read
#endif

#endif

!   
#if defined key_1DV_force_towfile || defined key_1DV_force_ufile
   CALL cas1DV_read_forcing
#endif


   ! River input at t
   IF (nbfl  > 0) CALL river(ssh)

#ifdef key_obstructions
   ! Updating obstruction parameters
   CALL obst_update
#endif

#ifdef key_substance
   ! Not used anymore. OBC data must be in netcdf format
!   CALL var_read_obc              ! Prescribe obc for substances

   ! interpolate data from files for cvrain and/or sub_flxatm 
   ! if one wants temporal variations of cvrain or of atm.deposition 
   ! 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
   
   ! forcing variables
   IF ( nv_dri>0 .AND. (.NOT.l_driv_cst) ) THEN
     IF (t >= t_driv) &
       CALL user_ftime_vdriv(t,cv_wat &
#ifdef key_sedim_MUSTANG
                              ,cv_sed &
#endif
                               )
   END IF

    ! intermediate variables
   IF (nv_int > 0) THEN
     CALL user_vintermediate(t,cv_wat &
#ifdef key_sedim_MUSTANG
                              ,cv_sed &
#endif
                                                    )
   END IF
#endif

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

   ! Nudging : correct the fields
   ! ----------------------------
#ifdef key_nudging
   !IF (MOD(nbouc,10)==0) CALL nud_correct_fields
   CALL nud_correct_fields(l_nud_ssh,l_nud_sal,l_nud_temp,l_nud_boz,l_nud_u3d,l_nud_ubt,l_nud_zonalmean)
   CALL second_mars(ztime2)
   xt_nudging = xt_nudging + ztime2 - ztime1
   ztime1 = ztime2
#endif

   ! simulate the dynamics and tracers
   ! ---------------------------------

   IF (num_testcase==18 .OR. num_testcase==19 .OR. num_testcase==20 .OR. num_testcase==24) THEN

   !************************************************************************
   ! Test advection scheme without any use of the dynamics equations
   !************************************************************************

     CALL tssub_eq_test

   ELSE


!!!#ifdef key_agrif
!!!     IF (.NOT.l_modele2d .AND. t >= tdeb3d) THEN
!!!#else
     IF (MOD(nbouc,2) == 1 .AND. .NOT.l_modele2d .AND. t >= tdeb3d) THEN
!!!#endif

#if ! defined key_offline

   ! Update the time step

       CALL dyn3d_dt(uz,vz)


   ! Test the stability of the water column and
   ! apply vertical mixing on temp and sal if required

       IF (l_stability) &
         CALL dyn3d_stabts(sal,temp,boz)

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

   ! key_offline
#endif

#ifdef key_obstructions
#ifndef key_sedim_MUSTANG
       CALL obst_comp_botstress
#endif
#endif
     END IF     ! Ends loop on odd nbouc and 3d run

#ifdef key_flocmod         
       IF (t .GE. tdebsubs) THEN
         CALL flocmod_main
       ENDIF
#endif

#if defined key_parsub_newdt && defined key_substance
         !  modif pour recalculer le ndt_part pour chaque groupe et optimiser
       IF (t .GE. tdebsubs+dt) THEN

         DO num_ws=1,nb_ws_max
            iv=iv_ws_max(num_ws,1)
#ifdef key_sedim_MUSTANG            
            IF(itypv(iv) == 1) THEN
              ! gravel 
               ndt_part(num_ws)=1
            ELSE
              IF (.NOT. l_subs2D(iv)) THEN
                CALL_MPI MAX_ALL_MPI_RSH(ws3max(iv))
#endif
#ifdef key_biolo
                CALL_MPI MAX_ALL_MPI_RSH(ws3max_bio(iv))
                ww=0.5_rsh*(ws3max_bio(iv)+abs(ws3max_bio(iv)))
#ifdef key_sedim_MUSTANG
                ws3max(iv)=max(ws3max(iv),ww)
#else
                ws3max(iv)=ww
#endif
#endif
#if defined key_sedim_MUSTANG || defined key_biolo
                ndt_part_anc=ndt_part(num_ws)
                IF(ws3max(iv) .NE. 0.0_rsh)ndt_part(num_ws)=NINT(dt/MIN((dhmin/ws3max(iv)),dt))
#else
                ndt_part_anc=ndt_part(num_ws)
                IF(sv_part(iv) .NE. 0.0_rsh)ndt_part(num_ws)=NINT(dt/MIN((dhmin/sv_part(iv)),dt))
#endif
#ifdef key_MPI_2D
                IF (MASTER .AND. ndt_part_anc /=ndt_part(num_ws)) WRITE(iscreen,*) 'ndt_part(',num_ws,')= ',ndt_part(num_ws)
#else               
                IF (ndt_part_anc /=ndt_part(num_ws)) WRITE(iscreen,*) 'ndt_part(',num_ws,')= ',ndt_part(num_ws)
#endif
#ifdef key_sedim_MUSTANG
              ENDIF
            ENDIF
#endif

         ENDDO
       ENDIF
#endif

   ! Contaminant - equilibrium reactions 
   ! ----------------------------------
#if defined key_contaminant && defined key_substance
   ! after ==> before
       IF ( nreacmax > 0 .AND. l_treat_re_after .AND. t >= tdebsubs ) &
         CALL conta_reac_after(cv_wat,sal,temp,ssh)
#endif

  !   END IF     ! Ends loop on odd nbouc and 3d run


#if defined key_sedim_MUSTANG && defined key_substance
    IF (t .GE. tdebsubs) THEN
#if defined key_biolo && ! defined key_sedim_partial
         dtbio=0.5_rlg*dt
         CALL biolo_sksc_sed(ssh)
#endif

!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!$OMP PARALLEL
         ! Forcings and settling evaluations (bottom shear stress..)
         ! and remodelling (consolidation, diffusion, bioturbation)
         ! variables en arguments definies dans coupleur_define_MUSTANG.h
         CALL sed_MUSTANG_first(IMIN_BOUCL,IMAX_BOUCL,JMINP1_BOUCL,     &
              JMINP2_BOUCL,JMAX_BOUCL,WATER_CONCENTRATION,              &
              SALINITY_MOD,TEMPERATURE_MOD ,                            &
              THICK_LAYER_WAT,VERT_VISC,RESIDUAL_THICKNESS_WAT,         &
#if ! defined key_MARS
               WATER_FLUX_INPUTS,SETTL_FLUX_w2s,EROS_FLUX_s2w,WAT_SETTL,&
               SETTL_FLUX_TEMP_w2s,SETTL_FLUX_SAL_w2s,                  &
               EROS_FLUX_TEMP_s2w,EROS_FLUX_SAL_s2w, &
               CORFLUX_SAND,CORFLUY_SAND, &
#endif

              BATHY_h0,WATER_ELEVATION,RHOREF,DENSWAT_BOTTOM,           &
              CELL_SURF,SALREF_LIN,TEMPREF_LIN,TRANSPORT_TIME_STEP,     &
              CURRENT_TIME,BAROTROP_ZONAL_VELOCITY_U,BAROTROP_MERIDIONAL_VELOCITY_V)
           !IF(t>1100.)write(*,*)t,'apres intro'
              

!$OMP END PARALLEL
!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    END IF
#endif

   !************************************************************************
   ! 1er demi-pas de temps
   ! Estimate ssh,u,uz,sal,temp
   !************************************************************************

     IF (MOD(nbouc,2) == 1) THEN

   ! Estimate of the dynamics
   ! ------------------------

       IF (l_modele2d .OR. t < tdeb3d) CALL step2dxy
       IF (.NOT.l_modele2d .AND. t < tdeb3d) THEN
         DO j=ljminm2,ljmaxp2
           DO i=liminm2,limaxp2
             uz(:,i,j)=u(i,j)
             vz(:,i,j)=v(i,j)
           END DO
         END DO
       END IF
       IF (.NOT.l_modele2d .AND. t >= tdeb3d) CALL step3dxy


   ! Estimate of diagnoses
   ! ---------------------

       CALL second_mars(ztime1)

       IF (l_diag) CALL diagnostic(u,v,uz,vz,ssh,sal,temp, &
#ifdef key_substance
                                   cv_wat,                 &
#endif
                                   ect,boz)
       CALL second_mars(ztime2)
       xt_diag = xt_diag + ztime2 - ztime1
       ztime1 = ztime2


!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if defined key_sedim_MUSTANG && defined key_substance
!$OMP PARALLEL
       IF (t .GE. tdebsubs) THEN
         ! SETTL_FLUXSUM_w2s = somme des depots effectifs (si implicite) pendant le transport et si plusieurs pas fractionnaires
         CALL sed_MUSTANG_deposition(IMIN_BOUCL,IMAX_BOUCL,              &
                     JMINP1_BOUCL,JMINP2_BOUCL,JMAX_BOUCL,              &
                     WATER_ELEVATION,CELL_SURF,BATHY_h0,CELL_DX,CELL_DY,&
                     RESIDUAL_THICKNESS_WAT,WATER_CONCENTRATION,          &
#if ! defined key_MARS
                     SETTL_FLUXSUM_w2s,                                 &
#endif
                     SALINITY_MOD,TEMPERATURE_MOD)
                     
       END IF
!$OMP END PARALLEL
#endif
!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


#if defined key_substance
       IF(l_bilan .AND. t .GE. tdebsubs) THEN
          IF (t >= t_out_budget) CALL sub_budget_out(t,cv_wat  &
#if defined key_sedim_MUSTANG
                                         ,cv_sed         &
#endif
                                         ,ssh)
       END IF
#endif

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

   ! Save variables at each time step for debugging purpose
   ! ------------------------------------------------------

       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

     END IF

!************************************************************************
!  2eme demi-pas de temps
!  Estimate ssh,v,vz,sal,temp...
!************************************************************************

     IF (MOD(nbouc,2) == 0) THEN

   ! Estimate of the dynamics
   ! ------------------------

       IF (l_modele2d .OR. t < tdeb3d) CALL step2dyx
       IF (.NOT.l_modele2d .AND. t < tdeb3d) THEN
         DO j=ljminm2,ljmaxp2
           DO i=liminm2,limaxp2
             uz(:,i,j)=u(i,j)
             vz(:,i,j)=v(i,j)
           END DO
         END DO
       END IF
       IF (.NOT.l_modele2d .AND. t >= tdeb3d) CALL step3dyx


   ! Estimate of diagnoses
   ! ---------------------

       CALL second_mars(ztime1)

       IF (l_diag) CALL diagnostic(u,v,uz,vz,ssh,sal,temp, &
#ifdef key_substance
                                   cv_wat,                 &
#endif
                                   ect,boz)
       CALL second_mars(ztime2)
       xt_diag = xt_diag + ztime2 - ztime1
       ztime1 = ztime2


   ! Sediment 
   ! ---------------------
!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#if defined key_sedim_MUSTANG && defined key_substance

       IF (t .GE. tdebsubs) THEN
!$OMP PARALLEL
          CALL sed_MUSTANG_deposition(IMIN_BOUCL,IMAX_BOUCL,              &
                     JMINP1_BOUCL,JMINP2_BOUCL,JMAX_BOUCL,              &
                     WATER_ELEVATION,CELL_SURF,BATHY_h0,CELL_DX,CELL_DY,&
                     RESIDUAL_THICKNESS_WAT,WATER_CONCENTRATION,          &
#if ! defined key_MARS
                     SETTL_FLUXSUM_w2s,                                 &
#endif
                     SALINITY_MOD,TEMPERATURE_MOD)            
                     
!$OMP END PARALLEL
         IF(l_morphocoupl) CALL sed_morpho

       END IF
#endif

   ! Balance  
   ! ---------------------
#if defined key_substance
       IF(l_bilan .AND. t .GE. tdebsubs) THEN
          IF (t >= t_out_budget) CALL sub_budget_out(t,cv_wat  &
#if defined key_sedim_MUSTANG
                                         ,cv_sed      &
#endif   
                                         ,ssh)
       END IF
#endif

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

   ! Save variables at each time step for debugging purpose
   ! ------------------------------------------------------

       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


   ! Fixed variables 
   !  A VOIR si valables aussi pour BIOLO ou si deja prises en compte  ???
   !  pour l instant mis seulement pour contaminant
   ! ----------------------------------------------
#if defined key_contaminant && defined key_substance

   !   fixed variable which are not advected but which are submitted to reactions
       cv_wat(nv_adv+1:nv_state,:,limin:limax,ljmin:ljmax)=          &
                 cv_wat(nv_adv+1:nv_state,:,limin:limax,ljmin:ljmax) &
                 +dt*0.5_rlg*dcdt(nv_adv+1:nv_state,:,limin:limax,ljmin:ljmax)
!cval ne devrait on pas echanger la ?
       CALL second_mars(ztime2)
       xt_subs = xt_subs + ztime2 - ztime1
       ztime1 = ztime2
#endif

     END IF

   END IF    ! end test on test cas casAdvCst or casAdvSpot

   CALL second_mars(ztime1)

#ifdef key_oasis_mars_atm
   il_oasis_time = int(t-tdeb-dt)
   IF (il_oasis_time .GE. 0) CALL snd_fields_to_atm(il_oasis_time)
#endif

#ifdef key_oasis_mars_ww3
   il_oasis_time = int(t-tdeb-dt)
   IF (il_oasis_time .GE. 0) CALL snd_fields_to_ww3(il_oasis_time)
#endif

   ! Wave - current interaction with PALM
   ! ------------------------------------

#if defined key_wave_qeuler_dynamics && defined key_palm_ww3
   ! Get the forcing from WW3 (PALM)
#include "palmww3.h"
   ! Estimate the wave forcing to the ocean (quasi-eulerian dynamics)
   CALL wave_qeuler ! Hum, not sure the new version is compatible with palm
#endif


   ! Outputs
   ! -------

   CALL output_mng

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

#ifdef key_diagdyn
   CALL diagdyn_output
   xt_diag = xt_diag + ztime2 - ztime1
   ztime1 = ztime2
#endif

#ifdef key_wave_qeuler_diag
   CALL wave_diag_qeuler_output
   CALL second_mars(ztime2)
   xt_diag = xt_diag + ztime2 - ztime1
   ztime1 = ztime2
#endif

   ! Statistics
   ! ----------

   IF (l_maxu .AND. l_modele2d) THEN
     umx = -100.0_rsh
     vmx = -100.0_rsh
!$OMP PARALLEL default(none) &
!$OMP& PRIVATE(i,j,liumx,ljumx,livmx,ljvmx,lumx,lvmx) &
#ifdef key_agrif
!$OMP& SHARED(ishift_obc_west,ishift_obc_east,jshift_obc_south,jshift_obc_north, &
!$OMP&        jmax,imax) &
#endif
!$OMP& SHARED(umx,vmx,u,v,iumx,jumx,ivmx,jvmx,limin,limax,ljmin,ljmax,hx,hy,h0fond)
     lumx = umx
     lvmx = vmx
!$OMP do SCHEDULE(static)
     !DO j = ljmin,ljmax
     DO j=MAX0(ljmin,jshift_obc_south), MIN0(ljmax,jmax-jshift_obc_north)
       !DO i = limin,limax
       DO i=MAX0(limin,ishift_obc_west), MIN0(limax,imax-ishift_obc_east)
         IF (ABS(u(i,j)) > lumx .AND. hx(i,j) /= -valmanq -h0fond) THEN
           lumx = ABS(u(i,j))
           liumx = i
           ljumx = j
         END IF
         IF (ABS(v(i,j)) > lvmx .AND. hy(i,j) /= -valmanq -h0fond ) THEN
           lvmx = ABS(v(i,j))
           livmx = i
           ljvmx = j
         END IF
       END DO
     END DO
!$OMP END DO nowait
!$OMP CRITICAL
     IF (lumx.GT.umx) THEN
       umx = lumx
       iumx = liumx
       jumx = ljumx
     END IF
     IF (lvmx.GT.vmx) THEN
       vmx = lvmx
       ivmx = livmx
       jvmx = ljvmx
     END IF
!$OMP END CRITICAL
!$OMP END PARALLEL
     CALL_MPI FINDMAX_MPI_3(umx,iumx,jumx)
     CALL_MPI FINDMAX_MPI_3(vmx,ivmx,jvmx)
     IF_MPI (MASTER) THEN
     WRITE(iscreen,*) 'DATE : ',cdate
     WRITE(iscreen,*) 'dt=',dt
     WRITE(iscreen,*)'umax = ',umx,'en',iumx,jumx,'hx= ',hx_g(iumx,jumx)
     WRITE(iscreen,*)'vmax = ',vmx,'en',ivmx,jvmx,'hy= ',hy_g(ivmx,jvmx)
     ENDIF_MPI

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

   END IF

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

   IF (MOD(nbouc,2) == 0) THEN
     l_boolsauv=.FALSE.
     l_boolsort=.FALSE.

     IF (l_saverestart_1file.OR.l_saverestart_bydate) THEN
       CALL out_saverestart(u,uz,v,vz,       &
                            ssh(limin:limax,ljmin:ljmax),        &
                            sal,temp,intz_exp_vstar,igdv_old,    &
#ifdef key_substance
                            cv_wat,                              &
#endif
                            ect   )
#if defined key_sedim_MUSTANG && defined key_substance
       CALL sed_outsaverestart
#endif
     END IF     

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

   END IF     ! Ends loop on odd nbouc

   IF_AGRIF ( t > tfin .AND. .NOT.Agrif_Root() .AND. .NOT.Agrif_GrandMother() &
   INST_AGRIF    .AND. Agrif_NbStepint()==(Agrif_IRhot()-1) ) CALL dealloc
   IF_AGRIF ( t > tfin .AND. .NOT.Agrif_Root() .AND. .NOT.Agrif_GrandMother() &
   INST_AGRIF    .AND. Agrif_NbStepint()==(Agrif_IRhot()-1) ) CALL dealloc_inter

   PRINT_DBG*, 'exit STEP routine'

  END SUBROUTINE step
