telemac3d.F

Go to the documentation of this file.

!                   ********************
                    SUBROUTINE telemac3d
!                   ********************
!
     &(pass, nit_ori)
!
!***********************************************************************
! TELEMAC3D   V8P2
!***********************************************************************
!
!brief
!
!history  JACEK A. JANKOWSKI PINXIT
!+        01/03/1999
!+
!+   FORTRAN95 VERSION
!
!history  J-M HERVOUET (LNHE)
!+        05/05/2010
!+        V6P0
!+   K-OMEGA MODEL BY HOLGER WEILBEER (ISEB/UHA)
!
!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  J-M HERVOUET (LNHE)
!+        02/08/2011
!+        V6P1
!+   CALL MITTIT(18,AT,LT) changed into CALL MITTIT(19,AT,LT)
!+   CALL MITTIT(19,AT,LT) changed into CALL MITTIT(20,AT,LT)
!+   2 fractional steps were not correctly labelled in the listing
!
!history  J-M HERVOUET (LNHE)
!+        12/08/2011
!+        V6P2
!+   Calls to CHECK and BIL3D changed
!
!history  J-M HERVOUET (LNHE)
!+        02/04/2012
!+        V6P2
!+   Clean restart implemented.
!
!history  J-M HERVOUET (LNHE)
!+        01/06/2012
!+        V6P2
!+   Call to vector before call to Tel4del corrected (GRAZCO)
!+   Initialisation of TAN after call to condim.
!
!history  J-M HERVOUET (LNHE)
!+        18/129/2012
!+        V6P3
!+   Call to IFAB3DT added, arguments of cstkep removed.
!
!history  J-M HERVOUET (LNHE)
!+        25/01/2013
!+        V6P3
!+   TAN renamed TRN, copy of TRN on TA moved from after CONDIM to
!+   after BIEF_SUITE, FLULIM set to 1 before first call to PREADV
!
!history  J-M HERVOUET (LNHE)
!+        11/03/2013
!+        V6P3
!+   Call to METEO modified. Stop if variables not found for a 2D
!+   continuation.
!
!history  R. KOPMANN (EDF R&D, LNHE)
!+        16/04/2013
!+        V6P3
!+   Adding the file format in calls to FIND_IN_SEL.
!
!history  J-M HERVOUET (LNHE)
!+        18/03/2013
!+        V6P3
!+   Dealing with the newly created FILE FOR 2D CONTINUATION.
!
!history  J-M HERVOUET (LNHE)
!+        25/04/2013
!+        V6P3
!+   AKN and EPN initialised in case of computation continued, for the
!+   first call to PREADV.
!+   Mesh better updated in case of coupling with Sisyphe.
!
!history  J-M HERVOUET (LNHE)
!+        20/09/2013
!+        V6P3
!+   CALL PLANE_BOTTOM added at the beginning of time loop (otherwise
!+   when calling kepcl3 IPBOT is done with ZPROP at the first iteration
!+   and with Z for the others, while ZPROP is always sent as argument.
!+   This could trigger unexpected divisions by 0.
!
!history  J-M HERVOUET (LNHE)
!+        15/11/2013
!+        V6P3
!+   After second call to bief_suite, checking that Z has been found,
!+   otherwise stop
!
!history  C. VILLARET & T. BENSON & D. KELLY (HR-WALLINGFORD)
!+        27/02/2014
!+        V7P0
!+   New developments in sediment merged on 25/02/2014.
!
!history  J-M HERVOUET (LNHE)
!+        14/03/2014
!+        V7P0
!+   CALL BIL3D put out of the IF(S3D_SEDI) test. Address of depth-averaged
!+   tracers from 38 to 37+NTRAC in ALIRE2D.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        19/03/2014
!+        V7P0
!+   Boundary segments have now their own numbering, independent of
!+   boundary points numbering. Differents calls changed accordingly.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        02/05/2014
!+        V7P0
!+   Argument ZR added to FONVAS. S3D_HDEPupdated differently after calling
!+   Sisyphe, to avoid truncation errors that would give S3D_HDEP<0.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        31/07/2014
!+        V7P0
!+   Call to METEO moved just before the first call to FSGRAD, not just
!+   after (atmospheric pressure gradients now systematically added
!+   to free surface gradients).
!
!history  C VILLARET (HRW+EDF) & J-M HERVOUET (EDF - LNHE)
!+        18/09/2014
!+        V7P0
!+   Calls to sisyphe and wac chenged.
!
!history  G. ANTOINE & M. JODEAU & J.M. HERVOUET (EDF - LNHE)
!+        13/10/2014
!+        V7P0
!+   New developments in sediment for mixed sediment transport
!
!history  R. ATA (EDF LAB, LNHE)
!+        05/11/2014
!+        V7P0
!+   add optional variables to meteo in a sake of harmonization
!+   with telemac-2d
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        31/03/2015
!+        V7P1
!+   Just a few extra debugger prints, up to CALL KEPINI, where was the
!+   last user bug I looked for.
!
!history Y AUDOUIN (LNHE)
!+       25/05/2015
!+       V7P0
!+       Modification to comply with the hermes module
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        26/06/2015
!+        V7P1
!+   2D and 3D RESULT FILE can be optional. Tests for writing them added.
!
!history  A. JOLY (EDF LAB, LNHE)
!+        27/08/2015
!+        V7P1
!+   Imposed flowrates on the bed.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        21/01/2016
!+        V7P1
!+   Initial conditions of results file in restart mode were forgotten
!+   they are useless but now requested by the Hermes module...
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        08/02/2016
!+        V7P2
!+   Adding the argument HPROP in the call to Sisyphe.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        24/03/2016
!+        V7P2
!+   Adapting to new CVDF3D, saving Z at time T(n) in ZN.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        27/05/2016
!+        V7P2
!+   Allowing k-epsilon model on a direction and not on the other.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        30/05/2016
!+        V7P2
!+   In a continued computation, K, Epsilon and the dynamic pressure
!+   must not be read if they are not necessary because their arrays
!+   are not allocated.
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        22/08/2016
!+        V7P2
!+   Adding TB2 in the call to CVDF3D.
!
!history  M.JODEAU (EDF LAB, LNHE)
!+        08/2016
!+        V7P3
!+   Water quality: AED2 coupling
!
!history  R. ATA (EDF LAB, LNHE)
!+        12/01/2017
!+        V7P3
!+   Bug fix in the call of DERIVE: now zchar, nplan, transf are well
!    introduced
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        11/09/2017
!+        V7P3
!+   Adding NELMAX2 in list of arguments of ifab3d, ifab3dt, TBORD,
!+   flux_ef_vf_3d and make_zconv.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_telemac3d
      USE meteo_telemac !, ONLY: WINDX,WINDY,PATMOS,TAIR
      USE declarations_waqtel, ONLY : waqprocess,ind_sf,ind_csf
      USE oilspill
      USE interface_telemac3d, ex_telemac3d => telemac3d
      USE interface_telemac2d
      USE interface_sisyphe, ONLY: sisyphe
      USE interface_gaia, ONLY: gaia_step
      USE interface_tomawac, ONLY: wac
      USE tel4del, ONLY: tel4delwaq
      USE declarations_gaia, ONLY: nsicla,fludp,fludpt,fluer,nsusp_tel,
     &                       xmvs0,num_isusp_icla,setdep
      USE declarations_waqtel,ONLY: waqprocess
      USE gotm_coupling
!
      USE declarations_special
      IMPLICIT NONE
!
!-----------------------------------------------------------------------
! DECLARES LOCAL VARIABLES FOR TELEMAC3D
!-----------------------------------------------------------------------
!
      INTEGER,          INTENT(IN) :: PASS
      INTEGER,          INTENT(IN) :: NIT_ORI
!
!-----------------------------------------------------------------------
!
      INTEGER ITRAC,ISOUSI,IPOIN
      INTEGER SCHDVI_HOR,SCHDVI_VER,SCHCVI_HOR,SCHCVI_VER
      INTEGER IBID,I,K,I3D,IP
!
!
      DOUBLE PRECISION TETADIVER
      DOUBLE PRECISION UMIN,  UMAX,  SIGMAU, VMIN,  VMAX, SIGMAV
      DOUBLE PRECISION WMIN,  WMAX,  SIGMAW
      DOUBLE PRECISION TAMIN, TAMAX, SIGMTA,TETATRA
!
      LOGICAL CLUMIN, CLUMAX, CLVMIN, CLVMAX, CLWMIN, CLWMAX
      LOGICAL CTAMIN, CTAMAX, YASEM3D,YAS0U,YAS1U
      LOGICAL CLKMIN, CLKMAX, CLEMIN, CLEMAX, CLNUMIN,CLNUMAX
      LOGICAL YAWCHU,NEWDIF,LBID,LBID2,BC
!
      CHARACTER(LEN=24), PARAMETER :: CODE1='TELEMAC3D               '
      CHARACTER(LEN=16) FORMUL
!
      INTRINSIC mod
!
      TYPE(slvcfg) :: SLVD
!
      DOUBLE PRECISION, POINTER, DIMENSION(:) :: SAVEZ
      DOUBLE PRECISION, POINTER, DIMENSION(:,:) :: FLUXBIL

      TYPE(bief_obj), POINTER :: MULTI1, MULTI2, MULTI3, MULTI4
!
!=======================================================================
!
!
! TODO: temporary !!! Remove also in deall_telemac3d
      IF(.NOT.inclus(coupling,'GAIA  ').AND.pass.LE.0) THEN
        isusp=1
        setdep=-1
        nsusp_tel=0
!
        CALL allblo(fludp, 'FLUDP ')
        CALL allblo(fludpt,'FLUDPT')
        CALL allblo(fluer, 'FLUER ')
!
        CALL bief_allvec_in_block(fludp,1,1,'FLUDP ',ielm2h,1,2,mesh2d)
        CALL bief_allvec_in_block(fludpt,1,1,'FLUDPT',ielm2h,1,2,mesh2d)
        CALL bief_allvec_in_block(fluer,1,1,'FLUER ',ielm2h,1,2,mesh2d)
!
        ALLOCATE(num_isusp_icla(1))
        num_isusp_icla(1)=1
!
      ENDIF
! temporary !!!
#if defined COMPAD
      CALL ad_telemac3d_begin
#endif
!
!     READS TRACERS IN PREVIOUS FILES
!
      IF(ntrac.GT.0) THEN
        DO i=adr_trac,adr_trac+ntrac-1
          alire3d(i)=1
        ENDDO
        DO i=adr_trac_2d,adr_trac_2d+ntrac-1
          alire2d(i)=1
        ENDDO
!
!       THIS IS S3D_ESOMT...
!       IF(S3D_SEDI) ALIRE2D(37)=1
!
      ENDIF
!
!     DO NOT READ K AND EPSILON IF NOT NECESSARY
!
      IF(iturbh.NE.3.AND.iturbv.NE.3.AND.iturbv.NE.6.AND.
     &   iturbh.NE.7.AND.iturbv.NE.7) THEN
        alire3d(8)=0
        alire3d(9)=0
      ENDIF
!
!     DO NOT READ DYNAMIC PRESSURE IF NOT NECESSARY
!
      IF(.NOT.nonhyd) alire3d(12)=0
!
!=======================================================================
! FOR DROGUES (CALLS TO FLOT3D WILL INCREASE OR DECREASE NFLOT)
!=======================================================================
!
      nflot=0
!
!=======================================================================
! FOR COMPUTING FLUXES OF ADVECTED VARIABLES
!=======================================================================
!
!     NO FLUX COMPUTED FOR U,V,W,K,EPSILON
      DO i=1,5
        calcflu(i)=.false.
      ENDDO
!     DEPENDING ON BILMAS FOR TRACERS
      IF(ntrac.GT.0) THEN
        DO i=6,5+ntrac
          calcflu(i)=bilmas
        ENDDO
      ENDIF
!
!=======================================================================
! FOR TAKING INTO ACCOUNT RAIN IN ADVECTION OF VARIOUS VARIABLES
!=======================================================================
!
!     NO RAIN FOR U,V,W,K,EPSILON
      DO i=1,5
        calcrain(i)=.false.
      ENDDO
!     DEPENDING OF RAIN FOR TRACERS
      IF(ntrac.GT.0) THEN
        DO i=6,5+ntrac
          calcrain(i)=rain
        ENDDO
      ENDIF
!
!=======================================================================
! INITIALISATION: READS, PREPARES AND CHECKS
!=======================================================================
!
      IF(pass.EQ.0) THEN
        WRITE(lu,*) 'INITIALISING TELEMAC3D FOR ',code1
        WRITE(lu,*) 'INITIALISING TELEMAC3D'
      ELSEIF(pass.EQ.1) THEN
        GO TO 700
      ELSEIF(pass.NE.-1) THEN
        WRITE(lu,*) 'WRONG ARGUMENT PASS: ',pass
        CALL plante(1)
        stop
      ENDIF

      CALL telemac3d_init
!
!=======================================================================
! THE TIME LOOP BEGINS HERE
!=======================================================================
!
      IF(pass.EQ.0) THEN
        WRITE(lu,*) 'TELEMAC3D INITIALISED'
        RETURN
      ENDIF
!
700   CONTINUE
!
      lt = lt+1

      at = at + dt
!
#if defined COMPAD
      CALL ad_telemac3d_timestep_begin
#endif
!
!     SAVING ORIGINAL ELEVATIONS (FOR DISTRIBUTIVE SCHEMES)
!
      CALL os('X=Y     ',x=zn,y=z3)
!
      IF(debug.GT.0) WRITE(lu,*) 'BOUCLE EN TEMPS LT=',lt
      infogr = .false.
      IF (mod(lt,lisprd) == 0) infogr = .true.
      infogr = listin .AND. infogr
      IF (infogr) CALL mittit(1,at,lt)
!
!=======================================================================
!
!     IPBOT HAS BEEN MODIFIED FOR CVDF3D IN THE PREVIOUS TIME STEP,
!     IT IS RESTORED HERE WITH ZPROP
!     NOTE: DIFFERENT IPBOT_Z AND IPBOT_ZPROP WOULD BE CLEARER....
      IF(lt.GT.1) THEN
        CALL plane_bottom(ipbot%I,zprop%R,npoin2,nplan,sigmag,optban)
      ENDIF
!
!=======================================================================
! SOURCES : COMPUTES INPUTS WHEN VARYING IN TIME
!           IF NO VARIATION IN TIME QSCE2=QSCE AND TASCE2=TASCE
!=======================================================================
!
      IF(nptsce.GT.0) THEN
        DO i=1,nptsce
          qsce2(i)=t3d_debsce(at,i,qsce)
        ENDDO
        IF(ntrac.GT.0) THEN
          DO i=1,nptsce
            DO itrac=1,ntrac
              ta_sce%ADR(itrac)%P%R(i)=t3d_trsce(at,i,itrac)
            ENDDO
          ENDDO
        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,infogr,
     &            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)
        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
!
!=======================================================================
! END OF CAMILLE LEQUETTE'S MODIFICATIONS
!=======================================================================
!
!     COUPLING WITH TOMAWAC
!
      IF(inclus(coupling,'TOMAWAC').AND.
     &   percou_wac*((lt-1)/percou_wac).EQ.lt-1) THEN
!
        CALL config_code(3)
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TOMAWAC'
        CALL t3d_wac_cpl_update(nit_ori)
        CALL wac(part=1)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TOMAWAC'
        CALL config_code(1)
!
      ENDIF
!
!=======================================================================
!
! SAVES H, TA, TP, AK, EP
! IN    HN,TRN,TPN,AKN,EPN
!
      CALL os ( 'X=Y     ', x=hn,    y=h     )
      CALL os ( 'X=Y     ', x=volun, y=volu  )
      IF(ncsize.GT.1) CALL os('X=Y     ',x=volunpar,y=volupar)
      CALL os ( 'X=Y     ', x=un,    y=u     )
      CALL os ( 'X=Y     ', x=vn,    y=v     )
      IF(nonhyd) CALL os ( 'X=Y     ' , x=wn, y=w)
      CALL os ( 'X=Y     ', x=gradzn,y=gradzs)
!     TRACERS (IF LT=1 DONE AFTER CALL CONDIM AND READ_DATASET)
      IF(ntrac.GT.0.AND.lt.GT.1) CALL os ('X=Y     ', x=trn, y=ta)
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.iturbv.EQ.7.OR.iturbh.EQ.7) THEN
        CALL os ( 'X=Y     ', x=akn, y=ak )
        CALL os ( 'X=Y     ', x=epn, y=ep )
      ENDIF
      IF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
        CALL os( 'X=Y      ', x=nun, y=nu )
      ENDIF
!
      IF(bilmas) THEN
        massen_water = masse_water
        CALL os ( 'X=Y     ', x=massen, y=masse )
      ENDIF
!
! COMPUTES MEAN UN AND VN IN THE VERTICAL
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE VERMOY'
      CALL vermoy(un2d%R,vn2d%R,un%R,vn%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'
!
! INTEGRATES USTOKES AND VSTOKES OVER DEPTH
!
      IF(inclus(coupling,'TOMAWACT3D')) THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE VERMOY POUR VIT STOKES'
          CALL vermoy(us2d%R,vs2d%R,ustokes%R,vstokes%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 POUR VIT STOKES'
      ENDIF
!
!-----------------------------------------------------------------------
!
! COMPUTES FRICTION COEFFICIENT
!
!     TIME VARIATIONS OF RUGOF (CORSTR IS IN TELEMAC-2D LIBRARY)
!     MUST BE USER-IMPLEMENTED - NOTHING DONE IN STANDARD
      CALL corstr
!
      IF(.NOT.inclus(coupling,'TOMAWACT3D')) THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE COEFRO'
        CALL coefro(cf,h,un2d,vn2d,karman,kfrot,rugof,grav,mesh2d,t2_01,
     &              .false.)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE COEFRO'
      ELSE
!       FRICTION COEFFICIENT MODIFIED TO TAKE INTO ACCOUNT WAVES+CURRENTS
        DO ip=1,npoin2
          cf%R(ip) = cfwc%R(ip)
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
! CHECKS AND HARMONISES THE BOUNDARY CONDITION TYPES
!
      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'
!
! BOUNDARY CONDITIONS FOR THE K-EPSILON MODEL
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.iturbv.EQ.7.OR.iturbh.EQ.7) THEN
        CALL kepicl(likbof%I,liebof%I,liubof%I,
     &              likbol%I,liebol%I,liubol%I,
     &              likbos%I,liebos%I,
     &              nptfr2,nplan,npoin2,kent,ksort,kentu)
      ENDIF
!
! BOUNDARY CONDITIONS FOR THE S-A MODEL
!
      IF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SAPICL'
        CALL sapicl(linubof%I, liubof%I,
     &              linubol%I, liubol%I,
     &              linubos%I,
     &              nptfr2, nplan, npoin2, kent, ksort)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SAPICL'
      ENDIF
!-----------------------------------------------------------------------
! FORCING AT THE BOUNDARIES
!
! METEOROLOGICAL CONDITIONS
!
      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
!
!-----------------------------------------------------------------------
!
!     SEDIMENT
!
      IF(s3d_sedi) THEN
!
!       COMPUTES THE SEDIMENT SETTLING VELOCITY
!
        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)
!
!       BOUNDARY CONDITIONS FOR THE SEDIMENTOLOGY
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CLSEDI'
!GA: CLSEDI HAS BEEN MODIFIED TO ALLOW TWO SEDIMENT CLASSES
        IF(s3d_mixte) THEN

          DO itrac = ntrac-1,ntrac

            CALL clsedi
     &       (atabof%ADR(itrac)%P%R,btabof%ADR(itrac)%P%R,
     &        s3d_wchu%R,
     &        z, h, deltar%R, t3_01, t3_02%R,
     &        s3d_epai, s3d_cfdep,s3d_conc, s3d_hdep%R, s3d_fluer%R,
     &        s3d_fludpt%R, litabf%ADR(itrac)%P%I,
     &        klog, npoin3, npoin2, nplan, s3d_ncouch,
     &        dt, rho0, s3d_rhos,
     &        s3d_tocd,s3d_mpart,s3d_toce,uetcar%R,
     &        grav,s3d_sedco,s3d_dmoy,s3d_cref,zref,cf,s3d_ac,
     &        s3d_kspratio,s3d_icr,s3d_icq,
     &        rugof,s3d_setdep,s3d_hsed,
     &        s3d_wcs%R, s3d_epaico%R, s3d_epainco%R, s3d_mixte,
     &        s3d_sednco, s3d_fludptc%R, s3d_fludptnc%R, s3d_fluerc%R,
     &        s3d_fluernc%R, ntrac, itrac)
!           ATABOF AND BTABOF ARE NO LONGER 0 FOLLOWING CLSEDI
            atabof%ADR(itrac)%P%TYPR='Q'
            btabof%ADR(itrac)%P%TYPR='Q'

          ENDDO

        ELSEIF (s3d_sedco.OR.s3d_sednco) THEN

        CALL clsedi
     &   (atabof%ADR(ntrac)%P%R,btabof%ADR(ntrac)%P%R,
     &    s3d_wchu%R,
     &    z, h, deltar%R, t3_01, t3_02%R,
     &    s3d_epai, s3d_cfdep,s3d_conc, s3d_hdep%R, s3d_fluer%R,
     &    s3d_fludpt%R, litabf%ADR(ntrac)%P%I,
     &    klog, npoin3, npoin2, nplan, s3d_ncouch,
     &    dt, rho0, s3d_rhos,
     &    s3d_tocd,s3d_mpart,s3d_toce,uetcar%R,
     &    grav,s3d_sedco,s3d_dmoy,s3d_cref,zref,cf,s3d_ac,s3d_kspratio,
     &    s3d_icr,s3d_icq,rugof,s3d_setdep,s3d_hsed,
     &    s3d_wcs%R, s3d_epaico%R, s3d_epainco%R, s3d_mixte, s3d_sednco,
     &    s3d_fludptc%R, s3d_fludptnc%R, s3d_fluerc%R, s3d_fluernc%R,
     &    ntrac, itrac)
!
!         ATABOF AND BTABOF ARE NO LONGER 0 FOLLOWING CLSEDI
          atabof%ADR(ntrac)%P%TYPR='Q'
          btabof%ADR(ntrac)%P%TYPR='Q'
!
        ELSE
!
          WRITE(lu,*) ' '
          WRITE(lu,*) 'SEDI3D : ERROR ON SEDIMENT KEY WORD'
          CALL plante(1)
          stop
!
        ENDIF

        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CLSEDI'
      ENDIF
!
!     PREPARING BOUNDARY CONDITIONS FOR THOMPSON METHOD
!
      IF(thomfr.AND.nfrliq.GT.0) THEN
!
        CALL cpstvc(h,t2_01)
        CALL prebor(hbor%R,ubor2d%R,vbor2d%R,taborl,u2d%R,v2d%R,h%R,
     &              t2_01%R,ta,mesh2d%NBOR%I,
     &              mesh2d%NPOIN,mesh2d%NPTFR,
!    &              NTRAC SET TO ZERO PROVISIONALLY
     &              0    ,nfrliq,frtype,numliq%I)
!       RESTORING USER BOUNDARY CONDITIONS BEFORE CALLING BORD3D
!       TO AVOID UNDUE CALLS TO SL3, ETC.
        CALL thomps_bc(3)
!
      ENDIF
!
! UPDATES BOUNDARY CONDITION VALUES
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BORD3D'
      CALL bord3d(nfrliq)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BORD3D'
!
      IF(thomfr.AND.nfrliq.GT.0) THEN
!
!       NOW THAT BORD3D HAS BEEN CALLED
!       CHANGING AGAIN BOUNDARY CONDITIONS FOR THOMPSON
        CALL thomps_bc(2)
!
!       UBORL/VBORL UPDATED IN BORD3D, SO ARE UBOR2D/VBOR2D
        DO i=1,nptfr2
          ubor2d%R(i) = uborl%R(i)
          vbor2d%R(i) = vborl%R(i)
        ENDDO
!
      ENDIF
!
! BOUNDARY CONDITIONS FOR THE VELOCITY ON LATERAL BOUNDARIES
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TBORD'
      CALL tbord(auborl%R,
     &           rugol%R,
     &           mesh2d%DISBOR%R,mesh2d%NELBOR%I,mesh2d%NULONE%I,
     &           mesh2d%IKLE%I,nelmax2,
     &           u%R,v%R,w%R,
     &           nbor2%I,npoin2,nplan,nptfr2,dnuvih,dnuviv,
     &           karman,lisrul,kfrotl,
     &           uetcal%R,nonhyd,
     &           t2_02%R,mesh2d)
      IF(kfrotl.EQ.0) THEN
        auborl%TYPR='0'
      ELSE
        auborl%TYPR='Q'
      ENDIF
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TBORD, APPEL DE TFOND'
!
! BOUNDARY CONDITIONS FOR THE VELOCITY ON THE BOTTOM
!
      CALL tfond(auborf%R,
     &           cf%R,un2d%R,vn2d%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)
      auborf%TYPR='Q'
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TFOND'
!
! BOUNDARY CONDITIONS FOR K-EPSILON MODEL + COMPUTES CONSTRAINTS
! AT THE BOTTOM AND LATERAL BOUNDARIES IF K-EPSILON IS REQUIRED
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3) THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE KEPCL3'
        CALL kepcl3(kborf%R,eborf%R,likbof%I,liebof%I,
     &              kborl%R,eborl%R,likbol%I,liebol%I,liubol%I,
     &              kbors%R,ebors%R,
     &              likbos%I,liebos%I,
     &              mesh2d%DISBOR%R,ak%R,h%R,zprop%R,
     &              nbor2%I,npoin2,nplan,nptfr2,
     &              karman,cmu,
     &              kmin,emin,
     &              kent,kentu,ksort,kadh,klog,
     &              uetcar%R,fict)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE KEPCL3'
!
      ELSEIF(iturbv.EQ.7.OR.iturbh.EQ.7) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE KOMCL3'
        CALL komcl3(kborf%R,eborf%R,likbof%I,liebof%I,liubof%I,
     &              kborl%R,eborl%R,likbol%I,liebol%I,liubol%I,
     &              ebors%R,liebos%I,
     &              mesh2d%DISBOR%R,ak%R,
     &              u%R,v%R,h%R,zprop%R,
     &              nbor2%I,npoin2,nplan,nptfr2,
     &              karman,betas,omstar,
     &              kmin,emin,
     &              kentu,kent,kadh,klog,uetcar%R,uetcal%R)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE KOMCL3'
!
      ELSEIF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SACL3'
        CALL sacl3(nuborf%R,linubof%I,
     &             nuborl%R,linubol%I,liubol%I,
     &             h%R,zprop%R,
     &             nbor2%I,npoin2,nplan,nptfr2,
     &             karman,uetcar%R,numin,
     &             kent,kentu,ksort,kadh,klog,fict)
!
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SACL3'
      ENDIF
!
! CLIPS HBOR
!
      IF(optban.EQ.2) THEN
        CALL clip(hbor,hmin,.true.,1.d6,.false.,0)
      ENDIF
!
!-----------------------------------------------------------------------
!
!     THOMPSON BOUNDARY CONDITIONS
!
      IF(thomfr) THEN
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE THOMPS'
!     T2_01 IS HERE A VALUE OF H SAVED BEFORE CALLING BORD3D
      CALL thomps(hbor%R,ubor2d%R,vbor2d%R,taborl,u2d,v2d,t2_01,
     &            ta,zf,mesh2d%X%R,mesh2d%Y%R,mesh2d%NBOR%I,
     &            frtype,t2_02,t2_03,t2_04,t2_06,t2_07,
     &            lihbor%I,liubol%I,livbol%I,it1%I,
     &            t2_08%R,t2_09%R,w1%R,t2_21,t2_22,t2_23,
     &            tac,t2_10,mesh2d%SURDET%R,mesh2d%IKLE%I,
     &            mesh2d%IFABOR%I,mesh2d%NELEM,mesh2d,
     &            mesh2d%XNEBOR%R,mesh2d%YNEBOR%R,
!                                                   NTRAC SET TO ZERO
     &            mesh2d%NPOIN,mesh2d%NPTFR,dt,grav,0    ,
     &            nfrliq,kent,kentu,msk,maskel,
     &            mesh2d%NELMAX,11,t2_11%R,numliq%I,mat2d%ADR(1)%P%X%R,
     &            t2_12%R,t2_13%R,t2_14%R,it3,it4,
     &            t2_15,t2_16,t2_17,t2_18,t2_19,t2_20,t3_01)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE THOMPS'
!
!     DUPLICATING ON THE VERTICAL
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE THOMPS_2DTO3D'
      CALL thomps_2dto3d
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE THOMPS_2DTO3D'
!
      ENDIF
!
!-----------------------------------------------------------------------
! SOURCE TERMS
!
      IF(nptsce.GT.0) THEN
        CALL findksce(npoin2,nplan,z3%R,nptsce,isce,zsce,ksce,infogr)
      ENDIF
      IF(nbuse.GT.0) THEN
        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
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TRISOU'
      CALL trisou
     & (s0u%R,s0v%R, s0u,s0v,un%R,vn%R,x,y,z,
     &  t3_01%R, deltar, mesh3d, fcor, coriol, ntrac,
     &  at, surfa2%R, t3_02%R, t3_02, w1%R,
     &  mesh3d%M%X%R(1:6*nelem3),mesh3d%M%X%R(6*nelem3+1:12*nelem3),
     &  grav, npoin3, nelem3, npoin2, nelem2, nplan, netage,
     &  ikle3%I, lv, msk, maskel%R, inchyd,
     &  volu,svide,ielm3,maskel,nrejeu,isce,ksce,qsce2,
     &  u_sce%R,v_sce%R,
     &  gradzs%ADR(1)%P,gradzs%ADR(2)%P,mesh2d,
     &  t3_03, t3_03%R, t3_04, t3_04%R, longit,
     &  yasem3d,schcvi,denlaw,fxh,fyh,courou,npth,t3d_files,t3dbi1)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TRISOU, APPEL DE SOURCE'
!
      CALL source(s0u, s0v, s0w, s1u, s1v, s1w,
     &            u, v, ws, w,
     &            volu, volun,t3_01,
     &            npoin3, ntrac, lt, at, dt, prive, nonhyd,
     &            npoin2, nsce,isce,ksce,qsce2,u_sce%R,v_sce%R,maxsce)
!
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOURCE'
!
!     SAVES BOUNDARY VALUES FOR TIME TN
!
      IF(nsousi.GT.1) THEN
        DO ip=1,nptfr3
          uborsave%R(ip)=un%R(nbor3%I(ip))
          vborsave%R(ip)=vn%R(nbor3%I(ip))
        ENDDO
        IF(nonhyd) THEN
          DO ip=1,nptfr3
            wborsave%R(ip)=wn%R(nbor3%I(ip))
          ENDDO
        ENDIF
        IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.
     &     iturbv.EQ.7.OR.iturbh.EQ.7) THEN
          DO ip=1,nptfr3
            kborsave%R(ip)=akn%R(nbor3%I(ip))
            eborsave%R(ip)=epn%R(nbor3%I(ip))
          ENDDO
        ELSEIF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
          DO ip=1, nptfr3
            nuborsave%R(ip)=nun%R(nbor3%I(ip))
          ENDDO
        ENDIF
        IF(ntrac.GT.0) THEN
          DO itrac=1,ntrac
            DO ip=1,nptfr3
              trborsave%ADR(itrac)%P%R(ip)=
     &        trn%ADR(itrac)%P%R(nbor3%I(ip))
            ENDDO
          ENDDO
        ENDIF
      ENDIF
!
!=======================================================================
! THE SUB-ITERATIONS LOOP BEGINS HERE
!=======================================================================
!
      subiter: DO isousi = 1,nsousi
!
#if defined COMPAD
      CALL ad_telemac3d_subiteration_begin
#endif
!
!     RESTORES BOUNDARY VALUES FOR TIME TN
!
      IF(isousi.GT.1) THEN
        DO ip=1,nptfr3
          un%R(nbor3%I(ip))=uborsave%R(ip)
          vn%R(nbor3%I(ip))=vborsave%R(ip)
        ENDDO
        IF(nonhyd) THEN
          DO ip=1,nptfr3
            wn%R(nbor3%I(ip))=wborsave%R(ip)
          ENDDO
        ENDIF
        IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.
     &     iturbv.EQ.7.OR.iturbh.EQ.7) THEN
          DO ip=1,nptfr3
            akn%R(nbor3%I(ip))=kborsave%R(ip)
            epn%R(nbor3%I(ip))=eborsave%R(ip)
          ENDDO
        ELSEIF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
          DO ip=1,nptfr3
            nun%R(nbor3%I(ip))=nuborsave%R(ip)
          ENDDO
        ENDIF
        IF(ntrac.GT.0) THEN
          DO itrac=1,ntrac
            DO ip=1,nptfr3
              trn%ADR(itrac)%P%R(nbor3%I(ip))=
     &        trborsave%ADR(itrac)%P%R(ip)
            ENDDO
          ENDDO
        ENDIF
      ENDIF
!
!     BUILDS THE MESH FOR PROPAGATION STEP
!
      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'
!
      IF(isousi.GT.1) THEN
!       REBUILDS THE INITIAL MESH
!       NOTE: EVOLUTION OF ZF IS NOT TAKEN INTO ACCOUNT HERE - INVESTIGATE
        CALL calcot(z,hn%R)
        CALL os('X=Y     ',x=volu,y=volun)
        IF(ncsize.GT.1) CALL os('X=Y     ',x=volupar,y=volunpar)
        CALL grad2d(gradzf%ADR(1)%P,gradzf%ADR(2)%P,zprop,nplan,svide,
     &              unsv2d,t2_02,t2_03,t2_04,
     &              ielm2h,mesh2d,msk,maskel)
        CALL fsgrad(gradzs,zflats,z(npoin3-npoin2+1:npoin3),
     &              zf,ielm2h,mesh2d,msk,maskel,
     &              unsv2d,t2_01,npoin2,optban,svide)
      ENDIF
!
!     SOURCES AND SINKS OF WATER
!
!     TEMPORARILY PUTS ZPROP IN MESH3D%Z
      savez     =>mesh3d%Z%R
      mesh3d%Z%R=>zprop%R
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOURCES_SINKS'
      CALL sources_sinks
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOURCES_SINKS'
!
      IF(inclus(coupling,'TOMAWACT3D')) THEN
!       ADD 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
!     RESTORES Z
      mesh3d%Z%R=>savez

      IF(bedbou)THEN
        IF(debug.GT.0) WRITE(lu,*) 'DEUXIEME APPEL DE BED_FLUXES'
        CALL bed_fluxes
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BED_FLUXES'
      ENDIF
!
!     SETS ADVECTION AND DIFFUSION PARAMETERS TO MONITOR CVDF3D
!     DIFFUSION AND SOURCE TERMS ARE DONE IN WAVE_EQUATION
!     IN CVDF3D (THIS IS DONE IN WAVE_EQUATION)
!
!     DIFFUSION OF U AND V IS DONE IN WAVE_EQUATION
      schdvi_hor = 0
      schdvi_ver = schdvi
!
      schcvi_hor = schcvi
      schcvi_ver = schcvi
!     ADVECTION IS NOT DONE AT THE FIRST TIME-STEP (THIS WAS VERSION 6.1)
!     IF(LT.EQ.1.AND.ISOUSI.EQ.1) THEN
!       SCHCVI_HOR = 0
!       SCHCVI_VER = 0
!     ENDIF
!
!     WHEN SCHCVI=ADV_SUP DIFF3D IS CALLED AND
!     SOURCE TERMS WOULD BE TREATED TWICE
      yas0u=.false.
      yas1u=.false.
!
!-----------------------------------------------------------------------
! ADVECTION-DIFFUSION STEP FOR VELOCITY COMPONENTS
!-----------------------------------------------------------------------
!
!     HERE DIFFUSION IS DONE IN MESH3D%Z, IT IS DIFFERENT FROM ZPROP IF
!     FROM THE SECOND SUB-ITERATION ON. SO IPBOT IS REDONE HERE
!
      IF(isousi.GT.1) THEN
        CALL plane_bottom(ipbot%I,z,npoin2,nplan,sigmag,optban)
      ENDIF
!
      IF(infogr) THEN
        IF (nonhyd) THEN
          CALL mittit(17,at,lt)
        ELSE
          CALL mittit(4,at,lt)
        ENDIF
      ENDIF
!
      sigmau = 1.d0
      umin   = 0.d0
      umax   = 1.d0
      clumin = .false.
      clumax = .false.
      yawchu = .false.
!     YASEM3D = DONE IN TRISOU
      newdif=.true.

      fluxbil(1:maxfro+nsce+1,1:5+ntrac) => fluxb%R
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR U'
      CALL cvdf3d
     & (ud,uc,un,viscvi,sigmau,s0u,yas0u,s1u,yas1u,
     &  uborl, uborf, ubors, auborl, auborf, aubors,
     &  buborl, buborf, bubors, liubol, liubof, liubos,
     &  fluxbil(1:maxfro+nsce+1,1),
     &  flux%R(1), fluext,fluextpar,umin, clumin, umax, clumax,
     &  schcvi_hor,schdvi_hor,slvdvi,trbavi,infogr,newdif,
     &  calcflu(1),t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04,mesh3d,ikle3,maskel,mtra1,
     &  nptfr3,mmurd,murd_tf,volu3d,volu3dpar,volu3d,volu3dpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,smu,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,tetadi,yawchu,wchu,s3d_wchu,
     &  agglod,nsce,sources,u_sce%R,numliq%I,dirflu,nfrliq,
     &  volut,zt,zprop,calcrain(1),pluie,parapluie,0.d0,flodel,flopar,
     &  sigmag,ipbot%I,maxadv,fludpt,fludp,fluer,volu2d,v2dpar,setdep,
     &  s3d_fludpt,s3d_fludp,s3d_fluer,s3d_setdep,optsou,zn%R,
     &  optadv_vi,nco_dist,nsp_dist,tb2)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR U'
!
      sigmav = 1.d0
      vmin   = 0.d0
      vmax   = 1.d0
      clvmin = .false.
      clvmax = .false.
      yawchu = .false.
!     YASEM3D = DONE IN TRISOU
!     MDIFF ALREADY COMPUTED FOR U
      newdif=.false.
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR V'
!     USE OF AUBORL,AUBORF,AUBORS IS NOT A MISTAKE
      CALL cvdf3d
     & (vd,vc,vn,viscvi,sigmav,s0v,yas0u,s1v,yas1u,
     &  vborl, vborf, vbors, auborl,auborf,aubors,
     &  bvborl, bvborf, bvbors, livbol, livbof, livbos,
     &  fluxbil(1:maxfro+nsce+1,2),
     &  flux%R(2), fluext,fluextpar,vmin, clvmin, vmax, clvmax,
     &  schcvi_hor,schdvi_hor,slvdvi,trbavi,infogr,newdif,
     &  calcflu(2),t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04, mesh3d , ikle3 , maskel , mtra1,
     &  nptfr3,mmurd,murd_tf,volu3d,volu3dpar,volu3d,volu3dpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,smv,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,tetadi,yawchu,wchu,s3d_wchu,
     &  agglod,nsce,sources,v_sce%R,numliq%I,dirflu,nfrliq,
     &  volut,zt,zprop,calcrain(2),pluie,parapluie,0.d0,flodel,flopar,
     &  sigmag,ipbot%I,maxadv,fludpt,fludp,fluer,volu2d,
     &  v2dpar,setdep,s3d_fludpt,s3d_fludp,s3d_fluer,s3d_setdep,optsou,
     &  zn%R,optadv_vi,nco_dist,nsp_dist,tb2)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR V'
!
      IF(nonhyd) THEN
!
        sigmaw = 1.d0
        wmin   = 0.d0
        wmax   = 1.d0
        clwmin = .false.
        clwmax = .false.
        yasem3d= .false.
        yawchu = .false.
        newdif=.true.
!       TETADI MAY BE EQUAL TO 2 FOR U AND V, WHEN THE WAVE EQUATION
!       IS USED - NOT DONE ON W SO FAR
        tetadiver = min(tetadi,1.d0)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR W'
!       USE OF AUBORL,AUBORF,AUBORS IS NOT A MISTAKE
        CALL cvdf3d
     & (wd,wc,wn,viscvi,sigmaw,s0w,.true.,s1w,.true.,
     &  wborl, wborf, wbors, auborl, auborf, aubors,
     &  bwborl, bwborf, bwbors, liwbol, liwbof, liwbos,
     &  fluxbil(1:maxfro+nsce+1,3),
     &  flux%R(3), fluext,fluextpar,wmin, clwmin, wmax, clwmax,
     &  schcvi_ver,schdvi_ver,slvdvi,trbavi,infogr,newdif,
     &  calcflu(3),t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04, mesh3d , ikle3 , maskel , mtra1,
     &  nptfr3,mmurd,murd_tf,volu3d,volu3dpar,volu3d,volu3dpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,sem3d,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,
     &  tetadiver,yawchu,wchu,s3d_wchu,agglod,nsce,sources,w_sce%R,
     &  numliq%I,dirflu,nfrliq,volut,zt,zprop,calcrain(3),pluie,
     &  parapluie,0.d0, flodel,flopar,
     &  sigmag,ipbot%I,maxadv,fludpt,fludp,fluer,volu2d,
     &  v2dpar,setdep,s3d_fludpt,s3d_fludp,s3d_fluer,s3d_setdep,optsou,
     &  zn%R,optadv_vi,nco_dist,nsp_dist,tb2)
!
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR W'
      ENDIF
!
!-----------------------------------------------------------------------
! DIFFUSION AND PROPAGATION STEP BY WAVE_EQUATION
!-----------------------------------------------------------------------
!
      IF(infogr) THEN
        CALL mittit(6,at,lt)
      ENDIF
!     TEMPORARILY PUTS ZPROP IN MESH3D%Z
      savez     =>mesh3d%Z%R
!     ALL PROPAGATION WILL BE DONE WITH ZPROP INSTEAD OF Z
      mesh3d%Z%R=>zprop%R
!     IPBOT HAS BEEN MODIFIED FOR CVDF3D, IT IS RESTORED HERE WITH ZPROP
      IF(isousi.GT.1) THEN
        CALL plane_bottom(ipbot%I,zprop%R,npoin2,nplan,sigmag,optban)
      ENDIF
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WAVE_EQUATION'
!
      CALL wave_equation(isousi)
!
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WAVE_EQUATION'
!
!     RESTORES Z
      mesh3d%Z%R=>savez
!
!-----------------------------------------------------------------------
! CLIPS NEGATIVE DEPTHS
!-----------------------------------------------------------------------
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CORRECTION_DEPTH_3D'
      CALL correction_depth_3d(mesh2d%GLOSEG%I,mesh2d%GLOSEG%DIM1)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CORRECTION_DEPTH_3D'
!
!-----------------------------------------------------------------------
! BUILDS NEW MESH WITH THE NEW FREE SURFACE
!-----------------------------------------------------------------------
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CALCOT'
      CALL calcot(z,h%R)
!     IPBOT UPDATED ACCORDINGLY, E.G. FOR CALLS TO PREDIV AND CVDF3D
      CALL plane_bottom(ipbot%I,z,npoin2,nplan,sigmag,optban)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CALCOT'
!
!----------------------------------------------------------------------
!
!     GENERATES DATA FOR DELWAQ
!
      IF(inclus(coupling,'DELWAQ')) THEN
!
!     COMPUTING FLODEL (POINT TO POINT FLUXES)
!
      formul = 'VGRADP       HOR'
      formul(8:8) = '2'
!     ADVECTION FLUXES PER NODE (STORED IN MESH3D%W%R)
!     THE ASSEMBLED RESULT IN T3_04 IS NOT USED HERE
      savez     =>mesh3d%Z%R
      mesh3d%Z%R=>zprop%R
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE VECTOR'
      CALL vector(t3_04,'=',formul,ielm3,-1.d0,dm1,svide,grazco,
     &            uconv,vconv,svide,mesh3d,msk,maskel)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE VECTOR'
      CALL flux_ef_vf_3d(flodel%R,mesh2d%W%R,mesh3d%W%R,
     &                   mesh2d%NSEG,nelem2,nelmax2,
     &                   mesh2d,.true.,
     &                   2,2,mesh3d%TYPELM+1,nplan,
!                          2: HORIZONTAL FLUXES FROM TOP TO BOTTOM
     &                   mesh3d%IKLE%I,mesh3d%NELMAX,mesh2d%KNOLG%I)
!     FLUX LIMITATION (FLULIM IS 2D, SO NUMBERING FROM TOP TO BOTTOM
!                      MAKES NO PROBLEM)
      IF(opt_hneg.EQ.2) THEN
        CALL flux3dlim(flodel%R,flulim%R,nplan,mesh2d%NSEG,npoin2,1)
      ENDIF
      mesh3d%Z%R=>savez
!
!     NOW CALLING TEL4DELWAQ WITH FLODEL COMPLETED
!
!     SENDS UCONV AND VCONV AS ADVECTING FIELD (SEE WAVE_EQUATION)
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TEL4DELWAQ'
      CALL tel4delwaq(npoin3,npoin2,mesh2d%NSEG,
     &  mesh2d%IKLE%I,mesh2d%ELTSEG%I,mesh2d%GLOSEG%I,
     &  mesh2d%GLOSEG%DIM1,x,y,mesh3d%NPTFR,lihbor%I,mesh3d%NBOR%I,
     &  nplan,at,dt,lt,nit_ori,h%R,hprop%R,mesh3d%Z%R,uconv%R,
     &  vconv%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
!
!----------------------------------------------------------------------
!
! MASKING
!
      IF(isousi.EQ.nsousi) THEN
        IF(msk) THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE MASK3D'
        IF(msk) 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,*) 'RETOUR DE MASK3D'
        ENDIF
      ENDIF
!
! COMPUTES SURFACE GRADIENTS AT TIME LEVEL N+1 AND DSSUDT
!
      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=Y-Z   ', x=dssudt, y=h, z=hn )
      CALL os( 'X=CX    ', x=dssudt, c=1.d0/dt )
!
! COMPUTES THE VOLUMES ASSOCIATED WITH 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(ncsize.GT.1) THEN
        CALL os('X=Y     ',x=volupar,y=volu)
        CALL parcom(volupar,2,mesh3d)
      ENDIF
!
! IN 2D, ONLY IF MASKING (OTHERWISE NOTHING CHANGED)
!
      IF(msk) CALL masbas2d(volu2d,v2dpar,unsv2d,
     &                      ielm2h,mesh2d,msk,maskel,t2_01,svide)
!
!-----------------------------------------------------------------------
! CONTINUITY STEP (NON-HYDROSTATIC OPTION) IN NEW MESH
!-----------------------------------------------------------------------
!
      IF(nonhyd.AND..NOT.dpwaveq) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'DANS NONHYDRO1'
        IF(infogr) CALL mittit(19,at,lt)
!
        CALL os ('X=Y     ', x=w , y=wd  )
!
!-----------------------------------------------------------------------
!
! COMPUTES THE DYNAMIC PRESSURE
!
!       WITH WAVE EQUATION, DYNAMIC PRESSURE HERE IS INCREMENTAL
!       THUS WITHOUT BOUNDARY CONDITIONS
        bc=.NOT.dpwaveq
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PREDIV'
        CALL prediv(dp,u,v,w,infogr,bc,1,.true.,.true.,.true.)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PREDIV'
!
!-----------------------------------------------------------------------
! VELOCITY PROJECTION STEP
!-----------------------------------------------------------------------
!
        IF(infogr) CALL mittit(20,at,lt)
!
        CALL velres(u,v,w,dp,
     &             t3_01,t3_02,t3_03,msk,maskel,mesh3d,
     &             svide,ielm3,optban,t3_04,.true.,npoin3,npoin2,
     &             sigmag,ipbot%I,aggloh,ksort,nptfr3,liubol,concor)
!
!       BOUNDARY CONDITIONS ON W AT THE BOTTOM AND FREE SURFACE
!
!       FREE SURFACE (NOT ALWAYS TO BE DONE, DSSUDT IS SOMETIMES TOO BIG)
!
        IF(cldyn) THEN
!
          CALL ov('X=Y     ',x=w%R(npoin3-npoin2+1:npoin3), y=dssudt%R,
     &                       dim1=npoin2)
          CALL ov('X=X+YZ  ',x=w%R(npoin3-npoin2+1:npoin3),
     &                       y=gradzs%ADR(1)%P%R,
     &                       z=u%R(npoin3-npoin2+1:npoin3), dim1=npoin2)
          CALL ov('X=X+YZ  ',x=w%R(npoin3-npoin2+1:npoin3),
     &                       y=gradzs%ADR(2)%P%R,
     &                       z=v%R(npoin3-npoin2+1:npoin3), dim1=npoin2)
!
        ENDIF
!
!       BOTTOM
!
        IF(velprobot) THEN
          IF(sigmag.OR.optban.EQ.1) THEN
            DO i=1,npoin2
              DO ip=0,ipbot%I(i)
                i3d=ip*npoin2+i
                w%R(i3d)=gradzf%ADR(1)%P%R(i)*u%R(i3d)
     &                  +gradzf%ADR(2)%P%R(i)*v%R(i3d)
              ENDDO
            ENDDO
          ELSE
            DO i=1,npoin2
              w%R(i)=gradzf%ADR(1)%P%R(i)*u%R(i)
     &              +gradzf%ADR(2)%P%R(i)*v%R(i)
            ENDDO
          ENDIF
        ENDIF
!
!       RE-ENSURES THE DIRICHLET BOUNDARY CONDITIONS AND U.N = 0
!
        CALL airwik2(lihbor%I, uborf%R, vborf%R, wborf%R,
     &               liubof%I, livbof%I, liwbof%I,
     &               uborl%R, vborl%R, wborl%R,
     &               liubol%I, livbol%I, liwbol%I,
     &               ubors%R, vbors%R, wbors%R,
     &               liubos%I, livbos%I, liwbos%I,
     &               u%R,v%R,w%R,mesh2d%XNEBOR%R,mesh2d%YNEBOR%R,
     &               nbor2%I,nptfr2,nplan,npoin2,kent,kadh,klog,kentu,
     &               velprolat)
!
      ENDIF ! IF NONHYD
!
!=======================================================================
! INTERNAL 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)
!       TFOND IS CALLED FOR AN UPDATE OF
!       UETCAR at TIME N+1
        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)
!
!       NOW RUNS ONE TURN OF GAIA'S TIME LOOP AND RETURNS CALL
        CALL config_code(5)
        IF(debug.GT.0) WRITE(lu,*) 'CALL GAIA_STEP'
!
        CALL gaia_step(
     &       lt,graprd,lisprd,nit_ori,u2d,v2d,h,zf,
     &       uetcar,deltar,cf,rugof,code1,
     &       u,v,at,viscvi,dt,charr_tel,susp_tel,
     &       rho0,grav,
!            VARIABLES TRANSMITTED FROM TOMAWAC
     &       dirmoy,hm0,tpr5,orbvel,.true.)
        IF(debug.GT.0) WRITE(lu,*) 'BACK GAIA_STEP'
        CALL config_code(1)
!
!       RETRIEVES ORIGINAL U AND V STRUCTURE
        CALL cpstvc(un,u)
        CALL cpstvc(vn,v)
!
      ENDIF
!
!=======================================================================
! INTERNAL COUPLING WITH SISYPHE
!=======================================================================
!
      IF( inclus(coupling,'SISYPHE')   .AND.
     &   (percou_sis*(lt/percou_sis).EQ.lt.OR.lt.EQ.1) ) THEN
!
!       U AND V WITH 2D STRUCTURE : BOTTOM VELOCITY AS A 2D VARIABLE
        CALL cpstvc(u2d,u)
        CALL cpstvc(v2d,v)
!
!       NOW RUNS ONE TURN OF SISYPHE'S TIME LOOP AND RETURNS
        CALL config_code(2)
        ibid=1
        lbid=.false.
        lbid2=.false.
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SISYPHE'
        CALL sisyphe(1,lt,graprd,lisprd,nit_ori,u2d,v2d,h,hn,zf,
     &               uetcar,cf,rugof,lbid,ibid,lbid2,code1,percou_sis,
     &               u,v,at,viscvi,dt*percou_sis,s3d_charr,s3d_susp,
!                          1 PRECLUDES THE USE OF THE 4 FOLLOWING ARGUMENTS
     &               flbor,1,dm1,uconv,vconv,zconv,
!                    VARIABLES TRANSMITTED FROM TOMAWAC
     &               dirmoy,hm0,tpr5,orbvel,.true.)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SISYPHE'
        CALL config_code(1)
!
!       S3D_HDEPMUST BE UPDATED BECAUSE SISYPHE CHANGED ZF
        IF(s3d_sedi) CALL os('X=Y-Z   ',x=s3d_hdep,y=zf,z=zr)
!
!       RETRIEVES ORIGINAL U AND V STRUCTURE
        CALL cpstvc(un,u)
        CALL cpstvc(vn,v)
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!     SEDIMENT
!
      IF(nsusp_tel.GT.0) THEN
!
!       COMPUTES THE SEDIMENT SETTLING VELOCITY
!
        CALL compute_settling_vel(wchu,u,v,
     &       t3_01,t3_02,t3_03,svide,mesh3d,ielm3,npoin2,npoin3,
     &       nplan,msk,maskel,uetcar,ta,hn)
!
!       BOUNDARY CONDITIONS FOR THE
!       SEDIMENTOLOGY
!
        CALL compute_bc_sedi
!
      ENDIF
!
!-----------------------------------------------------------------------
!     PREPARING SOURCE TERMS FOR ADVECTION-DIFFUSION STEP
!-----------------------------------------------------------------------
!
!     PREPARING SOURCE TERMS FOR K-EPSILON AND K-OMEGA MODELS
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.iturbv.EQ.7.OR.iturbh.EQ.7) THEN
!
        IF (infogr) CALL mittit(7,at,lt)
!
        s0ak%TYPR='Q'
        s0ep%TYPR='Q'
        s1ak%TYPR='Q'
        s1ep%TYPR='Q'
!
        IF(iturbv.EQ.3.OR.iturbh.EQ.3) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOUKEP'
        CALL soukep(s0ak%R,s0ep%R,s1ak%R,s1ep%R,
     &              u,v,w,deltar,ri%R,t3_01,t3_02,t3_03,t3_04,
     &              t3_05,t3_06,t3_07,t3_08,t3_09,
     &              t3_10,ak%R,ep%R,c1,c2,cmu,grav,
     &              npoin3,msk,maskel,mesh3d,ielm3,svide,
     &              vent,windx,windy,npoin2,kmin,prandtl)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOUKEP'
!
        ENDIF
!
        IF(iturbv.EQ.7.OR.iturbh.EQ.7) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOUKOM'
        CALL soukom(s0ak,s0ep,s1ak,s1ep,u,v,w,
     &              deltar,t3_01,t3_02,t3_03,
     &              t3_04,t3_05,t3_06,t3_07,t3_08,
     &              t3_09,t3_10,t3_12,t3_13,
     &              t3_14,t3_15,t3_16,t3_17,
     &              rotat,ak,ep,alpha,beta,betas,grav,
     &              t3_11,npoin3,msk,maskel,mesh3d,ielm3,svide)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOUKOM'
!
        ENDIF
!
      ELSEIF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
!
        IF (infogr) CALL mittit(33,at,lt)
        s0nu%TYPR='Q'
        s1nu%TYPR='Q'
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOUSA'

        CALL sousa(s0nu%R,s1nu%R,u,v,w,rotan,strain,t3_01,
     &              nu,npoin3,msk,maskel,mesh3d,
     &              ielm3,svide,wdist,npoin2, iturbv)
!
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOUSA'
      ENDIF
!
!     PREPARING SOURCE TERMS FOR TRACERS
!
      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
!         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
!
      ENDIF
!
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOURCE_TRAC'
      IF(ntrac.GT.0) CALL source_trac
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOURCE_TRAC'
!
!----------------------------------------------------------------------
! ADVECTION-DIFFUSION STEP FOR ALL ADVECTED VARIABLES
!----------------------------------------------------------------------
!
!     ALL ADVECTION SCHEMES EXCEPT SUPG
!
      IF (infogr .AND. (.NOT.nonhyd)) CALL mittit(9,at,lt)
      IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PREADV'
      CALL preadv(w,ws,zprop,isousi,lt,volu,volun)
      IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PREADV'
!
!----------------------------------------------------------------------
!     NOW CVDF3D WILL DO SUPG AND DIFFUSION
!----------------------------------------------------------------------
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3.OR.iturbv.EQ.7.OR.iturbh.EQ.7) THEN
!
        clkmin = .true.
        clkmax = .true.
        yasem3d = .false.
        yawchu = .false.
        newdif = .true.
        tetatra=min(tetadi,1.d0)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR AK'
        CALL cvdf3d
     & (ak,akc,akn,viscvi,sigmak,s0ak,.true.,s1ak,.true.,
     &  kborl, kborf, kbors, akborl, akborf, akbors,
     &  bkborl, bkborf, bkbors, likbol, likbof, likbos,
     &  fluxbil(1:maxfro+nsce+1,1),
     &  flux%R(1), fluext,fluextpar,kmin, clkmin, kmax, clkmax,
     &  schcke,schdke,slvdke,trbake,infogr,newdif,calcflu(4),
     &  t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04, mesh3d , ikle3 , maskel , mtra1,
     &  nptfr3,mmurd,murd_tf,volu,volupar,volun ,volunpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,sem3d,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,tetatra,
     &  yawchu,wchu,s3d_wchu,agglod,nsce,sources,ak_sce%R,
     &  numliq%I,dirflu,nfrliq,volut,zt,zprop,calcrain(4),
     &  pluie,parapluie,0.d0,flodel,flopar,sigmag,ipbot%I,maxadv,
     &  fludpt,fludp,fluer,volu2d,v2dpar,setdep,
     &  s3d_fludpt,s3d_fludp,s3d_fluer,s3d_setdep,optsou,
     &  zn%R,optadv_ke,nco_dist,nsp_dist,tb2)
!
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR AK'
!
        clemin  = .true.
        clemax  = .true.
        yasem3d = .false.
        yawchu  = .false.
!
!       NEGLECTS MOLECULAR DIFFUSIVITY...
!       DIFFUSION MATRIX NOT RECOMPUTED
        newdif = .false.
        CALL om('M=CM    ',m=mdiff,c=sigmae/sigmak,mesh=mesh3d)
!
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR EP'
        CALL cvdf3d
     & (ep,epc,epn,viscvi,sigmae,s0ep,.true.,s1ep,.true.,
     &  eborl, eborf, ebors, aeborl, aeborf, aebors,
     &  beborl, beborf, bebors, liebol, liebof, liebos,
     &  fluxbil(1:maxfro+nsce+1,1),
     &  flux%R(1), fluext,fluextpar,emin, clemin, emax, clemax,
     &  schcke,schdke,slvdke,trbake,infogr,newdif,calcflu(5),
     &  t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04, mesh3d , ikle3 , maskel , mtra1,
     &  nptfr3,mmurd,murd_tf,volu,volupar,volun,volunpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,sem3d,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,tetatra,
     &  yawchu,wchu,s3d_wchu,agglod,nsce,sources,ep_sce%R,
     &  numliq%I,dirflu,nfrliq,volut,zt,zprop,calcrain(5),
     &  pluie,parapluie,0.d0,flodel,flopar,sigmag,ipbot%I,
     &  maxadv,fludpt,fludp,fluer,volu2d,v2dpar,
     &  setdep,s3d_fludpt,s3d_fludp,s3d_fluer,s3d_setdep,
     &  optsou,zn%R,optadv_ke,nco_dist,nsp_dist,tb2)
!
!
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR EP'
!
      ELSEIF(iturbv.EQ.5.OR.iturbv.EQ.9) THEN
!
        clnumin = .true.
        clnumax = .false.
        yasem3d = .false.
        yawchu = .false.
        newdif = .true.
        tetatra=min(tetadi,1.d0)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CVDF3D POUR NU'
        CALL cvdf3d
     & (nu,nuc,nun,viscnu,sigmanu,s0nu,.true.,s1nu,.true.,
     &  nuborl, nuborf, nubors, anuborl, anuborf, anubors,
     &  bnuborl, bnuborf, bnubors, linubol, linubof, linubos,
     &  fluxbil(1:maxfro+nsce+1,1),
     &  flux%R(1),fluext,fluextpar,numin,clnumin,numax,clnumax,
     &  schcke,schdke,slvdke,trbake,infogr,newdif,calcflu(4),
     &  t2_01,t2_03,
     &  t3_01,t3_02,t3_03,t3_04, mesh3d , ikle3 , maskel , mtra1,
     &  nptfr3,mmurd,murd_tf,volu,volupar,volun ,volunpar,
     &  nbor3,npoin3,npoin2,dt,msk,nelem3,
     &  nplan,ielm3,msupg,ielm2h,ielm2v,mdiff,mtra2,
     &  inchyd,maskbr,maskpt,sem3d,yasem3d,svide,it1,
     &  trav3,mesh2d,optban,tetatra,
     &  yawchu,wchu,s3d_wchu,agglod,nsce,sources,nu_sce%R,
     &  numliq%I,dirflu,nfrliq,volut,zt,zprop,calcrain(4),
     &  pluie,parapluie,0.d0,
     &  flodel,flopar,sigmag,ipbot%I,maxadv,fludpt,fludp,
     &  fluer,volu2d,v2dpar,setdep,s3d_fludpt,s3d_fludp,s3d_fluer,
     &  s3d_setdep,optsou,zn%R,optadv_ke,nco_dist,
     &  nsp_dist,tb2)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CVDF3D POUR NU'
      ENDIF
!
!-----------------------------------------------------------------------
!
! COMPUTES THE VISCOSITIES VISCVI, VISCTA AND VISCTP
!
      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'
!
      ENDIF
!
      IF(iturbv.EQ.3.OR.iturbh.EQ.3) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE VISCKE'
        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'
!
      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.iturbv.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_STEP'
          CALL gotm_coupling_step
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE GOTM_COUPLING_STEP'
      ENDIF
!
      IF(inclus(coupling,'TOMAWACT3D')) THEN
!       ADDING VERTICAL MIXING DUE TO WAVES
        CALL os('X=X+Y   ',x=viscvi%ADR(3)%P,y=fdk)
      ENDIF
!
      IF(optban.EQ.1) THEN
!
        CALL visclip(viscvi,viscta,h,nplan,npoin2,ntrac,hlim)
!
      ENDIF
!
!=======================================================================
!     OIL SPILL MODEL (UNDER DEVELOPMENT IN MYGRHYCAR PROJECT)
!=======================================================================
!
      IF(spill_model) THEN
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE OIL_SPILL_3D'
          CALL oil_spill_3d(lt,ielm2h,mesh2d,nflot_max,t3d_files,
     &                     maxlu_t3d,npoin2,t3dmig,uconv,vconv,wconv,
     &                     nflot,nplan,mesh3d,at,dt,grav,cf,x,y,z,h,hn,
     &                     ielm3,npoin3,nelem2,xflot,yflot,zflot,shpflo,
     &                     shzflo,tagflo,clsflo,eltflo,etaflo,floprd,
     &                     t3dflo,it1,it2,t3_01,t3_02,t3_03,mtra1,mtra2,
     &                     viscvi,windx,windy,unsv3d,ntrac,trn,trav3,
     &                     atabos,t2_17,vent)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE OIL_SPILL_3D'
!
      ENDIF
!
!-----------------------------------------------------------------------
! ADVECTION-DIFFUSION OF TRACERS
!
!GA: POINTERS ARE USED TO ALLOW TWO SEDIMENT CLASSES (MULTI1,2,3,4)
      IF(ntrac.GT.0) THEN
!
        IF (infogr) CALL mittit(5,at,lt)
!
          sigmta = 1.d0
          tamin  = 0.d0
          tamax  = 1.d0
          ctamin = .false.
          ctamax = .false.
          yasem3d = .false.
          newdif = .true.
          tetatra=min(tetadi,1.d0)
!
!     IF MICROPOL MODULE IS ACTIVATED, TRACERS BOTTOM SEDIMENT (IND_SF)
!     AND MICROPOLLUTANT ADSORBED BY BOTTOM SEDIMENT (IND_CSF)
!     ARE NEITHER ADVECTED NOR DIFFUSED (NO ADVECTION DONE IN LECDON)
        IF(inclus(coupling,'WAQTEL').AND.
     &     7*int(waqprocess/7).EQ.waqprocess) THEN
          CALL os('X=0     ', x=viscta%ADR(ind_sf)%P)
          CALL os('X=0     ', x=viscta%ADR(ind_csf)%P)
        ENDIF
!
        DO itrac = 1,ntrac
!
          multi1 => s3d_fludpt
          multi2 => s3d_fludp
          multi3 => s3d_fluer
          multi4 => s3d_wchu
!
          IF(itrac.GE.ind_sed.AND.itrac.LE.ind_sed+nsusp_tel-1) THEN
            isusp=itrac-ind_sed+1
            yawchu=.true.
!           SOLVER STRUCTURE
            slvd=slvdta(itrac)
          ELSEIF(s3d_sedi.AND.itrac.EQ.ntrac) THEN
            yawchu=.true.
!           SOLVER STRUCTURE
            slvd=s3d_slvdse
          ELSE
            yawchu=.false.
!           SOLVER STRUCTURE
            slvd=slvdta(itrac)
          ENDIF
!
          IF(s3d_sedi.AND.s3d_mixte.AND.itrac.EQ.(ntrac-1)) THEN
            yawchu=.true.
            slvd=s3d_slvdse
            multi1 => s3d_fludptnc
            multi2 => s3d_fludpnc
            multi3 => s3d_fluernc
            multi4 => s3d_wcs
          ELSEIF(s3d_sedi.AND.s3d_mixte.AND.itrac.EQ.ntrac) THEN
            multi1 => s3d_fludptc
            multi2 => s3d_fludpc
            multi3 => s3d_fluerc
          ENDIF
!
          IF(debug.GT.0) THEN
            WRITE(lu,*) 'APPEL DE CVDF3D POUR TRACEUR ',itrac
          ENDIF
!
          CALL cvdf3d
     &   (ta%ADR(itrac)%P,tac%ADR(itrac)%P,trn%ADR(itrac)%P,
     &    viscta%ADR(itrac)%P,sigmta,
     &    s0ta%ADR(itrac)%P,.true.,s1ta%ADR(itrac)%P,.true.,
     &    taborl%ADR(itrac)%P,taborf%ADR(itrac)%P,tabors%ADR(itrac)%P,
     &    atabol%ADR(itrac)%P,atabof%ADR(itrac)%P,atabos%ADR(itrac)%P,
     &    btabol%ADR(itrac)%P,btabof%ADR(itrac)%P,btabos%ADR(itrac)%P,
     &    litabl%ADR(itrac)%P,litabf%ADR(itrac)%P,litabs%ADR(itrac)%P,
     &    fluxbil(1:maxfro+nsce+1,5+itrac),
     &    flux%R(5+itrac),fluext,fluextpar,
     &    tamin,ctamin,tamax,ctamax,schcta(itrac),
     &    schdta,slvd,trbata,infogr,newdif,calcflu(5+itrac),
     &    t2_01,t2_03,t3_01,t3_02,t3_03,t3_04,mesh3d,ikle3,maskel,
     &    mtra1,nptfr3,mmurd,murd_tf,volu,volupar,volun,volunpar,
     &    nbor3,npoin3,npoin2,dt,msk,nelem3,nplan,ielm3,msupg,
     &    ielm2h,ielm2v,mdiff,mtra2,inchyd,maskbr,maskpt,sem3d,yasem3d,
     &    svide,it1,trav3,mesh2d,optban,tetatra,
     &    yawchu,wchu%ADR(itrac)%P,multi4,agglod,nsce,sources,
     &    ta_sce%ADR(itrac)%P%R,
     &    numliq%I,dirflu,nfrliq,volut,zt,zprop,calcrain(5+itrac),
     &    pluie,parapluie,train(itrac),flodel,flopar,sigmag,ipbot%I,
     &    maxadv,fludpt%ADR(num_isusp_icla(isusp))%P,
     &    fludp%ADR(num_isusp_icla(isusp))%P,
     &    fluer%ADR(num_isusp_icla(isusp))%P,volu2d,v2dpar,setdep,
     &    multi1,multi2,multi3,s3d_setdep,optsou,
     &    zn%R,optadv_tr(itrac),nco_dist,nsp_dist,tb2)
!
!         NEWDIF=.FALSE. (POSSIBLE IF SIGMTA UNCHANGED)
!
          IF(debug.GT.0) THEN
            WRITE(lu,*) 'RETOUR DE CVDF3D POUR TRACEUR ',itrac
          ENDIF
!
        ENDDO
!
!-----------------------------------------------------------------------
! COMPUTES DELRA RHO / RHO FOR THE BUOYANCY TERMS
!
        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(debug.GT.0) WRITE(lu,*) 'RETOUR DE DRSURR'
!
        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
!
      ENDIF
!
#if defined COMPAD
      CALL ad_telemac3d_subiteration_end
#endif
!
      END DO subiter
!
!-----------------------------------------------------------------------
!
! SEDIMENTOLOGY : FLUDP SENT TO GAIA
!
      IF(nsusp_tel.GT.0) THEN
        DO itrac=ind_sed, ind_sed+nsusp_tel-1
          isusp=itrac-ind_sed+1
!
          IF(setdep.NE.1) THEN
            IF(optban.EQ.1) THEN
              DO ipoin=1,npoin2
!               correction for tidal flats: take the first point above crushed planes
!               IPBOT=0: no tidal flats
!               IPBOT=NPLAN-1: dry element
                IF(ipbot%I(ipoin).NE.nplan-1) THEN
                  fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin)=
     &                  fludpt%ADR(num_isusp_icla(isusp))%P%R(ipoin)*
     &                  ta%ADR(itrac)%P%R(ipbot%I(ipoin)*npoin2+ipoin)
                  fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin)=
     &                  max(fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin),
     &                      0.d0)
                ELSE
                  fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin)=0.d0
                ENDIF
              ENDDO
            ELSE
              DO ipoin=1,npoin2
                fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin)=
     &                fludpt%ADR(num_isusp_icla(isusp))%P%R(ipoin)*
     &                ta%ADR(itrac)%P%R(ipoin)
!               FLUDP MUST BE POSITIVE, EVEN IF TA<0 DUE TO TRUNCATION ERRORS PROBLEM
!               SEEN WITH TA=-1.D-87 !!!!!
                fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin)=
     &                max(fludp%ADR(num_isusp_icla(isusp))%P%R(ipoin),
     &                0.d0)
              ENDDO
            ENDIF
          ENDIF
!
        ENDDO
      ENDIF
!
! SEDIMENTOLOGY SEDI3D: BOTTOM TREATMENT
!
      IF(s3d_sedi) THEN
!
!       FONVAS DOES ZF=ZR+S3D_HDEP, THUS S3D_HDEPMUST INCLUDE BEDLOAD
!       EROSION, HAS BEEN TAKEN INTO ACCOUNT INTO CLSEDI ABOVE
!   GA: FONVAS HAS BEEN MODIFIED TO ALLOW TWO SEDIMENT CLASSES
        IF (s3d_mixte) THEN
          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FONVAS S3D_MIXTE'
          CALL fonvas
     &       (s3d_epai,s3d_conc,
     &        s3d_hdep%R,s3d_fludp%R,s3d_fludpt%R,s3d_fluer%R,zf%R,
     &        ta%ADR(ntrac)%P%R,
     &        npoin2,npoin3,s3d_ncouch,dt,
     &        s3d_zf_s%R,s3d_esomt%R,volu2d,s3d_masdep,s3d_setdep,zr%R,
     &        ta%ADR(ntrac-1)%P%R,s3d_fludptc%R,s3d_fludptnc%R,
     &        s3d_fluerc%R,s3d_fluernc%R,s3d_mixte,s3d_fludpc%R,
     &        s3d_fludpnc%R,s3d_pvsco%R,s3d_pvsnco%R,s3d_cfdep,
     &        s3d_epaico%R,s3d_epainco%R)

          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FONVAS S3D_MIXTE'
!
        ELSEIF (s3d_sedco.OR.s3d_sednco) THEN

          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FONVAS '
          CALL fonvas
     &       (s3d_epai,s3d_conc,
     &        s3d_hdep%R,s3d_fludp%R,s3d_fludpt%R,s3d_fluer%R,zf%R,
     &        ta%ADR(ntrac)%P%R,
     &        npoin2,npoin3,s3d_ncouch,dt,
     &        s3d_zf_s%R,s3d_esomt%R,volu2d,s3d_masdep,s3d_setdep,zr%R,
     &        ta%ADR(ntrac)%P%R,s3d_fludptc%R,s3d_fludptnc%R,
     &        s3d_fluerc%R,s3d_fluernc%R,s3d_mixte,s3d_fludpc%R,
     &        s3d_fludpnc%R,s3d_pvsco%R,s3d_pvsnco%R,s3d_cfdep,
     &        s3d_epaico%R,s3d_epainco%R)

          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FONVAS '
        ELSE
!
          WRITE(lu,*) ' '
          WRITE(lu,*) 'SEDI3D : ERROR ON SEDIMENT KEY WORD'
          CALL plante(1)
          stop
        ENDIF
      ENDIF
!
!
      IF(inclus(coupling,'GAIA  ').AND.nsusp_tel.GT.0) THEN
!
!       U AND V WITH 2D STRUCTURE : BOTTOM VELOCITY AS A 2D
!       VARIABLE
        CALL cpstvc(u2d,u)
        CALL cpstvc(v2d,v)
!
!       NOW RUNS ONE TURN OF GAIA'S TIME LOOP
!       AND RETURNS
        CALL config_code(5)
        IF(debug.GT.0) WRITE(lu,*) 'CALLING GAIA_STEP 2'
        CALL gaia_step(
     &       lt,graprd,lisprd,nit_ori,u2d,v2d,h,zf,
     &       uetcar,deltar,cf,rugof,code1,
     &       u,v,at,viscvi,dt,charr_tel,susp_tel,
     &       rho0,grav,
!            VARIABLES TRANSMITTED FROM TOMAWAC
     &       dirmoy,hm0,tpr5,orbvel,.true.)
        IF(debug.GT.0) WRITE(lu,*) 'BACK GAIA_STEP 2'
        CALL config_code(1)
!
!       check if this is necessary
!       RETRIEVES ORIGINAL U AND V STRUCTURE
        CALL cpstvc(un,u)
        CALL cpstvc(vn,v)
!       end check
      ENDIF
!
! UPDATES GEOMETRY IF THE BOTTOM HAS EVOLVED
!
      IF(inclus(coupling,'SISYPHE').OR.inclus(coupling,'GAIA  ')
     &   .OR.s3d_sedi) THEN
!
!       COPIES MODIFIED BOTTOM TOPOGRAPHY INTO Z AND ZPROP
        CALL ov('X=Y     ', x=z(1:npoin2), y=zf%R, dim1=npoin2)
        CALL ov('X=Y     ', x=zprop%R(1:npoin2), y=zf%R, dim1=npoin2)
!       COMPUTES NEW BOTTOM GRADIENTS AFTER SEDIMENTATION
        CALL grad2d(gradzf%ADR(1)%P,gradzf%ADR(2)%P,zprop,nplan,svide,
     &              unsv2d,t2_02,t2_03,t2_04,
     &              ielm2h,mesh2d,msk,maskel)
!       COMPUTES NEW Z COORDINATES
        CALL calcot(z,h%R)
!       USEFUL ? NOT SURE, IS DONE AT EACH TIMESTEP ELSEWHERE, SO..
!       CALL CALCOT(ZPROP%R,HPROP%R)
        CALL fsgrad(gradzs,zflats,z(npoin3-npoin2+1:npoin3),
     &              zf,ielm2h,mesh2d,msk,maskel,
     &              unsv2d,t2_01,npoin2,optban,svide)
        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(ncsize.GT.1) THEN
          CALL os('X=Y     ',x=volupar,y=volu)
          CALL parcom(volupar,2,mesh3d)
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
! PREPARES 2D AND 3D OUTPUT
!
!     ALWAYS CALLED (SEE E.G. OUTPUT MAXZ)
!
      CALL preres_telemac3d
!
!     CALLED IF OUTPUTS REQUESTED
!
      IF(mod(lt,graprd).EQ.0.AND.lt.GE.gradeb) THEN
!
!       3D OUTPUT
!
        IF(t3d_files(t3dres)%NAME(1:1).NE.' ') THEN
          CALL bief_desimp(t3d_files(t3dres)%FMT,varso3,
     &                     npoin3,t3d_files(t3dres)%LU,at,lt,
     &                     lisprd,graprd,
     &                     sorg3d,sorim3,maxva3,text3,gradeb,lisdeb)
        ENDIF
!
!       3D OUTPUT FOR RESTART
!
        IF(lt.EQ.nit_ori.AND.restart_mode
     &              .AND.t3d_files(t3drst)%NAME(1:1).NE.' ') THEN
          CALL bief_desimp(t3d_files(t3drst)%FMT,varso3,npoin3,
     &                     t3d_files(t3drst)%LU,at,lt,
     &                     1,nit_ori,sorest,soris3,maxva3,text3,0,0)
        ENDIF
!
!       2D OUTPUT
!
        IF(t3d_files(t3dhyd)%NAME(1:1).NE.' ') THEN
          CALL bief_desimp(t3d_files(t3dhyd)%FMT,varsor,
     &                     npoin2,t3d_files(t3dhyd)%LU,at,lt,
     &                     lisprd,graprd,
     &                     sorg2d,sorimp,maxvar,texte,gradeb,lisdeb)
        ENDIF
!
      ENDIF
!
! SEDIMENTOLOGY OUTPUT FOR SEDI
!
      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
!
! OPTIONAL USER OUTPUT
!
      CALL utimp(gradeb,graprd,lisdeb,lisprd)
!
! SEDIMENT OUTPUT FOR SEDI
!
      IF(s3d_sedi) THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE IMPSED'
        CALL impsed(s3d_ivide%R,s3d_epai,s3d_conc,s3d_temp%R,s3d_hdep%R,
     &              s3d_pdepo%R,s3d_fluer%R,
     &              zr%R,zf%R,ta%ADR(ntrac)%P%R,s3d_wchu%R,x,y,
     &              npoin2,npoin3,s3d_npfmax,s3d_ncouch,npf%I,lt,
     &              s3d_rhos,s3d_cfmax,
     &              s3d_cfdep,s3d_epai0,s3d_tasse,s3d_gibson,
     &              prive,lisprd)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE IMPSED'
      ENDIF
!
!=======================================================================
! DROGUES/FLOATS/BUOYS
!=======================================================================
!
      IF(nflot_max.GT.0.AND..NOT.spill_model) THEN
!
        IF(spheri) THEN
          CALL os('X=Y/Z   ',x=uconv,y=uconv,z=mesh3d%COSLAT)
          CALL os('X=Y/Z   ',x=vconv,y=vconv,z=mesh3d%COSLAT)
        ENDIF
!
!       ADDING AND REMOVING DROGUES
!
        IF(debug.GT.0) WRITE(lu,*) 'CALLING FLO3D'
        CALL flot3d(xflot%R,yflot%R,zflot%R,nflot,nflot_max,x,y,z,
     &              mesh3d%IKLE%I,
     &              mesh3d%NELEM,mesh3d%NELMAX,npoin3,nplan,
     &          tagflo%I,clsflo%I,shpflo%R,shzflo%R,eltflo%I,etaflo%I,
     &              mesh3d,lt,nit_ori,at)
        IF(debug.GT.0) WRITE(lu,*) 'BACK FROM FLO3D'
!
        IF(infogr) CALL mittit(12,at,lt)
!
!       MOVING THEM
!
        IF(debug.GT.0) WRITE(lu,*) 'CALLING DERIVE'
        CALL derive(uconv%R,vconv%R,wconv%R,dt,at,
     &              x,y,z,
     &              mesh2d%IKLE%I,mesh3d%IFABOR%I,lt,ielm3,uconv%ELM,
     &              3,3,
     &              npoin3,npoin2,nelem2,mesh2d%NELMAX,
     &              mesh2d%SURDET%R,xflot%R,yflot%R,zflot%R,
     &          shpflo%R,shzflo%R,tagflo%I,clsflo%I,eltflo%I,etaflo%I,
     &          nflot,nflot_max,mesh3d,
     &              it1%I,t3_01%R,t3_02%R,t3_03%R,it2%I,
!                                              NO STOCHASTIC DIFFUSION
     &              mtra1%X%R,mtra2%X%R,npoin3,0,svide,
     &              nplan,zchar%R,transf)
        IF(debug.GT.0) WRITE(lu,*) 'BACK FROM DERIVE'
!
!       PRINTING THEM
!
        IF(debug.GT.0) WRITE(lu,*) 'CALLING UTIMP_DROGUES'
        CALL utimp_drogues( lt,at,npoin2,npoin3,
     &    xflot%R,yflot%R,zflot%R,tagflo%I,clsflo%I,
     &    nflot,nflot_max,floprd,deja_drogues_utimp,
     &    t3d_files(t3dflo),t3d_files(t3dblo), mardat,martim )
        IF(debug.GT.0) WRITE(lu,*) 'BACK FROM UTIMP_DROGUES'
!
        IF(spheri) THEN
          CALL os('X=XY    ',x=uconv,y=mesh3d%COSLAT)
          CALL os('X=XY    ',x=vconv,y=mesh3d%COSLAT)
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
! MASS BALANCE FOR THE CURRENT TIMESTEP
!
      IF (bilmas) THEN
!
        IF (.NOT.infmas) infogr = .false.
        infogr = infogr .AND. listin
        IF (infogr) CALL mittit(10,at,lt)
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE MASS3D'
        CALL mass3d(infogr)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE MASS3D'
!
!
        IF(s3d_sedi) THEN
!
!         DETERMINE S3D_MASSUSP: MASS IN SUSPENSION
!         S3D_MASBED: MASS OF SEDIMENT BED
!         S3D_MASDEP: DEPOSITED MASS
          s3d_massusp= masse%R(5+ntrac)
!         INITIALISATION
          IF(lt.EQ.1) s3d_massusp0= s3d_massusp
!
!GA: SED3D HAS BEEN MODIFIED TO ALLOW TWO SEDIMENT CLASSES
          IF(infogr) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SED3D'
            CALL sed3d(s3d_masbed,s3d_masbed0, s3d_masdep,
     &                 s3d_epai,s3d_conc,
     &                 t2_01%R,
     &                 npoin2,s3d_ncouch,
     &                 at,volu2d%R,
     &                 s3d_cfdep,s3d_epaico%R,s3d_epainco%R,s3d_mixte)
!
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SED3D'
          ENDIF
!
        ENDIF
!
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BIL3D'
        CALL bil3d(lt,nit_ori,mesh3d%IKLBOR%I,nptfr2,netage)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIL3D'
!
      ENDIF
!
! COMPARES WITH REFERENCE FILE
!
      IF(valid.AND.t3d_files(t3dres)%NAME(1:1).NE.' ') THEN
        IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BIEF_VALIDA'
        CALL bief_valida(trav3,textp3,
     &                   t3d_files(t3dref)%LU,t3d_files(t3dref)%FMT,
     &                   varso3,text3,
     &                   t3d_files(t3dres)%LU,t3d_files(t3dres)%FMT,
     &                   maxva3,npoin3,lt,nit_ori,alire3d)
        IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIEF_VALIDA'
      ELSEIF(valid.AND.t3d_files(t3dres)%NAME(1:1).EQ.' '.AND.
     &       lt.EQ.nit_ori) THEN
        WRITE(lu,*) ' '
        WRITE(lu,*) ' '
        WRITE(lu,*) 'NO VALIDATION IF NO'
        WRITE(lu,*) '3D RESULT FILE!!!'
      ENDIF
!
!
! CHECKS VALUES SHARED BETWEEN SUBDOMAINS
!
!     CALL CHECK_DIGITS(H ,T2_01,MESH2D)
!     CALL CHECK_DIGITS(U ,T3_01,MESH3D)
!     CALL CHECK_DIGITS(V ,T3_01,MESH3D)
!     CALL CHECK_DIGITS(W ,T3_01,MESH3D)
!     IF(NTRAC.GT.0) THEN
!       DO ITRAC=1,NTRAC
!         CALL CHECK_DIGITS(TA%ADR(ITRAC)%P,T3_01,MESH3D)
!       ENDDO
!     ENDIF
!
#if defined COMPAD
      CALL ad_telemac3d_timestep_end
#endif
!
      IF(lt.LT.nit) GO TO 700
! END OF TIME LOOP
!
!=======================================================================
! THE TIME LOOP ENDS HERE
!=======================================================================
!
      IF(pass.NE.1) THEN
        IF(listin) WRITE(lu,19)
19      FORMAT(/,1x,'END OF TIME LOOP',////)
      ENDIF
!
!-----------------------------------------------------------------------
!
#if defined COMPAD
      CALL ad_telemac3d_end
#endif
!
      RETURN
      END