sisyphe.F

Go to the documentation of this file.

!                   ******************
                    SUBROUTINE sisyphe
!                   ******************
!
     &(part,loopcount,grafcount,listcount,telnit,
     & u_tel,v_tel,h_tel,hn_tel,zf_tel,uetcar,cf_tel,ks_tel,
     & constflow,nsis_cfd,sisyphe_cfd,code,pericou,
     & u3d,v3d,t_tel,visc_tel,dt_tel,charr_tel,susp_tel,
     & flbor_tel,solsys,dm1,uconv_tel,vconv_tel,zconv,
     & thetaw_tel,hw_tel,tw_tel,uw_tel,yagout,api_iter,grcomp)
!
!***********************************************************************
! SISYPHE   V7P3
!***********************************************************************
!
!brief  The real main program of Sisyphe, with the time loop.
!
!history  R-S MOURADI (LNHE)
!+        24/03/2016
!+        V7P1
!+  Adding 2 possible values for PART, to run SISYPHE from the API without coupling
!+  Adding an additional variable API_ITER, current iteration number for the API.
!
!history  C. LENORMANT; J-M HERVOUET; C. MACHET; C. VILLARET; U. MERKEL; R. KOPMANN
!+        20/03/2011
!+        V6P1
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  C. VILLARET (LNHE)
!+        01/03/2012
!+        V6P2
!+   Second call to condim moved upwards out of test for calling BIEF_SUITE.
!
!history  MAK (HRW)
!+        31/05/2012
!+        V6P2
!+   modifications to include CSRATIO
!
!history  JWI (HRW)
!+        14/06/2012
!+        V6P2
!+   added (several) lines to use wave orbital velocities directly if found in hydro file
!
!history  PAT (LNHE)
!+        18/06/2012
!+        V6P2
!+   updated version with HRW's development (wave orbital velocities) + Soulsby-van Rijn's concentration
!
!history  CV (LNHE)
!+        12/06/2012
!+        V6P2
!+   added AT0 to CALL CONLIT(MESH%NBOR%I,AT0)
!
!history  CV (LNHE)
!+        02/07/2012
!+        V6P2
!+   DT changed to DTS in CALL FLUSEC_SISYPHE
!
!history  CV (LNHE)
!+        28/08/2012
!+        V6P2
!+   Modification call to init_sediment and suspension_main
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        08/03/2013
!+        V6P3
!+   Adding multiclass treatment in slides.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        22/03/2013
!+        V6P3
!+   Adding arguments THETAW_TEL, HW_TEL, TW_TEL for variables
!+   transmitted from Tomawac to Sisyphe through Telemac-2D or 3D in the
!+   of triple coupling.
!
!history  R. KOPMANN (EDF R&D, LNHE)
!+        16/04/2013
!+        V6P3
!+   Adding the file format in the call to FONSTR.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        02/01/2014
!+        V7P0
!+   KNOGL removed from call to flusec_sisyphe.
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        28/04/2014
!+        V7P0
!+   Use of KP1BOR removed, replaced by IKLBOR.
!+   TPREC replaced by (TPREC-AT0) in formulas giving NUMEN0.
!+   OPTSUP replaced by OPTADV in the call to suspension_main
!+   (see keyword SCHME OPTION FOR ADVECTION)
!
!history  J-M HERVOUET (EDF R&D, LNHE)
!+        02/05/2014
!+        V7P0
!+   ZF_SIS renamed ZF_TEL (it is ZF in the memory of Telemac-2D or 3D
!+   so clearer like this isn't it?
!+   When a new ZF is received from Telemac, ELAY updated consequently,
!+   this was a long time overlooked bug, that can be seen when coupling
!+   with Telemac-3D with non-erodable beds.
!
!history  C VILLARET (HRW+EDF) & J-M HERVOUET (EDF - LNHE)
!+        18/09/2014
!+        V7P0
!+   Adding the variable UW_TEL in argument (orbital velocity)
!+   for getting it back to the calling program.
!
!history Y AUDOUIN (LNHE)
!+       25/05/2015
!+       V7P0
!+       Modification to comply with the hermes module
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        08/02/2016
!+        V7P2
!+   Adding the variable HPROP.
!
!history  R KOPMANN (BAW)
!+        10/05/2016
!+        V7P2
!+ CALFA,SALFA dependent of grain classes
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        21/06/2016
!+        V7P2
!+   Changing the formula for VALNIT when coupling. For some reason the
!+   old formula is no longer valid.
!
!history S.E.BOURBAN (HRW))
!+       01/01/2017
!+       V7P2
!+   Adding differentiated variables
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!
!history  N. HUYBRECHTS & P. TASSI
!+        December 2017
!+        V7P3
!+   Correction of input argument for the subroutine TASSEMENT_2
!+   ES replaced by ES_VASE
!
!history  B.GLANDER (BAW)
!+        06/12/2018
!+        V7P2
!+  NEW VARIABLE: ZRL  REFERENCE LEVEL FOR NESTOR
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CF_TEL         |<->| QUADRATIC FRICTION COEFFICIENT FROM TELEMAC
!| CHARR_TEL      |<->| LOGICAL, BED LOAD OR NOT: Sent to TELEMAC-2D
!| CODE           |-->| NAME OF CALLING PROGRAMME (TELEMAC2D OR 3D)
!| CONSTFLOW      |<->| LOGICAL, CONSTANT FLOW DISCHARGE OR NOT (A SUPPRIMER)
!| CSRATIO        |<->| EQUILIBRIUM CONCENTRATION FOR SOULSBY-VAN RIJN EQ.
!| DM1            |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV)
!| DT_TEL         |-->| TIME STEP FROM TELEMAC
!| FLBOR_TEL      |-->| FLOW FLUXES AT BOUNDARIES
!| GRAFCOUNT      |-->| PERIOD OF GRAPHICAL OUTPUTS
!| GRCOMP         |-->| COUNTER FOR GRAPHICAL OUTPUTS
!| HN_TEL         |-->| WATER DEPTH FROM TEL HN
!| H_TEL          |-->| WATER DEPTH FROM TEL H (N+1)
!| H_PROP         |-->| WATER DEPTH FROM DIVERGENCE TERM IN CONTINUITY.
!| KS_TEL         |-->| BED ROUGHNESS SENT TO TELEMAC
!| LISTCOUNT      |-->| PERIODE DE SORTIE LISTING
!| LOOPCOUNT      |-->| NUMERO DE L'ITERATION
!| UW_TEL         |-->| ORBITAL VELOCITY
!| NSIS_CFD       |---| NUMBER OF ITERATIONS FOR TELEMAC (CONSTANT FLOW DISCHARGE OPTION-TO BE SUPRESSED)
!| PART           |-->| IF -1, NOT COUPLING : ON PASSE TOUTE LA
!|                |   | SUBROUTINE. SINON, INDIQUE LA PARTIE DE LA
!|                |   | SUBROUTINE DANS LAQUELLE ON PASSE
!| PERICOU        |-->| COUPLING PERIOD FOR BEDLOAD
!| SISYPHE_CFD    |<->| LOGICAL, CONSTANT FLOW DISCHARGE OR NOT
!| SOLSYS         |-->|1 OR 2. IF 2 ADVECTION FIELD IS UCONV + DM1*GRAD(ZCONV)
!| SUSP_TEL       |<->|LOGICAL, SUSPENDED LOAD OR NOT: Sent to TELEMAC-2D
!| TELNIT         |-->| NUMBER OF TELEMAC ITERATIONS
!| T_TEL          |-->| CURRENT TIME IN CALLING PROGRAMME
!| U3D,V3D        |-->| 3D VELOCITY SENT BY TELEMAC 3D
!| UCONV_TEL      |-->| ADVECTION VELOCITY FROM TELEMAC2D (X-DIRECTION)
!| UETCAR         |<->| SQUARE OF THE FRICTION VELOCITY (COUPLING WITH TEL 3D)
!| U_TEL          |-->| U VELOCITY FROM TELEMAC
!| VCONV_TEL      |-->| ADVECTION VELOCITY FROM TELEMAC2D (Y-DIRECTION)
!| VISC_TEL       |-->| VELOCITY DIFFUSIVITY (TELEMAC-2D)
!| V_TEL          |-->| V VELOCITY FROM TELEMAC
!| ZCONV          |-->| THE PIECE-WISE CONSTANT PART OF ADVECTION FIELD
!|                |   | IS DM1*GRAD(ZCONV), SEE SOLSYS.
!| ZF_TEL         |<->| BOTTOM ELEVATION OF THE CALLING TELEMAC
!| API_ITER       |<->| OPTIONAL - CURRENT ITERATION NUMBER IN THE API_SISYPHE
!| YAGOUT         |-->| LOGICAL: IF YES GRAPHIC OUTPUT (STEERED BY T2D)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE interface_sisyphe, ex_sisyphe => sisyphe
      USE bief
      USE declarations_telemac
      USE declarations_sisyphe
      USE interface_hermes
      USE declarations_special
!
      USE interface_parallel, ONLY : p_max,p_min,p_max
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER,           INTENT(IN)    :: PART,LOOPCOUNT,GRAFCOUNT
      INTEGER,           INTENT(IN)    :: LISTCOUNT,TELNIT,PERICOU
      CHARACTER(LEN=24), INTENT(IN)    :: CODE
      TYPE(bief_obj),    INTENT(IN)    :: U_TEL,V_TEL,H_TEL,HN_TEL
      TYPE(bief_obj),    INTENT(INOUT) :: ZF_TEL,UETCAR,KS_TEL
      INTEGER,           INTENT(INOUT) :: NSIS_CFD
      LOGICAL,           INTENT(INOUT) :: CONSTFLOW,SISYPHE_CFD
      TYPE(bief_obj),    INTENT(IN)    :: U3D,V3D,VISC_TEL
      TYPE(bief_obj),    INTENT(INOUT) :: CF_TEL
      DOUBLE PRECISION,  INTENT(IN)    :: T_TEL
      LOGICAL,           INTENT(INOUT) :: CHARR_TEL,SUSP_TEL
      DOUBLE PRECISION,  INTENT(IN)    :: DT_TEL
      INTEGER,           INTENT(IN)    :: SOLSYS
      TYPE(bief_obj),    INTENT(IN)    :: FLBOR_TEL,DM1,ZCONV
      TYPE(bief_obj),    INTENT(IN)    :: UCONV_TEL,VCONV_TEL
      TYPE(bief_obj),    INTENT(IN)    :: THETAW_TEL,HW_TEL,TW_TEL
      TYPE(bief_obj),    INTENT(IN)    :: UW_TEL
      LOGICAL,           INTENT(IN)    :: YAGOUT
      INTEGER, OPTIONAL, INTENT(IN)    :: API_ITER,GRCOMP
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, PARAMETER :: NHIST = 0
      INTEGER, PARAMETER :: NSOR = 100
      INTEGER            :: NLISS
      INTEGER            :: I,J,K,MN,MT,ISOUS,IMA,IMI,IELEB,KP1,IVAR
      INTEGER            :: IERR
      INTEGER            :: IMIN,IMAX,NUMENX,NUMDEB
      INTEGER            :: TROUVE(maxvar+10)
      DOUBLE PRECISION   :: DTS,BID,XMA,XMI
      DOUBLE PRECISION   :: XMIN,XMAX
      DOUBLE PRECISION   :: AT,AT2
      LOGICAL            :: ENTETS,YAZR
!
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEZF,SAVEQU,SAVEQV
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEZ
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEUW
!
      INTEGER, ALLOCATABLE :: NULLT(:)
      DOUBLE PRECISION, ALLOCATABLE :: NULLD(:)
!
!     SAVES LOCAL VARIABLES
!     NUMEN0 : 1ST RECORD TO READ
      INTEGER :: NUMEN0
      INTEGER :: DEBUT,FIN
!
!     VARIABLES TO READ IF COMPUTATION IS CONTINUED
!     --------------------------------
!     0 : DISCARD
!     1 : READ  (SEE SUBROUTINE NOMVAR)
!
!     HYDRO + EVOLUTION
      INTEGER :: ALIRE(maxvar) =
     &         (/ (1,i=1,9),(0,i=10,maxvar) /)
!   WAVES ONLY
      INTEGER :: ALIRH(maxvar) =
     &         (/ (0,i=1,11),(1,i=12,14),(0,i=15,maxvar) /)
!
!    FOR VALIDATION, EACH VARIABLE IN THE FILE IS COMPARED
!
      INTEGER :: ALIRV(maxvar) = (/ (1,i=1,maxvar) /)
!
!======================================================================!
!======================================================================!
!                               PROGRAM                                !
!======================================================================!
!======================================================================!
!
#if defined COMPAD
      CALL ad_sisyphe_begin
      CALL ad_sisyphe_initialisation_begin
#endif
!
!------------------------------------------------------------------
!     PART 1 : INITIALISATION
!------------------------------------------------------------------
!
      percou = pericou !UHM!!!!!!!
!
      IF(part==0.OR.part==-1.OR.part==2) THEN
        IF(debug.GT.0) WRITE(lu,*) 'INITIALIZATION'
!
        WRITE(lu,*) 'PART 0 : INITIALISING SISYPHE'
!
!       INITIALISES THE CONSTANTS
!
        CALL init_constant(karim_holly_yang,karman,pi)
!
!      READS THE WAVE DATA IN THE HYDRODYNAMIC FILE
!
        IF(houle.AND.sis_files(siscou)%NAME(1:1).EQ.' ') THEN
          alire(12)=1
          alire(13)=1
          alire(14)=1
          alire(22)=1
        ENDIF
!
!       READS THE SEDIMENTOLOGICAL DATA IN THE CONTINUATION FILE
!
        IF(debu) THEN
          alire(15)=1
          alire(16)=1
          alire(17)=1
          alire(18)=1
          alire(19)=1
          alire(20)=1
          alire(21)=1
!         TO READ REFERENCE LEVEL ZRL FOR NESTOR IN THE RESTART FILE
          IF(nestor)  alire(23)=1
!         READS AVAI FROM THE PREVIOUS COMPUTATION FILE
          DO i=1,nsicla*nomblay
            alire(23+i)=1
          ENDDO
!         READS CS (CONCENTRATION) FROM THE PREVIOUS COMPUTATION FILE
          IF(susp) THEN
            DO i=1,nsicla
              alire(23+(nomblay+1)*nsicla+i)=1
            ENDDO
          ENDIF
!         READS THE LAYER THICKNESSES
          DO i=1,nomblay
            alire(29+(nomblay+4)*nsicla+i)=1
          ENDDO
        ENDIF
!       V6P2 CV lecture concentration des couches
        IF(tass) THEN
          DO i=1,nomblay
            alire(29+(nomblay+4)*nsicla+i+nomblay)=1
          ENDDO
        ENDIF
!
!       INITIALISING DEPOSITS ON THE BOTTOM
!       THEY MUST BE USED IN BILAN_SISYPHE EVEN IF SUSP=.FALSE.
!
        DO i=1,nsicla
          masdep(i)=0.d0
          masdept(i)=0.d0
        ENDDO
!
!       DIFFERENTIATED VARIABLES
!       FOR READING GRADIENTS IN SELAFIN FILES
!
        k = 29+max(4,npriv)+nsicla*(nomblay+4)+2*nomblay+varcl%N
        IF(nadvar.GT.0) THEN
          DO ivar=1,nadvar
!           SEE NOMVAR_SISYPHE
            alire(k+ivar) = 1
          ENDDO
        ENDIF
!
! --------  INITIALISES (SETS TO 0) THE ARRAYS
!
        CALL init_zero
!
!       -- INITIALISES THE REFERENCE LEVEL ZRL FOR NESTOR
!          IT MUST BE EXACLY 123456789.0D0
        zrl%R(:) = 123456789.0d0
!
! --------  END OF INITIALISATION
!
!       DISCRETISATION : LINEAR FOR THE TIME BEING
!
!       IELMT HARD-CODED IN LECDON
!
        ielmx  = max(ielmt,ielmh_sis,ielmu_sis)
        nelmax = nelem
!
!=======================================================================
!
! : 1          READS, PREPARES AND CONTROLS THE DATA
!
!=======================================================================
!
!       READS THE BOUNDARY CONDITIONS AND INDICES FOR THE BOUNDARY NODES
!       EBOR IS READ HERE FOR THE FIRST CLASS ONLY
!       SEE CONLIT FOR OTHER CLASSES
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE LECLIM'
        ALLOCATE(nullt(nptfr),stat=ierr)
        CALL check_allocate(ierr,'NULLT')
        ALLOCATE(nulld(nptfr*2),stat=ierr)
        CALL check_allocate(ierr,'NULLD')
        CALL leclim(lihbor%I,liqbor%I,nullt,liebor%I,q2bor%R,nulld,
     &              nulld,ebor%ADR(1)%P%R,nulld,afbor%R,bfbor%R,
     &              mesh%NPTFR,'SIS',.true.,
     &              sis_files(sisgeo)%FMT,sis_files(sisgeo)%LU,
     &              kent,kent,-1,-1,-1,-1,
     &              numliq%I,mesh,boundary_colour=boundary_colour%I)
        DEALLOCATE(nullt)
        DEALLOCATE(nulld)
        IF(debug.GT.0) WRITE(lu,*) 'END_LECLIM'
!
!-----------------------------------------------------------------------
!
!       COMPLEMENTS THE DATA STRUCTURE FOR BIEF
!
        IF(debug.GT.0) WRITE(lu,*) 'INBIEF'
        CALL inbief(it1%I,klog,it2,it3,it4,lvmac,ielmx,
     &                 0.d0,spheri,mesh,t1,t2,optass,produc,equa)
!
        IF(debug.GT.0) WRITE(lu,*) 'END_INBIEF'
!
!-----------------------------------------------------------------------
!
!       LOCATES THE BOUNDARIES
!
        IF(ncsize.GT.1) THEN
          nfrliq=0
          DO i=1,nptfr
            nfrliq=max(nfrliq,numliq%I(i))
          ENDDO
          nfrliq=p_max(nfrliq)
          WRITE(lu,*) ' '
          WRITE(lu,*) 'LIQUID BOUNDARIES:',nfrliq
!         NFRLIQ CHECKED (NFRSOL NOT USED IN PARALLEL)
          IF(nfrliq.GT.maxfro) THEN
            WRITE(lu,*) 'INCREASE THE MAXIMUM NUMBER OF BOUNDARIES'
            WRITE(lu,*) 'CURRENTLY AT ',maxfro
            WRITE(lu,*) 'TO THE VALUE ',nfrliq
            CALL plante(1)
            stop
          ENDIF
        ELSE
          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FRONT2'
          CALL front2(nfrliq,
     &                liebor%I,liebor%I,
     &                mesh%X%R,mesh%Y%R,mesh%NBOR%I,mesh%KP1BOR%I,
     &                it1%I,npoin,nptfr,klog,.true.,numliq%I,maxfro)
          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FRONT2'
        ENDIF
        IF(nfrliq.GT.maxfro) THEN
          WRITE(lu,*) 'FRONT2: SIZE OF ARRAYS EXCEEDED'
          WRITE(lu,*) '        INCREASE THE KEYWORD'
          WRITE(lu,*) '        MAXIMUM NUMBER OF BOUNDARIES'
          WRITE(lu,*) '        IN THE CALLING PROGRAM'
          WRITE(lu,*) '        THE CURRENT VALUE IS ',maxfro
          WRITE(lu,*) '        THE VALUE SHOULD BE ',nfrliq
          CALL plante(1)
          stop
        ENDIF
!
!-----------------------------------------------------------------------
!       LOOKS FOR BOTTOM AND BOTTOM FRICTION IN THE GEOMETRY FILE :
!-----------------------------------------------------------------------
!
        IF(debug.GT.0) WRITE(lu,*) 'FONSTR'
        CALL fonstr(t1,zf,t2,chestr,sis_files(sisgeo)%LU,
     &              sis_files(sisgeo)%FMT,
     &              sis_files(sisfon)%LU,sis_files(sisfon)%NAME,
     &              mesh,sfon,.true.,
     &              n_names_priv,names_prive,prive)
        IF(debug.GT.0) WRITE(lu,*) 'END_FONSTR'
!
!
!-----------------------------------------------------------------------
!       LOOKS FOR RADSEC FOR SECONDARY CURRENTS IN THE GEOMETRY FILE :
!-----------------------------------------------------------------------
!
!
!       RADSEC IS LANGUAGE INDEPENDENT
!        CALL FIND_VARIABLE(RADSEC,'RADSEC          ',SIS_FILES(SISGEO)%LU,
!     &                   SIS_FILES(SISGEO)%FMT,W,OKRADSEC,TIME=BID)
! Initialising
        radsec%R=0.d0
        CALL find_variable(sis_files(sisgeo)%FMT,sis_files(sisgeo)%LU,
     &   'RADSEC          ', radsec%R,npoin,ierr)
!
!-----------------------------------------------------------------------
!       LOOKS FOR RFERENCE LEVEL FOR NESTOR IN THE GEOMETRY FILE :
!-----------------------------------------------------------------------
!       REFERENCE LEVEL IS LANGUAGE INDEPENDENT
!     Initialising
        CALL find_variable(sis_files(sisgeo)%FMT,sis_files(sisgeo)%LU,
     &   'REFERENCE LEVEL ', zrl%R,npoin,ierr)
        IF( ierr == 0 ) THEN
          WRITE(6,*)'-------------------------------------------'
          WRITE(6,*)
     &    'Found REFERENCE LEVEL (ZRL) for Nestor in Sis Geometry file'
          WRITE(6,*)'          max value of ZRL = ',maxval(zrl%R)
          WRITE(6,*)'          min value of ZRL = ',minval(zrl%R)
          WRITE(6,*)'-------------------------------------------'
        ENDIF
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!       PREPARES THE RESULTS FILE (OPTIONAL)
!
!       STANDARD SELAFIN FORMAT
!
        IF(debug.GT.0) WRITE(lu,*) 'WRITE_HEADER'
!       CREATES DATA FILE USING A GIVEN FILE FORMAT : FORMAT_RES
!       THE DATA ARE CREATED IN THE FILE: SISRES, AND ARE
!       CHARACTERISED BY A TITLE AND NAME OF OUTPUT VARIABLES
!       CONTAINED IN THE FILE.
        CALL write_header(sis_files(sisres)%FMT, ! RESULTS FILE FORMAT
     &                    sis_files(sisres)%LU,  ! LU FOR RESULTS FILE
     &                    titca,      ! TITLE
     &                    maxvar,     ! MAX NUMBER OF OUTPUT VARIABLES
     &                    texte,      ! NAMES OF OUTPUT VARIABLES
     &                    sorleo)     ! PRINT TO FILE OR NOT
        IF(debug.GT.0) WRITE(lu,*) 'END_WRITE_HEADER'
        IF(debug.GT.0) WRITE(lu,*) 'WRITE_MESH'
!       WRITES THE MESH IN THE OUTPUT FILE :
!       IN PARALLEL, REQUIRES NCSIZE AND NPTIR.
!       THE REST OF THE INFORMATION IS IN MESH.
!       ALSO WRITES : START DATE/TIME AND COORDINATES OF THE
!       ORIGIN.
        CALL write_mesh(sis_files(sisres)%FMT, ! RESULTS FILE FORMAT
     &                  sis_files(sisres)%LU,  ! LU FOR RESULTS FILE
     &                  mesh,
     &                  1,             ! NUMBER OF PLANES /NA/
     &                  mardat,        ! START DATE
     &                  martim,        ! START TIME
     &                  t1,t2,
     &                  ncsize.GT.1, nptir,
     &                  ngeo=sis_files(sisgeo)%LU,
     &                  geoformat=sis_files(sisgeo)%FMT)
        IF(debug.GT.0) WRITE(lu,*) 'END_WRITE_MESH'
!
!   --- FILLS IN MASKEL BY DEFAULT
!
        IF(msk) CALL os ('X=C     ', x=maskel, c=1.d0)
!
!       BUILDING THE MASK FOR LIQUID BOUNDARIES
!       A SEGMENT IS LIQUID IF BOTH ENDS ARE NOT SOLID
!
        DO ieleb = 1, mesh%NELEB
          k=mesh%IKLBOR%I(ieleb)
          kp1=mesh%IKLBOR%I(ieleb+mesh%NELEBX)
          IF(liebor%I(k).NE.2.AND.liebor%I(kp1).NE.2) THEN
            mask%R(ieleb) = 1.d0
          ELSE
            mask%R(ieleb) = 0.d0
          ENDIF
        ENDDO
!
!=======================================================================
!
! : 2                  INITIALISES
!
!=======================================================================
!
        pass      = .true.
!
        pass_susp = .true.
        vcumu     = 0.d0
        mass_gf   = 0.d0
!
!
!----   DETERMINES THE NUMBER OF EVENTS (1ST LOOP)       : NCALCU
!                      NUMBER OF TIMESTEPS (2ND LOOP)    : NIDT
!                      TOTAL NUMBER OF TIMESTEPS         : NIT
!                      INITIAL TIME                      : AT0
!                      TIMESTEP                          : DT
!
!
!       INITIALIZING IN CASE OF COUPLING
        IF(part.EQ.0.OR.(part.EQ.2.AND.code(1:7).EQ.'TELEMAC')) THEN
          at0=t_tel
          dt = dt_tel
          ncalcu = 1
          nidt   = 1
          nit=telnit
        ELSE
          at0=max(tprec,0.d0)
          dt=delt
          IF(perma) THEN
            ncalcu=1
            nidt=npas
            nit=nidt
            nsous=1
          ELSE
            ncalcu = nmaree
! COMPUTES DT AFTER READING THE HYDRO FILE
!           NIDT =  NINT ( PMAREE / DT + 0.1D0 )
!           NIT=NIDT*NCALCU
! ELSE
            nidt=npas
            nit=nidt*ncalcu
          ENDIF
        ENDIF
!
! UNSTEADY MODE : DT IS COMPUTED FROM THE HYDRO FILE
!                 NUMEN: TOTAL NUMBER OF RECORDS
        IF(sis_files(sishyd)%NAME(1:1).NE.' ')  THEN
!         JUST TO GET NUMEN AND DT
          ! GET THE RECORD NUMBER OF THE LAST TIMESTEP
          ! AND THE TIME STEP BETWEEN TWO RECORDS
          CALL get_data_ntimestep(sis_files(sishyd)%FMT,
     &                  sis_files(sishyd)%LU,numen,ierr)
          CALL check_call(ierr,'SISYPHE:GET_DATA_NTIMESTEP')
          IF(perma) THEN
            dt = delt
          ELSE
            ! GET TIME FOR ONE BEFORE LAST RECORD
            CALL get_data_time(sis_files(sishyd)%FMT,
     &         sis_files(sishyd)%LU,numen-2,at2,ierr)
            CALL check_call(ierr,'SISYPHE:GET_DATA_TIME:AT2')
            ! GET TIME FOR LAST RECORD
            CALL get_data_time(sis_files(sishyd)%FMT,
     &         sis_files(sishyd)%LU,numen-1,at,ierr)
            CALL check_call(ierr,'SISYPHE:GET_DATA_TIME:AT')
            dt = at - at2
            nidt =  nint( pmaree / dt + 0.1d0 )
            IF(abs(nidt*dt-pmaree) > 1.d-3) THEN
              WRITE(lu,102) nidt*dt
            ENDIF
            nit  = ncalcu * nidt
          ENDIF
        ENDIF
!
!       VALIDATES AGAINST THE LAST TIMESTEP
!
        valnit = nit
!
102     FORMAT(/,
     &         'CAUTION : THE PERIOD OF COMPUTATION IS NOT A MULTIPLE',
     &         /,'OF THE HYDRODYNAMIC FILE PRINTOUT PERIOD.',/,
     &       'THE LENGTH OF COMPUTATION WILL THEREFORE BE',g16.7,/,
     &       'SECONDS')
!
!  SISYPHE ONLY
!  -----------------------------------------------------------------------
!
!  NUMEN : NUMBER OF RECORDS IN THE HYDRODYNAMIC FILE
!  DT    : TIMESTEP OF THE HYDRODYNAMIC RECORDS
!  NUMEN0: 1ST RECORD TO READ FROM HYDRODYNAMIC FILE
!  TPREC : START TIME
!
        IF(part.EQ.-1.OR.(part.EQ.2.AND.code(1:7).NE.'TELEMAC')) THEN
          IF(.NOT.perma) THEN
            IF(tprec.GE.0.d0) THEN
              numen0 = nint((tprec-at0)/dt)+1
            ELSE
              numen0 = numen-int(pmaree/dt+1.1d0)
            ENDIF
          ELSE
            IF(tprec.GE.0.d0) THEN
              numen0 = nint((tprec-at0)/dt)+1
            ELSE
              numen0 = numen
            ENDIF
          ENDIF
!
          IF(sis_files(sishyd)%NAME(1:1).NE.' ')  THEN
            IF(debug.GT.0) WRITE(lu,*) 'READ_DATASET'
            ! The time step are numbered from 0
            numen0 = numen0 - 1
            CALL read_dataset(sis_files(sishyd)%FMT,
     &                sis_files(sishyd)%LU,varsor,npoin,numen0,at,
     &                textpr,trouve,alire,.true.,perma,maxvar)
!
!           TRACES IF WAVE DATA HAVE BEEN FOUND
!
            IF(houle) THEN
!             NOTE: BIEF_ALLVEC SETS %TYPR TO '?'
              IF(trouve(12).EQ.1) hw%TYPR='Q'
              IF(trouve(13).EQ.1) tw%TYPR='Q'
              IF(trouve(14).EQ.1) thetaw%TYPR='Q'
              IF(trouve(22).EQ.1) uw%TYPR='Q'
              IF(uw%TYPR=='Q') THEN
                WRITE(lu,*)
                WRITE(lu,*) 'WAVE RESULTS IN :',sis_files(sishyd)%NAME
                WRITE(lu,*)
                WRITE(lu,*) 'BOTTOM ORBITAL VELOCITY FOUND'
                WRITE(lu,*) 'THESE WILL BE USED DIRECTLY'
                WRITE(lu,*)
              ELSEIF(hw%TYPR=='Q'.AND.tw%TYPR=='Q') THEN
                WRITE(lu,*)
                WRITE(lu,*) 'WAVE RESULTS IN :',sis_files(sishyd)%NAME
                WRITE(lu,*)
                WRITE(lu,*) 'WAVE HEIGHT AND PERIOD FOUND'
                WRITE(lu,*) 'BOTTOM VELOCITY WILL BE COMPUTED IN CALCUW'
                WRITE(lu,*)
              ENDIF
            ENDIF
            IF(debug.GT.0) WRITE(lu,*) 'END_READ_DATASET'
            IF(debug.GT.0) WRITE(lu,*) 'RESCUE_SISYPHE'
            CALL rescue_sisyphe(hn%R,z%R,zf%R,
     &                        zr%R,es,hw%R,tw%R,thetaw%R,npoin,nomblay,
     &                        nsicla,trouve,alire,pass,icf,.true.,
     &                        maxvar)
            IF(debug.GT.0) WRITE(lu,*) 'END_RESCUE_SISYPHE'
!
!       DIFFERENTIATED VARIABLES
!
            IF( nadvar.GT.0 ) THEN
              DO ivar = 1,nadvar
                CALL ad_set_sisyphe(ivar,advar%ADR(ivar)%P)
              ENDDO
            ENDIF
!
          ENDIF
!
        ENDIF
!
!---- RESUMES SISYPHE COMPUTATION
!
        yazr=.false.
        IF(sis_files(sispre)%NAME(1:1).NE.' ')  THEN
!
!         READS THE HYDRO AND SEDIMENTOLOGICAL VARIABLES
!
          IF(debug.GT.0) WRITE(lu,*) 'READ_DATASET'
          CALL read_dataset(sis_files(sispre)%FMT,sis_files(sispre)%LU,
     &                      varsor,npoin,numenx,at0,textpr,trouve,alire,
     &                      .true.,.true.,maxvar)
!
          IF( trouve(23) == 1 ) THEN  ! found REFERENCE LEVEL for Nestor in sis previous comp. file
            WRITE(6,*)'-------------------------------------------'
            WRITE(6,*)
     &      'Found REFERENCE LEVEL (ZRL) for Nestor in the '
            WRITE(6,*)
     &      'Sisyphe restart file.'
            WRITE(6,*)'          max value of ZRL = ',maxval(zrl%R)
            WRITE(6,*)'          min value of ZRL = ',minval(zrl%R)
            WRITE(6,*)'-------------------------------------------'
          ENDIF
!
          IF(debug.GT.0) WRITE(lu,*) 'END_READ_DATASET'
!
          IF(debug.GT.0) WRITE(lu,*) 'RESCUE_SISYPHE'
          CALL rescue_sisyphe(hn%R,z%R,zf%R,
     &                        zr%R,es,hw%R,tw%R,thetaw%R,npoin,nomblay,
     &                        nsicla,trouve,alire,pass,icf,.true.,
     &                        maxvar)
          IF(trouve(9).EQ.1) yazr=.true.
          IF(debug.GT.0) WRITE(lu,*) 'SORTIE DE RESCUE_SISYPHE'
!
!         CHANGES THE UNITS OF CONCENTRATIONS
!
          IF(susp.AND.unit) THEN
            DO i=1,nsicla
              IF(trouve(23+(nomblay+1)*nsicla+i).EQ.1) THEN
                CALL os('X=CX    ',x=cs%ADR(i)%P,c=1.d0/xmvs)
              ENDIF
            ENDDO
          ENDIF
!
        ENDIF
!
!----   READS THE LAST RECORD : WAVE FILE
!
!       NOTE : SIS_FILES(SISCOU)%NAME SET TO ' ' IF HOULE=NO
!
        IF(sis_files(siscou)%NAME(1:1).NE.' ')  THEN
!
          IF(debug.GT.0) WRITE(lu,*) 'SUITE_HOULE'
          WRITE(lu,*) ' LECTURE HOULE :',sis_files(siscou)%NAME
          CALL read_dataset(sis_files(siscou)%FMT,sis_files(siscou)%LU,
     &                      varsor,npoin,numenx,at,textpr,trouve,alirh,
     &                      .true.,.true.,maxvar)
          IF(debug.GT.0) WRITE(lu,*) 'END_SUITE_HOULE'
!         TRACES IF WAVE DATA HAVE BEEN FOUND
          IF(trouve(12).EQ.1) hw%TYPR='Q'
          IF(trouve(13).EQ.1) tw%TYPR='Q'
          IF(trouve(14).EQ.1) thetaw%TYPR='Q'
!
        ENDIF
!
        IF(code(1:7) == 'TELEMAC'.AND.(part==0.OR.part==2)) THEN
!
          at0=t_tel
!
          WRITE(lu,*) 'INITIALISATION EN CAS DE COUPLAGE : PART=',part
!         INFORMATION ON SUSPENSION SENT BACK
          charr_tel = charr
          susp_tel = susp
!
!         OV INSTEAD OF OS IN ORDER TO AVOID PROBLEMS WITH QUASI-BUBBLE ELEMENTS
!         OPERATES ONLY ON THE (1:NPOIN) RANGE OF THE TELEMAC FIELDS
!         IT IS A HIDDEN DISCRETISATION CHANGE
!
          CALL ov('X=Y     ', x=u2d%R, y=u_tel%R, dim1=npoin)
          CALL ov('X=Y     ', x=v2d%R, y=v_tel%R, dim1=npoin)
          CALL ov('X=Y     ', x=hn%R, y=h_tel%R, dim1=npoin)
!         BOTTOM GIVEN BY CALLING PROGRAMME
          CALL os('X=Y     ', x=zf, y=zf_tel)
!
!         CLIPS NEGATIVE DEPTHS
!
          IF(optban.GT.0) THEN
            DO i = 1,npoin
              IF(hn%R(i).LT.hmin) THEN
                u2d%R(i)=0.d0
                v2d%R(i)=0.d0
                hn%R(i) = hmin
              ENDIF
            ENDDO
          ENDIF
!
!         CASE OF TRIPLE COUPLING
!
          IF(inclus(coupling,'TOMAWAC')) THEN
!           incident wave direction
            CALL os('X=Y     ',x=thetaw,y=thetaw_tel)
!           Wave period
            CALL os( 'X=Y     ', x=tw, y=tw_tel)
!           significant wave height
            CALL os( 'X=Y     ', x=hw , y=hw_tel)
!           Bottom orbital velocity
            CALL os( 'X=Y     ', x=uw , y=uw_tel)
            hw%TYPR='Q'
            tw%TYPR='Q'
            thetaw%TYPR='Q'
            uw%TYPR='Q'
          ENDIF
!
        ENDIF
!
!  ---- END COUPLING  -------------
!
        IF(debug.GT.0) WRITE(lu,*) 'CONDIM_SISYPHE'
        IF(.NOT.debu) THEN
        CALL condim_sisyphe
     &        (u2d%R,v2d%R,qu%R,qv%R,hn%R,zf%R,z%R,esomt%R,thetaw%R,
     &         q%R,hw%R,tw%R,mesh%X%R,mesh%Y%R,npoin,at0,pmaree)
        ENDIF
        IF(debug.GT.0) WRITE(lu,*) 'END_CONDIM_SISYPHE'
!
!       AT THIS LEVEL U2D,V2D,H AND ZF MUST HAVE BEEN DEFINED
!       EITHER BY READ_DATASET, CONDIM_SISYPHE OR CALLING PROGRAM
!
!       NOW COMPUTES FUNCTIONS OF U2D,V2D,HN AND ZF
!
!       FREE SURFACE
        IF(debug.GT.0) WRITE(lu,*) 'FREE SURFACE'
        CALL os('X=Y+Z   ', x=z, y=zf, z=hn)
        IF(debug.GT.0) WRITE(lu,*) 'END FREE SURFACE'
!
        IF(code(1:7).NE.'TELEMAC') THEN
!         PRODUCT H*
          CALL os('X=YZ    ', x=qu, y=u2d, z=hn)
!         PRODUCT H*V
          CALL os('X=YZ    ', x=qv, y=v2d, z=hn)
!         DISCHARGE
          CALL os('X=N(Y,Z)', x=q, y=qu, z=qv)
        ENDIF
!
!       CHECKS THE WAVE DATA
!
        IF(houle) THEN
          IF(uw%TYPR.EQ.'Q') THEN
!           IF(HW%TYPR    .NE.'Q'.OR.
          ELSEIF(hw%TYPR    .NE.'Q'.OR.
     &      tw%TYPR    .NE.'Q'.OR.
     &      thetaw%TYPR.NE.'Q') THEN
            WRITE(lu,*) ' '
            WRITE(lu,*) ' '
            WRITE(lu,*) 'MISSING WAVE DATA'
            IF(lng.EQ.lng_fr) THEN
              IF(hw%TYPR.NE.'Q') WRITE(lu,*) 'HAUTEUR HM0'
              IF(tw%TYPR.NE.'Q') WRITE(lu,*) 'PERIODE PIC TPR5'
              IF(thetaw%TYPR.NE.'Q') WRITE(lu,*) 'DIRECTION MOY'
            ENDIF
            IF(lng.EQ.lng_en) THEN
              IF(hw%TYPR.NE.'Q') WRITE(lu,*) 'WAVE HEIGHT HM0'
              IF(tw%TYPR.NE.'Q') WRITE(lu,*) 'PEAK PERIOD TPR5'
              IF(thetaw%TYPR.NE.'Q') WRITE(lu,*) 'MEAN DIRECTION'
            ENDIF
            CALL plante(1)
            stop
          ENDIF
        ENDIF
!
! END OF HYDRODYNAMIC INITIALISATION
!
!
!        COMPUTES AREAS (WITHOUT MASKING)
!
        IF(debug.GT.0) WRITE(lu,*) 'VECTOR FOR VOLU2D'
        CALL vector(volu2d,'=','MASBAS          ',
     &              ielmh_sis,1.d0,
     &              t1,t1,t1,t1,t1,t1,mesh,.false.,maskel)
        IF(debug.GT.0) WRITE(lu,*) 'END VECTOR FOR VOLU2D'
!       V2DPAR : LIKE VOLU2D BUT IN PARALLEL VALUES COMPLETED AT
!                INTERFACES BETWEEN SUBDOMAINS
        CALL os('X=Y     ',x=v2dpar,y=volu2d)
        IF(ncsize.GT.1) CALL parcom(v2dpar,2,mesh)
!       INVERSE OF VOLUMES (DONE WITHOUT MASKING)
        CALL os('X=1/Y   ',x=unsv2d,y=v2dpar,
     &          iopt=2,infini=0.d0,zero=1.d-12)
!
!       START OF MODIFICATIONS FOR MIXED SEDIMENTS
!
!       SETTING THE NON-ERODABLE BED (IT CAN BE SET BEFORE
!                                     IF COMPUTATION CONTINUED, I.E. DEBU)
!
        IF(.NOT.debu.OR..NOT.yazr) THEN
          IF(debug.GT.0) WRITE(lu,*) 'NOEROD'
          CALL noerod(hn%R,zf%R,zr%R,z%R,mesh%X%R,
     &                mesh%Y%R,npoin,choix,nliss)
          IF(debug.GT.0) WRITE(lu,*) 'END NOEROD'
        ENDIF
!
!       INITIALISATION FOR SEDIMENT
!
        IF(debug.GT.0) WRITE(lu,*) 'INIT_SEDIMENT'
        CALL init_sediment(nsicla,elay,zf,zr,npoin,
     &                     avail,fracsed_gf,ava0,lgrafed,calwc,
     &                     xmvs,xmve,grav,vce,xwc,fdm,calac,ac,
     &                     sedco,es,es_sable,es_vase,
     &                     ncouch_tass,conc_vase,
     &                     ms_sable, ms_vase,acladm,
     &                     unladm,toce_sable,
     &                     conc,debu,mixte)
        IF(debug.GT.0) WRITE(lu,*) 'END INIT_SEDIMENT'
!
!
! END OF MODIFICATIONS CV
!
!
! MEAN VELOCITY
!======================================================================
        CALL os('X=N(Y,Z)',x=unorm,y=u2d,z=v2d)
! =====================================================================
!  WAVE ORBITAL VELOCITY
! =====================================================================
        IF(houle) THEN
!         JWI 31/05/2012 - added lines to use wave orbital velocities
!         directly if found in hydro file; otherwise compute with CALCUW
          IF(uw%TYPR.NE.'Q') THEN
            CALL calcuw(uw%R,hn%R,hw%R,tw%R,grav,npoin)
          ENDIF
        ENDIF
!
! =====================================================================
!
        IF(debug.GT.0) WRITE(lu,*) 'TOB_SISYPHE'
        CALL tob_sisyphe(tob,tobw,mu,ks,ksp,ksr,cf,fw,
     &                   chestr,uetcar,cf_tel,ks_tel,code ,
     &                   kfrot,icr,kspratio,houle,
     &                   grav,xmve,xmvs,vce,karman,zero,
     &                   hmin,hn,acladm,unorm,uw,tw,npoin,kspred,iks)
        IF(debug.GT.0) WRITE(lu,*) 'END TOB_SISYPHE'
!
!       INITIALISATION FOR TRANSPORT
!
        IF(debug.GT.0) WRITE(lu,*) 'INIT_TRANSPORT'
        CALL init_transport(trouve,debu,hiding,nsicla,npoin,
     &     t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,
     &     charr,qs_c,qsxc,qsyc,calfa_cl,salfa_cl,coefpn,slopeff,
     &     susp,qs_s,qs,qscl,qscl_c,qscl_s,qsclxs,qsclys,
     &     unorm,u2d,v2d,hn,cf,mu,tob,tobw,uw,tw,thetaw,fw,houle,
     &     avail,acladm,unladm,ksp,ksr,ks,
     &     icf,hidfac,xmvs,xmve,grav,vce,hmin,karman,
     &     zero,pi,ac,imp_inflow_c,zref,icq,cstaeq,csratio,
     &     cmax,cs,cs0,seccurrent,bijk,
     &     ielmt,fdm,xwc,fd90,sedco,vitce,partheniades,vitcd,
     &     u3d,v3d,code)
        IF(debug.GT.0) WRITE(lu,*) 'END INIT_TRANSPORT'
!
! ---------- DEBUT IMPRESSIION INITIALISATION =================
!
        CALL entete_sisyphe(1,at0,0)
!       PREPARES RESULTS
!
! CONCENTRATION OUTPUT IN G/L
!
        IF(unit) CALL os('X=CX    ',x=cs,c=xmvs)
        CALL predes(0,at0,.true.,code)
!
!       PRINTS OUT THE RESULTS
!
        IF(debug.GT.0) WRITE(lu,*) 'BIEF_DESIMP'
        CALL bief_desimp(sis_files(sisres)%FMT,varsor,
     &                   npoin,sis_files(sisres)%LU,
     &                   at0,0,lispr,leopr,sorleo,sorimp,maxvar,
     &                   texte,ptinig,ptinil,
     &                   ileo=yagout)
        IF(debug.GT.0) WRITE(lu,*) 'END BIEF_DESIMP'
!
        IF(unit) CALL os('X=CX    ',x=cs,c=1.d0/xmvs)
!
!===============FIN IMPRESSION CONDITIONS INITIALES =================
!
!      COUPLING
!
        IF(nestor) THEN
          leopr = grafcount
          CALL nestor_interface(1)
          WRITE(lu,*) 'SISYPHE COUPLED WITH NESTOR'
        ENDIF
!
        IF(code(1:7).NE.'SISYPHE') THEN
          WRITE(lu,*) 'SISYPHE COUPLED WITH: ',code
        ENDIF
!
!      COUPLING WITH TELEMAC-2D OR 3D
!
        IF(code(1:7).EQ.'TELEMAC') ncalcu = 1
!
!=======================================================================
!
!     INITIAL CONDITION FOR CONSTANT FLOW DISCHARGE
!
        IF(lcondis) THEN
          sisyphe_cfd = lcondis
          nsis_cfd    = ncondis
          constflow   = .false.
        ELSE
          sisyphe_cfd = .false.
          nsis_cfd    = 1
          constflow   = .false.
        ENDIF
!
!------------------------------------------------------------------
!
#if defined COMPAD
      CALL ad_sisyphe_initialisation_end
#endif
!
!=======================================================================
!
!     END OF INITIALISATIONS
        IF(debug.GT.0) WRITE(lu,*) 'END_INITIALIZATION'
      ENDIF ! IF (PART==0 OR PART = -1 OR PART=2)!
!=======================================================================
!
      IF(part==1.OR.part==-1.OR.part==3) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'TIME_LOOP'
!
!=======================================================================
!
! : 3                    /* LOOP ON TIME */
!
!=======================================================================
!
!----   STOPS THE COMPUTATION WHEN THE REQUIRED NUMBER OF ITERATIONS IS 0
!
        IF(nit == 0) THEN
          WRITE(lu,201)
201       FORMAT(' STOP IN SISYPHE, NUMBER OF ITERATIONS EQ.0')
          CALL plante(1)
          stop
        ENDIF
!
!---------------------------------------------------------------------
!       STARTS THE COMPUTATIONS
!---------------------------------------------------------------------
!       LOOP ON THE NUMBER OF EVENTS
!       (IN STEADY STATE: LOOP ON THE TIMESTEPS)
!---------------------------------------------------------------------
!
        IF(code(1:7) == 'TELEMAC') THEN
!         VALNIT WILL BE USED FOR CALLING VALIDA
          valnit = (telnit/pericou)*pericou
!         TO AVOID 2 SUCCESSIVE CALLS TO VALIDA
!         WHEN BEDLOAD AND SUSPENSION
          IF(grafcount.GT.telnit) valnit=nit+1
!         CHARR, SUSP AND TIME STEP MONITORED BY CALLING PROGRAM
          charr = charr_tel
          susp= susp_tel
          at0=t_tel
          dt=mofac*dt
!
        ENDIF
!
        DO mn = 1, ncalcu
!
          IF(part.EQ.3) THEN
            ! API_ITER optional argument is only present when ipart.eq.3
            IF(api_iter.GT.1) GOTO 111
          ENDIF
!------------------------------------------------------------------
!
#if defined COMPAD
      CALL ad_sisyphe_timestep_begin
#endif
!
!------------------------------------------------------------------
!
          IF(.NOT.perma.AND.code(1:7).NE.'TELEMAC') THEN
!           DETERMINES THE FIRST RECORD TO BE READ :
!           NUMDEB IS THE FIRST RECORD TO BE READ FROM THE HYDRO
!           FILE
            numdeb=numen0
            IF(numdeb+nidt > numen) THEN
              WRITE(lu,203)
203           FORMAT(1x,'THE HYDRODYNAMIC FILE IS NOT LONG ENOUGH')
              CALL plante(1)
            ENDIF
          ENDIF
!
!         LOOP ON THE RECORDS (IF PERMA NIDT=1)
!         ------------------------------
111       IF(part.EQ.3.AND.code(1:7).NE.'TELEMAC') THEN
            debut=api_iter
            fin=min(api_iter,nidt) !LIMITED BY HYDRO FILE LAST ITERATION
          ELSE
            debut=1
            fin=nidt
          ENDIF
          DO mt = debut, fin
!          DO MT = 1, NIDT
!
!------------------------------------------------------------------
!
!     ALGORITHMIC DIFFERENTIATION
#if defined COMPAD
      CALL ad_sisyphe_looprecords_begin
#endif
!
!------------------------------------------------------------------
!
!  ----   DETERMINES THE TIMESTEP NUMBER :
!
            IF(part.EQ.3.AND.code(1:7).NE.'TELEMAC') THEN
              lt = api_iter
            ELSE
              lt = (mn-1)*nidt + mt
            ENDIF
!
            IF(code(1:7) == 'TELEMAC') THEN
              dt=dt_tel*mofac
              lt    = loopcount
              leopr = grafcount
              lispr = listcount
              nsous=1
            ENDIF
!
!  ----     PRINTOUTS TO LISTING :
!
            entets = .false.
            IF(lispr*(lt/lispr).EQ.lt) THEN
              entet = .true.
            ELSE
              entet = .false.
            ENDIF
!
!  ----     PRINTOUTS TO C-VSP !!! UHM
!
            IF(cvsmpperiod.EQ.0) cvsmpperiod = leopr
            cvsm_out  = .false.
            IF(cvsmpperiod*(lt/cvsmpperiod).EQ.lt) cvsm_out = .true.
!
!----       READS AND UPDATES H AND ZF
!----       IF 1ST PASS OR UNSTEADY AND NO COUPLING
!
            dts = dt / nsous
!
            isous = 0
!
            IF(.NOT.perma.OR.pass) THEN
!
              IF(debug.GT.0) WRITE(lu,*) 'CONDIM_SISYPHE'
              CALL condim_sisyphe
     &        (u2d%R,v2d%R,qu%R,qv%R,hn%R,zf%R,z%R,esomt%R,thetaw%R,
     &        q%R,hw%R,tw%R,mesh%X%R,mesh%Y%R,npoin,at0,pmaree)
              IF(debug.GT.0) WRITE(lu,*) 'END_CONDIM_SISYPHE'
!
!             NOTE : NAME FOR SISHYD SET TO ' ' IF COUPLING
!
              IF(sis_files(sishyd)%NAME(1:1).NE.' ')  THEN
!
!               WORK ON ZF,QU,QV,Z WILL BE IN FACT DONE ON:
!               T4,DEL_QU,DEL_QV AND DEL_Z
!               BY PLAYING WITH POINTERS
!               THEN THE INCREMENT IN ONE TIME STEP WILL BE COMPUTED
!               AND THE INCREMENT FOR ONE SUB TIME-STEP DEDUCED
!               THEN AT EVERY SUB TIME-STEP THIS INCREMENt WILL BE ADDED
!               AND WILL PROGRESSIVELy UPDATE THE VARIABLES FROM TIME N
!               TO TIME N+1
                savezf=>zf%R
                savequ=>qu%R
                saveqv=>qv%R
                savez =>z%R
                IF(houle.AND.uw%TYPR=='Q') saveuw=>uw%R
                zf%R  =>t4%R
                qu%R  =>del_qu%R
                qv%R  =>del_qv%R
                z%R   =>del_z%R
                IF(houle.AND.uw%TYPR=='Q') uw%R  =>del_uw%R
!
                numdeb=numdeb+1
!
                IF(entet) WRITE(lu,*) 'DEFINITION INITIALE DES VITESSES'
!
                CALL read_dataset(sis_files(sishyd)%FMT,
     &                            sis_files(sishyd)%LU,varsor,npoin,
     &                            numdeb,bid,textpr,trouve,alire,
     &                            entet,perma,maxvar)
!
                IF(debug.GT.0) WRITE(lu,*) 'RESCUE_SISYPHE_NOTPERMA'
                CALL rescue_sisyphe_notperma
     &               (qu%R,qv%R,q%R,u2d%R,v2d%R,hn%R,z%R,t4%R,
     &                hw%R,tw%R,thetaw%R,npoin,trouve,alire,icf,entet,
     &                maxvar)
                IF(debug.GT.0) WRITE(lu,*) 'END_RESCUE_SISYPHE_NOTPERMA'
!
!               BACK TO ORIGINAL ADDRESSES
                zf%R=>savezf
                qu%R=>savequ
                qv%R=>saveqv
                z%R=>savez
                IF(houle.AND.uw%TYPR=='Q') uw%R=>saveuw
!
!               INCREMENT OF QU, QV AND Z PER SUB-TIME-STEP
                DO i = 1,npoin
                  del_qu%R(i) = (del_qu%R(i)-qu%R(i))/nsous
                  del_qv%R(i) = (del_qv%R(i)-qv%R(i))/nsous
                  del_z%R(i)  = (del_z%R(i) -z%R(i)) /nsous
                  IF(houle.AND.uw%TYPR=='Q') THEN
                    del_uw%R(i) = (del_uw%R(i)-uw%R(i))/nsous
                  ENDIF
                ENDDO
!
!               UPDATES UNSTEADY HYDRO
!               (TO BE MOVED TO RESCUE_SISYPHE_NOTPERMA)
!               -----------------------------------
!               CLIPS NEGATIVE DEPTHS
!               COMPUTES U2D AND V2D
!
                CALL os('X=Y-Z   ', x=hn, y=z, z=zf)
!
                IF(optban.GT.0) THEN
                  DO i = 1,npoin
                    IF(hn%R(i).LT.hmin) THEN
                      u2d%R(i)=0.d0
                      v2d%R(i)=0.d0
                      hn%R(i) = max(hn%R(i),hmin)
                    ELSE
                      u2d%R(i)=qu%R(i)/hn%R(i)
                      v2d%R(i)=qv%R(i)/hn%R(i)
                    ENDIF
                  ENDDO
                ELSE
                  CALL os('X=Y/Z   ', x=u2d, y=qu,   z=hn)
                  CALL os('X=Y/Z   ', x=v2d, y=qv,   z=hn)
                ENDIF
!
              ENDIF ! (SIS_FILES(SISHYD)%NAME(1:1) /=' ')
            ENDIF ! (NOT.PERMA.OR.PASS)
!
        IF(pass) THEN
!         IN STEADY STATE LOGICAL FOR READING SET TO FALSE
          IF (perma) pass = .false.
        ELSE
!         COMPUTES THE WATER DEPTH
          CALL os('X=Y-Z   ', x=hn, y=z, z=zf)
        ENDIF
!
!       COUPLING
!
        IF(code(1:7).EQ.'TELEMAC') THEN
!
!         OV INSTEAD OF OS IN ORDER TO AVOID PROBLEMS WITH QUASI-BUBBLE ELEMENTS
!         OPERATES ONLY ON THE (1:NPOIN) RANGE OF THE TELEMAC FIELDS
!         IT IS A HIDDEN DISCRETISATION CHANGE
!
          CALL ov('X=Y     ', x=u2d%R, y=u_tel%R, dim1=npoin)
          CALL ov('X=Y     ', x=v2d%R, y=v_tel%R, dim1=npoin)
          CALL ov('X=Y     ', x=hn%R, y=h_tel%R, dim1=npoin)
!         ZF MAY BE MODIFIED BY CALLING PROGRAM
          CALL os('X=Y     ', x=zf, y=zf_tel)
!         ELAY MUST BE UPDATED CONSEQUENTLY
!         IF COUPLED WITH TELEMAC3D, NO BOTTOM CHANGE FROM TELEMAC2D
          IF(code(1:9).EQ.'TELEMAC3D') THEN
            CALL os('X=Y-Z   ',x=elay,y=zf,z=zr)
          ENDIF
!         CLIPS NEGATIVE DEPTHS
          IF(optban.GT.0) THEN
            DO i = 1,hn%DIM1
              IF(hn%R(i).LT.hmin) THEN
                u2d%R(i)=0.d0
                v2d%R(i)=0.d0
                hn%R(i)=hmin
              ENDIF
            ENDDO
          ENDIF
!         FREE SURFACE
          CALL os('X=Y+Z   ', x=z, y=zf, z=hn)
!
!         COPY OF TOMAWAC VARIABLES
!
          IF(inclus(coupling,'TOMAWAC')) THEN
!           INCIDENT WAVE DIRECTION
            CALL os( 'X=Y     ', x=thetaw, y=thetaw_tel)
!           Wave period
            CALL os( 'X=Y     ', x=tw, y=tw_tel)
!           SIGNIFICANT HEIGHT
            CALL os( 'X=Y     ', x=hw, y=hw_tel)
!           SIGNIFICANT HEIGHT
            CALL os( 'X=Y     ', x=uw, y=uw_tel)
            hw%TYPR='Q'
            tw%TYPR='Q'
            thetaw%TYPR='Q'
            uw%TYPR='Q'
          ENDIF
!
        ENDIF
!
!       END OF COUPLING
! =========================================================================
! TREATMENT OF TIDAL FLATS, DEFINITION OF THE MASKS
! =====================================================================!
!
        IF(optban.EQ.2) THEN
!
! ----    BUILDS MASKING BY ELEMENTS
!
          CALL os ('X=Y     ', x=msktmp, y=maskel)
          CALL os ('X=C     ', x=maskel, c=1.d0)
          IF(code(1:7) == 'TELEMAC') THEN
!           MASKS ARE DERIVED FROM THE NON-CLIPPED VALUES OF H
!           PROVIDED BY TELEMAC
            CALL masktf(maskel%R,h_tel%R,hmin,mesh%IKLE%I,
     &                  nelem,npoin)
          ELSE
            CALL masktf(maskel%R,hn%R,hmin,mesh%IKLE%I,
     &                  nelem,npoin)
          ENDIF
!
        ENDIF
!
! ----   BUILDS THE MASK OF THE POINTS FROM THE MASK OF THE ELEMENTS
! ----   AND CHANGES IFAMAS (IFABOR WITH MASKING)
!
        IF(msk) CALL maskto(maskel%R,maskpt,ifamas%I,
     &                      mesh%IKLE%I,
     &                      mesh%IFABOR%I,mesh%ELTSEG%I,mesh%NSEG,
     &                      nelem,ielmt,mesh)
!
! ------------------------------------------------------------------
!  START OF SUB-ITERATIONS IN UNSTEADY STATE
!
! ------------------------------------------------------------------
!
702     CONTINUE
!
        isous = isous + 1
!
!------------------------------------------------------------------
!
#if defined COMPAD
        CALL ad_sisyphe_subiteration_begin
#endif
!
!------------------------------------------------------------------
!
        IF(code(1:7).NE.'TELEMAC') at0=at0+dts
        IF(entet.AND.isous.EQ.1) CALL entete_sisyphe(2,at0,lt)
        IF(entet.AND.isous.EQ.nsous) entets=.true.
!
!---------------------------------------------------------------------
!        FRICTION COEFFICIENT VARIABLE IN TIME
!---------------------------------------------------------------------
!
        CALL corstr_sisyphe
!
! ----  TREATS THE BOUNDARY CONDITIONS
!
        IF(debug.GT.0) WRITE(lu,*) 'CONLIT'
        CALL conlit(mesh%NBOR%I,at0)
        IF(debug.GT.0) WRITE(lu,*) 'END CONLIT'
!
! =======================================================================
!
!       IF 'VARIABLE TIME-STEP = YES' NSOUS WILL BE COMPUTED FURTHER DOWN
!       THE CONPUTATION OF THE TIMESTEP SIS HAS BEEN MOVED BEFORE READING
!       THE HYDRO CONDITIONS
!
!  ---  MEAN DIAMETER FOR THE ACTIVE-LAYER AND UNDER-LAYER
!
        IF(.NOT.mixte.AND.nsicla.GT.1) CALL mean_grain_size
!
!  ---  MEAN VELOCITY UNORM
!
        CALL os('X=N(Y,Z)',x=unorm,y=u2d,z=v2d)
!
!  ---  WAVE ORBITAL VELOCITY --> UW
!
        IF(houle) THEN
          IF(uw%TYPR.NE.'Q') THEN
            CALL calcuw(uw%R,hn%R,hw%R,tw%R,grav,npoin)
          ENDIF
        ENDIF
!
        CALL tob_sisyphe
     &    (tob,tobw, mu, ks, ksp,ksr,cf, fw,
     &     chestr, uetcar, cf_tel,ks_tel, code ,
     &     kfrot, icr, kspratio,houle,
     &     grav,xmve,  xmvs, vce, karman,zero,
     &     hmin,hn, acladm, unorm,uw, tw, npoin,kspred,iks)
!
!  END OF INITIALISATION
!
      ! ******************** !
      ! BEDLOAD COMPUTATION  !
      ! ******************** !
        IF(charr) THEN
          IF(debug.GT.0) WRITE(lu,*) 'BEDLOAD_MAIN'
          CALL bedload_main
     &        (acladm,ksp,ksr,volu2d,unsv2d,
     &         cf,ebor,fw,hn,liqbor,mask,maskel,
     &         maskpt,qbor,u2d,v2d,s,unladm,uw,thetaw,
     &         mu,tob,tobw,tw,zf,debug,hidfac,icf,
     &         ielmt,kddl,kdir,kent,klog,kneu,ksort,
     &         npoin,nptfr,nsicla,
     &         optban,beta,fd90,fdm,grav,hidi,hmin,
     &         vce,csf_sable,xmve,xmvs,xwc,pi,karman,zero,
     &         karim_holly_yang,msk,susp,vf,
     &         mesh,
     &         elay,liebor,limtec,masktr,
     &         it1,t1,t2,t3,t4,t5,t6,t7,t8,t9,
     &         t10,t11,t12,t13,unorm,ac,dts,
     &         avail,breach,calfa_cl,coefpn,hiding,
     &         qscl_c,qscl_s,qs_c,qsclxc,qsxc,qsclyc,
     &         qsyc,salfa_cl,zf_c,zfcl_c,entets,
     &         seccurrent,slopeff,phised,devia,beta2,bijk,
     &         sedco,houle,u3d,v3d,code,flbcla,maxadv)
!
          IF(debug.GT.0) WRITE(lu,*) 'END_BEDLOAD_MAIN'
!
!         UPDATES THE BOTTOM
!
          IF(.NOT.stat_mode) THEN
            CALL os('X=X+Y   ',x=zf,y=zf_c)
          ENDIF
!
!         UPDATES THE LAYERS  --> ELAY
!
          IF(.NOT.mixte.AND.nsicla.GT.1) THEN
            IF(debug.GT.0) WRITE(lu,*) 'LAYER AFTER BEDLOAD'
            IF(vsmtype.EQ.0) THEN
              CALL layer(zfcl_c,nlayer,zr,zf,estrat,elay,volu2d,
     &                   acladm,nsicla,npoin,elay0,voltot,es,
     &                   avail,const_alayer,t2%R,it1%I)
            ELSE
              CALL cvsp_main(zfcl_c,zf,nsicla,npoin)
            ENDIF
            IF(debug.GT.0) WRITE(lu,*) 'END_LAYER'
          ELSE
            CALL os('X=Y-Z   ',x=elay,y=zf,z=zr)
          ENDIF
! END OF BEDLOAD
        ENDIF
      ! ********************** !
      ! SUSPENSION COMPUTATION !
      ! ********************** !
        IF(susp) THEN
!
          IF(debug.GT.0) WRITE(lu,*) 'SUSPENSION_MAIN'
          CALL suspension_main
     &(slvtra,hn,hn_tel,mu,tob,fdm,fd90,ksp,ksr,ks,
     & volu2d,v2dpar,unsv2d,afbor,bfbor,zf,licbor,
     & ifamas,maskel,maskpt,u2d,v2d,nsicla,
     & npoin,nptfr,ielmt,optdif,resol,lt,nit,optban,optadv,
     & opdtra,kent,ksort,klog,kneu,kdir,kddl,
     & debug,dts,csf_sable,zero,grav,xkx,xky,
     & karman,xmve,xmvs,vce,hmin,xwc,vitcd,partheniades,
     & bilma,msk,charr,imp_inflow_c,mesh,zf_s,cs,
     & cst,ctild,cbor,disp,it1,it2,it3,it4,tb,t1,t2,t3,t4,
     & t8,t9,t10,t11,t12,t14,te1,clt,te2,te3,s,am1_s,am2_s,mbor,
     & elay,limdif,masktr,teta_susp,ac,mased0,masini,masten,mastou,
     & es,es_sable, es_vase,avail,entets,pass_susp,
     & zfcl_s,hprop,fludpt,fludp,fluer,disp_c,kx,ky,kz,uconv,
     & vconv,qsxs,qsys,qsclxs,qsclys,qscl_s,qs_s,qs_c,
     & cstaeq,csratio,icq,mastcp,masfin,masdept,masdep,massou,corr_conv,
     & zref,sedco,visc_tel,code,dift,dm1,uconv_tel,vconv_tel,
     & zconv,solsys,flbor_tel,flbor_sis,flbortra,numliq%I,nfrliq,
     & mixte,ncouch_tass,conc,toce_vase,toce_sable,
     & fluer_vase,toce_mixte,ms_sable,ms_vase,dirflu,maxadv)
          IF(debug.GT.0) WRITE(lu,*) 'END_SUSPENSION_MAIN'
!
!      UPDATES THE BOTTOM
!
! mak:
          IF(.NOT.stat_mode) THEN
            CALL os('X=X+Y   ',x=zf,y=zf_s)
          ENDIF
! end mak
!        CALL OS('X=X+Y   ',X=ZF,Y=ZF_S)
!
!      UPDATES THE LAYERS
!      REDEFINES THE LAYER OF ERODABLE SEDIMENT
!      EXTENDED GRANULOMETRY (TO BE REPLACED WITH NOMBLAY>1
!
          IF(.NOT.mixte.AND.nsicla.GT.1) THEN
            IF(debug.GT.0) WRITE(lu,*) 'LAYER AFTER SUSPENSION'
            IF(vsmtype.EQ.0) THEN
              CALL layer(zfcl_s,nlayer,zr,zf,estrat,elay,volu2d,
     &                   acladm,nsicla,npoin,elay0,voltot,es,
     &                   avail,const_alayer,t2%R,it1%I)
            ELSE
              CALL cvsp_main(zfcl_s,zf,nsicla,npoin)
            ENDIF
            IF(debug.GT.0) WRITE(lu,*) 'END_LAYER'
          ELSE
            CALL os('X=Y-Z   ',x=elay,y=zf,z=zr)
          ENDIF
! END OF SUSPENSION
        ENDIF
!
! RECONSTITUTES THE BEDLOAD AND/OR SUSPENSION DATA
! -----------------------------------------------------
!
        IF(debug.GT.0) WRITE(lu,*) 'QS_RESULT'
!
        DO i = 1, nsicla
          CALL os('X=Y+Z   ',x=t1,y=qsclxc%ADR(i)%P,z=qsclxs%ADR(i)%P)
          CALL os('X=Y+Z   ',x=t2,y=qsclyc%ADR(i)%P,z=qsclys%ADR(i)%P)
          CALL os('X=N(Y,Z)', x=qscl%ADR(i)%P,y=t1,z=t2)
          IF(i.EQ.1) THEN
            CALL os('X=Y     ', x=qsx, y=t1)
            CALL os('X=Y     ', x=qsy, y=t2)
          ELSE
            CALL os('X=X+Y   ', x=qsx, y=t1)
            CALL os('X=X+Y   ', x=qsy, y=t2)
          ENDIF
        ENDDO
        CALL os('X=N(Y,Z)', x=qs, y=qsx, z=qsy)
        IF(debug.GT.0) WRITE(lu,*) 'END_QS_RESULT'
!
!=======================================================================
!
!     MAXIMUM BOTTOM SLOPE : EVOL IN T1
!
      IF(slide) THEN
!
        IF(entet) CALL entete_sisyphe(14,at0,lt)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE MAXSLOPE'
        CALL maxslope(phised,zf%R,zr%R,mesh%XEL%R,mesh%YEL%R,
     &                mesh%NELEM,
     &                mesh%NELMAX,npoin,mesh%IKLE%I,t1,unsv2d,mesh,
     &                zfcl_ms,avail,nomblay,nsicla)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE MAXSLOPE'
        CALL os('X=X+Y   ',x=zf,y=t1)
!
        IF(.NOT.mixte.AND.nsicla.GT.1) THEN
          IF(debug.GT.0) WRITE(lu,*) 'LAYER AFTER SLIDE'
          IF(vsmtype.EQ.0) THEN
            CALL layer(zfcl_ms,nlayer,zr,zf,estrat,elay,volu2d,
     &                 acladm,nsicla,npoin,elay0,voltot,es,
     &                 avail,const_alayer,t2%R,it1%I)
          ELSE
            CALL cvsp_main(zfcl_ms,zf,nsicla,npoin)
          ENDIF
          IF(debug.GT.0) WRITE(lu,*) 'END_LAYER'
        ELSE
          CALL os('X=Y-Z   ',x=elay,y=zf,z=zr)
        ENDIF
!
      ENDIF
!
!========================================================================
!
!     SETTLING: EVOLUTION COMPUTED IN T3
!
      IF(tass) THEN
!
        IF(entet) THEN
          IF(.NOT.charr.AND..NOT.susp.AND..NOT.slide) THEN
            CALL entete_sisyphe(2,at0,lt)
          ENDIF
          CALL entete_sisyphe(15,at0,lt)
        ENDIF
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TASSEMENT'
        IF(itass.EQ.1) THEN
          CALL tassement(npoin,dts,elay,t3,t2,avail,nsicla,
     &                 es,xmvs,xkv,trans_mass,conc_vase,ncouch_tass,
     &                 ms_sable%R,ms_vase%R)
        ELSEIF(itass.EQ.2.OR.itass.EQ.3) THEN
!!!! ONLY FOR ONE CLASS
          CALL tassement_2(npoin,dts,elay,t3,t2,lt,xmvs,xmve,
     &        nomblay,es_vase,conc_vase,ms_vase%R,xwc(1),
     &        coef_n,conc_gel,conc_max)
!        ELSEIF(ITASS.EQ.3) THEN
!!!! ONLY FOR ONE CLASS
!          CALL TASSEMENT_3(ZF,NPOIN,DTS,ELAY,
!     &               T3,T2,LT,XMVS,XMVE,GRAV,NOMBLAY,
!     &               ES,CONC_VASE,CONC,IVIDE,MS_VASE%R,XWC(1),
!     &                 TRA01,TRA02,TRA03,CONC_GEL,COEF_N,CONC_MAX)
        ENDIF
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TASSEMENT'
!
!       UPDATES ZF (ELAY HAS BEEN UPDATED IN TASSEMENT)
!
        IF(.NOT.stat_mode) THEN
          CALL os('X=X+Y   ',x=zf,y=t3)
        ENDIF
!
      ENDIF
!
!=======================================================================
! : 5        COMPUTES THE EVOLUTIONS FOR THIS CYCLE OF TIMESTEP
!            AND UPDATES AFTER THIS COMPUTATION
!=======================================================================
!
! ----  COMPUTES  THE EVOLUTIONS FOR THIS (SUB) TIMESTEP
!
      IF(charr) THEN
        CALL os('X=Y     ',x=e,y=zf_c)
      ELSE
        CALL os('X=0     ',x=e)
      ENDIF
      IF(susp)  CALL os('X=X+Y   ',x=e,y=zf_s)
      IF(slide) CALL os('X=X+Y   ',x=e,y=t1)
      IF(tass)  CALL os('X=X+Y   ',x=e,y=t3)
!
      CALL os('X=X+Y   ', x=esomt, y=e)
!
!  UPDATES
!
      IF(part.EQ.-1.OR.(part.EQ.3.AND.code(1:7).NE.'TELEMAC')) THEN
!
! mak:
        IF(.NOT.stat_mode) THEN
          CALL os('X=X-Y   ',x=hn,y=e)
        ENDIF
! end mak
!      CALL OS('X=X-Y   ',X=HN,Y=E)
        IF(optban.GT.0) THEN
          DO i = 1,hn%DIM1
            IF(hn%R(i).LT.hmin) THEN
              u2d%R(i)=0.d0
              v2d%R(i)=0.d0
              hn%R(i) =hmin
            ELSE
              u2d%R(i)= qu%R(i)/hn%R(i)
              v2d%R(i)= qv%R(i)/hn%R(i)
            ENDIF
          ENDDO
        ELSE
          CALL os('X=Y/Z   ', x=u2d, y=qu,   z=hn)
          CALL os('X=Y/Z   ', x=v2d, y=qv,   z=hn)
        ENDIF
!
!=======================================================================
! : 6     STOPS IF EVOLUTIONS GREATER THAN EMAX = RC*(INITIAL DEPTH)
!=======================================================================
!
!       DETERMINES THE MAXIMUM EVOLUTION THRESHOLD
        DO i = 1, npoin
          emax%R(i) = rc*max(hn%R(i),hmin)
        ENDDO
!
! ----  STOPS WHEN THE EVOLUTIONS ARE GREATER THAN A CERTAIN THRESHOLD
!      THIS TEST IS ONLY CALLED IN 'SISYPHE ONLY' MODE
!
        IF(debug.GT.0) WRITE(lu,*) 'ARRET'
          CALL sis_arret(esomt,emax,hn,varsor,npoin,mn,
     &                   sis_files(sisres)%LU,sis_files(sisres)%FMT,
     &                   maxvar,at0,rc,texte,
     &                   sorleo,sorimp,t1,t2)
          IF(debug.GT.0) WRITE(lu,*) 'END_ARRET'
        ENDIF
!
! ----     CONSTANT FLOW DISCHARGE
!
        IF(lcondis) THEN
          CALL condis_sisyphe(constflow)
        ELSE
          constflow =.false.
        ENDIF
!
!=======================================================================
! : 8     MASS BALANCE
!=======================================================================
!       COMPUTES THE COMPONENTS OF SAND TRANSPORT FOR THE MASS BALANCE,
!       GRAPHIC OUTPUTS AND VALIDATION STAGE
!
        IF(bilma) THEN
! CVL Bilan in volume
! not valid for cohesif or mixte
          IF(.NOT.mixte.AND.(.NOT.sedco(1))) THEN
          IF(debug.GT.0) WRITE(lu,*) 'BILAN_SISYPHE'
          CALL bilan_sisyphe(e,esomt,t1,
     &                       vcumu,dts,nptfr,entets,
     &                       zfcl_c,zfcl_s,zfcl_ms,
     &                       nsicla,voltot,
     &                       numliq%I,nfrliq,flbcla,lt,nit,
     &                       npoin,volu2d,csf_sable,masdep,masdept,
     &                       charr,susp,slide)
          ENDIF
          IF(debug.GT.0) WRITE(lu,*) 'END_BILAN_SISYPHE'
        ENDIF
!
!       CONTROL SECTIONS
!
        IF(ncp.GT.0) THEN
          IF(debug.GT.0) WRITE(lu,*) 'FLUSEC_SISYPHE'
          CALL flusec_sisyphe(u2d,v2d,hn,
     &                        qsxc,qsyc,charr,qsxs,qsys,susp,
     &                        mesh%IKLE%I,
     &                        mesh%NELMAX,mesh%NELEM,
     &                        mesh%X%R,mesh%Y%R,
     &                        dts,ncp,ctrlsc,entets,at0)
          IF(debug.GT.0) WRITE(lu,*) 'END_FLUSEC_SISYPHE'
        ENDIF
!
!-----------------------------------------------------------------------
!
        IF(.NOT.perma.AND.sis_files(sishyd)%NAME(1:1).NE.' ') THEN
!
!         UPDATES THE HYDRO
!
!         IF READING ON HYDRODYNAMIC FILE, INCREMENTS QU, QV AND Z
          IF(sis_files(sishyd)%NAME(1:1).NE.' ')  THEN
            CALL os('X=X+Y   ', x=qu, y=del_qu)
            CALL os('X=X+Y   ', x=qv, y=del_qv)
            CALL os('X=X+Y   ', x=z , y=del_z)
! JWI 31/05/2012 - added line to use wave orbital velocities directly if found in hydro file
            IF(houle.AND.uw%TYPR=='Q') CALL os('X=X+Y   ',x=uw,y=del_uw)
! JWI END
          ENDIF
          CALL os('X=Y-Z   ', x=hn, y=z, z=zf)
!         CLIPS NEGATIVE DEPTHS
          IF(optban.GT.0) THEN
            DO i = 1, npoin
              IF(hn%R(i).LT.hmin) THEN
                u2d%R(i)=0.d0
                v2d%R(i)=0.d0
                hn%R(i) = max(hn%R(i),hmin)
              ELSE
                u2d%R(i)= qu%R(i)/hn%R(i)
                v2d%R(i)= qv%R(i)/hn%R(i)
              ENDIF
            ENDDO
          ELSE
            CALL os('X=Y/Z   ', x=u2d, y=qu,   z=hn)
            CALL os('X=Y/Z   ', x=v2d, y=qv,   z=hn)
          ENDIF
!
        ENDIF
!
!------------------------------------------------------------------
!
#if defined COMPAD
        CALL ad_sisyphe_subiteration_end
#endif
!
!------------------------------------------------------------------
!
!       END OF THE LOOP ON SUB-TIMESTEPS NSOUS
!
        IF(isous.LT.nsous) GO TO 702
!
!=======================================================================
! : 9        PRINTS OUT EXTREME VALUES
!=======================================================================
!
        IF(entet.AND.charr) THEN
          WRITE(lu,*)
          CALL maxi(xmax,imax,e%R,npoin)
          IF(ncsize.GT.1) THEN
            xma=p_max(xmax)
            IF(xmax.EQ.xma) THEN
              ima=mesh%KNOLG%I(imax)
            ELSE
              ima=0
            ENDIF
            ima=p_max(ima)
          ELSE
            ima=imax
            xma=xmax
          ENDIF
          WRITE(lu,372) xma,ima
372       FORMAT(' MAXIMAL EVOLUTION        : ',g16.7,' NODE  : ',i6)
          CALL mini(xmin,imin,e%R,npoin)
          IF(ncsize.GT.1) THEN
            xmi=p_min(xmin)
            IF(xmin.EQ.xmi) THEN
              imi=mesh%KNOLG%I(imin)
            ELSE
              imi=0
            ENDIF
            imi=p_max(imi)
          ELSE
            imi=imin
            xmi=xmin
          ENDIF
          WRITE(lu,374) xmi,imi
374       FORMAT(' MINIMAL EVOLUTION        : ',g16.7,' NODE  : ',i6)
!
          IF(const_alayer) THEN
            IF(nsicla.GT.1.AND.xmi.LT.-0.5d0*elay0) THEN
              WRITE(lu,886)
886           FORMAT(' EROSION GREATER THAN ONE LAYER THICKNESS !')
            ENDIF
            IF(nsicla.GT.1.AND.xma.GT.elay0) THEN
              WRITE(lu,888)
888           FORMAT(' DEPOSITION MORE THAN ONE LAYER THICKNESS !')
            ENDIF
          ELSE
            DO j=1,npoin
              elay0 = 3.d0*acladm%R(j)
              IF(nsicla.GT.1.AND.e%R(j).LT.-0.5d0*elay0) THEN
                WRITE(lu,886)
              ENDIF
              IF(nsicla.GT.1.AND.e%R(j).GT.elay0) THEN
                WRITE(lu,888)
              ENDIF
            ENDDO
          ENDIF
        ENDIF
        IF(entet) THEN
          CALL maxi(xmax,imax,esomt%R,npoin)
          IF(ncsize.GT.1) THEN
            xma=p_max(xmax)
            IF(xmax.EQ.xma) THEN
              ima=mesh%KNOLG%I(imax)
            ELSE
              ima=0
            ENDIF
            ima=p_max(ima)
          ELSE
            ima=imax
            xma=xmax
          ENDIF
          WRITE(lu,882) xma,ima
882       FORMAT(' TOTAL MAXIMAL EVOLUTION  : ',g16.7,' NODE  : ',i6)
          CALL mini(xmin,imin,esomt%R,npoin)
          IF(ncsize.GT.1) THEN
            xmi=p_min(xmin)
            IF(xmin.EQ.xmi) THEN
              imi=mesh%KNOLG%I(imin)
            ELSE
              imi=0
            ENDIF
            imi=p_max(imi)
          ELSE
            imi=imin
            xmi=xmin
          ENDIF
          WRITE(lu,884) xmi,imi
884       FORMAT(' TOTAL MINIMAL EVOLUTION  : ',g16.7,' NODE  : ',i6)
        ENDIF
!
!=======================================================================
! : 10         PRINTS OUT RESULTS AT THIS TIMESTEP
!              AND COMPARES AGAINST A REFERENCE FILE
!=======================================================================
!
        IF(unit) CALL os('X=CX    ', x=cs, c= xmvs)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PREDES'
        CALL predes(lt,at0,yagout,code)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PREDES'
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BIEF_DESIMP'
        IF(PRESENT(grcomp))THEN
          CALL bief_desimp(sis_files(sisres)%FMT,varsor,
     &                     npoin,sis_files(sisres)%LU,
     &                     at0,lt,lispr,leopr,
     &                     sorleo,sorimp,maxvar,texte,ptinig,ptinil,
     &                     ileo=yagout,compgraph=grcomp+1)
        ELSE
          CALL bief_desimp(sis_files(sisres)%FMT,varsor,
     &                     npoin,sis_files(sisres)%LU,
     &                     at0,lt,lispr,leopr,
     &                     sorleo,sorimp,maxvar,texte,ptinig,ptinil)
        ENDIF
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIEF_DESIMP'
!
        IF(unit) CALL os('X=CX    ', x=cs,  c= 1.d0/xmvs)
!
!       SENDS THE NEW ZF TO TELEMAC-2D OR 3D
!
        IF(code(1:7).EQ.'TELEMAC') THEN
          CALL ov ('X=Y     ', x=zf_tel%R, y=zf%R, dim1=npoin)
        ENDIF
!
!       THE SUBROUTINE VALIDA FROM THE LIBRARY IS STANDARD
!       IT CAN BE MODIFIED FOR EACH PARTICULAR CASE
!       BUT ITS CALL MUST BE LEFT IN THE LOOP ON TIME
!
!------------------------------------------------------------------
!
#if defined COMPAD
          CALL ad_sisyphe_looprecords_end
#endif
!
!------------------------------------------------------------------
!
        IF(valid) THEN
          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BIEF_VALIDA'
          CALL bief_valida(tb,textpr,sis_files(sisref)%LU,
     &                     sis_files(sisref)%FMT,
     &                     varsor,texte,sis_files(sisres)%LU,
     &                     sis_files(sisres)%FMT,
     &                     maxvar,npoin,lt,valnit,alirv)
          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIEF_VALIDA'
        ENDIF
!
!       END OF THE LOOP ON THE RECORDS
        ENDDO !MT
!
!------------------------------------------------------------------
!
#if defined COMPAD
          CALL ad_sisyphe_timestep_end
#endif
!
!=======================================================================
!
!       END OF THE LOOP ON THE NUMBER OF EVENTS
!
        ENDDO !MN
!
!-----------------------------------------------------------------------
!
        IF(debug.GT.0) WRITE(lu,*) 'END_TIME_LOOP'
      ENDIF
!
!-----------------------------------------------------------------------
!
#if defined COMPAD
      CALL ad_sisyphe_end
#endif
!
!-----------------------------------------------------------------------
!
      RETURN
      END