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

  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---------------------------------------------------------------------





























































!!! definitions de dimensions, variables, parametres pour MUSTANG
!!! et couplage avec autre modele




!&E==========================================================================
!&E                   ***    coupleur_MUSTANG  ***
!&E==========================================================================

!/* Spatial MPI limits grid
!*/





!/* Spatial Grid limits definition  of loops
!*/







!/* dimensions table definition
!*/







!/* general variable hydro , bathy, time ... defined in hydro model but using by MUSTANG
!*/


























!/* number of advected and total substances
!*/
! nombre de variables transportees

! nombre de variables totales (variable transportees,+variables fixees, +variables forcantes etc..)


!/* name of variables in MUSTANG (dimensions order could changed)
!*/
! layers thickness (ordre des dimensions differentes)

! vertical eddy viscosity  (ordre des dimensions differentes)


!/* name of fluxes exchange between MUSTANG and hydro model
!*/











!/* dimensions of variables in hydro modele
!*/!! dimensions des variables 3D ou 4D dans le modele hydro !! (si pas MARS)

!!!#define DIM_MODELE_HYDRO_cw_bottom nv_tot,liminm1:limaxp2,ljminm1:ljmaxp2
!#define DIM_COORD_SIG 0:kmax+1





!/* Laterale Erosion
!*/
! Erosion laterale dans les mailles voisines (pourrait dependre architecture grille)













!/* sliding proces of fluid mud
!*/
! glissement depot : pente du fond dans les mailles voisines  (pourrait dependre architecture grille)










!! * 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 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


   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

   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

   USE sed_MUSTANG_MARS,    ONLY : sed_outsaverestart










   USE wave,            ONLY : wave_read








   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



   REAL(KIND=rlg),DIMENSION(2)         :: ztime1,ztime2


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

   !PRINT*, 'enter STEP routine'
!mh because t might be equal to tfin + epsilon
   !IF ( 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



! 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)




! Waves at t


   CALL wave_read




!



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




! 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 &

                              ,cv_sed &

                               )
   END IF

! intermediate variables
   IF (nv_int > 0) THEN
     CALL user_vintermediate(t,cv_wat &

                              ,cv_sed &

                                                    )
   END IF


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

! Nudging : correct the fields
! ----------------------------


! 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



! 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



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





! Contaminant - equilibrium reactions
! ----------------------------------


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



    IF (t .GE. tdebsubs) THEN


!!!!!!!!!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(limin,limax,ljminp1,     &
              ljminp2,ljmaxm1,cv_wat,              &
              sal,temp ,                            &
              eps,nz,h0fond,         &


              h0,ssh,rhoref,roswat_bot,           &
              surf,saliref_lin,tetaref_lin,halfdt,     &
              t,u,v)
!IF(t>1100.)write(*,*)t,'apres intro'
              

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


!************************************************************************
! 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, &

                                   cv_wat,                 &

                                   ect,boz)
       CALL second_mars(ztime2)
       xt_diag = xt_diag + ztime2 - ztime1
       ztime1 = ztime2


!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!$OMP PARALLEL
       IF (t .GE. tdebsubs) THEN
! flx_w2s_sum = somme des depots effectifs (si implicite) pendant le transport et si plusieurs pas fractionnaires
         CALL sed_MUSTANG_deposition(limin,limax,              &
                     ljminp1,ljminp2,ljmaxm1,              &
                     ssh,surf,h0,dx,dy,&
                     h0fond,cv_wat,          &

                     sal,temp)
                     
       END IF
!$OMP END PARALLEL

!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



       IF(l_bilan .AND. t .GE. tdebsubs) THEN
          IF (t >= t_out_budget) CALL sub_budget_out(t,cv_wat  &

                                         ,cv_sed         &

                                         ,ssh)
       END IF


       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, &

                                       cv_wat,                 &

                                       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, &

                                   cv_wat,                 &

                                   ect,boz)
       CALL second_mars(ztime2)
       xt_diag = xt_diag + ztime2 - ztime1
       ztime1 = ztime2


! Sediment
! ---------------------
!!!!!!!!!Module MUSTANG !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


       IF (t .GE. tdebsubs) THEN
!$OMP PARALLEL
          CALL sed_MUSTANG_deposition(limin,limax,              &
                     ljminp1,ljminp2,ljmaxm1,              &
                     ssh,surf,h0,dx,dy,&
                     h0fond,cv_wat,          &

                     sal,temp)            
                     
!$OMP END PARALLEL
         IF(l_morphocoupl) CALL sed_morpho

       END IF


! Balance
! ---------------------

       IF(l_bilan .AND. t .GE. tdebsubs) THEN
          IF (t >= t_out_budget) CALL sub_budget_out(t,cv_wat  &

                                         ,cv_sed      &

                                         ,ssh)
       END IF


       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, &

                                       cv_wat,                 &

                                       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
! ----------------------------------------------


     END IF

   END IF    ! end test on test cas casAdvCst or casAdvSpot

   CALL second_mars(ztime1)





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




! Outputs
! -------

   CALL output_mng

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





! 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) &

!$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 FINDMAX_MPI_3(umx,iumx,jumx)
     CALL FINDMAX_MPI_3(vmx,ivmx,jvmx)
     IF (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

     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,    &

                            cv_wat,                              &

                            ect   )

       CALL sed_outsaverestart

     END IF     

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

   END IF     ! Ends loop on odd nbouc

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

   !PRINT*, 'exit STEP routine'

  END SUBROUTINE step
