semimp.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE semimp
!                   *****************
!
     & (f    , cf   ,  xk   ,  nf   , ndire, npoin2,
     &  iangnl,tstot,  tsder,  told , tnew,  z0new, twnew,
     &  taux1, taux2,  taux3,  taux4, taux5, taux6, taux7,
     &  mdia,  ianmdi, coemdi, fbor)
!
!***********************************************************************
! TOMAWAC   V7P3
!***********************************************************************
!
!brief    SOLVES THE INTEGRATION STEP OF THE SOURCE TERMS USING
!+                A SCHEME WITH VARIABLE DEGREE OF IMPLICITATION.
!
!history  M. BENOIT
!+        26/03/95
!+        V1P0
!+   CREATED
!
!history  M. BENOIT
!+        07/11/96
!+        V1P2
!+   MODIFIED
!
!history
!+        25/08/2000
!+        V5P0
!+   MODIFIED
!
!history  JMH
!+        16/12/2008
!+        V5P9
!+   BETA HAS BEEN ADDED TO THE LIST OF ARGUMENTS AND
!
!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  E.G.RENOU (EDF), G.MATTAROLO (EDF)
!+        12/05/2011
!+        V6P1
!+   MODIFIED: integration of new source terms, developed by
!+   E.G. Renou.
!+     - modification of the variables in argument list
!+     - modification of the local variable declarations
!+     - modification concerning friction velocity and roughness
!+       length calculation
!+     - calls to subroutines QWINDL, QWIND3, QMOUT2, QNLIN2,
!+       QNLIN3
!
!history  G.MATTAROLO (EDF - LNHE)
!+        27/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!history  U.H.MERKEL
!+        27/06/2012
!+        V6P2
!+   Renamed SUM to SOMME, due to NAG compiler
!
!history  E. GAGNAIRE-RENOU (EDF - LNHE)
!+        12/03/2013
!+        V6P3
!+   HF diagnostic tail is not necessarily imposed
!
!history  VITO BACCHI (EDF - LNHE)
!+        12/09/2014
!+        V7P0
!+   Friction due to vegetation added.
!
!history THIERRY FOUQUET (EDF-LNHE)
!+       19/11/2014
!+       V7P0
!+   BAJ MODELING
!
!history A JOLY (EDF-LNHE)
!+       18/05/2017
!+       V7P3
!+   New condition to stop source terms being added to open
!+   boundaries.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CF             |-->| ADVECTION FIELD ALONG FREQUENCY
!| COEMDI         |-->| COEFFICIENTS USED FOR MDIA METHOD
!| F              |<->| DIRECTIONAL SPECTRUM
!| FBOR           |-->| SPECTRAL VARIANCE DENSITY AT THE BOUNDARIES
!| FMOY           |<--| MEAN FREQUENCIES F-10
!| FREQ           |-->| DISCRETIZED FREQUENCIES
!| IANGNL         |-->| ANGULAR INDICES TABLE
!| IANMDI         |-->| ANGULAR INDICES TABLE FOR MDIA
!| MDIA           |-->| NUMBER OF CONFIGURATIONS FOR MDIA METHOD
!| NBOR           |-->| GLOBAL NUMBER OF BOUNDARY POINTS
!| NF             |-->| NUMBER OF FREQUENCIES
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| NPTFR          |-->| NUMBER OF BOUNDARY POINTS
!| TAUX1          |<->| WORK TABLE
!| TAUX2          |<->| WORK TABLE
!| TAUX3          |<->| WORK TABLE
!| TAUX4          |<->| WORK TABLE
!| TAUX5          |<->| WORK TABLE
!| TAUX6          |<->| WORK TABLE
!| TAUX7          |<->| WORK TABLE
!| TNEW           |<->| WORK TABLE
!| TOLD           |<->| WORK TABLE
!| TSDER          |<--| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
!| TSTOT          |<--| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
!| TWNEW          |<->| WIND DIRECTION AT TIME N+1
!| USNEW          |<->| FRICTION VELOCITY AT TIME N+1
!| VARIAN         |-->| SPECTRUM VARIANCE
!| XK             |-->| DISCRETIZED WAVE NUMBER
!| XKMOY          |<--| AVERAGE WAVE NUMBER
!| Z0NEW          |<->| SURFACE ROUGHNESS LENGTH AT TIME N+1
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac, ONLY : deupi,t3_01,t3_02,texveb,
     &                                 namewx,namewy,unitveb,phasveb,
     &                                 source_on_bnd,debug,cg,porous,
     &                                 amorp,lt, gravit, cimpli,
     &                                 ifrbj, ifrtg , ifrro, ifrih,
     &                                 diaghf, coefhs, nsits,  smout,
     &  fmoy, varian, xkmoy, sfrot, svent, lvent, strif, vent, vensta,
     &  luvef,  luveb, namvef, namveb, fmtveb, fmtvef,
     &  sbrek,  xdtbrk, ndtbrk, stria, proinf, df_lim, limit,
     &  vegetation, cbaj, sdscu,freq, dfreq, depth,
     &  teta, tailf, raisf, nbor, nptfr, lifbor,
     &  tv1   ,tv2   , uv1, uv2, vv1, vv2, uv, vv,
     &  tauwav,usold, usnew, twold, z0old, at, dtsi, indiv, nvwin
! FROM DECLARATION_TOMAWAC :
! DIAGHF          OPTION FOR SPECTRUM DIAGNOSTIC TAIL
! CBAJ            CHOICE OF THE CENTRAL FREQUENCY CALCULUS
! COEFHS          MAXIMUM VALUE OF THE RATIO HM0 ON D
! DF_LIM          WORK TABLE
! DTSI            INTEGRATION TIME STEP (SECONDS)
! F1              MINIMAL DISCRETIZED FREQUENCY
! INDIV           WIND FILE FORMAT
! LIMIT           TYPE OF WAVE GROWTH LIMITER MODEL SELECTED
! LVENT           LINEAR WAVE GROWTH MODEL SELECTION
! NDTBRK          NUMBER OF TIME STEPS FOR BREAKING SOURCE TERM
! NSITS           NUMBER OF ITERATIONS FOR THE SOURCE TERMS
! PROINF          LOGICAL INDICATING INFINITE DEPTH ASSUMEPTION
! RAISF           FREQUENTIAL RATIO
! SBREK           DEPTH-INDUCED BREAKING DISSIPATION MODEL
! SDSCU           DISSIPATION BY STRONG CURRENT
! SMOUT           SELECTION OF WHITE CAPPING SOURCE TERM MODEL
! STRIA           SELECTION OF THE TRIAD INTERACTION MODEL
! STRIF           SELECTION OF QUADRUPLET INTERACTION MODEL
! SFROT           SELECTION OF THE BOTTOM FRICTION DISSIPATION
! SVENT           SELECTION OF THE WIND GENERATION MODEL
! TAILF           SPECTRUM TAIL FACTOR
! TAUWAV          STRESS DUE TO THE WAVES
! TETA            DISCRETIZED DIRECTIONS
! TWOLD           WIND DIRECTION AT TIME N
! TNEW            WORK TABLE
! TOLD            WORK TABLE
! TV1             TIME T1 IN THE WIND FILE
! TV2             TIME T2 IN THE WIND FILE
! USOLD           FRICTION VELOCITY AT TIME N
! U1,V1           WIND SPEED AT TIME T1 IN THE WIND FILE
! U2,V2           WIND SPEED AT TIME T2 IN THE WIND FILE
! UV,VV           WIND DATA INTERPOLATED OVER 2D MESH
! VEGETATION      IF YES, VEGETATION TAKEN INTO ACCOUNT
! VENSTA          INDICATES IF THE WIND IS STATIONARY
! VENT            INDICATES IF WIND IS TAKEN INTO ACC
! XDTBRK          COEFFICIENT OF TIME SUB-INCREMENTS FOR BREAKING
! Y               ORDINATES OF POINTS IN THE MESH
! Z0OLD           SURFACE ROUGHNESS LENGTH AT TIME N
      USE interface_tomawac, ex_semimp => semimp
!
      USE declarations_telemac
      USE declarations_special
      IMPLICIT NONE
!
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NPOIN2,NDIRE,NF, IANGNL(*)
      DOUBLE PRECISION, INTENT(INOUT) :: Z0NEW(npoin2), TWNEW(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: TAUX1(npoin2), TAUX2(npoin2),
     &                                   taux3(npoin2), taux4(npoin2),
     &                                   taux5(npoin2), taux6(npoin2),
     &                                   taux7(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: F(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TOLD(npoin2,ndire)
      DOUBLE PRECISION, INTENT(INOUT) :: TNEW(npoin2,ndire)
      DOUBLE PRECISION, INTENT(INOUT) :: XK(npoin2,nf)
      DOUBLE PRECISION, INTENT(IN) :: CF(*)
!....MDIA method declarations
      INTEGER, INTENT(IN) :: MDIA, IANMDI(*)
      DOUBLE PRECISION, INTENT(IN) ::  COEMDI(*)
      DOUBLE PRECISION, INTENT(IN)   :: FBOR(nptfr,ndire,nf)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER ISITS,IFF,IP,JP,K,NVENT,IFCAR,MF1,MF2,MFMAX,IDT
      INTEGER IPTFR
      DOUBLE PRECISION AUX1,AUX2,AUX3,AUX4,COEF
      DOUBLE PRECISION FM1,FM2,TDEB,TFIN,VITVEN
      DOUBLE PRECISION VITMIN,HM0,HM0MAX,DTN,SOMME,AUXI,USMIN
!                              MDIA, HERE HARDCODED
      DOUBLE PRECISION  XCCMDI(4)
!
      LOGICAL TROUVE(3)
!
      DOUBLE PRECISION  CPHAS , SEUILF
!
      CHARACTER(LEN=8) FMTVEN
!
!-----------------------------------------------------------------------
!
      vitmin=1.d-3
!
!     ------------------------------------------------------------------
!     CHOPS THE SPECTRUM IN ACCORDANCE WITH THE BATHYMETRY
!     -----------------------------------------------------------------
!
      IF(.NOT.proinf) THEN
!
!       0.1 COMPUTES THE TOTAL VARIANCE OF THE SPECTRUM
!       -----------------------------------------------
!
        IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE TOTNRJ'
        CALL totnrj(varian, f, nf, ndire, npoin2)
        IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE TOTNRJ'
!
!       0.2 COMPUTES THE CORRECTION COEFFICIENT ON THE SPECTRUM
!       -------------------------------------------------------
!
        DO ip=1,npoin2
          hm0max=coefhs*depth(ip)
          hm0 =max(4.d0*sqrt(varian(ip)),1.d-20)
          taux1(ip)=min((hm0max/hm0)**2,1.d0)
        ENDDO
!
!       0.3 CORRECTS THE SPECTRUM
!       --------------------------
!
        DO iff=1,nf
          DO jp=1,ndire
            DO ip=1,npoin2
              f(ip,jp,iff)=f(ip,jp,iff)*taux1(ip)
            ENDDO
          ENDDO
        ENDDO
!
      ENDIF
!
!     ----------------------------------------------------------------
!     IF THE COMPUTATION INCLUDES STATIONARY WINDS, DUPLICATES THE
!     CONDITIONS AT THE START OF THE TIME STEP TO THE END OF THE TIME
!     STEP. (THIS IS BECAUSE ARRAYS TWNEW, USNEW AND Z0NEW ARE WORKING
!     ARRAYS USED IN DUMP2D BETWEEN 2 CALLS TO SEMIMP).
!
!     TODO: QUESTION JMH 16/09/2014 : WHAT IF VENT=.FALSE. AND USNEW USED
!     IN BETWEEN ???? IT WOULD BE SAFER TO HAVE USOLD AND USNEW
!     WITH THIS NAME AND ONLY USED FOR THIS...
!
!     ----------------------------------------------------------------
!
      IF(vent.AND.vensta) THEN
        DO ip=1,npoin2
          twnew(ip)=twold(ip)
        ENDDO
!
        IF (svent.GE.2.OR.(lvent.EQ.1.AND.svent.NE.1).OR.
     &                                (smout.EQ.2.AND.svent.NE.1)) THEN
          DO ip=1,npoin2
            usnew(ip)=usold(ip)
            z0new(ip)=z0old(ip)
          ENDDO
        ENDIF
      ENDIF
!
!     -----------------------------------------------------------------
!     START OF THE MAIN LOOP ON THE NUMBER OF TIME STEPS (NSITS)
!     FOR INTEGRATION OF THE SOURCE TERMS, BY PROPAGATION TIME STEP
!     -----------------------------------------------------------------
!
      DO isits=1,nsits
!
!       1. ASSIGNS THE START AND END DATES OF TIME STEP
!       ===============================================
!
        tdeb=at-(nsits-isits+1)*dtsi
        tfin=tdeb+dtsi
!
!
!       2. UPDATES (IF HAS TO) THE WIND ARRAYS
!       ======================================
!
        IF(vent.AND..NOT.vensta) THEN
!
!         2.1 UPDATES THE WIND FIELD FOR DATE TFIN
!         ---------------------------------------------------
!
          IF(namveb(1:1).NE.' '.OR.namvef(1:1).NE.' ') THEN
            IF(namvef(1:1).NE.' ') THEN
              nvent=luvef
              fmtven=fmtvef
            ELSE
              nvent=luveb
              fmtven=fmtveb
            ENDIF
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE NOUDON'
            CALL noudon(uv,namewx,2, vv,namewy,2,
     &                  vv,'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX',0,
     &                  npoin2, nvent,fmtven,tfin,tv1,tv2,
     &                  uv1,uv2,vv1,vv2,vv1,vv2,indiv,
     &                  'WIND   ',nvwin,texveb,trouve,unitveb,phasveb)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE NOUDON'
          ELSE
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE ANAVEN'
            CALL anaven
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE ANAVEN'
          ENDIF
!
!         2.2 COMPUTES THE WIND DIRECTION
!         -----------------------------------
!
          DO ip=1,npoin2
            vitven=sqrt(uv(ip)**2+vv(ip)**2)
            IF(vitven.GT.vitmin) THEN
              twnew(ip)=atan2(uv(ip),vv(ip))
            ELSE
              twnew(ip)=0.d0
            ENDIF
          ENDDO
!
!         2.3 COMPUTES THE FRICTION VELOCITIES AND ROUGHNESS LENGTHS
!         ------------------------------------------------------------
!
          IF(svent.GE.2.OR.(lvent.EQ.1.AND.svent.NE.1).OR.
     &        (smout.EQ.2.AND.svent.NE.1)) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE USTAR2'
            CALL ustar2( usnew, npoin2)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE USTAR2'
          ENDIF
!
        ENDIF
!
        IF(vent.AND.svent.EQ.1) THEN
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE USTAR1'
          CALL ustar1(usnew, z0new, tauwav, npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE USTAR1'
        ENDIF
!
!
!       3. COMPUTES MEAN PARAMETERS FOR THE DIRECTIONAL SPECTRUM
!       =========================================================
!
!       3.1 COMPUTES THE TOTAL VARIANCE OF THE SPECTRUM
!       -----------------------------------------------
!
        IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE TOTNRJ'
        CALL totnrj(varian, f, nf, ndire, npoin2)
        IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE TOTNRJ'
!
        IF (cbaj.EQ.0) THEN
!
!       3.2 COMPUTES THE MEAN FREQUENCY OF THE SPECTRUM
!       -----------------------------------------------
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREMOY'
          CALL fremoy(fmoy, f, nf, ndire, npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREMOY'
!
!       3.3 COMPUTES THE MEAN WAVE NUMBER OF THE SPECTRUM
!       -------------------------------------------------
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE KMOYEN'
          CALL kmoyen (xkmoy, xk , f, nf, ndire, npoin2,
     &         taux1 , taux2 , taux3 )
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE KMOYEN'
        ELSEIF (cbaj.EQ.1) THEN
!
!       3.2 COMPUTES THE MEAN FREQUENCY OF THE SPECTRUM
!       -----------------------------------------------
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREM01'
          CALL frem01 (fmoy, f, nf, ndire, npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREM01'
!
!       3.3 COMPUTES THE MEAN WAVE NUMBER OF THE SPECTRUM
!       -------------------------------------------------
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE KMOYE2'
          CALL kmoye2(xkmoy, xk, f, nf, ndire , npoin2,
     &         taux1 , taux2 , taux3 )
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE KMOYE2'
        ELSE
          WRITE(lu,*) 'UNKNOWN VALUE OF BAJ:',cbaj
          CALL plante(1)
          stop
        ENDIF
!
!       4. COMPUTES THE CONTRIBUTIONS OF THE SOURCE TERMS FOR GENERATION,
!          WHITECAPPING AND INTERACTIONS BETWEEN QUADRUPLETS
!       =============================================================
!
!       4.1 INITIALISES THE ARRAYS FOR THE SOURCE TERMS
!       ----------------------------------------------------
        DO iff=1,nf
          DO jp=1,ndire
            DO ip=1,npoin2
              tstot(ip,jp,iff)=0.d0
              tsder(ip,jp,iff)=0.d0
            ENDDO
          ENDDO
        ENDDO
!
!       4.2 GENERATION BY WIND
!       ----------------------
!
        IF(vent) THEN
          IF(svent.EQ.1) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QWIND1'
            CALL qwind1
     &( tstot , tsder , f     , xk    , usold , usnew , twold , twnew ,
     &  z0old , z0new , nf    , ndire , npoin2, told  , tnew  ,
     &  taux2 , taux3 , taux4 , taux5 , taux6 , taux7 )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QWIND1'
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE STRESS'
            CALL stress
     &( tauwav, tstot , f     , usnew , twnew , z0new ,
     &  npoin2, ndire , nf    , taux1 , taux2 , taux3 )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE STRESS'
          ELSEIF(svent.EQ.2) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QWIND2'
            CALL qwind2
     &( tstot , tsder , f     , xk    , usold , usnew , twold , twnew ,
     &  nf    , ndire , npoin2, t3_01%R,t3_02%R )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QWIND2'
          ELSEIF(svent.EQ.3) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QWIND3'
            CALL qwind3
     &( tstot , tsder , f     , xk    , usold , usnew , twold , twnew ,
     &  nf    , ndire , npoin2, taux1 , taux2 , taux3 , taux4 )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QWIND3'
!
          ENDIF
!
!       ADDS THE LINEAR WIND GROWTH SOURCE TERME
!       """""""""""""""""""""""""""""""""""""""
!
          IF(lvent.EQ.1) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QWINDL'
            CALL qwindl(tstot, usold, usnew, twold, twnew, nf, ndire,
     &                  npoin2, t3_01%R, t3_02%R, taux5, taux6)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QWINDL'
          ENDIF
!
        ELSE
!
          DO ip=1,npoin2
            usnew(ip)=0.d0
          ENDDO
!
        ENDIF

        IF (cbaj.EQ.1) THEN
! Calculation pf mean frequency fmean_WS put in the  tabular TAUX7
          seuilf = 1.d-20
          DO ip=1,npoin2
            taux1(ip) = 0.0d0
            taux2(ip) = 0.0d0
          ENDDO
!
          DO iff=1,nf
            aux3=deupi/dble(ndire)*dfreq(iff)
            aux4=aux3/freq(iff)
            DO jp=1,ndire
              DO ip=1,npoin2
                cphas=deupi*freq(iff)/xk(ip,iff)
                auxi=28.0d0/cphas*usnew(ip)*cos(teta(jp)-twnew(ip))
                IF ((tstot(ip,jp,iff).GT.0).OR.(auxi.GE.1.0d0)) THEN
                  taux1(ip) = taux1(ip) + f(ip,jp,iff)*aux3
                  taux2(ip) = taux2(ip) + f(ip,jp,iff)*aux4
                ENDIF
              ENDDO
            ENDDO
          ENDDO
!
          DO ip=1,npoin2
            IF (taux1(ip).LT.seuilf) THEN
              taux7(ip) = seuilf
            ELSE
              taux7(ip) = taux1(ip)/taux2(ip)
            ENDIF
          ENDDO
        ENDIF

!       4.3 NON-LINEAR INTERACTIONS BETWEEN QUADRUPLETS
!       -----------------------------------------------
!
        IF(strif.EQ.1) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QNLIN1'
          CALL qnlin1(tstot, tsder, iangnl, nf, ndire, npoin2,
     &                f, xkmoy, taux1, taux2, taux3, taux4, taux5,taux6)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QNLIN1'
!
        ELSEIF (strif.EQ.2) THEN
!
!         sets XCCMDI values for MDIA method
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QNLIN2'
          xccmdi(1)=8.360d7
          xccmdi(2)=7.280d7
          xccmdi(3)=3.340d7
          xccmdi(4)=2.570d6
          DO k=1,mdia
            xccmdi(k)=xccmdi(k)/dble(mdia)
          ENDDO
!         alls MDIA method
          DO k=1,mdia
            CALL qnlin2
     &( tstot , tsder , ianmdi((k-1)*ndire*16+1:k*ndire*16) ,
     &  coemdi((k-1)*32+1:k*32) , nf    , ndire,
     &  npoin2, f   , xkmoy , taux1 , taux2 , xccmdi(k))
          ENDDO
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QNLIN2'
!....calls GQM method
        ELSEIF (strif.EQ.3) THEN
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QNLIN3'
          CALL qnlin3(tstot , tsder , f     ,npoin2, ndire , nf )
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QNLIN3'
        ENDIF
!
!       4.4 WHITE-CAPPING DISSIPATION
!       -------------------------------------------------
!
        IF(smout.EQ.1) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QMOUT1'
          CALL qmout1
     &( tstot, tsder, f     , xk    , varian, fmoy, xkmoy,
     &  nf    , ndire , npoin2, taux1 )
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QMOUT1'
!
        ELSEIF(smout.EQ.2) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QMOUT2'
          CALL qmout2
     &( tstot , tsder , f     , xk    , varian, fmoy  , xkmoy , usold ,
     &  usnew , nf    , ndire , npoin2, taux1 , taux2 , taux5 , taux6 )
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QMOUT2'
!
        ENDIF
!
!       4.5 BOTTOM FRICTION DISSIPATION
!       -------------------------------
!
        IF(sfrot.EQ.1.AND..NOT.proinf) THEN
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QFROT1'
          CALL qfrot1
     &( tstot , tsder , f     , xk    , nf    ,ndire , npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QFROT1'
        ELSEIF(sfrot.NE.0) THEN
          WRITE(lu,*) 'OPTION FOR BOTTOM FRICTION DISSIPATION'
          WRITE(lu,*) 'UNKNOWN: SFROT=',sfrot
          CALL plante(1)
          stop
        ENDIF
!
!       5. UPDATES THE SPECTRUM - TAKES THE SOURCE TERMS INTO ACCOUNT
!         (GENERATION, WHITECAPPING AND QUADRUPLET INTERACTIONS)
!       ==============================================================
!
!
!       COMPUTATION OF LIMITING FACTOR INSERTED
!       IN THE LOOP TO HAVE DF_LIM(NPOIN2) INSTEAD
!       OF DF_LIM(NPOIN2,NF)
!
!     if nf is growing, inverse if limit and loop on nf
        DO iff=1,nf
!         LIMITING FACTOR TAKEN FROM WAM-CYCLE 4
          IF(limit.EQ.1) THEN
            coef=0.62d-4*dtsi/1200.d0
            auxi=coef/freq(iff)**5
            DO ip=1,npoin2
              df_lim(ip)=auxi
            ENDDO
!         LIMITING FACTOR FROM HERSBACH AND JANSSEN (1999)
          ELSEIF(limit.EQ.2) THEN
            coef=3.d-7*gravit*freq(nf)*dtsi
            auxi=coef/freq(iff)**4
            usmin=gravit*5.6d-3/freq(iff)
            DO ip=1,npoin2
              df_lim(ip)=auxi*max(usnew(ip),usmin)
            ENDDO
!
!           LIMITING FACTOR FROM LAUGEL ou BAJ
!
          ELSEIF (limit.EQ.3) THEN
            coef=3.0d-7*gravit*dtsi
            auxi=coef/freq(iff)**4
            usmin=gravit*5.6d-3/freq(iff)
            DO ip=1,npoin2
              df_lim(ip)=auxi*max(usnew(ip),usmin)*taux7(ip)
            ENDDO
!
          ELSEIF(limit.NE.0) THEN
            WRITE(lu,*) 'UNKNOWN LIMITING FACTOR:',limit
            CALL plante(1)
            stop
          ENDIF
          IF(limit.NE.0) THEN
            DO jp=1,ndire
              DO ip=1,npoin2
                aux1=max(1.d0-dtsi*tsder(ip,jp,iff)*cimpli,1.d0)
                aux2=dtsi*tstot(ip,jp,iff)/aux1
                aux3=min(abs(aux2),df_lim(ip))
                aux4=sign(aux3,aux2)
                f(ip,jp,iff)=max(f(ip,jp,iff)+aux4,0.d0)
              ENDDO
            ENDDO
          ELSE
            DO jp=1,ndire
              DO ip=1,npoin2
                aux1=max(1.d0-dtsi*tsder(ip,jp,iff)*cimpli,1.d0)
                aux2=dtsi*tstot(ip,jp,iff)/aux1
                f(ip,jp,iff)=max(f(ip,jp,iff)+aux2,0.d0)
              ENDDO
            ENDDO
          ENDIF
        ENDDO
!
        IF(diaghf.EQ.1) THEN
!
!       6. TREATS THE HIGH FREQUENCIES DIFFERENTLY
!       =======================================================
!
!       6.1 COMPUTES THE MEAN FREQUENCY OF THE SPECTRUM
!       ----------------------------------------------
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREMOY'
          CALL fremoy(fmoy, f, nf, ndire, npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREMOY'
!
          aux1=gravit/(7.d0*deupi*freq(1))
          aux2=2.5d0/freq(1)
          aux3=1.d0/log10(raisf)
!
!     IF CBAJ and loop inversed 6.2 and 6.3 written twice with a different formula for FM2
          IF(cbaj.EQ.0) THEN
            DO ip=1,npoin2
!
!       6.2 COMPUTES THE LAST FREQUENCY OF THE DISCRETISED SPECTRUM.
!           THIS FREQUENCY IS THE MAXIMUM OF (FM1=4.*FPM ; FM2=2.5*FMOY).
!           ITS INDEX IS MFMAX.
!       -------------------------------------------------------------
!
              fm1 =aux1/max(usnew(ip),1.d-90)
              fm2 =aux2*fmoy(ip)
              mf1=int(aux3*log10(fm1)+1.d0)
              mf2=int(aux3*log10(fm2)+1.d0)
              mfmax=max(min(max(mf1,mf2),nf),1)
!
!       6.3 MODIFIES THE HIGH FREQUENCY PART OF THE SPECTRUM
!           A DECREASE IN F**(-TAILF) IS IMPOSED BEYOND
!           FREQ(MFMAX).  (TAILF=5 IN WAM-CYCLE 4)
!       -------------------------------------------------------------
!
              DO iff=mfmax+1,nf
                aux4=(freq(mfmax)/freq(iff))**tailf
                DO jp=1,ndire
                  f(ip,jp,iff)=aux4*f(ip,jp,mfmax)
                ENDDO
              ENDDO
            ENDDO

          ELSE
            DO ip=1,npoin2
              fm1 =aux1/max(usnew(ip),1.d-90)
              fm2 =aux2*taux7(ip)
              mf1=int(aux3*log10(fm1)+1.d0)
              mf2=int(aux3*log10(fm2)+1.d0)
              mfmax=max(min(max(mf1,mf2),nf),1)
              DO iff=mfmax+1,nf
                aux4=(freq(mfmax)/freq(iff))**tailf
                DO jp=1,ndire
                  f(ip,jp,iff)=aux4*f(ip,jp,mfmax)
                ENDDO
              ENDDO
!
            ENDDO
          ENDIF
        ELSEIF(diaghf.GE.2) THEN
          WRITE(lu,*) 'OPTION FOR DIAGNOSTIC TAIL'
          WRITE(lu,*) 'UNKNOWN: DIAGHF=',diaghf
          CALL plante(1)
          stop
        ENDIF
!
!
!       7. TAKES THE BREAKING SOURCE TERM INTO ACCOUNT
!       =================================================
!
!        IF((SBREK.GT.0.OR.STRIA.GT.0).AND..NOT.PROINF) THEN
!VB mofid
        IF(((sbrek.GT.0.OR.stria.GT.0.OR.vegetation.OR.porous).AND.
     &     .NOT.proinf).OR.sdscu.EQ.2) THEN
!VB fin modif
!
!         7.1 COMPUTES A REPRESENTATIVE FREQUENCY
!         ------------------------------------------
          IF (sbrek.GT.0.AND.sbrek.LT.5) THEN
            IF (sbrek.EQ.1) ifcar = ifrbj
            IF (sbrek.EQ.2) ifcar = ifrtg
            IF (sbrek.GE.3) ifcar = ifrro
            IF (sbrek.GE.4) ifcar = ifrih
!
            IF (ifcar.EQ.1) THEN
!
!             MEAN FREQUENCY FMOY
!             - - - - - - - - - - - -
              IF (cbaj.EQ.1) THEN
                IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREMOY'
                CALL fremoy(taux3, f, nf, ndire, npoin2 )
                IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREMOY'
              ELSE
                DO ip=1,npoin2
                  taux3(ip)=fmoy(ip)
                ENDDO
              ENDIF

            ELSE IF (ifcar.EQ.2) THEN
!
!             MEAN FREQUENCY F01
!             - - - - - - - - - - -
              IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREM01'
              CALL frem01( taux3, f, nf, ndire, npoin2)
              IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREM01'

            ELSE IF (ifcar.EQ.3) THEN
!
!             MEAN FREQUENCY F02
!             - - - - - - - - - - -
              IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREM02'
              CALL frem02( taux3, f, nf, ndire, npoin2)
              IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREM02'

            ELSE IF (ifcar.EQ.4) THEN
!
!             PEAK FREQUENCY (DISCRETE FREQUENCY WITH MAX VARIANCE)
!             - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
              IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREPIC'
              CALL frepic( taux3, f, nf, ndire, npoin2)
              IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREPIC'

            ELSE IF (ifcar.EQ.5) THEN
!
!             PEAK FREQUENCY (READ WITH EXPONENT 5)
!             - - - - - - - - - - - - - - - - - - - - - - - - - -
              IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FPREAD'
              CALL fpread( taux3, f, nf, ndire, npoin2, 5.d0)
              IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FPREAD'

            ELSE IF (ifcar.EQ.6) THEN
!
!             PEAK FREQUENCY (READ WITH EXPONENT 8)
!             - - - - - - - - - - - - - - - - - - - - - - - - - -
              IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FPREAD'
              CALL fpread( taux3, f, nf, ndire, npoin2, 8.d0)
              IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FPREAD'

            ELSE

              WRITE(lu,*) 'WAVE FREQUENCY NOT EXPECTED......IFCAR=',
     &                       ifcar
              CALL plante(1)
              stop
            ENDIF
!
        ENDIF
!
!.........LOOP ON SUB-TIME STEPS FOR BREAKING
!         = = = = = = = = = = = = = = = = = = = = = = = = = = =
        somme=(xdtbrk**ndtbrk-1.d0)/(xdtbrk-1.d0)
        dtn=dtsi/somme
!
        DO idt=1,ndtbrk
!         7.2 INITIALISES THE ARRAYS FOR THE SOURCE-TERMS
!         ----------------------------------------------------
          DO iff=1,nf
            DO jp=1,ndire
              DO ip=1,npoin2
                tstot(ip,jp,iff)=0.d0
              ENDDO
            ENDDO
          ENDDO
!
!         7.3 COMPUTES THE TOTAL VARIANCE OF THE SPECTRUM
!         --------------------------------------------
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE TOTNRJ'
          CALL totnrj(varian, f, nf, ndire, npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE TOTNRJ'
!
!
!         7.4 COMPUTES THE WAVE BREAKING CONTRIBUTION
!         --------------------------------------
!
!         7.4.1 BREAKING ACCORDING TO BATTJES AND JANSSEN (1978)
!         - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
          IF(sbrek.EQ.1) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QBREK1'
          CALL qbrek1
     & ( tstot , f     , taux3 , varian, nf    , ndire , npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QBREK1'
!
!
!         7.4.2 BREAKING ACCORDING TO THORNTON AND GUZA (1983)
!         - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
          ELSEIF(sbrek.EQ.2) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QBREK2'
          CALL qbrek2
     & ( tstot , f     , taux3 , varian, nf    , ndire ,
     &   npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QBREK2'
!
!
!         7.4.3 BREAKING ACCORDING TO ROELVINK (1993)
!         - - - - - - - - - - - - - - - - - - - - - -
!
          ELSEIF(sbrek.EQ.3) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QBREK3'
          CALL qbrek3
     &( tstot , f     , taux3 , varian, nf    , ndire ,
     &  npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QBREK3'
!
!
!         7.4.4 BREAKING ACCORDING TO IZUMIYA AND HORIKAWA (1984)
!         - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
          ELSEIF(sbrek.EQ.4) THEN
!
          IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QBREK4'
          CALL qbrek4
     &( tstot , f     ,taux3,varian, nf    , ndire , npoin2)
          IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QBREK4'
!
          ELSEIF(sbrek.NE.0) THEN
            WRITE(lu,*) 'BREAKING FORMULATION NOT PROGRAMMED: ',
     &                     sbrek
            CALL plante(1)
            stop
          ENDIF
!
!       7.5 NON-LINEAR INTERACTIONS BETWEEN FREQUENCY TRIPLETS
!       -----------------------------------------------------------
          IF(stria.EQ.1) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE FREMOY'
            CALL fremoy( fmoy, f, nf, ndire, npoin2)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE FREMOY'
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QTRIA1'
            CALL qtria1
     &( f     , xk    , nf    , ndire , npoin2, tstot , varian, fmoy  )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QTRIA1'
!
          ELSEIF(stria.EQ.2) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QTRIA2'
            CALL qtria2
     &( f     , xk    , nf    , ndire , npoin2, tstot )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QTRIA2'
          ENDIF
!
!
!         7.6 WAVE BLOCKING DISSIPATION
!         -----------------------------
          IF(sdscu.EQ.2) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QDSCUR'
            CALL qdscur
     &( tstot , tsder , f     , cf    , xk    , usold , usnew ,
     &  nf    , ndire , npoin2, taux2 ,t3_01%R,t3_02%R)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QDSCUR'
          ENDIF
!
!======================================================================
!         7.7 VEGETATION
!VBA PRISE EN COMPTE VEGETATION
!======================================================================
!
          IF(vegetation) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QVEG'
            CALL qveg( tstot, tsder, f, varian, fmoy, xkmoy, nf,
     &                   ndire  ,npoin2   )
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QVEG'
          ENDIF
!======================================================================
!         4.6 POROUS MEDIA
!TF PRISE EN COMPTE POROSITE
!======================================================================
!
          IF(porous) THEN
            IF(debug.EQ.2) WRITE(lu,*) '     APPEL DE QPOROS'
            CALL qporos( tstot , tsder , f , cg, lt, xk,
     &                   nf    , ndire  , npoin2, amorp)
            IF(debug.EQ.2) WRITE(lu,*) '     RETOUR DE QPOROS'
          ENDIF
!
!======================================================================
!VBA PRISE EN COMPTE VEGETATION
!======================================================================
!
!
!         7.9 UPDATES THE SPECTRUM - TAKES THE BREAKING SOURCE TERM
!             INTO ACCOUNT (EXPLICIT EULER SCHEME)
!         ---------------------------------------------------------
!
        DO iff=1,nf
          DO jp=1,ndire
            DO ip=1,npoin2
              f(ip,jp,iff)=max(f(ip,jp,iff)+dtn*tstot(ip,jp,iff),0.d0)
            ENDDO
          ENDDO
        ENDDO
!
        dtn=dtn*xdtbrk
!
        ENDDO
!
        ENDIF
!
!
!       8. TRANSFERS DATA FROM NEW TO OLD FOR THE NEXT TIME STEP
!       ==============================================================
        IF(vent) THEN
          DO ip=1,npoin2
            usold(ip)=usnew(ip)
            z0old(ip)=z0new(ip)
            twold(ip)=twnew(ip)
          ENDDO
        ENDIF
!
!======================================================================
!        9 IGNORE BOUNDARY CONDITIONS ON IMPOSED BOUNDARIES
!======================================================================
! ALONG THE IMPOSED BOUNDARIES, REWRITE THE CORRECT SPECTRUM
! IF THE SOURCE TERMS SHOULD NOT BE TAKEN INTO ACCOUNT
!
        IF(.NOT.source_on_bnd) THEN
          DO iptfr=1,nptfr
            IF(lifbor(iptfr).EQ.kent) THEN
              DO iff=1,nf
                DO ip=1,ndire
                  f(nbor(iptfr),ip,iff)=fbor(iptfr,ip,iff)
                ENDDO
              ENDDO
            ENDIF
          ENDDO
        ENDIF
      ENDDO

!
!     -----------------------------------------------------------------
!     END OF THE MAIN LOOP ON THE NUMBER OF TIME STEPS (NSITS)
!     FOR INTEGRATION OF THE SOURCE TERMS, BY PROPAGATION TIME STEP
!     -----------------------------------------------------------------
!
      RETURN
      END