  SUBROUTINE dyn3d_dt(uz,vz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dyn3d_dt  ***
   !&E
   !&E ** Purpose : adaptative computation of the time-step by examining
   !&E              by examining the surface current number.
   !&E
   !&E ** Description : The time-step increase is applied after a lap of tobs,
   !&E                  if during this observation period the stability criterion
   !&E                  has been exceeded.
   !&E
   !&E ** Called by : stepu
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Used ij-arrays : uz, vz, dy_v, htx, hty
   !&E
   !&E ** Modified variables : dt, cflumaxmax, cflvmaxmax, tdebobs, tobs
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2004-12-29
   !&E       !  2007-09    (V. Garnier) homogeneneity of dx dy...
   !&E       !  2015-09    (S. Petton) detection of cfl over all layers 
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters,ONLY : limin,limax,ljmin,kmax,ljminp1,ljmaxm1,rsh,rlg, &
                         liminm2,limaxp2,ljminm2,ljmaxp2,epsilon,        &
                         jshift_obc_south,jshift_obc_north,jmax,valmanq,dxmin
   USE comvars2d, ONLY : ig,id,dx_u,dy_v,htx,hty,h0
   USE comvars2d, ONLY : cflumaxmax,cflvmaxmax,cflcrt ,cdate
   USE comvars2d, ONLY : dt,halfdt,tobs,dtmin,t,tdebobs,dtmax,iscreen
#if defined key_agrif
   USE comvars2d, ONLY : dtmaxi
#endif
   USE_MPI toolmpi, ONLY : FINDMAX_MPI,MASTER

#ifdef key_suivi_traceur3d
   USE comvars2d, ONLY : h0_g
   USE comvarp2d, ONLY : sal,temp,ssh 
   USE parameters,ONLY : valmanq
   USE_MPI toolmpi, ONLY : FINDMAX_MPI_5
#endif
#ifdef key_suivi_vitesse3d
   USE_MPI toolmpi, ONLY : FINDMAX_MPI_6
#endif
   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(in) :: uz,vz

   !! * Local declarations
   INTEGER        :: i,j,iumax,jumax,ivmax,jvmax,k
   INTEGER        :: liumx,ljumx,livmx,ljvmx
   REAL(KIND=rsh) :: flumx,flvmx
   REAL(KIND=rsh) :: uzmax,htxmax,dxumax,vzmax,htymax,dyvmax
   REAL(KIND=rsh) :: cflumax,cflvmax,cflmax,cflu,cflv
   REAL_AGRIF(KIND=rlg) :: dtmaxp

#ifdef key_suivi_traceur3d
   REAL(kind=rsh) :: tmax,smax,ssht,sshs,ltmax,lsmax
   INTEGER        :: itmax,jtmax,ktmax,ismax,jsmax,ksmax
   INTEGER        :: litmax,ljtmax,lktmax,lismax,ljsmax,lksmax
#endif
#ifdef key_suivi_vitesse3d
   INTEGER        :: kumax,kvmax,lkumx,lkvmx
#endif

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

   PRINT_DBG*,'enter DYN3DDT'

   ! on examine le nb de courant max en surface
   !-------------------------------------------

   cflumax=-100.0_rsh
   cflvmax=-100.0_rsh
   iumax=limin
   jumax=ljmin
   ivmax=limin
   jvmax=ljmin

#ifdef key_suivi_vitesse3d
   kumax=1
   kvmax=1
#endif

#ifdef key_suivi_traceur3d
   tmax=-100.0_rsh
   smax=-100.0_rsh
   itmax=limin
   jtmax=ljmin
   ktmax=1
   ismax=limin
   jsmax=ljmin
   ksmax=1
#endif
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP& PRIVATE(i,j,liumx,ljumx,livmx,ljvmx,flumx,flvmx,cflu,cflv, &
#ifdef key_suivi_vitesse3d
!$OMP&        lkumx,lkvmx,&
#endif
#ifdef key_suivi_traceur3d
!$OMP&        ltmax,litmax,ljtmax,lktmax,lsmax,lismax,ljsmax,lksmax,&
#endif
!$OMP&        k) &
!$OMP& SHARED(cflumax,cflvmax,uz,vz,dt,dx_u,iumax,jumax,ivmax,jvmax,ig,id,dy_v, &
#ifdef key_suivi_vitesse3d
!$OMP&        kumax,kvmax, &
#endif
#ifdef key_suivi_traceur3d
!$OMP&        itmax,jtmax,ktmax,tmax,ismax,jsmax,ksmax,smax,temp,sal, &
#endif
#ifdef key_agrif
!$OMP&        jmax,jshift_obc_south,jshift_obc_north, &
#endif
!$OMP&        ljminp1,ljmaxm1,ljmin,limin,limax)

   flumx = cflumax
   flvmx = cflvmax

#ifdef key_suivi_traceur3d
   ltmax = tmax
   lsmax = smax
#endif

#ifdef key_MPI_2D
!$OMP do SCHEDULE(runtime)
#else
!$OMP do SCHEDULE(static)
#endif

   !DO j=ljminp1,ljmaxm1
   DO j=MAX0(ljminp1,jshift_obc_south+1), MIN0(ljmaxm1,jmax-1-jshift_obc_north)
     DO i=MAX0(limin,ig(j)),MIN0(limax,id(j)-1)
      IF_CURV (dx_u(i,j)/=dxmin .and. dy_v(i,j)/=dxmin) THEN

#ifdef key_suivi_vitesse3d
       DO k=1,kmax
         cflu=ABS(uz(k,i,j))*dt/dx_u(i,j)
         IF(cflu > flumx) THEN
           flumx=cflu
           liumx=i
           ljumx=j
           lkumx=k
         ENDIF

         cflv=ABS(vz(k,i,j))*dt/dy_v(i,j)
         IF(cflv > flvmx) THEN
           flvmx=cflv
           livmx=i
           ljvmx=j
           lkvmx=k
         ENDIF
       END DO
#else
       cflu=ABS(uz(kmax,i,j))*dt/dx_u(i,j)
       IF(cflu > flumx) THEN
         flumx=cflu
         liumx=i
         ljumx=j
       ENDIF

       cflv=ABS(vz(kmax,i,j))*dt/dy_v(i,j)
       IF(cflv > flvmx) THEN
         flvmx=cflv
         livmx=i
         ljvmx=j
       ENDIF
#endif

#ifdef key_suivi_traceur3d
       do k=1,kmax
         if ( temp(k,i,j) > ltmax .and. temp(k,i,j)/=valmanq) then
           ltmax=temp(k,i,j)
           litmax=i
           ljtmax=j
           lktmax=k
         end if
         if ( sal(k,i,j) > lsmax .and. sal(k,i,j)/=valmanq ) then
           lsmax=sal(k,i,j)
           lismax=i
           ljsmax=j
           lksmax=k
         end if
       end do
#endif
      ENDIF_CURV
     ENDDO
   ENDDO
!$OMP end do nowait

!$OMP CRITICAL
   IF(flumx.GT.cflumax) THEN
     cflumax = flumx
     iumax = liumx
     jumax = ljumx
#ifdef key_suivi_vitesse3d
     kumax = lkumx
#endif
   ENDIF
   IF(flvmx.GT.cflvmax) THEN
     cflvmax = flvmx
     ivmax = livmx
     jvmax = ljvmx
#ifdef key_suivi_vitesse3d
     kvmax = lkvmx
#endif
   ENDIF
#ifdef key_suivi_traceur3d
   IF(ltmax.GT.tmax) THEN
     tmax  = ltmax
     itmax = litmax
     jtmax = ljtmax
     ktmax = lktmax
   ENDIF
   IF(lsmax.GT.smax) THEN
     smax  = lsmax
     ismax = lismax
     jsmax = ljsmax
     ksmax = lksmax
   ENDIF
#endif
!$OMP END CRITICAL
!$OMP END PARALLEL
#if ! defined key_suivi_vitesse3d
   uzmax=uz(kmax,iumax,jumax)
   htxmax=htx(iumax,jumax)
   dxumax=dx_u(iumax,jumax)
   vzmax=vz(kmax,ivmax,jvmax)
   htymax=hty(ivmax,jvmax)
   dyvmax=dy_v(ivmax,jvmax)
   CALL_MPI FINDMAX_MPI(cflumax,iumax,jumax,uzmax,htxmax,dxumax)
   CALL_MPI FINDMAX_MPI(cflvmax,ivmax,jvmax,vzmax,htymax,dyvmax)

   IF_MPI (MASTER) THEN 
   WRITE(iscreen,*) 'DATE : ',cdate
   WRITE(iscreen,*) 'dt=',dt
!   WRITE(iscreen,*) 'cflumax=',cflumax
!   WRITE(iscreen,*) 'cflvmax=',cflvmax
   WRITE(iscreen,*) 'max en u = ',uzmax,iumax,jumax,htxmax
   WRITE(iscreen,*) 'max en v = ',vzmax,ivmax,jvmax,htymax
   ENDIF_MPI
#else
   uzmax=uz(kumax,iumax,jumax)
   htxmax=htx(iumax,jumax)
   dxumax=dx_u(iumax,jumax)
   vzmax=vz(kvmax,ivmax,jvmax)
   htymax=hty(ivmax,jvmax)
   dyvmax=dy_v(ivmax,jvmax)
   CALL_MPI FINDMAX_MPI_6(cflumax,iumax,jumax,kumax,uzmax,htxmax,dxumax)
   CALL_MPI FINDMAX_MPI_6(cflvmax,ivmax,jvmax,kvmax,vzmax,htymax,dyvmax)

   IF_MPI (MASTER) THEN 
   WRITE(iscreen,*) 'DATE : ',cdate
   WRITE(iscreen,*) 'dt=',dt
   WRITE(iscreen,'(a11,f17.8,2x,2(i4,3x),f17.8,2x,i3)') 'max en u = ',uzmax,iumax,jumax,htxmax,kumax
   WRITE(iscreen,'(a11,f17.8,2x,2(i4,3x),f17.8,2x,i3)') 'max en v = ',vzmax,ivmax,jvmax,htymax,kvmax
   ENDIF_MPI
#endif

#ifdef key_suivi_traceur3d
   ssht=ssh(itmax,jtmax)
   sshs=ssh(ismax,jsmax)
   CALL_MPI FINDMAX_MPI_5(tmax,itmax,jtmax,ktmax,ssht)
   CALL_MPI FINDMAX_MPI_5(smax,ismax,jsmax,ksmax,sshs)
   IF_MPI (MASTER) THEN 
   WRITE(iscreen,'(a11,f17.8,2x,2(i4,3x),f17.8,2x,i3,2x,f17.8)') 'max temp = ', &
       tmax,itmax,jtmax,h0_g(itmax,jtmax)+ssht,ktmax,h0_g(itmax,jtmax)
   WRITE(iscreen,'(a11,f17.8,2x,2(i4,3x),f17.8,2x,i3,2x,f17.8)') 'max sal  = ', &
       smax,ismax,jsmax,h0_g(ismax,jsmax)+sshs,ksmax,h0_g(ismax,jsmax)
   ENDIF_MPI
#endif

   IF(MAX(cflumax,cflvmax) > cflcrt) THEN
     IF_MPI (MASTER) THEN 
     WRITE(iscreen,*) 'ALERTE ADS!!!!!!!!!!!!!'
     ENDIF_MPI
     IF_AGRIF (agrif_root()) THEN
       dt=MIN(0.83_rsh*cflcrt*dyvmax/MAX(REAL(ABS(vzmax),rlg),epsilon), &
              0.83_rsh*cflcrt*dxumax/MAX(REAL(ABS(uzmax),rlg),epsilon))
       dt=MAX(dt,dtmin)
       halfdt=dt*0.5_rlg
       IF_MPI (MASTER) THEN 
       WRITE(iscreen,*) 'HIIIIIIIII!!!! nouveau dt=',dt
       ENDIF_MPI
     ELSE_AGRIF
       INST_AGRIF dtmaxp=MAX(MIN(0.83_rsh*cflcrt*dyvmax/MAX(REAL(ABS(vzmax),rlg),epsilon), &
       INST_AGRIF                0.83_rsh*cflcrt*dxumax/MAX(REAL(ABS(uzmax),rlg),epsilon)),dtmin)
! Set the next time step of the parent grid
       CALL_AGRIF Agrif_Set_Parent(dtmax, MAX(Agrif_parent(dtmin),dtmaxp*Agrif_Rhot()))
       IF_AGRIF (cflumax > cflcrt) CALL Agrif_Set_Parent(cflumaxmax, cflumax)
       IF_AGRIF (cflvmax > cflcrt) CALL Agrif_Set_Parent(cflvmaxmax, cflvmax)
       IF_MPI (MASTER) THEN 
       INST_AGRIF WRITE(iscreen,*) 'HIIIIIIIII!!!! desired dt=',dtmaxp
       ENDIF_MPI
     ENDIF_AGRIF
   ENDIF

#if defined key_agrif
! Maximum time step may have been restricted by a child grid
! Done whaterver the rank because a grand-daughter must slow down a daughter which will slow down the mother
! The slackening is done with a time step late
! dtmaxi = initialized through paraspec.txt (=dtmax) !!!
   IF (dt > dtmax) THEN
     dt = dtmax
     halfdt=dt*0.5_rlg
     IF_MPI (MASTER) THEN
     WRITE(iscreen,*) 'HIIIIIIIII!!!! nouveau dt=',dt,' due to AGRIF'
     ENDIF_MPI
     dtmax = dtmaxi
   ELSE IF (dt==dtmax) THEN
     dtmax = dtmaxi ! In case of dtmax (from child grid) equals dtmin (of root grid) and to allow the speed up later
   ENDIF

! Maintain fine grid time in agreement with coarse grid time step and time refinement factor
   IF (.NOT.Agrif_Root()) THEN
     dt = Agrif_Parent(dt)/Agrif_Rhot()
     IF (Agrif_Grandmother()) THEN
         dt = Agrif_Parent(dt) * Agrif_irhot() ! For the grandmother grid, Agrif_Parent is Agrif_Root
                                               ! A multiplication is required (instead of a division)
     END IF
     halfdt=dt*0.5_rlg
   ENDIF
#endif


   cflumaxmax=MAX(cflumaxmax,cflumax)
   cflvmaxmax=MAX(cflvmaxmax,cflvmax)

   !***************************************
   ! fin de la periode d observation
   !***************************************

   IF (t > tdebobs+tobs) THEN
     cflmax=MAX(cflumaxmax,cflvmaxmax)
     IF_MPI (MASTER) THEN 
     WRITE(iscreen,*) 'ATTENTION  le ',cdate
     WRITE(iscreen,*) 'cflumaxmax=',cflumaxmax
     WRITE(iscreen,*) 'cflvmaxmax=',cflvmaxmax
     ENDIF_MPI

     ! Reinitialisation for the next observation period 
     cflumaxmax=0.0_rsh
     cflvmaxmax=0.0_rsh
     tdebobs=t

     IF_AGRIF (.NOT.Agrif_Root()) RETURN ! Only RootGrid can accelerate

     IF(cflmax < cflcrt*0.84_rsh) THEN
       IF_MPI (MASTER) THEN 
       WRITE(iscreen,*) 'TURBOOOOOOOOOOOOOOOOOO.......!!!!!! old dt=',dt
       ENDIF_MPI
       dt=MIN(dt*1.2_rlg,dtmax)
       halfdt=dt*0.5_rlg
       IF_MPI (MASTER) THEN 
       WRITE(iscreen,*) 'Nouveau dt : ',dt
       ENDIF_MPI
     ENDIF
   ENDIF


   PRINT_DBG*,'exit DYN3DDT'

  END SUBROUTINE dyn3d_dt
