!                   *************************
                    SUBROUTINE TELEMAC3D_INIT
!                   *************************
!
!***********************************************************************
! TELEMAC3D
!***********************************************************************
!
! Pour permettre la màj bathy automatique en cas de reprise de calcul.
! T3D_CORFON est appelée plus tard, lors du traitement de la reprise
! de calcul 3D.
!
!***********************************************************************
!
!brief Initialisation step of telemac3d
!
!history  Y. AUDOUIN (LNHE)
!+        23/10/2017
!+        V8P0
!+   Creating of the file
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE BIEF
      USE DECLARATIONS_TELEMAC
      USE DECLARATIONS_TELEMAC3D
      USE DECLARATIONS_GAIA, ONLY:XMVS0,NSUSP_TEL,NSICLA,NUM_ISUSP_ICLA
      USE METEO_TELEMAC !, ONLY: WINDX,WINDY,PATMOS,TAIR
      USE OILSPILL
      USE INTERFACE_TELEMAC3D, EX_TELEMAC3D_INIT => TELEMAC3D_INIT
      USE INTERFACE_TELEMAC2D
      USE INTERFACE_SISYPHE, ONLY: SISYPHE
      USE INTERFACE_GAIA, ONLY: GAIA_INIT
      USE INTERFACE_TOMAWAC, ONLY: WAC
      USE TEL4DEL, ONLY: TEL4DELWAQ
      USE DECLARATIONS_WAQTEL,ONLY: WAQPROCESS,RAYAED2
      USE EXCHANGE_WITH_ATMOSPHERE, ONLY: SOLRAD
      USE GOTM_COUPLING
!
      USE DECLARATIONS_SPECIAL
      USE INTERFACE_PARALLEL, ONLY : P_MAX
      IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
      INTEGER DATE(3),TIME(3)
      INTEGER :: ISOUSI
      LOGICAL :: TRAC
      INTEGER TROUVE(MAXVAR+10)
      INTEGER IBID,I,II,K,NFRLIQ0
      CHARACTER(LEN=24), PARAMETER :: CODE1='TELEMAC3D               '
      INTEGER ITRAC,IERR
      LOGICAL LBID, LBID2
      DOUBLE PRECISION RAY_SOL,NEBU
#if defined COMPAD
      CALL AD_TELEMAC3D_INITIALISATION_BEGIN
#endif
!
      LT     = 0       ! INITIALISES TIMESTEP
!     INITIALISES NUMBER OF SUB-ITERATIONS, LOOK IN PREADV
      ISOUSI = 0
      IF(NTRAC.GT.0) THEN
        TRAC=.TRUE.
      ELSE
        TRAC=.FALSE.
      ENDIF
!
!     DATE AND TIME (NOW TAKEN FROM KEYWORDS MARTIM AND MARDAT)
!
      DATE(1) = MARDAT(1)
      DATE(2) = MARDAT(2)
      DATE(3) = MARDAT(3)
      TIME(1) = MARTIM(1)
      TIME(2) = MARTIM(2)
      TIME(3) = MARTIM(3)
!
      INFOGR = LISTIN
!
!-----------------------------------------------------------------------
!
! 2D BOUNDARY CONDITIONS:
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE LECLIM'
      CALL LECLIM
     & (LIHBOR%I,LIUBOL%I,LIVBOL%I,IT4%I,HBOR%R,UBOR2D%R,VBOR2D%R,
     &  T2_01%R,T2_02%R,T2_03%R,T2_04%R,NPTFR2,'T3D',.FALSE.,
     &  T3D_FILES(T3DGEO)%FMT,T3D_FILES(T3DGEO)%LU,
     &  KENT,KENTU,KSORT,KADH,KLOG,KINC,NUMLIQ%I,MESH2D)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE LECLIM'
!
!-----------------------------------------------------------------------
!
! STORE HBOR, UBOR2D, VBOR2D IN CASE OF COUPLING WITH TOMAWACT3D
!
      IF(INCLUS(COUPLING,'TOMAWACT3D')) THEN
        DO K=1,NPTFR2
          HBORCLI%R(K) = HBOR%R(K)
          UBORCLI%R(K) = UBOR2D%R(K)
          VBORCLI%R(K) = VBOR2D%R(K)
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
! MESH ORGANISATION - 2D LEVEL
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE INBIEF POUR MESH2D'
      CALL INBIEF(LIHBOR%I,KLOG,IT1,IT2,IT3,
     &            LVMAC,IELMX,LATIT,SPHERI,MESH2D,
     &            T2_01,T2_02,OPTASS2D,PRODUC,EQUA)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE INBIEF'
!
! CORRECTS THE NORMAL VECTORS AT THE POINTS
! WHERE LIQUID AND SOLID BOUNDARIES MEET
!
      CALL CORNOR(MESH2D%XNEBOR%R,MESH2D%YNEBOR%R,
     &            MESH2D%XSGBOR%R,MESH2D%YSGBOR%R,
     &            NPTFR2,KLOG,LIHBOR%I,
     &            T2_01,T2_02,MESH2D,MESH2D%IKLBOR%I,
     &            MESH2D%NELEB,MESH2D%NELEBX)
!
! 3D BOUNDARY CONDITIONS (SO FAR SAME FILE AS 2D)
! T2_02 IS AUBOR IN T2D, COULD BE KEPT
! THIS TIME BOUNDARY COLOURS ARE READ
!
      CALL LECLIM
     & (LIHBOR%I,LIUBOL%I,LIVBOL%I,LITABL%ADR(1)%P%I,
     &  HBOR%R,UBORL%R,VBORL%R,TABORL%ADR(1)%P%R,AUBORL%R,
     &  ATABOL%ADR(1)%P%R,BTABOL%ADR(1)%P%R,
     &  NPTFR2,'T3D',TRAC,T3D_FILES(T3DGEO)%FMT,T3D_FILES(T3DGEO)%LU,
     &  KENT,KENTU,KSORT,KADH,KLOG,KINC,
     &  NUMLIQ%I,MESH3D,BOUNDARY_COLOUR%I)
!
! MESH ORGANISATION - 3D LEVEL
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE INBIEF POUR MESH3D'
      CALL INBIEF(LIHBOR%I,KLOG,IT1,IT2,IT3,
     &            LVMAC,IELM3,LATIT,SPHERI,MESH3D,
     &            T3_01,T3_02,OPTASS,PRODUC,EQUA,MESH2D)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE INBIEF'
!
! INITIALISES 3D BOUNDARY CONDITION ATTRIBUTES FOR BOUNDARY NODES
! DUPLICATES 2D CONDITIONS ON THE VERTICAL
!
      CALL LIMI3D
!
! COMPLETES IFABOR IN 3D
!
      IF(IELM3.EQ.41) THEN
        CALL IFAB3D
     &   (MESH3D%IFABOR%I,LIUBOF%I,LIUBOL%I,LIUBOS%I,
     &    MESH2D%KP1BOR%I,MESH2D%NELBOR%I,
     &    MESH2D%NULONE%I,IKLE2%I,
     &    NELEM2,NELMAX2,NPOIN2,NPTFR2,NPLAN,NETAGE,
     &    KLOG,TRANSF)
      ELSEIF(IELM3.EQ.51) THEN
        CALL IFAB3DT
     &   (MESH3D%IFABOR%I,
     &    LIUBOF%I,LIUBOS%I,
     &    IKLE2%I,IKLE3%I,
     &    NELEM2,NELMAX2,MESH3D%NELMAX,NPOIN2,NETAGE,
     &    KLOG)
      ENDIF
!
! CONTROLS MESH
!
      CALL CHECK(IKLE2%I,NBOR2%I,MESH2D%NELBOR%I,MESH3D%IKLBOR%I,
     &           MESH3D%NELEB,MESH3D%NELEBX,
     &           IKLE3%I,MESH3D%NELBOR%I,MESH3D%NULONE%I,
     &           MESH3D%NULONE%DIM1,MESH3D%NULONE%DIM2,NBOR3%I,
     &           NELEM2,NPTFR2,MESH3D%NELMAX,NPTFR3,LISTIN)
!
! INITIALISES PRIVATE VECTOR BLOCK
!
      IF(NPRIV.GT.0)   CALL OS('X=0     ',X=PRIVE)
      IF(NPRIV2D.GT.0) CALL OS('X=0     ',X=PRIVE2D)
!
! LOOKS FOR THE BOTTOM AND BOTTOM FRICTION VARIABLES IN THE GEOMETRY FILE
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE FONSTR'
      CALL FONSTR(T2_01,ZF,T2_02,RUGOF,T3D_FILES(T3DGEO)%LU,
     &            T3D_FILES(T3DGEO)%FMT,
     &            T3D_FILES(T3DFON)%LU,T3D_FILES(T3DFON)%NAME,
     &            MESH2D,RUGOF0,LISTIN,
     &            N_NAMES_PRIV2D,NAMES_PRIVE2D,PRIVE2D)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE FONSTR'
      I=11
      IF(RUGOF%ELM.NE.I) CALL CHGDIS(RUGOF,I,RUGOF%ELM,MESH2D)
!
!*** Modif EPO *********************************************************
!
! CORRECTS THE BOTTOM
!
!  - SMOOTHES ACCORDING TO THE LISFON VALUE
!  - CHANGES THE BOTTOM TOPOGRAPHY (FORTRAN)
!
!      IF(DEBUG.GT.0) WRITE(LU,*) 'CALLING T3D_CORFON'
!      CALL T3D_CORFON(ZF,T2_01,T2_02,ZF%R,T2_01%R,T2_02%R,
!     &                X,Y,PRIVE,NPOIN2,LISFON,LISFON_AFTER,.FALSE.,
!     &                MASKEL,MATR2H,MESH2D,SVIDE)
!      IF(DEBUG.GT.0) WRITE(LU,*) 'RETURN FROM T3D_CORFON'	  
!
! FINELY ANALYSES THE BATHYMETRY
! IN T2D CALLED IF (OPTBAN == 2)
!
!      IF(MSK) CALL TOPOGR(ZF%R,T2_01%R,ZFE%R,IKLE2%I,MESH2D%IFABOR%I,
!     &  MESH2D%NBOR%I, MESH2D%NELBOR%I, MESH2D%NULONE%I,
!     &  IT1%I, IT2%I, IT3%I, NELEM2, NPTFR2, NPOIN2, MXPTVS2)
!
!*** Fin Modif EPO *****************************************************
!
!=======================================================================
! VARIOUS INITIALISATIONS
!=======================================================================
!
!     COUNTS THE LIQUID BOUNDARIES
!
      NFRLIQ0=NFRLIQ
      IF(NCSIZE.GT.1) THEN
        NFRLIQ=0
        DO I=1,NPTFR2
          NFRLIQ=MAX(NFRLIQ,NUMLIQ%I(I))
        ENDDO
        NFRLIQ=P_MAX(NFRLIQ)
        WRITE(LU,*) ' '
        WRITE(LU,*) 'NUMBER OF LIQUID BOUNDARIES:',NFRLIQ
        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,
     &              LIHBOR%I,LIUBOL%I,X,Y,NBOR2%I,MESH2D%KP1BOR%I,
     &              IT1%I,NPOIN2,NPTFR2,KLOG,LISTIN,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
      IF(NFRLIQ0.NE.0.AND.NFRLIQ0.NE.NFRLIQ) THEN
        WRITE(LU,*) ' '
        WRITE(LU,*) 'OPTION FOR TIDAL BOUNDARY CONDITIONS:'
        WRITE(LU,*) 'GIVE THE SAME NUMBER OF VALUES AS THE NUMBER '
        WRITE(LU,*) 'OF LIQUID BOUNDARIES, I.E.:',NFRLIQ
        CALL PLANTE(1)
        STOP
      ENDIF
!
!     3D EXTENSION OF NUMLIQ
!
      CALL NUMLIQ_3D(NUMLIQ%I,NUMLIQ_ELM,NPLAN,MESH2D%NPTFR,
     &               MESH2D%IKLBOR%I,MESH2D%NELEB,MESH2D%NELEBX)
!
!     SAVING AND MODIFYING BOUNDARY CONDITIONS FOR THOMPSON METHOD
!
      IF(THOMFR.AND.NFRLIQ.GT.0) THEN
        CALL THOMPS_BC(1)
        CALL THOMPS_BC(2)
      ENDIF
!
!=======================================================================
!
!     READS THE STAGE-DISCHARGE CURVES FILE
!
      IF(T3D_FILES(T3DPAR)%NAME(1:1).NE.' ') THEN
        CALL T3D_READ_FIC_CURVES(T3D_FILES(T3DPAR)%LU,NFRLIQ,
     &                           STA_DIS_CURVES,PTS_CURVES)
      ENDIF
!
!     SETS TURBULENCE CONSTANTS (ALL MODELS)
!
      CALL CSTSA
      CALL CSTKEP
!
!-----------------------------------------------------------------------
! READS INITIAL CONDITIONS FROM A PREVIOUS 3D COMPUTATION FILE
! OR SETS THEM IN FORTRAN
!
      AKEP = .TRUE.
      AKOM = .TRUE.
      AKSA = .TRUE.
!
!     STARTS FROM A 2D FILE (U,V AND H ARE READ TO BE USED IN CONDIM)
!                            AT IS ALSO INITIALISED
!     IF A COMPUTATION CONTINUED IS ALSO ASKED, WHAT IS READ HERE MAY
!     NOT BE USED, BUT AS CONDIM IS ALWAYS CALLED, WILL AVOID A CRASH
!     IN CONDIM.
!
      IF(SUIT2) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE SUITE AVEC UN FICHIER 2D'
          CALL READ_DATASET(T3D_FILES(T3DS2D)%FMT,T3D_FILES(T3DS2D)%LU,
     &                   VARSOR,MESH2D%NPOIN,IBID,AT,TEXTPR,TROUVE,
     &                   ALIRE2D,LISTIN,.TRUE.,MAXVAR)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE SUITE'
        IF(TROUVE(1).NE.1.OR.TROUVE(2).NE.1.OR.
     &     TROUVE(4).NE.1) THEN
          WRITE(LU,*)
          WRITE(LU,*) 'TELEMAC3D: VARIABLES U2D, V2D OR H NOT FOUND'
          WRITE(LU,*) '           2D CONTINUATION IMPOSSIBLE'
          WRITE(LU,*) '           IF IT IS A LANGUAGE PROBLEM'
          WRITE(LU,*) '           YOU CAN MODIFY NOMVAR_2D_IN_3D'
          CALL PLANTE(1)
          STOP
        ENDIF
      ENDIF
!	  
!     COPIES THE BOTTOM TOPOGRAPHY INTO Z (= MESH3D%Z%R)
!    (IF IT IS A CONTINUATION, Z WILL BE ALSO FOUND
!     IN THE PREVIOUS RESULTS FILE. ANYWAY THE COPY IS USEFUL HERE
!     TO AVOID A CRASH IN CONDIM)
!
      CALL OV('X=Y     ', X=Z(1:NPOIN2), Y=ZF%R, DIM1=NPOIN2)
!
!
!     GET OFFSET I_ORIG AND J_ORIG IN THE MESH
!
      CALL GET_MESH_ORIG(T3D_FILES(T3DGEO)%FMT,T3D_FILES(T3DGEO)%LU,
     &                   I_ORIG,J_ORIG,IERR)
!
      IF(I_ORIG.NE.0.OR.J_ORIG.NE.0) THEN
!       IF MODELLING TIDE WITH GEOGRAPHIC SYSTEM DIFFERENT FROM UTM OR
!       LAMBERT PROJECTIONS
        IF(TIDALTYPE.NE.0.AND.
     &     GEOSYST.NE.2.AND.GEOSYST.NE.3.AND.GEOSYST.NE.4) THEN
!         WARNING MESSAGE FOR EVERY TIDAL SOLUTION (I.E. WITH JMJ, 
!         POSSIBLE IF HARMONIC CONSTANTS FILE ALREADY GENERATED)
          IF(GEOSYST.EQ.1.OR.(SPHERI.AND..NOT.(PROTYP.EQ.2))) THEN
            WRITE(LU,*) 'WHEN MODELLING TIDE, THE OFFSET _ORIG IS'
            WRITE(LU,*) 'NOT IMPLEMENTED FOR GEOGRAPHIC SYSTEM OTHER'
            WRITE(LU,*) 'THAN UTM OR LAMBERT PROJECTIONS'
!           BUT ERROR MESSAGE IF USING OSU SOLUTIONS LIKE TPXO
            IF(TIDALDB.EQ.2) THEN
              WRITE(LU,*) 'IF USING OSU TIDAL SOLUTIONS LIKE TPXO'
              WRITE(LU,*) 'YOU HAVE TO IMPLEMENT YOUR TRANSFORMATION'
              WRITE(LU,*) 'WITH GEOGRAPHIC SYSTEM = 0 (USER DEFINED)'
              WRITE(LU,*) 'GOOD LUCK!'
              CALL PLANTE(1)
              STOP
            ENDIF
          ENDIF
        ENDIF
      ENDIF
!
!     NOW CALLS CONDIM EVEN IF A COMPUTATION IS CONTINUED
!    (DONE TO RETRIEVE ZSTAR)
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE CONDIM'
      CALL CONDIM
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE CONDIM'
!
!     COMPUTES TRANSF AND ZCHAR
!
      CALL TRANSF_ZCHAR(TRANSF,ZCHAR,ZSTAR,TRANSF_PLANE,NPLAN)
!
!     CLIPS POSSIBLE NEGATIVE DEPTHS SET BY USER
!
      CALL OS('X=+(Y,C)',X=H,Y=H,C=0.D0)
!
!     IF COMPUTATION CONTINUED, RETRIEVES SOME VARIABLES + Z + DEPTH
!
      IF(.NOT.DEBU) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE SUITE'
        CALL READ_DATASET(T3D_FILES(T3DPRE)%FMT,T3D_FILES(T3DPRE)%LU,
     &                  VARSO3,MESH3D%NPOIN,START_RECORD,AT,TEXTP3,
     &                  TROUVE,ALIRE3D,LISTIN,START_RECORD.EQ.0,MAXVAR)
        IF(TROUVE(1).NE.1) THEN
          WRITE(LU,*) ' '
          WRITE(LU,*) 'Z COORDINATE MISSING IN THE'
          WRITE(LU,*) 'PREVIOUS COMPUTATION FILE'
          CALL PLANTE(1)
          STOP
        ENDIF
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE SUITE'
!
!*** Modif EPO *********************************************************
!
! Appel de T3D_CORFON ici: conservation des Z (3D) passés et calcul des 
! nouvelles hauteurs d'eau avec ZF mis à jour:
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'CALLING T3D_CORFON'
        CALL T3D_CORFON(ZF,T2_01,T2_02,ZF%R,T2_01%R,T2_02%R,
     &                X,Y,PRIVE,NPOIN2,LISFON,LISFON_AFTER,.FALSE.,
     &                MASKEL,MATR2H,MESH2D,SVIDE)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETURN FROM T3D_CORFON'
	  
        DO K=1,NPOIN2
          H%R(K)=Z(K+NPOIN2*(NPLAN-1))-ZF%R(K)
          !ZF%R(K)=Z(K)
        ENDDO
!
!*** Fin Modif EPO *****************************************************
!
!       SEE VARSO3 IN POINT FOR INDICES 8 AND 9 (K AND EPSILON)
        IF(TROUVE(8).EQ.1.AND.TROUVE(9).EQ.1) THEN
          AKEP=.FALSE.
          AKOM=.FALSE.
          AKSA=.FALSE.
          IF(ITURBV.NE.6) THEN
!           WILL BE USED BY FIRST CALL TO PREADV
            CALL OS('X=Y     ',X=AKN,Y=AK)
            CALL OS('X=Y     ',X=EPN,Y=EP)
          ENDIF
          CALL OS('X=Y     ',X=NUN,Y=NU)
        ENDIF
!
        IF(TROUVE(19).EQ.1) THEN
!         RETRIEVING DH AND HN
          DO I=1,NPOIN2
            DH%R(I)=DHHN%R(I)
            HN%R(I)=DHHN%R(I+NPOIN2)
          ENDDO
        ELSE
          DO I=1,NPOIN2
            HN%R(I)=H%R(I)
          ENDDO
        ENDIF
!
      ENDIF
!
!     TIME OPTIONALLY RESET TO ZERO
!
      IF(RAZTIM) THEN
        AT=0.D0
        WRITE(LU,*) 'ELAPSED TIME RESET TO ZERO'
      ENDIF
!
!     ONLY TA IS INITIALISED IN CONDIM OR READ_DATASET
!
      IF(NTRAC.GT.0) CALL OS ('X=Y     ', X=TRN, Y=TA)
!
!     INITIALISES SEDIMENT PROPERTIES
!
      IF(S3D_SEDI) THEN
!
!    THE CALL TO CONDIS SHOULD ALWAYS BE DONE FIRST
!    VALUES THEN MAY BE MODIFIED FROM SEDIMENTOLOGICAL 2D RESULT
!
          IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE CONDIS'
          CALL CONDIS(S3D_IVIDE%R,S3D_EPAI,S3D_TREST,S3D_CONC,
     &                S3D_TEMP%R,S3D_HDEP%R,
     &                ZR%R,ZF%R,X,Y,NPOIN2,NPF%I,
     &                S3D_NCOUCH,S3D_TASSE,S3D_ITASS, S3D_RHOS,S3D_XKV,
     &                S3D_CFDEP,S3D_ESOMT, S3D_TOCE, S3D_SEDCO,
     &                S3D_CONC_LAYER,S3D_TOCE_LAYER, S3D_ES_LAYER,
     &                S3D_SEDNCO, S3D_MIXTE,S3D_EPAICO%R, S3D_EPAINCO%R,
     &                S3D_PVSCO%R, S3D_PVSNCO%R,S3D_PVSNCO0)
!
          IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE CONDIS'
!
!        ELSE
!
        IF(T3D_FILES(T3DSUS)%NAME(1:1).NE.' ') THEN
!
! IMPORTANT NOTE: SUISED NOW READS A 2D FILE PREVIOUS SEDIMENTOLOGICAL
!                 FILE
!
          CALL SUISED(S3D_EPAI,S3D_HDEP%R,ZR%R,ZF%R,T2_05,
     &                NPOIN2,S3D_NCOUCH,S3D_ITASS,
     &                T3D_FILES(T3DSUS)%LU,T3D_FILES(T3DSUS)%FMT)
!
        ENDIF
!
! LOOK FOR CRITICAL SHEAR STRESS FOR EROSION IN GEOMETRY FILE
        IF (S3D_READ_TOCE) THEN
          IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MUDSTRESS3D'
          CALL MUDSTRESS3D(T3D_FILES(T3DGEO)%LU,T3D_FILES(T3DGEO)%FMT,
     &                 S3D_LAYTOCE,S3D_NCOUCH,MESH2D)
          IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE MUDSTRESS3D'
        ENDIF
!
!       SO FAR CONSTANT MEAN DIAMETER=S3D_D50
        CALL OS('X=C     ',X=S3D_DMOY,C=S3D_D50)
      ENDIF
!
! CLIPS H AND COMPUTES Z, HPROP AND ZPROP
! NOTE : HMIN = -1000.0 IN DICTIONARY BUT HMIN IS AT LEAST 0.0
!        IF OPTBAN=2
!
      IF(OPTBAN.EQ.2) THEN
        CALL CLIP (H, HMIN, .TRUE., 1.D6, .FALSE., 0)
      ENDIF
!
      CALL CALCOT(Z,H%R)
!     BOTTOM OF ZPROP UPDATED
      CALL OV('X=Y     ', X=ZPROP%R(1:NPOIN2), Y=Z3%R(1:NPOIN2),
     &        DIM1=NPOIN2)
!
!-----------------------------------------------------------------------
! MASKING:
!
      IF(MSK) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MASK3D'
        CALL MASK3D(MESH3D%IFABOR%I,MASKEL%R,MASKPT,MASKBR%R,
     &              X2%R,Y2%R,ZF%R,ZFE%R,H%R,HMIN,AT,LT,IT1%I,
     &              MESH3D%NELBOR%I,NELMAX2,NELEM2,NPOIN2,MESH2D%NELEB,
     &              NPLAN,NETAGE,IELM3,MESH2D)
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MASK3D'
      ENDIF
!
! MESH FOR PROPAGATION (IF NOT FOUND IN PREVIOUS RESULTS FILE)
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MESH_PROP'
      CALL MESH_PROP(HPROP,HN,H,PROLIN,HAULIN,TETAH,NSOUSI,ZPROP,
     &               IPBOT,NPOIN2,NPLAN,OPTBAN,SIGMAG,OPT_HNEG,
     &               MESH3D,VOLU3D,VOLU3DPAR,
     &               UNSV3D,MASKEL,IELM3)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE MESH_PROP'
!
! INITIALISES THE MEAN VELOCITY IN 2D
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VERMOY'
      CALL VERMOY(U2D%R,V2D%R,U%R,V%R,2,Z,
     &            T3_01%R,T3_02%R,T3_03%R,1,NPLAN,NPOIN2,NPLAN,OPTBAN)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VERMOY'
!
! INITIALISES THE VARIABLES COMING FROM TOMAWAC
! AND OUTPUT GRAPHICAL VARIABLES
!
      IF(INCLUS(COUPLING,'TOMAWACT3D')) THEN
        CALL OS('X=0     ',X=WIP)
        CALL OS('X=0     ',X=WIPDX)
        CALL OS('X=0     ',X=WIPDY)
        CALL OS('X=0     ',X=USTOKES)
        CALL OS('X=0     ',X=VSTOKES)
        CALL OS('X=0     ',X=WSTOKES)
        CALL OS('X=0     ',X=WST1)
        CALL OS('X=0     ',X=WST2)
        CALL OS('X=0     ',X=US2D)
        CALL OS('X=0     ',X=VS2D)
        CALL OS('X=0     ',X=DUS2D)
        CALL OS('X=0     ',X=DVS2D)
        CALL OS('X=0     ',X=DUS2DX)
        CALL OS('X=0     ',X=DVS2DY)
        CALL OS('X=0     ',X=USTX)
        CALL OS('X=0     ',X=VSTY)
        CALL OS('X=0     ',X=FDX)
        CALL OS('X=0     ',X=FDY)
        CALL OS('X=0     ',X=FBX)
        CALL OS('X=0     ',X=FBY)
        CALL OS('X=0     ',X=CFWC)
        CALL OS('X=0     ',X=FDK)
        CALL OS('X=0     ',X=FWX)
        CALL OS('X=0     ',X=FWY)
      ENDIF

!-----------------------------------------------------------------------
! HARMONISES  BOUNDARY CONDITIONS
! INITIALISES BOUNDARY CONDITIONS FOR TELEMAC
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE LICHEK'
      CALL LICHEK(LIMPRO%I,NPTFR2,
     &            MESH2D%IKLBOR%I,MESH2D%NELEB,MESH2D%NELEBX)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE LICHEK'
!
!-----------------------------------------------------------------------
! INITIALISES THE VOLUMES ASSOCIATED WITH THE NODES
!
      CALL VECTOR(VOLU, '=', 'MASBAS          ',IELM3,1.D0-AGGLOH,
     &  SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,.FALSE.,MASKEL)
      IF(AGGLOH.GT.1.D-6) THEN
        CALL VECTOR(VOLU, '+', 'MASBAS2         ',IELM3,AGGLOH,
     &  SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,.FALSE.,MASKEL)
      ENDIF
!
!     IF NEW COMPUTATION OR DHHN NOT IN PREVIOUS RESULT FILE
!
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(19).NE.1)) THEN
        CALL OS('X=Y     ',X=VOLUN,Y=VOLU)
      ELSE
!       COMPUTING THE CORRECT VOLUN (ASSUMING AGGLOH HAS NOT CHANGED..)
        CALL CALCOT(Z,HN%R)
        CALL VECTOR(VOLUN, '=', 'MASBAS          ',IELM3,1.D0-AGGLOH,
     &    SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,.FALSE.,MASKEL)
        IF(AGGLOH.GT.1.D-6) THEN
          CALL VECTOR(VOLUN, '+', 'MASBAS2         ',IELM3,AGGLOH,
     &    SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH3D,.FALSE.,MASKEL)
        ENDIF
!       RESTORING Z DONE WITH H
        CALL CALCOT(Z,H%R)
      ENDIF
!
      IF(NCSIZE.GT.1) THEN
        CALL OS('X=Y     ',X=VOLUPAR,Y=VOLU)
        CALL PARCOM(VOLUPAR,2,MESH3D)
        CALL OS('X=Y     ',X=VOLUNPAR,Y=VOLUPAR)
      ENDIF
!
!     IN 2D
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MASBAS2D'
      CALL MASBAS2D(VOLU2D,V2DPAR,UNSV2D,
     &              IELM2H,MESH2D,MSK,MASKEL,T2_01,SVIDE)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE MASBAS2D'
!
!-----------------------------------------------------------------------
!
! INITIALISES THE METEOROLOGICAL VARIABLES
!
      IF(VENT.OR.ATMOS.OR.INCLUS(COUPLING,'WAQTEL')) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE METEO'
        CALL METEO(PATMOS%R,WINDX%R,WINDY%R,
     &             CST_WINDX,CST_WINDY,AT,LT,NPOIN2,VENT,ATMOS,
     &             T3ATMA,T3ATMB,
     &             T3D_FILES,LISTIN,CST_PATMOS,
     &             INCLUS(COUPLING,'WAQTEL'),PLUIE,OPTWIND)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE METEO'
!       RAIN TEMPERATURE EQUAL TO THE AIR TEMPERATURE OR 1.D0
!       IS BETTER THAN 0.D0
        IF(IND_T.NE.0.AND.INCLUS(COUPLING,'WAQTEL')) THEN
          TRAIN(IND_T) = MAX(TAIR%R(1),1.D0)
        ENDIF
        IF( (INCLUS(COUPLING,'WAQTEL') ).AND.
     &      (13*INT(WAQPROCESS/13).EQ.WAQPROCESS) ) THEN
          DO I=1,NPOIN2
            WINDSPD%R(I) = SQRT(WINDX%R(I)**2+WINDY%R(I)**2)
          ENDDO
        ENDIF
      ENDIF
!
!-----------------------------------------------------------------------
! FREE SURFACE AND BOTTOM GRADIENTS
! INITIALISES DSSUDT = 0
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE GRAD2D'
      CALL GRAD2D(GRADZF%ADR(1)%P,GRADZF%ADR(2)%P,ZPROP,NPLAN,SVIDE,
     &            UNSV2D,T2_02,T2_03,T2_04,IELM2H,MESH2D,MSK,MASKEL)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE GRAD2D'
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE FSGRAD'
      CALL FSGRAD(GRADZS,ZFLATS,Z(NPOIN3-NPOIN2+1:NPOIN3),
     &            ZF,IELM2H,MESH2D,MSK,MASKEL,
     &            UNSV2D,T2_01,NPOIN2,OPTBAN,SVIDE)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE FSGRAD'
!
      CALL OS('X=C     ',X=DSSUDT,C=0.D0)
!
!-----------------------------------------------------------------------
! INITIALISES THE CULVERT VARIABLES
      IF(NBUSE.GT.0) THEN
        ! DEFINE POINTERS TO FILL ISCE AND KSCE BEYOND NPTSCE+1
        ! UP TO NPTSCE+2*NBUSE+1
        ENTBUS  => ISCE(NPTSCE+1:NPTSCE+NBUSE)
        SORBUS  => ISCE(NPTSCE+NBUSE+1:NPTSCE+2*NBUSE)
        KENTBUS => KSCE(NPTSCE+1:NPTSCE+NBUSE)
        KSORBUS => KSCE(NPTSCE+NBUSE+1:NPTSCE+2*NBUSE)
        UBUS    => U_SCE%R(NPTSCE+1:NPTSCE+2*NBUSE)
        VBUS    => V_SCE%R(NPTSCE+1:NPTSCE+2*NBUSE)
        IF(DEBUG.GT.0) WRITE(LU,*) 'CALLING LECBUS'
        CALL LECBUS(RELAXB,NBUSE,ENTBUS,SORBUS,LRGBUS%R,
     &              HAUBUS%R,CLPBUS%I,ALTBUS%R,CSBUS%R,CEBUS%R,
     &              ANGBUS%R,LBUS%R,T3D_FILES(T3DBUS)%LU,MESH2D,
     &              CV%R,C56%R,CV5%R,C5%R,CTRASH%R,FRICBUS%R,
     &              LONGBUS%R,CIRC%I,DELBUS%R,AABUS%I)
        IF(DEBUG.GT.0) WRITE(LU,*) 'BACK FROM LECBUS'
      ENDIF
!-----------------------------------------------------------------------
! INITIALISES K AND EPSILON
! IF AKEP = .FALSE. K AND EPSILON HAVE BEEN GIVEN IN LECSUI OR CONDIM
!
      IF((ITURBV.EQ.3.OR.ITURBH.EQ.3).AND.AKEP) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE KEPINI'
        CALL KEPINI(AK%R,EP%R,U%R,V%R,Z,
     &             ZF%R,NPOIN2,NPLAN,DNUVIH,DNUVIV,KARMAN,CMU,KMIN,EMIN)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE KEPINI'
!       WILL BE USED BY FIRST CALL TO PREADV
        CALL OS('X=Y     ',X=AKN,Y=AK)
        CALL OS('X=Y     ',X=EPN,Y=EP)
      ENDIF
!
      IF((ITURBV.EQ.7.OR.ITURBH.EQ.7).AND.AKOM) THEN
        CALL OS('X=C     ',X=AK,C=KMIN)
        CALL OS('X=C     ',X=EP,C=EMIN)
        CALL OS('X=0     ',X=ROTAT)
      ENDIF
!
!-----------------------------------------------------------------------
!
! COMPUTES (DELTA RHO)/RHO FOR THE INITIAL OUTPUT
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE DRSURR'
      CALL DRSURR(DELTAR,TA,BETAC,T0AC,RHO,RHO0,XMVS0,S3D_RHOS,DENLAW,
     &            S3D_SEDI,NTRAC,IND_T,IND_S, IND_SED, NSUSP_TEL,
     &            S3D_MIXTE,NUM_ISUSP_ICLA,NSICLA)
!
      IF( (INCLUS(COUPLING,'WAQTEL') ).AND.
     &    (13*INT(WAQPROCESS/13).EQ.WAQPROCESS) ) THEN
!     ORDER: NUMBER OF THE PLANES FROM TOP TO BOTTOM
!     THEN NUMBER OF THE 2D NODES
        DO I=1,NPOIN2
          DO K=1,NPLAN
            RHO4AED2%R(K+(I-1)*NPLAN) = RHO%R(I+(NPLAN-K)*NPOIN2)
          ENDDO
        ENDDO
      ENDIF
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE DRSURR'
!
!-----------------------------------------------------------------------
!
! INITIALISES U* FOR OUTPUT OF INITIAL CONDITIONS AND SISYPHE
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE COEFRO'
      CALL COEFRO(CF,H,U2D,V2D,KARMAN,KFROT,RUGOF,GRAV,MESH2D,T2_01,
     &            .FALSE.)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE COEFRO'
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE TFOND'
      CALL TFOND(AUBORF%R,CF%R,U2D%R,V2D%R,U%R,V%R,W%R,KARMAN,
     &           LISRUF,DNUVIV,Z,NPOIN2,KFROT,RUGOF%R,UETCAR%R,
     &           NONHYD,OPTBAN,HN%R,GRAV,IPBOT%I,NPLAN)
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE TFOND'

!-----------------------------------------------------------------------
! INITIALISES NU FOR S-A MODEL
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE SAINI'
      IF((ITURBV.EQ.5.OR.ITURBV.EQ.9).AND.AKSA) THEN
        CALL SAINI(NU%R,U,V,W,Z,NPOIN2,NPLAN,KARMAN,NUMIN,UETCAR%R,
     &             H%R,KFROT,ITURBV,STRAIN,MSK,MASKEL,SVIDE,MESH3D)
        CALL OS('X=Y     ',X=NUN,Y=NU)
      ENDIF
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE SAINI'
!
!-----------------------------------------------------------------------
!
! COMPUTES THE VISCOSITIES VISCVI AND VISCTA
!
      IF(ITURBH.EQ.1.OR.ITURBV.EQ.1) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISCOS'
        CALL VISCOS(VISCVI,VISCTA,DNUTAV,DNUTAH,
     &              DNUVIV,DNUVIH,NTRAC,ITURBH,ITURBV)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISCOS'
!
      ENDIF
!
      IF(ITURBV.EQ.2) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISCLM'
        CALL VISCLM(VISCVI,VISCTA,RI,U,V,DELTAR,Z3,H,
     &              T3_01, T3_02, T3_03, T3_04, T3_05, T3_07,
     &              MESH3D, IELM3, GRAV, NPLAN,
     &              NPOIN3, NPOIN2, NTRAC, MSK, MASKEL,
     &              MIXING,DAMPING,DNUVIV,DNUTAV,KARMAN,
     &              PRANDTL,KFROT,RUGOF,ZF,LINLOG,IPBOT%I)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISCLM'
!
        IF(INCLUS(COUPLING,'TOMAWACT3D')) THEN
!         ADDING VERTICAL MIXING DUE TO WAVES
          CALL OS('X=X+Y   ',X=VISCVI%ADR(3)%P,Y=FDK)
        ENDIF
!
      ENDIF
!
      IF(ITURBV.EQ.3.OR.ITURBH.EQ.3) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISCKE'
!       FOR FIRST PRINTOUT (RI ONLY DONE IN SOUKEP LATER)
        CALL OS('X=0     ',X=RI)
        CALL VISCKE(VISCVI,VISCTA,AK,EP,NTRAC,CMU,
     &              DNUVIH,DNUVIV,DNUTAH,DNUTAV,EMIN,
     &              ITURBH,ITURBV,PRANDTL)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISCKE'
!
        IF(INCLUS(COUPLING,'TOMAWACT3D').AND.ITURBV.EQ.3) THEN
!         ADDING VERTICAL MIXING DUE TO WAVES
          CALL OS('X=X+Y   ',X=VISCVI%ADR(3)%P,Y=FDK)
        ENDIF
!
      ENDIF
!
      IF(ITURBH.EQ.4) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISSMA'
        CALL VISSMA(VISCVI,VISCTA,DNUTAH,DNUVIH,DNUVIV,DNUTAV,
     &              U,V,W,T3_01,T3_02,T3_03,T3_04,T3_05,T3_06,
     &              SVIDE,MESH3D,
     &              IELM3,NTRAC,MSK,MASKEL,ITURBV)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISSMA'
!
      ENDIF
!
      IF(ITURBV.EQ.7.OR.ITURBH.EQ.7) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISCKO'
        CALL VISCKO(VISCVI,VISCTA,ROTAT,AK,EP,NTRAC,
     &              DNUVIH,DNUVIV,DNUTAH,DNUTAV,ITURBH,ITURBV,
     &              T3_01,T3_02)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISCKO'
!
      ENDIF
!
      IF(ITURBV.EQ.5.OR.ITURBH.EQ.5.OR.ITURBV.EQ.9.OR.ITURBH.EQ.9) THEN

        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VISCSA'
        CALL VISCSA(VISCVI,VISCNU,NU,
     & DNUVIH,DNUVIV,AK,EP,STRAIN)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VISCSA'
      ENDIF
!
      IF(ITURBV.EQ.6) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE GOTM_COUPLING_INIT'
          CALL GOTM_COUPLING_INIT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE GOTM_COUPLING_INIT'
      ENDIF

      IF(OPTBAN.EQ.1) THEN
!
        CALL VISCLIP(VISCVI,VISCTA,H,NPLAN,NPOIN2,NTRAC,HLIM)
!
      ENDIF

!     ----------------------------------------
!      COMPUTES WALL DISTANCES-2D COMPUTATION
!     ---------------------------------------
      IF(ITURBV.EQ.5.OR.ITURBH.EQ.5.OR.ITURBV.EQ.9.OR.ITURBH.EQ.9) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WALLDIST'
        CALL WALLDIST
     &      (WDIST,T2_11,T2_12,T2_13,T2_14,T2_15,FLBOR,TRAV2,
     &      MAT2D%ADR(1)%P,MAT2D%ADR(2)%P,SVIDE,
     &      LIUBOL,IELM2H,NPTFR2,MESH2D)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WALLDIST'
      ENDIF

!     ------------------------------------
!     PREPARES THE 3D OUTPUT FILE :
!     ------------------------------------
!
!     OUTPUT FOR THE INITIAL CONDITIONS
!
      IF (INFOGR) CALL MITTIT(1,AT,LT)
!
!     PREPARES THE 2D AND 3D OUTPUT FOR INITIAL CONDITIONS
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE PRERES_TELEMAC3D'
      CALL PRERES_TELEMAC3D
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE PRERES_TELEMAC3D'
!
!-----------------------------------------------------------------------
!
!     COUPLING WITH DELWAQ
!
      IF(INCLUS(COUPLING,'DELWAQ')) THEN
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE TEL4DELWAQ'
        CALL TEL4DELWAQ(NPOIN3,NPOIN2,MESH2D%NSEG,IKLE2%I,
     &  MESH2D%ELTSEG%I,MESH2D%GLOSEG%I,
     &  MESH2D%GLOSEG%DIM1,X,Y,MESH3D%NPTFR,LIHBOR%I,MESH3D%NBOR%I,
     &  NPLAN,AT,DT,LT,NIT,H%R,HPROP%R,MESH3D%Z%R,U%R,V%R,
     &  TA%ADR(MAX(IND_S,1))%P%R,
     &  TA%ADR(MAX(IND_T,1))%P%R,VISCVI%ADR(3)%P%R,TITCAS,
     &  T3D_FILES(T3DGEO)%NAME,T3D_FILES(T3DCLI)%NAME,WAQPRD,
     &  T3DDL1,T3DDL2,T3DDL3,T3DDL5,T3DDL6,T3DDL7,T3DL11,T3DDL4,
     &  T3DDL8,T3DDL9,T3DL10,INFOGR,NELEM2,SALI_DEL,TEMP_DEL,VELO_DEL,
     &  DIFF_DEL,MARDAT,MARTIM,FLODEL%R,V2DPAR%R,MESH2D%KNOLG%I,
     &  T3D_FILES)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE TEL4DELWAQ'
!
      ENDIF
!
!     3D OUTPUT (OPTIONAL)
!
      IF(T3D_FILES(T3DRES)%NAME(1:1).NE.' ') THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_HEADER'
        CALL WRITE_HEADER(T3D_FILES(T3DRES)%FMT, ! RESULT FILE FORMAT
     &                    T3D_FILES(T3DRES)%LU,  ! RESULT FILE LU
     &                    TITCAS,     ! TITLE
     &                    MAXVA3,     ! MAX NUMBER OF OUTPUT VARIABLES
     &                    TEXT3,      ! NAMES OF OUTPUT VARIABLES
     &                    SORG3D)     ! OUTPUT OR NOT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_HEADER'
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_MESH'
        CALL WRITE_MESH(T3D_FILES(T3DRES)%FMT, ! RESULT FILE FORMAT
     &                  T3D_FILES(T3DRES)%LU,  ! RESULT FILE LU
     &                  MESH3D,
     &                  NPLAN,           ! NUMBER OF PLANE /NA/
     &                  DATE,            ! START DATE
     &                  TIME,            ! START HOUR
     &                  T3_01,T3_02,     ! WORKING ARRAYS
     &                  NCSIZE.GT.1,NPTIR,
     &                  NGEO=T3D_FILES(T3DGEO)%LU,
     &                  GEOFORMAT=T3D_FILES(T3DGEO)%FMT,
     &                  LATLONG=PROTYP0.EQ.3.AND.KEEP_LONLAT)! IF LONG-LAT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_MESH'
      ENDIF
!
!     3D OUTPUT FOR RESTART
!
      IF(RESTART_MODE.AND.T3D_FILES(T3DRST)%NAME(1:1).NE.' ') THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_HEADER EN RESTART'
        CALL WRITE_HEADER(T3D_FILES(T3DRST)%FMT, ! RESULT FILE FORMAT
     &                    T3D_FILES(T3DRST)%LU,  ! RESULT FILE LU
     &                    TITCAS,     ! TITLE
     &                    MAXVA3,     ! MAX NUMBER OF OUTPUT VARIABLES
     &                    TEXT3,      ! NAMES OF OUTPUT VARIABLES
     &                    SOREST)     ! OUTPUT OR NOT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_HEADER EN RESTART'
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_MESH EN RESTART'
        CALL WRITE_MESH(T3D_FILES(T3DRST)%FMT, ! RESULT FILE FORMAT
     &                  T3D_FILES(T3DRST)%LU,  ! RESULT FILE LU
     &                  MESH3D,
     &                  NPLAN,           ! NUMBER OF PLANE /NA/
     &                  DATE,            ! START DATE
     &                  TIME,            ! START HOUR
     &                  T3_01,T3_02,     ! WORKING ARRAYS
     &                  NCSIZE.GT.1, NPTIR,
     &                  NGEO=T3D_FILES(T3DGEO)%LU,
     &                  GEOFORMAT=T3D_FILES(T3DGEO)%FMT,
     &                  LATLONG=PROTYP0.EQ.3.AND.KEEP_LONLAT)!IF LONG-LAT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_MESH EN RESTART'
      ENDIF
!
!     THESE VARIABLES ARE INITIALISED FOR THE FIRST CALL TO BIEF_DESIMP
!
      IF(SORG3D(14)) CALL OS('X=0     ',X=UCONV)
      IF(SORG3D(15)) CALL OS('X=0     ',X=VCONV)
      IF(SORG3D(16)) CALL OS('X=0     ',X=WCONV)
      IF(SORG3D(18)) CALL OS('X=0     ',X=DM1)
      IF(SORG3D(19)) CALL OS('X=0     ',X=DHHN)
      IF(SORG3D(20)) CALL OS('X=0     ',X=UCONVC)
      IF(SORG3D(21)) CALL OS('X=0     ',X=VCONVC)
      IF(SORG3D(22)) CALL OS('X=0     ',X=UD)
      IF(SORG3D(23)) CALL OS('X=0     ',X=VD)
      IF(SORG3D(24)) CALL OS('X=0     ',X=WD)
      IF(SORG3D(25)) CALL OS('X=0     ',X=PRIVE%ADR(1)%P)
      IF(SORG3D(26)) CALL OS('X=0     ',X=PRIVE%ADR(2)%P)
      IF(SORG3D(27)) CALL OS('X=0     ',X=PRIVE%ADR(3)%P)
      IF(SORG3D(28)) CALL OS('X=0     ',X=PRIVE%ADR(4)%P)
!
!     INITIAL CONDITIONS OF 3D RESULTS FILE
!
      IF(T3D_FILES(T3DRES)%NAME(1:1).NE.' ') THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE BIEF_DESIMP'
        CALL BIEF_DESIMP(T3D_FILES(T3DRES)%FMT,VARSO3,
     &                   NPOIN3,T3D_FILES(T3DRES)%LU,AT,LT,
     &                   LISPRD,GRAPRD,
     &                   SORG3D,SORIM3,MAXVA3,TEXT3,GRADEB,LISDEB)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE BIEF DESIMP'
      ENDIF
!
!     INITIAL CONDITIONS OF 3D OUTPUT FOR RESTART (NOT USED BY NOW MANDATORY)
!
      IF(RESTART_MODE.AND.T3D_FILES(T3DRST)%NAME(1:1).NE.' ') THEN
        CALL OS('X=0     ',X=UCONV)
        CALL OS('X=0     ',X=VCONV)
        CALL OS('X=0     ',X=DM1)
        CALL OS('X=0     ',X=DHHN)
        CALL OS('X=0     ',X=UCONVC)
        CALL OS('X=0     ',X=VCONVC)
        CALL OS('X=0     ',X=UD)
        CALL OS('X=0     ',X=VD)
        IF(NONHYD) THEN
          CALL OS('X=0     ',X=WCONV)
          CALL OS('X=0     ',X=WD)
        ENDIF
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE BIEF_DESIMP RESTART MODE'
        CALL BIEF_DESIMP(T3D_FILES(T3DRST)%FMT,VARSO3,NPOIN3,
     &                   T3D_FILES(T3DRST)%LU,AT,LT,
     &                   1,NIT,SOREST,SORIS3,MAXVA3,TEXT3,0,0)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE BIEF DESIMP RESTART MODE'
      ENDIF
!
!     SEDIMENTOLOGY OUTPUT
!
      IF(S3D_SEDI.AND.T3D_FILES(S3D_T3DSED)%NAME(1:1).NE.' ') THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE DESSED'
        CALL DESSED(NPF%I,S3D_IVIDE%R,S3D_EPAI,S3D_HDEP%R,
     &              S3D_TEMP%R,ZR%R,NPOIN2,S3D_NPFMAX,
     &              S3D_NCOUCH,GRAPRD,LT,S3D_DTC,S3D_TASSE,
     &              S3D_GIBSON,
     &              T3D_FILES(S3D_T3DSED)%LU,TITCAS,S3D_BIRSED,0)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE DESSED'
      ENDIF
!
! PREPARES THE 2D OUTPUT FILE
!
      IF(T3D_FILES(T3DHYD)%NAME(1:1).NE.' ') THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_HEADER'
        CALL WRITE_HEADER(T3D_FILES(T3DHYD)%FMT, ! FORMAT FICHIER RESULTAT
     &                    T3D_FILES(T3DHYD)%LU,  ! LU FICHIER RESULTAT
     &                    TITCAS,     ! TITRE DE L'ETUDE
     &                    MAXVAR,     ! MAX VARIABLES SORTIE
     &                    TEXTE,      ! NOMS VARIABLES SORTIE
     &                    SORG2D)     ! SORTIE OU PAS DES VARIABLES
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_HEADER'
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE WRITE_MESH'
        CALL WRITE_MESH(T3D_FILES(T3DHYD)%FMT, ! RESULT FILE FORMAT
     &                  T3D_FILES(T3DHYD)%LU,  ! RESULT FILE LU
     &                  MESH2D,
     &                  1,               ! NUMBER OF PLANE /NA/
     &                  DATE,            ! START DATE
     &                  TIME,            ! START HOUR
     &                  T2_01,T2_02,     ! WORKING ARRAYS
     &                  NCSIZE.GT.1, NPTIR,
     &                  NGEO=T3D_FILES(T3DGEO)%LU,
     &                  GEOFORMAT=T3D_FILES(T3DGEO)%FMT,
     &                  LATLONG=PROTYP0.EQ.3.AND.KEEP_LONLAT) ! IF LONG-LAT
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE WRITE_MESH'
!
!       2D OUTPUT
!
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE BIEF_DESIMP POUR 2D'
        CALL BIEF_DESIMP(T3D_FILES(T3DHYD)%FMT,VARSOR,
     &                   NPOIN2,T3D_FILES(T3DHYD)%LU,AT,LT,
     &                   LISPRD,GRAPRD,
     &                   SORG2D,SORIMP,MAXVAR,TEXTE,GRADEB,LISDEB)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE BIEF_DESIMP POUR 2D'
      ENDIF
!
!-----------------------------------------------------------------------
! INITIALISES MASS BALANCE AND CUMULATIVE FLUXES
!
      IF(BILMAS) THEN
!
        CALL MITTIT(10,AT,LT)
!
        CALL MASS3D(.TRUE.)
!
        CALL OS ( 'X=Y     ', X=MASINI, Y=MASSE)
        CALL OS ( 'X=0     ', X=FLUCUM         )
        MASINI_WATER=MASSE_WATER
        FLUXTOTCUM=0.D0
!       MAYBE NOT USEFUL
        CALL OS ( 'X=0     ', X=FLUX           )
!
!       INITIALIZATION OF BED PROPERTIES
!       S3D_MASBED0INITIAL MASS
!       S3D_ESOMT: CUMULATED BED EVOLUTION
!       --> Set to Zero but could be read from 2D file
!       S3D_MASDEP: Set to Zero
!
!GA: MASSED HAS BEEN MODIFIED TO ALLOW TWO SEDIMENT CLASSES
        IF(S3D_SEDI) THEN
          IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE MASSED'
          CALL MASSED(S3D_MASBED0,S3D_EPAI,S3D_CONC,T2_01%R,
     &                NPOIN2,S3D_NCOUCH,
     &                VOLU2D%R,
     &                S3D_CFDEP, S3D_EPAICO%R, S3D_EPAINCO%R, S3D_MIXTE)
          IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE MASSED'
          S3D_MASDEP= 0.D0
          CALL OS('X=0     ',X=S3D_ESOMT)
!         PRINT INITIAL MASS
          WRITE(LU,*) 'INITIAL MASS OF SEDIMENT BED :',S3D_MASBED0
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
! RETURNS WHEN THE NUMBER OF REQUIRED TIMESTEPS IS 0
!
      IF(NIT.EQ.0) THEN
        WRITE(LU,12)
        RETURN
      ENDIF
!
12    FORMAT(' BREAK IN TELEMAC-3D, NUMBER OF ITERATIONS ASKED NULL')
!
!-----------------------------------------------------------------------
!
! INITIALISES THE HORIZONTAL VELOCITY AFTER DIFFUSION
! IN ORDER TO ACCELERATE THE SOLVER CONVERGENCE
!
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(22).NE.1)) THEN
        CALL OS ( 'X=Y     ', X=UD, Y=U)
      ENDIF
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(23).NE.1)) THEN
        CALL OS ( 'X=Y     ', X=VD, Y=V)
      ENDIF
!
! INITIALISES THE FREE SURFACE AND DIFFERENT VERTICAL VELOCITIES
!
      IF(NONHYD) THEN
        IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(24).NE.1)) THEN
          CALL OS ( 'X=Y     ', X=WD,  Y=W  )
        ENDIF
      ENDIF
!
! SOURCE TERMS : FINDS LOCATION OF ALL THE SOURCES (USED IN PREADV HEREAFTER)
!                WILL SUBSEQUENTLY BE DONE AT EACH TIMESTEP
!
      IF(NPTSCE.GT.0) THEN
        IF(NREJETI.GT.0) THEN
!  SOURCE NODES GIVEN IN THE STEERING FILE
          DO I=1,NREJETI
            II = GLOBAL_TO_LOCAL_POINT(ISCE(I),MESH2D)
            ISCE(I)= II
            IF(II.NE.0) THEN
              XSCE(I)=MESH2D%X%R(II)
              YSCE(I)=MESH2D%Y%R(II)
            ENDIF
          ENDDO
        ELSE
!  SOURCE NODES COMPUTED IN PROXIM
!         LOCATION IN THE 2D MESH -> ISCE
          CALL PROXIM(ISCE,XSCE,YSCE,MESH2D%X%R,MESH2D%Y%R,NPTSCE,
     &                NPOIN2,IKLE2%I,NELEM2,NELMAX2)
        ENDIF
!       LOCATION ON THE VERTICAL -> KSCE
        CALL FINDKSCE(NPOIN2,NPLAN,Z3%R,NPTSCE,ISCE,ZSCE,KSCE,INFOGR)
      ENDIF
      IF(NBUSE.GT.0) THEN
!       LOCATION ON THE VERTICAL -> KENTBUS AND KSORBUS
        CALL FINDKSCE(NPOIN2,NPLAN,Z3%R,NBUSE,ENTBUS,ALTBUS%R,
     &                KENTBUS,INFOGR)
        CALL FINDKSCE(NPOIN2,NPLAN,Z3%R,NBUSE,SORBUS,ALTBUS%R,
     &                KSORBUS,INFOGR)
      ENDIF
!
!=======================================================================
!     PREPARATION OF ADVECTION FOR THE FIRST TIME STEP
!=======================================================================
!
!     WSCONV IS NOT INITIALISED BEFORE GOING INTO TRIDW2
      CALL OS('X=0     ',X=WSCONV)
!
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(14).NE.1)) THEN
        CALL OS('X=Y     ',X=UCONV,Y=U)
      ENDIF
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(15).NE.1)) THEN
        CALL OS('X=Y     ',X=VCONV,Y=V)
      ENDIF
!     USED ONLY FOR TRIDW3 IN PREADV
      IF(NONHYD) THEN
        IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(16).NE.1)) THEN
          CALL OS('X=Y     ',X=WCONV,Y=W)
        ENDIF
      ENDIF
!
      IF(.NOT.DEBU.AND.TROUVE(19).EQ.1) THEN
!       DH AND HN HAVE BEEN RECOVERED
!       OLD Z TEMPORARILY REDONE IN T3_01
        CALL OV('X=Y     ', X=T3_01%R(1:NPOIN2), Y=Z(1:NPOIN2),
     &          DIM1=NPOIN2)
        CALL CALCOT(T3_01%R,HN%R)
!       COMPUTING ZFLATS WITH THE OLD Z IN T3_01
        CALL FSGRAD(GRADZS,ZFLATS,T3_01%R(NPOIN3-NPOIN2+1:NPOIN3),
     &              ZF,IELM2H,MESH2D,MSK,MASKEL,
     &              UNSV2D,T2_01,NPOIN2,OPTBAN,SVIDE)
!       COMPUTING OLD ZCONV AND GRAZCO
        CALL MAKE_ZCONV(ZCONV,GRAZCO,ZFLATS,DH,HN,ZF,
     &                  TETAZCOMP,TETAH,NELEM2,NELMAX2,
     &                  OPTBAN,MESH2D%IKLE%I,MESH2D)
!       RESETTING NEW ZFLATS AND GRADZS
        CALL FSGRAD(GRADZS,ZFLATS,Z(NPOIN3-NPOIN2+1:NPOIN3),
     &              ZF,IELM2H,MESH2D,MSK,MASKEL,
     &              UNSV2D,T2_01,NPOIN2,OPTBAN,SVIDE)
      ELSE
        CALL OS('X=0     ',X=DH)
        CALL OS('X=0     ',X=ZCONV)
        CALL OS('X=0     ',X=GRAZCO)
      ENDIF
!
      IF(N_ADV(ADV_CAR).GT.0) THEN
        IF(DEBU.OR.TROUVE(20).NE.1.OR.TROUVE(21).NE.1) THEN
          CALL OS('X=Y     ',X=UCONVC,Y=U)
          CALL OS('X=Y     ',X=VCONVC,Y=V)
        ENDIF
      ENDIF
!
      IF(DEBU.OR.(.NOT.DEBU.AND.TROUVE(18).NE.1)) THEN
        CALL OS('X=0     ',X=DM1)
      ENDIF
!
!     INITIALISING SOURCES AND SMH
      CALL OS ('X=0     ',X=SMH)
!     SOURCES : COMPUTATION OF INITIAL INPUTS WHEN VARYING IN TIME
      IF(NPTSCE.GT.0) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPELS DE T3D_DEBSCE'
        DO I=1,NPTSCE
          QSCE2(I)=T3D_DEBSCE(AT,I,QSCE)
        ENDDO
        IF(DEBUG.GT.0) WRITE(LU,*) 'FIN DES APPELS DE T3D_DEBSCE'
        IF(NTRAC.GT.0) THEN
          IF(DEBUG.GT.0) WRITE(LU,*) 'APPELS DE T3D_TRSCE'
          DO I=1,NPTSCE
            DO ITRAC=1,NTRAC
              TA_SCE%ADR(ITRAC)%P%R(I)=T3D_TRSCE(AT,I,ITRAC)
            ENDDO
          ENDDO
          IF(DEBUG.GT.0) WRITE(LU,*) 'FIN DES APPELS DE T3D_TRSCE'
        ENDIF
      ENDIF
      IF(NBUSE.GT.0) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'CALLING BUSE'
        CALL BUSE(RELAXB,NBUSE,ENTBUS,SORBUS,GRAV,
     &            H%R,ZF%R,DBUS%R,LRGBUS%R,HAUBUS%R,CLPBUS%I,
     &            ALTBUS%R,CSBUS%R,CEBUS%R,ANGBUS%R,LBUS%R,
     &            NTRAC,TA,TBUS,UBUS,VBUS,U%R,V%R,LISTIN,
     &            CV%R,C56%R,CV5%R,C5%R,CTRASH%R,FRICBUS%R,
     &            LONGBUS%R,CIRC%I,DELBUS%R,OPTBUSE,V2DPAR,DT,
     &            SECBUS%R,MAXSCE,NPTSCE,NPOIN2,KSCE)
        IF(DEBUG.GT.0) WRITE(LU,*) 'BACK FROM BUSE'
        DO I=1,NBUSE
          QSCE2(NPTSCE+I)      =-DBUS%R(I)
          QSCE2(NPTSCE+NBUSE+I)= DBUS%R(I)
          !TODO: PRESCRIBE VERTICAL VELOCITIES, K AND EPSILON FOR
          !CULVERTS
          W_SCE%R(NPTSCE+I)= 0.D0
          W_SCE%R(NPTSCE+NBUSE+I)= 0.D0
          AK_SCE%R(NPTSCE+I)= 0.D0
          AK_SCE%R(NPTSCE+NBUSE+I)= 0.D0
          EP_SCE%R(NPTSCE+I)= 1.D-6
          EP_SCE%R(NPTSCE+NBUSE+I)= 1.D-6
        ENDDO
        IF(NTRAC.GT.0) THEN
          DO I=1,NBUSE
            DO ITRAC=1,NTRAC
              TA_SCE%ADR(ITRAC)%P%R(NPTSCE+I) = TBUS%ADR(ITRAC)%P%R(I)
              TA_SCE%ADR(ITRAC)%P%R(NPTSCE+NBUSE+I)=
     &                                      TBUS%ADR(ITRAC)%P%R(NBUSE+I)
            ENDDO
          ENDDO
        ENDIF
      ENDIF
      IF(DEBUG.GT.0) WRITE(LU,*) 'PREMIER APPEL DE SOURCES_SINKS'
      CALL SOURCES_SINKS
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE SOURCES_SINKS'
!
      IF(INCLUS(COUPLING,'TOMAWACT3D')) THEN
!       ADDING THE NEW TERMS ASSOCIATED WITH THE USTOKES AND VSTOKES
!       INTEGRATED OVER DEPTH
        CALL OS ('X=YZ     ',X=DUS2D,Y=US2D,Z=H)
        CALL VECTOR(DUS2DX,'=','GRADF          X',IELMH,1.D0,DUS2D,
     &              SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,.FALSE.,MASKEL,
     &              ASSPAR=.TRUE.)
!
        CALL OS ('X=YZ     ',X=DVS2D,Y=VS2D,Z=H)
        CALL VECTOR(DVS2DY,'=','GRADF          Y',IELMH,1.D0,DVS2D,
     &              SVIDE,SVIDE,SVIDE,SVIDE,SVIDE,MESH2D,.FALSE.,MASKEL,
     &              ASSPAR=.TRUE.)
!
        CALL OS('X=X-Y      ',X=SMH,Y=DUS2DX)
        CALL OS('X=X-Y      ',X=SMH,Y=DVS2DY)
      ENDIF
!
      IF(BEDBOU)THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'PREMIER APPEL DE BED_FLUXES'
        CALL BED_FLUXES
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE BED_FLUXES'
      ENDIF
!
!     FLULIM NOT INITIALISED AND USED IN PREADV, THROUGH FLUX3D
!
      IF(OPT_HNEG.EQ.2) THEN
        DO I=1,MESH2D%NSEG
          FLULIM%R(I)=1.D0
        ENDDO
      ENDIF
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'PREMIER APPEL DE PREADV'
!
      CALL PREADV(W,WS,ZPROP,ISOUSI,LT,VOLU,VOLUN)
!
      IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DU PREMIER APPEL DE PREADV'
!
!     NOW SETTING VOLUN=VOLU (IN CASE OF COMPUTATION CONTINUED IT HAS
!     BEEN RETRIEVED FROM THE PREVIOUS COMPUTATION)
!
      IF(.NOT.DEBU.AND.TROUVE(17).EQ.1) THEN
        CALL OS('X=Y     ',X=VOLUN,Y=VOLU)
      ENDIF
!
#if defined COMPAD
      CALL AD_TELEMAC3D_INITIALISATION_END
#endif
!
!=======================================================================
!     COUPLING WITH TOMAWAC
!=======================================================================
!
      IF(INCLUS(COUPLING,'TOMAWAC')) THEN
!
        WRITE (LU,*) 'TELEMAC-3D: INTERNAL COUPLING WITH TOMAWAC'
        CALL CONFIG_CODE(3)
        IF(DEBUG.GT.0) WRITE(LU,*) 'PREMIER APPEL DE TOMAWAC'
        CALL T3D_WAC_CPL_UPDATE(0)
        CALL WAC(PART=0)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE TOMAWAC'
        CALL CONFIG_CODE(1)
!
      ENDIF
!
!=======================================================================
!     COUPLING WITH GAIA OR SISYPHE
!=======================================================================
!
!     WRITES THE INITIAL CONDITIONS FOR U(Z=0), V(Z=0) AND H
!
      IF(COUPLING.NE.' ') THEN
        WRITE(LU,*) 'TELEMAC3D COUPLED WITH: ',COUPLING
      ENDIF
!
!     COUPLING WITH GAIA
!
      IF(INCLUS(COUPLING,'GAIA   ')) THEN
!
!       U AND V WITH 2D STRUCTURE : BOTTOM VELOCITY AS A 2D VARIABLE
        CALL CPSTVC(U2D,U)
        CALL CPSTVC(V2D,V)
!
        CALL CONFIG_CODE(5)
!       INPUT VARIABLES IN GAIA CANNOT BE HARDCODED
        IF(DEBUG.GT.0) WRITE(LU,*) 'CALL TO GAIA_INIT'
        CALL GAIA_INIT(
     &   GRAPRD,LISPRD,NIT,U2D,V2D,H,ZF,UETCAR,DELTAR,CF,RUGOF,CODE1,
     &   U,V,AT,DT,CHARR_TEL,SUSP_TEL,
     &   DIRMOY,HM0,TPR5,ORBVEL,.TRUE.,
     &   RHO0,GRAV)
        IF(DEBUG.GT.0) WRITE(LU,*) 'BACK FROM GAIA_INIT'
        CALL CONFIG_CODE(1)
!
!       RETRIEVES ORIGINAL U AND V STRUCTURE
!
        CALL CPSTVC(UN,U)
        CALL CPSTVC(VN,V)
!
      ENDIF
!
!     COUPLING WITH SISYPHE
!
      IF(INCLUS(COUPLING,'SISYPHE')) THEN
!
!       U AND V WITH 2D STRUCTURE : BOTTOM VELOCITY AS A 2D VARIABLE
        CALL CPSTVC(U2D,U)
        CALL CPSTVC(V2D,V)
!
        CALL CONFIG_CODE(2)
!       INOUT VARIABLES IN SISYPHE CANNOT BE HARDCODED
        IBID=1
        LBID=.FALSE.
        LBID2=.FALSE.
        IF(DEBUG.GT.0) WRITE(LU,*) 'PREMIER APPEL DE SISYPHE'
        CALL SISYPHE(0,LT,GRAPRD,LISPRD,NIT,U2D,V2D,H,H,ZF,UETCAR,
     &               CF,RUGOF,
     &               LBID,IBID,LBID2,CODE1,PERCOU_SIS,
     &               U,V,AT,VISCVI,DT,S3D_CHARR,S3D_SUSP,
!                          1 PRECLUDES THE USE OF THE FOLLOWING ARGUMENTS
     &               FLBOR,1,DM1,UCONV,VCONV,ZCONV,DIRMOY,HM0,TPR5,
     &               ORBVEL,.TRUE.)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DU PREMIER APPEL DE SISYPHE'
        CALL CONFIG_CODE(1)
!
!       RETRIEVES ORIGINAL U AND V STRUCTURE
!
        CALL CPSTVC(UN,U)
        CALL CPSTVC(VN,V)
!
      ENDIF
!
!=======================================================================
!
!     INITIALISES THE SEDIMENT SETTLING VELOCITY
!     NEGLECTS TURBULENCE HERE
!     S3D_WCHUCOMPUTED HERE IS USED IN BORD3D FOR ROUSE PROFILES
!     CV:Floculation and  hindered settling
!        Soulsby flocculation and hindered settling is now coded in VITCHU
!        along with all the previous telemac methods
!        NOTE: moved settling calculation to after CLSEDI since
!        soulsby floc requires bed shear stress
!
      IF(NSUSP_TEL.GT.0) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE COMPUTE_SETTLING_VEL'
        CALL COMPUTE_SETTLING_VEL(WCHU,U,V,
     &  T3_01,T3_02,T3_03,SVIDE,MESH3D,IELM3,NPOIN2,NPOIN3,
     &  NPLAN,MSK,MASKEL,UETCAR,TA,HN)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE COMPUTE_SETTLING_VEL'
      ENDIF
!
      IF(S3D_SEDI) THEN
        IF(DEBUG.GT.0) WRITE(LU,*) 'APPEL DE VITCHU'
        CALL VITCHU(S3D_WCHU,S3D_WCHU0,U,V,
     &              S3D_TURBA,S3D_TURBB,T3_01,T3_02,T3_03,SVIDE,MESH3D,
     &              IELM3,NPOIN2,NPOIN3,NPLAN,NTRAC,MSK,MASKEL,
     &              UETCAR,TA,HN,S3D_FLOC, S3D_FLOC_TYPE,
     &              S3D_HINDER,S3D_HIND_TYPE,S3D_CGEL,S3D_CINI)
        IF(DEBUG.GT.0) WRITE(LU,*) 'RETOUR DE VITCHU'
!
        IF(S3D_MIXTE) THEN
          CALL OS( 'X=C     ' , X=S3D_WCS, C=S3D_WCS0)
        ENDIF
      ENDIF
!
!=======================================================================
!
      IF( (INCLUS(COUPLING,'WAQTEL') ).AND.
     &    (13*INT(WAQPROCESS/13).EQ.WAQPROCESS) ) THEN
#if defined HAVE_AED2
        CALL OS('X=0     ',X=THICK4AED2)
!       ORDER: NUMBER OF THE PLANES FROM TOP TO BOTTOM
!       THEN NUMBER OF THE 2D NODES
        DO I=1,NPOIN2
!         K=1
!         THICK4AED2%R(1+(I-1)*NPLAN) = (Z(I+(NPLAN-1)*NPOIN2)
!    &                                - Z(I+(NPLAN-2)*NPOIN2))*0.5D0
          THICK4AED2%R(1+(I-1)*NPLAN) = (Z(I+NPOIN3-NPOIN2)
     &                                - Z(I+NPOIN3-2*NPOIN2))*0.5D0
          DO K=2,NPLAN-1
            THICK4AED2%R(K+(I-1)*NPLAN) = (Z(I+(NPLAN-K+1)*NPOIN2)
     &                                  - Z(I+(NPLAN-K-1)*NPOIN2))/2.D0
          ENDDO
!         K=NPLAN
          THICK4AED2%R(I*NPLAN) = (Z(I+NPOIN2)-Z(I))*0.5D0
        ENDDO
!
        DO I=1,NPOIN2
          BENTH4AED2%I(I) = NPLAN
          DO K=NPLAN,2,-1
            IF(THICK4AED2%R(K+(I-1)*NPLAN).LT.EPS_AED2) THEN
              BENTH4AED2%I(I) = K-1
            ENDIF
          ENDDO
        ENDDO
!
!       ORDER: NUMBER OF THE PLANES FROM TOP TO BOTTOM
!       THEN NUMBER OF THE 2D NODES
        DO I=1,NPOIN2
          DO K=1,NPLAN
            TEMP4AED2%R(K+(I-1)*NPLAN) =
     &      TA%ADR(IND_T)%P%R(I+(NPLAN-K)*NPOIN2)
            SALI4AED2%R(K+(I-1)*NPLAN) =
     &      TA%ADR(IND_S)%P%R(I+(NPLAN-K)*NPOIN2)
          ENDDO
        ENDDO
!
!       INITIALISES RAYAED2 BEFORE THE FIRST USE IN T3D_CONSTRUCT_AED2
        CALL SOLRAD(RAY_SOL,NEBU,MARDAT,MARTIM,AT,LATIT,LONGIT)
        CALL OS('X=C     ', X=RAYAED2, C=RAY_SOL)
!
        CALL T3D_CONSTRUCT_AED2(DT,TA,RHO4AED2,V2DPAR,NPLAN,
     &                          SALI4AED2,TEMP4AED2,WINDSPD,NPOIN2,
     &                          NPOIN3,THICK4AED2,BENTH4AED2,PLUIE)
#else
        WRITE(LU,*) 'ERROR: AED2 LIBRARY NOT COMPILED'
        CALL PLANTE(1)
#endif
      ENDIF
      END SUBROUTINE
