wac.F

Go to the documentation of this file.

!                   **************
                    SUBROUTINE wac
!                   **************
!
     &     (part, nit_ori)
!
!**********************************************************************
!     TOMAWAC   V7P3
!**********************************************************************
!
!     brief    MAIN SUBROUTINE OF TOMAWAC
!     +               SOLVES THE EQUATION FOR THE
!     +               DIRECTIONAL WAVE SPECTRUM
!
!     history  J-M HERVOUET (EDF - LNHE)
!     +        29/01/2013
!     +        V6P3
!     +Radiation stresses for Telemac now computed independently of the
!     +   printouts on results file.
!     +   Call to tomawac_constants moved to lecdon_tomawac.
!
!     history  J-M HERVOUET (EDF - LNHE)
!     +        22/03/2013
!     +        V6P3
!     +   New arguments DIRMOY_TEL,HM0_TEL, TPR5_TEL for transmission to
!     +   Sisyphe through Telemac-2D or 3D. Values computed in case of triple
!     +   coupling.
!
!     history  E. GAGNAIRE-RENOU (EDF - LNHE)
!     +        12/09/2014
!     +        V7P0
!     +   New arguments for SEMIMP for wave-blocking effects
!     +   Call the new routine LIMITE
!
!     history  J-M HERVOUET (EDF - LNHE)
!     +        15/09/2014
!     +        V7P0
!     +   Initialisation of SUSOLD and STRA34 if .NOT.VENT
!     +   (a former overlooked bug)
!     +   PERCOU_WAC now used (correction by C. Villaret (HR+EDF)
!     +   and Pablo*(Tassi+Santoro)).
!
!     history  C VILLARET (HRW+EDF) & J-M HERVOUET (EDF - LNHE)
!     +        18/09/2014
!     +        V7P0
!     +   Adding the variable ORBVEL_TEL in argument (orbital velocity)
!     +   for sending it back to the calling program.
!
!     history Y AUDOUIN (LNHE)
!     +       25/05/2015
!     +       V7P0
!     +       Modification to comply with the hermes module
!
!     history A JOLY (LNHE)
!     +       18/05/2017
!     +       V7P3
!     + Possibility to read a text file with coordinates of the nodes
!     +       for which a spectra will be written.
!
!     history A JOLY (LNHE)
!     +       18/05/2017
!     +       V7P3
!     +    Moved LIMWAC to the end of the time stepping operations to
!     +       better impose dirac boundary conditions.
!
!     history A JOLY (EDF-LNHE)
!     +       18/05/2017
!     +       V7P3
!     +  In semimp, new condition to stop source terms being added to
!     +       open boundaries.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!     | PART     |-->| -1: NO COUPLING
!     |          |   |  0: COUPLING WITH TELEMAC (INITIALISATION)
!     |          |   |  1: COUPLING WITH TELEMAC (LOOP OVER TIME STEPS)
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE couple_mod
      USE declarations_telemac
      USE declarations_tomawac
      USE interface_tomawac, ex_wac => wac
!
      USE declarations_special
      IMPLICIT NONE

!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!
      INTEGER,           INTENT(IN)      :: PART
      INTEGER,           INTENT(IN),OPTIONAL :: NIT_ORI
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!
      DOUBLE PRECISION DT_MIN,DT_MAX
      DOUBLE PRECISION, PARAMETER :: DEG2RAD = acos(0.0d0)/90.0d0
      INTEGER          DUMMY,LT_WAC, IREC, IERR, I
!
      INTEGER LT1
      INTEGER DATE(3),TIME(3),IP
!
      INTEGER ADC,MDC,JDC,HDC
      LOGICAL IMPRES, DEBRES
      INTEGER TOTAL_ITER
!-----------------------------------------------------------------------
!
      IF(part.EQ.wac_full_run.OR.
     &     part.EQ.wac_cpl_init.OR.
     &     part.EQ.wac_cpl_run.OR.
     &     part.EQ.wac_api_init) THEN
#if defined COMPAD
        CALL ad_tomawac_begin
        CALL ad_tomawac_initialisation_begin
#endif
!=====C
!     1  C INITIALISATIONS DES VARIABLES LOCALES
!=====C======================================
!     COUPLAGE : verification des conditions pour le couplage
!     TELEMAC-TOMAWAC
!
        IF(part.EQ.wac_cpl_init.OR.part.EQ.wac_cpl_run) THEN
          IF(debug.GT.0) WRITE(lu,*) 'ENTERING TOMAWAC, WITH T_TEL=',
     &         cpl_wac_data%AT_TEL,
     &         'NIT_TEL',cpl_wac_data%NIT_TEL
          IF(maree.OR.cousta.OR.dontel) THEN
            WRITE(lu,*) ''
            WRITE(lu,*) '***************************************'
            WRITE(lu,*) ' ATTENTION : COUPLING TELEMAC-TOMAWAC :'
            WRITE(lu,*) ' CURRENT/WATER LEVEL FILE CANNOT BE    '
            WRITE(lu,*) ' USED AS INPUT FILE.                   '
            WRITE(lu,*) '         END OF THE COMPUTATION        '
            WRITE(lu,*) '***************************************'
            CALL plante(1)
            stop
          ENDIF
!
          dt_max=max(dt,cpl_wac_data%DT_TEL*cpl_wac_data%PERCOU_WAC)
          dt_min=min(dt,cpl_wac_data%DT_TEL*cpl_wac_data%PERCOU_WAC)
!
          IF(abs(nint(dt_max/dt_min)-dt_max/dt_min).GT.1.d-6) THEN
            WRITE(lu,*) ''
            WRITE(lu,*) '***************************************'
            WRITE(lu,*) ' ATTENTION : COUPLING TELEMAC-TOMAWAC :'
            WRITE(lu,*) ' THE CHOSEN TIME STEPS ARE NOT MULTIPLE'
            WRITE(lu,*) ' OF EACH OTHER.                        '
            WRITE(lu,*) '         END OF THE COMPUTATION        '
            WRITE(lu,*) '***************************************'
            CALL plante(1)
            stop
          ENDIF
        ENDIF
!
!.....1.4 INITIALISATION DES TABLEAUX DATE ET TIME
!
        adc=int(ddc*1.d-8)
        mdc=int(ddc*1.d-6)
        jdc=int(ddc*1.d-4)
        hdc=int(ddc*1.d-2)
        date(1)=adc
        date(2)=mdc-100*adc
        date(3)=jdc-100*mdc
        time(1)=hdc-100*jdc
        time(2)=int(ddc-100.d0*hdc)
        time(3)=0
!
!.....1.5 READ ASCII FILES WITH COORDINATES OF POINTS FOR WHICH
!     SPECTRA WILL NEED TO BE READ
!
        IF((namwxy(1:1).NE.' ').AND.
     &       (part.EQ.wac_full_run.OR.part.EQ.wac_cpl_init.OR.
     &       part.EQ.wac_api_init))THEN
          IF(debug.GT.0)
     &         WRITE(lu,*)'CALLING READ_SPECTRA_COORDS FOR PRINTOUT'
          CALL read_spectra_coords(luwxy, npleo, xleo, yleo)
          IF(debug.GT.0)WRITE(lu,*)'BACK FROM  READ_SPECTRA_COORDS'
!     ALLOCATE THE MEMORY FOR THE NODE NUMBERS OF THE SPECTRA
          IF(ALLOCATED(noleo))DEALLOCATE(noleo)
          ALLOCATE(noleo(npleo))
        ENDIF
!
!.....1.5 READ ASCII FILES WITH COORDINATES OF POINTS FOR WHICH
!     SPECTRA WILL NEED TO BE READ
!
        IF((namixy(1:1).NE.' ').AND.
     &       (part.EQ.wac_full_run.OR.part.EQ.wac_cpl_init.OR.
     &       part.EQ.wac_api_init))THEN
          IF(debug.GT.0)
     &    WRITE(lu,*)'CALLING READ_SPECTRA_COORDS FOR IMPOSED SPECTRA'
          CALL read_spectra_coords(luixy, npspe,xspe,yspe)
          IF(debug.GT.0)WRITE(lu,*)'BACK FROM  READ_SPECTRA_COORDS'
        ENDIF
!
        IF (inclus(coupling,'TOMAWAC2')) THEN
          CALL get_mesh_nptir(wac_files(wacgeo)%FMT,
     &       wac_files(wacgeo)%LU,nptir,ierr)
        ENDIF
!=====C
!     2 LECTURE DES CONDITIONS LIMITES ET INDICES DES POINTS FRONTIERES
!=====C================================================================
!
!     COUPLAGE TELEMAC-TOMAWAC : initialisation
!
        IF(part.EQ.wac_full_run.OR.part.EQ.wac_cpl_init.OR.
     &       part.EQ.wac_api_init) THEN
          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WAC_INIT'
          CALL wac_init (part,impres,debres,date,time)
          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WAC_INIT'
        ENDIF
      ENDIF ! PART.LE.2
!
!     COUPLAGE TELEMAC-TOMAWAC
      IF(part.NE.wac_cpl_init.AND.part.NE.wac_api_init) THEN
!
!     Preparation au couplage : on calcule le nombre de cycle
!     de la boucle en temps, selon la duree du pas de temps dans
!     Telemac (CPL_WAC_DATA%DT_TEL) et dans Tomawac (DT).
!     On assigne a DUMMY la valeur originaire du nombre de pas
!     de temps specifie dans le fichier .cas de Tomawac.
!
        IF(part.EQ.wac_cpl_run) THEN
          at=cpl_wac_data%AT_TEL
          dummy=nit
          IF(dt.GT.cpl_wac_data%DT_TEL*cpl_wac_data%PERCOU_WAC) THEN
            WRITE(lu,*) ''
            WRITE(lu,*) '***************************************'
            WRITE(lu,*) ' ATTENTION : COUPLING TELEMAC-TOMAWAC :'
            WRITE(lu,*) ' TOMAWAC TIME STEP CAN NOT BE GREATER  '
            WRITE(lu,*) ' THAN TELEMAC TIME STEP                '
            WRITE(lu,*) '         END OF THE COMPUTATION        '
            WRITE(lu,*) '***************************************'
            CALL plante(1)
            stop
          ENDIF
          IF(dt.LE.cpl_wac_data%DT_TEL*cpl_wac_data%PERCOU_WAC) THEN
            nit = nint(cpl_wac_data%DT_TEL*cpl_wac_data%PERCOU_WAC/dt)
          ENDIF

!       RECEIVE DATA FROM TOMAWAC
!                                                       WAC2
          IF(inclus(coupling,'TOMAWAC2')) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE RECEIVE_COUPLE'
            CALL receive_couple(1,npoin2,nvartel2tom,tel2tom)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE RECEIVE_COUPLE'
            cpl_wac_data%U_TEL => u_tel
            cpl_wac_data%V_TEL => v_tel
            cpl_wac_data%H_TEL => h_tel
            cpl_wac_data%FX_WAC => fx_wac
            cpl_wac_data%FY_WAC => fy_wac
            cpl_wac_data%UV_WAC => uv_wac
            cpl_wac_data%VV_WAC => vv_wac
            cpl_wac_data%DIRMOY_TEL => dirmoy_tel
            cpl_wac_data%HM0_TEL => hm0_tel
            cpl_wac_data%TPR5_TEL => tpr5_tel
            cpl_wac_data%ORBVEL_TEL => orbvel_tel
          ENDIF
        ENDif!     Fin COUPLAGE
!
!     COUPLAGE TELEMAC-TOMAWAC : LT est defini comme le numero de pas
!     de temps effectifs de TOMAWAC. La variable LT_WAC
!     compte les pas de temps de chaque boucle, meme
!     quand TOMAWAC est appele par TELEMAC
!
!=====C
!     11 C BOUCLE EN TEMPS PRINCIPALE.
!=====C============================
!
        IF(debug.GT.0) WRITE(lu,*) 'TIME LOOP BEGINNING'
        DO lt_wac=1,nit
!
!------------------------------------------------------------------
!
#if defined COMPAD
          CALL ad_tomawac_timestep_begin
#endif
!
!------------------------------------------------------------------
!
!.....11.1 AFFECTATION DE LA DATE DE FIN DU PAS DE TEMPS COURANT.
!
          at=at+dt
!     TODO: Calcul de LT (NOTE JMH: WHY NOT LT=LT+1 ?)
          lt=nint((at-at0)/dt)
!
!     Fin COUPLAGE
!
          CALL impr(lisprd,lt,at,lt,3)
!
!     11.2 AFFECTATION DES CONDITIONS AUX LIMITES.
!
          IF (avant) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE LIMWAC'
            CALL limwac(f,       fbor, nptfr, ndire,  nf, npoin2,
     &      kent, prive, npriv, wac_files(impspe))
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE LIMWAC'
          ENDIF
!
!     11.2 MISE A ZERO DU SPECTRE SUR LES POINTS OU PROF < PROMIN
!     I GUESS I DO NOT HAVE THOSE ANYMORE....SINCE I CHANGED DEPTH
          IF (.NOT.proinf) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE ECRETE'
            IF (ecret) CALL ecrete(sf%R,sdepth%R,npoin2,ndire,nf,
     &                             promin)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE ECRETE'
          ENDIF
!
          IF (maree.OR.couran) lt1=(lt/lam)*lam
!
!......11.3 UPDATING DEPTH AND CURRENTS
!
!     COUPLING TELEMAC-TOMAWAC OR CURRENTS AND/OR DEPTH IN A FILE
!     THEY ARE UPDATED HERE.
!
          IF(maree.AND.lt.EQ.lt1.OR.
     &         (part.EQ.wac_cpl_run.AND.lt_wac.EQ.1)) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE CORMAR'
            CALL cormar(part, cpl_wac_data%U_TEL,
     &           cpl_wac_data%V_TEL,cpl_wac_data%H_TEL)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE CORMAR'
            DO ip=1,npoin2
              IF(depth(ip).LT.promin) depth(ip)=0.9d0*promin
            ENDDO
!
!......11.3.1 PREPARING PROPAGATION (METHOD OF CHARACTERISTICS).
!
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE INIPHY'
            CALL iniphy(xk, cg, b, npoin2, nf)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE INIPHY'
!
            IF(prop) THEN
!
              CALL impr(lisprd,lt,at,lt,1)
              CALL impr(lisprd,lt,at,lt,2)
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PREPRO 2'
!
              CALL prepro( stsder, ststot, ikle2, ibor, elt, eta, fre,
     &             xk,     cg,     itr01, npoin3, npoin2, nelem2,
     &             ndire,  nf,    couran.OR.part.EQ.wac_cpl_run)
!     Fin COUPLAGE
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PREPRO 2'
            ENDIF
          ENDIF
!     Fin cycle IF((MAREE.AND.LT.EQ.LT1).OR.
!     (PART.EQ.WAC_CPL_RUN.AND.LT_WAC.EQ.1))
!------------------------------------------------------------------
!     V6P2 Diffraction : diffraction term calculation
          IF(diffra.GT.0) THEN
            IF(lt.EQ.nptdif)THEN
              WRITE(lu,*)'*********************************'
              WRITE(lu,*)'DIFFRACTION IS TAKEN INTO ACCOUNT'
              WRITE(lu,*)'*********************************'
            ENDIF
            IF(lt.GE.nptdif) THEN
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PREDIF'
              CALL predif( stsder, ststot, ikle2, ibor,   elt,    eta,
     &        xk,     cg,     itr01, npoin3, npoin2, nelem2, ndire,
     &        nf, couran.OR.part.EQ.wac_cpl_run, f,
     &        rx, ry,     rxx,    ryy,    neigb)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PREDIF'
            ENDIF
          ENDIF
!     V6P2 End diffraction
!-------------------------------------------------------------------
!
!.....11.3 PROPAGATION (INTERPOLATION AU PIED DES CARACTERISTIQUES).
!
          IF(prop) THEN
            CALL impr(lisprd,lt,at,lt,5)
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE PROPA',lt
            CALL propa(f, b,    elt, eta, fre ,npoin3, npoin2,
     &           ndire,nf, couran.OR.part.EQ.wac_cpl_run, tsder)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE PROPA'
          ENDIF
!
!.....11.4 ECRETE PAR RAPPORT A UN SPECTRE LIMITE APRES PROPAGATION.
!
          IF(sdscu.EQ.1) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE LIMITE'
            CALL limite(f, freq, npoin2, ndire, nf)
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE LIMITE'
          ELSEIF(sdscu.GT.2) THEN
            WRITE(lu,*) ''
            WRITE(lu,*) '**************************'
            WRITE(lu,*) '   UNKNOWN OPTION FOR     '
            WRITE(lu,*) '    STRONG CURRENTS       '
            WRITE(lu,*) '**************************'
            CALL plante(1)
            stop
          ELSEIF(sdscu.EQ.2.AND..NOT.tsou) THEN
            WRITE(lu,*) ''
            WRITE(lu,*) '****************************'
            WRITE(lu,*) '  CONSIDERATION OF SOURCE   '
            WRITE(lu,*) '   TERMS MANDATORY FOR      '
            WRITE(lu,*) '    OPTION 2 FOR STRONG     '
            WRITE(lu,*) '         CURRENTS           '
            WRITE(lu,*) '****************************'
            CALL plante(1)
            stop
          ENDIF
!
!.....11.5 INTEGRATION DES TERMES SOURCES.
!
          IF(tsou) THEN
            CALL impr(lisprd,lt,at,nsits,4)
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SEMIMP'
            CALL semimp( f,    cf,     xk, nf, ndire, npoin2,
     &      iangnl, tstot, tsder, told,  tnew, tra35, tra36,
     &      tra37,  tra38, tra39, t1,    t2,   t3,     t4,
     &      mdia,  ianmdi, coemdi, fbor)

            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SEMIMP'
          ENDIF
!
!     11.8 AFFECTATION DES CONDITIONS AUX LIMITES.
!
          IF(.NOT.avant) THEN
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE LIMWAC'
            CALL limwac(f,   fbor, nptfr, ndire,  nf, npoin2, kent,
     &      prive, npriv,  wac_files(impspe))
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE LIMWAC'
          ENDIF
!
!.....10.6 TEST POUR SAVOIR SI ON IMPRIME OU PAS.
!     """""""""""""""""""""""""""""""""""""""""""
          impres=.false.
          debres=.false.
          IF(lt.GE.gradeb.AND.mod(lt-gradeb,graprd).EQ.0) impres=.true.
          IF(lt.EQ.gradeb) debres=.true.
!
          IF(impres) THEN
!
!.....11.7 PASSAGE EN FREQUENCE ABSOLUE.
!     """""""""""""""""""""""""""""""""""""""""""""
!
            IF(couran.OR.part.EQ.wac_cpl_run) THEN
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TRANSF'
              IF(lt.GE.gradeb.AND.mod(lt-gradeb,graprd).EQ.0) THEN
                CALL transf(tstot, f,   xk, itr11, itr12,itr13, tra31,
!     WORK ARRAY UNTIL CALL ECRSPE BELOW
     &               tra32, npoin2, ndire,nf)
              ENDIF
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TRANSF'
            ENDIF
!
!.....11.9 IMPRESSION (EVENTUELLE) DES VARIABLES SUR LE MAILLAGE 2D.
!
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE DUMP2D'
            IF(couran.OR.part.EQ.wac_cpl_run) THEN
              CALL dump2d(tstot, npoin3*nf)
            ELSE
              CALL dump2d(f, npoin3*nf)
            ENDIF
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE DUMP2D'
!
            IF(namres.NE.' ') THEN
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE BIEF_DESIMP'
              CALL bief_desimp(fmtres,varsor, npoin2,lures,at, lt,
     &             lisprd,graprd, sorleo,sorimp,maxvar,texte,gradeb,
     &         gradeb)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIEF_DESIMP'
            ENDIF
!
!.....11.10 IMPRESSION (EVENTUELLE) DES SPECTRES DIRECTIONNELS.
!
            IF(debug.GT.0) WRITE(lu,*) 'APPEL DE ECRSPE'
            irec = (lt - gradeb)/graprd
            IF(couran.OR.part.EQ.wac_cpl_run) THEN
              CALL ecrspe(tstot, ndire, nf, npoin2, irec,  tsder,
     &        noleo, npleo, debres,date,  time, mesh%KNOLG%I, mesh)
            ELSE
              CALL ecrspe(f,     ndire, nf, npoin2, irec,  tsder,
     &        noleo, npleo, debres,date,  time, mesh%KNOLG%I, mesh)
            ENDIF
            IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE ECRSPE'
!
          ENDIF
          IF(part.EQ.wac_cpl_run.AND.lt_wac.EQ.nit) THEN
            IF(cpl_wac_data%COUPL3D) THEN
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE TRANSF'
              CALL transf(tstot, f,    xk,    itr11, itr12,itr13,
!     WORK ARRAY UNTIL CALL ECRSPE BELOW TRA32 IS USED IN FRIC3D
     &        tra31, tra32, npoin2, ndire, nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE TRANSF'
!     DISSIPATION DUE TO DEPTH-INDUCED WAVE BREAKING
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FDISS3D'
              CALL fdiss3d ( cpl_wac_data%FDXW%R, cpl_wac_data%FDYW%R,
     &             npoin2, xk, ndire, f, nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FDISS3D'

!     DISSIPATION DUE TO BOTTOM FRICTION
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FBOTT3D'
              CALL fbott3d
     &             ( cpl_wac_data%FBXW%R, cpl_wac_data%FBYW%R, f,
     &             npoin2, xk, ndire,nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FBOTT3D'
!     The velocity STOKES DRIFT HORIZONTAL COMPONENTS
!     ONE SHOULD CALL TRANSF BEFORE !!!
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE UVSTOKES'
              CALL uvstokes( cpl_wac_data%USTW%R, cpl_wac_data%VSTW%R,
     &             cpl_wac_data%WSTW%R, tstot,  npoin2,xk,zf, ndire,
     &             cpl_wac_data%ZTELW%R, cpl_wac_data%NZW,nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE UVSTOKES'

!     The wave induced pressure J
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WIPJ'
              CALL wipj (cpl_wac_data%WIPW%R, tstot, npoin2, xk,
     &        cpl_wac_data%WIPDXW%R,cpl_wac_data%WIPDYW%R, ndire, nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WIPJ'
!
!     WAVE ENHANCED MIXING
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FDISSK'
              CALL fdissk( cpl_wac_data%FDKW%R, npoin2, ndire, f,
     &             cpl_wac_data%ZTELW%R,cpl_wac_data%NZW, tra38, tra37,
     &             nf)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FDISSK'

!     FRICTION COEFFICIENT - WAVES+CURRENTS

              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE VITFON'
              CALL vitfon(vifond, tstot, xk, nf,
     &             npoin2,ndire)
              IF(debug.GT.0) WRITE(lu,*) 'APPEL DE FRIC3D'
!     MODIFIED COEFFICIENT TO TAKE INTO ACCOUNT WAVES+CURRENTS
              CALL fric3d
     &             ( cpl_wac_data%CFWCW, npoin2, tra32,
     &             cpl_wac_data%U_TEL, cpl_wac_data%V_TEL, vifond)
              IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE FRIC3D'
!     SURFACE STRESS DUE TO WIND INPUT ENERGY AND WHITECAPPING
              IF(vent) THEN
                IF(svent.EQ.1) THEN
                  IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WINDISS1'
                  CALL windiss1( cpl_wac_data%FWX%R, cpl_wac_data%FWY%R,
     &                 npoin2, xk, ndire, f, nf)
                  IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WINDISS1'
                ELSEIF(svent.EQ.2) THEN
                  IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WINDISS2'
                  CALL windiss2( cpl_wac_data%FWX%R, cpl_wac_data%FWY%R,
     &                 npoin2, xk, ndire, f, nf)
                  IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WINDISS2'
                ELSEIF(svent.EQ.3) THEN
                  IF(debug.GT.0) WRITE(lu,*) 'APPEL DE WINDISS3'
                  CALL windiss3( cpl_wac_data%FWX%R, cpl_wac_data%FWY%R,
     &                 npoin2, xk, ndire, f, nf)
                  IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE WINDISS3'
                ENDIF
                IF(smout.EQ.1) THEN
                  CALL moudiss1( cpl_wac_data%FWX%R, cpl_wac_data%FWY%R,
     &                 npoin2, xk, ndire, f, nf)
                ELSEIF(smout.EQ.2)THEN
                  CALL moudiss2( cpl_wac_data%FWX%R, cpl_wac_data%FWY%R,
     &                 npoin2, xk, ndire, f, nf, t1, t2)
                ENDIF
              ENDIF

            ENDIF               !IF (CPL_WAC_DATA%COUPL3D)
!
!     RADIATION STRESSES COMPUTED HERE FOR TELEMAC INDEPENDENTLY
!     OF THE PRINTS TO RESULT FILE, WHICH IS NOT MANDATORY
!
!     STSTOT WORK ARRAY (ABSOLUTE FREQUENCY) IN ALL THIS IF BLOCK
            CALL transf(tstot, f, xk, itr11, itr12,itr13, tra31, tra32,
     &           npoin2,ndire, nf)
            CALL radiat(cpl_wac_data%FX_WAC%R, cpl_wac_data%FY_WAC%R,
     &           xk, tstot, cg, tsder, tra36, tra37, tra38, tra39)
            IF(vent) THEN
              CALL ov('X=Y     ', x=cpl_wac_data%UV_WAC%R, y=suv%R,
     &             dim1=npoin2)
              CALL ov('X=Y     ', x=cpl_wac_data%VV_WAC%R, y=svv%R,
     &             dim1=npoin2)
            ENDIF
            IF(inclus(coupling,'SISYPHE').OR.inclus(coupling,'GAIA'))
     &      THEN
!     3 VARIABLES THAT WILL BE TRANSMITTED TO SISYPHE OR GAIA
!     ALL THIS IF BLOCK ADAPTED FROM DUMP2D
!     MEAN DIRECTION
              CALL tetmoy(cpl_wac_data%DIRMOY_TEL%R, tstot, ndire,nf,
     &         npoin2)
            IF(trigo) THEN
              DO ip=1,npoin2
                cpl_wac_data%DIRMOY_TEL%R(ip) =
     &               (pisur2-cpl_wac_data%DIRMOY_TEL%R(ip))*gradeg
              ENDDO
            ELSE
              DO ip=1,npoin2
                cpl_wac_data%DIRMOY_TEL%R(ip) =
     &               cpl_wac_data%DIRMOY_TEL%R(ip)*gradeg
              ENDDO
            ENDIF
!     SIGNIFICANT WAVE HEIGHT
            CALL totnrj(tra37, tstot, nf, ndire, npoin2)
            DO ip=1,npoin2
              cpl_wac_data%HM0_TEL%R(ip)=4.d0*sqrt(tra37(ip))
            ENDDO
!     TPR5
            CALL fpread(cpl_wac_data%TPR5_TEL%R, tstot,
     &           nf, ndire, npoin2, 5.d0)
            DO ip=1,npoin2
              cpl_wac_data%TPR5_TEL%R(ip)=
     &      1.d0/min(max(cpl_wac_data%TPR5_TEL%R(ip),freq(1)),freq(nf))
            ENDDO
!     NEAR BED ORBITAL VELOCITY
            CALL vitfon(cpl_wac_data%ORBVEL_TEL%R, tstot, xk,  nf,
     &           npoin2,ndire)
          ENDIF
          IF (alg_dislodge) THEN
            CALL vitfon(cpl_wac_data%ORBVEL_TEL%R, tstot, xk,  nf,
     &           npoin2,ndire)
          ENDIF
!
        ENDIF
!
!------------------------------------------------------------------
!
#if defined COMPAD
        CALL ad_tomawac_timestep_end
#endif
!
!------------------------------------------------------------------
!
      ENDDO                     ! LT_WAC
!
!     Fin COUPLAGE : fin du cycle IF(PART.NE.0)
      ENDIF
!
!=====C
!     12 C IMPRESSIONS GLOBALES (EVENTUELLES) EN FIN DE CALCUL.
!=====C=====================================================
!
!     COUPLAGE TELEMAC-TOMAWAC : reset de la variable NIT (nombre
!     de pas de temps de TOMAWAC indique dans la fichier .cas)
      IF(part.EQ.wac_cpl_run) nit=dummy
      IF(PRESENT(nit_ori)) THEN
        total_iter = nit_ori
      ELSE
        total_iter = nit
      ENDIF
!     Fin COUPLAGE
!
      IF(abs(at-at0-total_iter*dt).LT.1.d-6) THEN
!
        IF(glob) THEN
          CALL impr(lisprd,nit,at,nit,6)
          IF(debug.GT.0) WRITE(lu,*) 'APPEL DE SOR3D'
          CALL sor3d(f,  ndire, nf, npoin2,  vent,
     &         couran.OR.part.EQ.wac_cpl_run,
     &         maree.OR.part.EQ.wac_cpl_run, titcas, tsder, mesh3d)
          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE SOR3D'
        ENDIF
!
!----------------------------------------------------------------------
!
!     VALIDATION DES RESULTATS SUR LE FICHIER DE REFERENCES
!
        IF(valid) THEN
          IF(debug.GT.0) WRITE(lu,*) 'CALLING BIEF_VALIDA'
          CALL bief_valida(bst1,  texte,  luref, fmtref,
     &         varsor, texte, lures, fmtres,
     &         maxvar,npoin2,   nit,    nit,alire)
          IF(debug.GT.0) WRITE(lu,*) 'RETOUR DE BIEF_VALIDA'
        ENDIF
!
      ENDIF
!
!------------------------------------------------------------------
!
#if defined COMPAD
      CALL ad_tomawac_end
#endif
!
!------------------------------------------------------------------
!------------------------------------------------------------------------
!
!   SEND DATA BACK TO TELEMAC
!
!                                                       WAC2
      IF (inclus(coupling,'TOMAWAC2')) THEN
        !CONVERT DIRECTIONS
        DO i=1,npoin2
          cosdir_tel%R(i) = cos(deg2rad*dirmoy_tel%R(i))
          sindir_tel%R(i) = sin(deg2rad*dirmoy_tel%R(i))
        ENDDO
        CALL send_couple(2,npoin2,nvartom2tel,tom2tel)
      ENDIF
!
      RETURN
      END