qmout2.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE qmout2
!                       *****************
!
     &( tstot, tsder, f    , xk    , enrj , fmoy , xkmoy , usold,  
     &  usnew, nf   , ndire, npoin2, taux1, f_int, betoto, betotn)
!
!**********************************************************************
! TOMAWAC   V6P3                                  23/06/2011
!**********************************************************************
!
!brief   COMPUTES THE CONTRIBUTION OF THE WHITECAPPING SINK TERM USING
!+          THE  PARAMETRISATION of VAN DER WESTHUYSEN (2007).
!
!reference    VAN DER WESTHUYSEN (2007): ADVANCES IN THE SPECTRAL
!+              MODELLING OF WIND WAVES IN THE NEARSHORE, PHD THESID,
!+              DELFT UNIVERSITY OF TECHNOLOGY
!
!history  E. GAGNAIRE-RENOU (EDF/LNHE)
!+        09/2010
!+        V6P0
!+
!
!history  G.MATTAROLO (EDF - LNHE)
!+        23/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!history  J-M HERVOUET (EDF/LNHE)
!+        23/12/2012
!+        V6P3
!+   A first optimisation + limitation of SINH argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| BETA           |<->| WORK TABLE
!| BETOTN         |<->| WORK TABLE
!| BETOTO         |<->| WORK TABLE
!| ENRJ           |-->| SPECTRUM VARIANCE
!| F              |-->| DIRECTIONAL SPECTRUM
!| FMOY           |-->| MEAN SPECTRAL FRQUENCY FMOY
!| XK             |-->| DISCRETIZED WAVE NUMBER
!| XKMOY          |-->| AVERAGE WAVE NUMBER
!| NF             |-->| NUMBER OF FREQUENCIES
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| TAUX1          |<->| WORK TABLE
!| TSDER          |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
!| TSTOT          |-->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
!| USNEW          |-->| FRICTION VELOCITY AT TIME N+1
!| USOLD          |<->| FRICTION VELOCITY AT TIME N
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
!  APPELS :    - PROGRAMME(S) APPELANT  : SEMIMP
!  ********    - PROGRAMME(S) APPELE(S) :    -
!
!  REMARKS:
!  ********
!
!  - THE CONSTANT CMOUT3 (Cdis,break) UTILISED IN WESTHUYSEN (2007)
!                                    IS EQUAL TO 5.0*10^(-5)
!  - THE CONSTANT CMOUT4 (Br) UTILISED IN WESTHUYSEN (2007)
!                                    IS EQUAL TO 1.75*10^(-3)
!  - THE CONSTANT CMOUT5 (Cdis,non-break) UTILISED IN WESTHUYSEN
!                                    (2007) IS EQUAL TO 3.29
!  - THE CONSTANT CMOUT6 (Delta) UTILISED IN WESTHUYSEN (2007)
!                                    IS EQUAL TO 0
!
      USE declarations_tomawac, ONLY : deupi,gravit, cmout3,cmout4, 
     &           cmout5, cmout6, cimpli, proinf, depth, freq
! Variables in TOMAWAC MODULE
! CMOUT3         WESTHUYSEN WHITE CAPPING DISSIPATION COEFFICIENT
! CMOUT4         WESTHUYSEN SATURATION THRES. FOR THE DISSIPATION
! CMOUT5         WESTHUYSEN WHITE CAPPING DISSIPATION COEFFICIENT
! CMOUT6         WESTHUYSEN WHITE CAPPING WEIGHTING COEFFICIENT
! CIMPLI         IMPLICITATION COEFFICIENT FOR SOURCE TERM INTEG.
!
      USE declarations_special
      USE interface_tomawac, ex_qmout2 => qmout2
      IMPLICIT NONE
!
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)             :: NF,NDIRE,NPOIN2
      DOUBLE PRECISION, INTENT(IN)    :: USNEW(npoin2),USOLD(npoin2)
      DOUBLE PRECISION, INTENT(IN)    :: FMOY(npoin2),XK(npoin2,nf)
      DOUBLE PRECISION, INTENT(IN)    :: ENRJ(npoin2),XKMOY(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: F_INT(npoin2),TAUX1(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: BETOTO(npoin2),BETOTN(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: F(npoin2,ndire,nf)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER JP,IFF,IP
      DOUBLE PRECISION PO,AUX,C1,C2,C3,P0O,P0N,W,SURDEUPIFREQ,B,DTETAR
      DOUBLE PRECISION BETAO,BETAN,CPHAS,CG1,SQBSCMOUT4,BETA,DEUKD,KD
      DOUBLE PRECISION SURCMOUT4
!
!-----------------------------------------------------------------------
!
!     DTETAR = DEUPI/DBLE(NDIRE)
!     F_INT WAS DIVIDED BY DEUPI AFTER IN FORMULAS, DIVISION REMOVED
      dtetar = 1.d0/dble(ndire)
      c1     = - cmout5*deupi**9/gravit**4
      c2     = - cmout5*deupi
      w = 25.d0
      surcmout4 = 1.d0/cmout4
!
      IF(proinf) THEN
!       DEEP WATER CASE, ARRAY DEPENDING ONLY ON THE SPATIAL MESH NODE
        DO ip = 1,npoin2
          taux1(ip) = c1 * enrj(ip)**2 * fmoy(ip)**9
        ENDDO
      ELSE
!       FINITE DEPTH CASE
        DO ip=1,npoin2
          taux1(ip) = c2 * enrj(ip)**2 * fmoy(ip) * xkmoy(ip)**4
        ENDDO
      ENDIF
!
!     LOOP ON THE DISCRETISED FREQUENCIES
!
      DO iff=1,nf
!
        surdeupifreq=1.d0/(deupi*freq(iff))
!
        DO ip=1,npoin2
          f_int(ip)=f(ip,1,iff)
        ENDDO
        DO jp=2,ndire
          DO ip=1,npoin2
            f_int(ip)=f_int(ip)+f(ip,jp,iff)
          ENDDO
        ENDDO
        DO ip=1,npoin2
          f_int(ip)=f_int(ip)*dtetar
        ENDDO
!
        IF(proinf) THEN
!
          DO ip = 1,npoin2
!
            cphas = xk(ip,iff)*surdeupifreq
            p0o=3.d0+tanh(w*(usold(ip)*cphas-0.1d0))
            p0n=3.d0+tanh(w*(usnew(ip)*cphas-0.1d0))
            cg1 = 0.5d0*gravit*surdeupifreq
            b   = cg1*f_int(ip)*xk(ip,iff)**3
            sqbscmout4=sqrt(b*surcmout4)
!           COMPUTES THE BREAK/NON-BREAK TRANSITION
            po = 0.5d0*(1.d0+tanh(10.d0*(sqbscmout4-1.d0)))
!           COMPUTES THE BREAK BETA
            c3=-cmout3*sqrt(gravit*xk(ip,iff))
            betao=c3*sqbscmout4**p0o
            betan=c3*sqbscmout4**p0n
!           COMPUTES THE NON-BREAK BETA
            aux = (freq(iff)/fmoy(ip))**2
            beta=taux1(ip)*aux*(1.d0-cmout6+cmout6*aux)
!           COMPUTES THE TOTAL BETA
            betoto(ip)=beta+po*(betao-beta)
            betotn(ip)=beta+po*(betan-beta)
!
          ENDDO
!
        ELSE
!
          DO ip = 1,npoin2
!
            cphas = xk(ip,iff)*surdeupifreq
            kd=min(xk(ip,iff)*depth(ip),350.d0)
            deukd=kd+kd
            cg1=( 0.5d0+xk(ip,iff)*depth(ip)/sinh(deukd) )/cphas
            b = cg1*f_int(ip)*xk(ip,iff)**3
            sqbscmout4=sqrt(b*surcmout4)
!           COMPUTES THE BREAK BETA
            c3=-cmout3*sqrt(gravit*xk(ip,iff))
            aux=tanh(kd)
            p0o=3.d0+tanh(w*(usold(ip)*cphas-0.1d0))
            p0n=3.d0+tanh(w*(usnew(ip)*cphas-0.1d0))
            betao=c3*sqbscmout4**p0o*aux**((2.d0-p0o)*0.25d0)
            betan=c3*sqbscmout4**p0n*aux**((2.d0-p0n)*0.25d0)
!           COMPUTES THE NON-BREAK BETA
            aux = xk(ip,iff) / xkmoy(ip)
!           COMPUTES THE TOTAL BETA
            beta=taux1(ip)*aux*(1.d0-cmout6+cmout6*aux)
!           COMPUTES THE BREAK/NON-BREAK TRANSITION
            po = 0.5d0*(1.d0+tanh(10.d0*(sqbscmout4-1.d0)))
            betoto(ip)=beta+po*(betao-beta)
            betotn(ip)=beta+po*(betan-beta)
!
          ENDDO
!
        ENDIF
!
!       TAKES THE SOURCE TERM INTO ACCOUNT
!
        DO jp = 1,ndire
          DO ip = 1,npoin2
            tstot(ip,jp,iff)=tstot(ip,jp,iff)
     &      +(betoto(ip)+cimpli*(betotn(ip)-betoto(ip)))*f(ip,jp,iff)
            tsder(ip,jp,iff)=tsder(ip,jp,iff) + betotn(ip)
          ENDDO
        ENDDO
!
      ENDDO
!
!-----------------------------------------------------------------------
!
      RETURN
      END