moudiss2.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE moudiss2
!                   *****************
     &  (fwx, fwy, npoin2, xk, ndire, fs,nf, taux1, f_int)
!  SURFACE STRESS DUE TO WIND INPUT ENERGY AND WHITECAPPING
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FS             |-->| VARIANCE DENSITY DIRECTIONAL SPECTRUM
!| FWX            |<--| SURFACE STRESS DUE TO WIND ALONG X
!| FWY            |<--| SURFACE STRESS DUE TO WIND ALONG Y
!| NF             |-->| NUMBER OF FREQUENCIES
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| XK             |-->| DISCRETIZED WAVE NUMBER
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      USE declarations_tomawac, ONLY : freq, dfreq, sintet, costet,
     & deupi, cmout3,cmout4, cmout5, cmout6, varian, fmoy, xkmoy,
     & depth, usold, proinf, gravit
!
      IMPLICIT NONE
!
      INTEGER, INTENT(IN) :: NPOIN2, NDIRE,NF
      DOUBLE PRECISION, INTENT(IN) :: FS(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(IN) :: XK(npoin2,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: FWX(npoin2), FWY(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: F_INT(npoin2),TAUX1(npoin2)
!
      DOUBLE PRECISION DTETAR, SIGMA, AUX, BETAMOU, AUX1
      DOUBLE PRECISION W, SURDEUPIFREQ, SQBSCMOUT4, SURCMOUT4
      DOUBLE PRECISION PO, P0O, KD, DEUKD
      DOUBLE PRECISION CG1, CPHAS, C3, C2, C1, BETAO, BETA, B
      INTEGER IP, JF, JP
      dtetar=deupi/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 * varian(ip)**2 * fmoy(ip)**9
        ENDDO
      ELSE
!       FINITE DEPTH CASE
        DO ip=1,npoin2
          taux1(ip) = c2 * varian(ip)**2 * fmoy(ip) * xkmoy(ip)**4
        ENDDO
      ENDIF
      DO jf=1,nf
        sigma=deupi*freq(jf)
        surdeupifreq=1.d0/(deupi*freq(jf))
        aux1=dfreq(jf)*dtetar
        DO ip=1,npoin2
          f_int(ip)=fs(ip,1,jf)
        ENDDO
        DO jp=2,ndire
          DO ip=1,npoin2
            f_int(ip)=f_int(ip)+fs(ip,jp,jf)
          ENDDO
        ENDDO
        DO ip=1,npoin2
          f_int(ip)=f_int(ip)*dtetar
        ENDDO
 !
        IF(proinf) THEN
!
          DO ip = 1,npoin2
!
            cphas = xk(ip,jf)*surdeupifreq
            p0o = 3.d0+tanh(w*(usold(ip)*cphas-0.1d0))
            cg1 = 0.5d0*gravit*surdeupifreq
            b   = cg1*f_int(ip)*xk(ip,jf)**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,jf))
            betao = c3*sqbscmout4**p0o
            betamou = beta+po*(betao-beta)

            DO jp=1,ndire
              fwx(ip)=fwx(ip)+((xk(ip,jf)/sigma)*sintet(jp)
     &             *betamou*fs(ip,jp,jf))*aux1
              fwy(ip)=fwy(ip)+((xk(ip,jf)/sigma)*costet(jp)
     &             *betamou*fs(ip,jp,jf))*aux1
            ENDDO
          ENDDO
        ELSE
!       FINITE WATER DEPTH (USES K).
          DO ip=1,npoin2
            cphas = xk(ip,jf)*surdeupifreq
            kd=min(xk(ip,jf)*depth(ip),350.d0)
            deukd=kd+kd
            cg1=( 0.5d0+xk(ip,jf)*depth(ip)/sinh(deukd) )/cphas
            b = cg1*f_int(ip)*xk(ip,jf)**3
            sqbscmout4=sqrt(b*surcmout4)
!           COMPUTES THE BREAK BETA
            c3=-cmout3*sqrt(gravit*xk(ip,jf))
            aux=tanh(kd)
            p0o=3.d0+tanh(w*(usold(ip)*cphas-0.1d0))
            betao=c3*sqbscmout4**p0o*aux**((2.d0-p0o)*0.25d0)
!           COMPUTES THE NON-BREAK BETA
            aux = xk(ip,jf) / 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)))
            betamou=beta+po*(betao-beta)

            DO jp=1,ndire
              fwx(ip)=fwx(ip)+((xk(ip,jf)/sigma)*sintet(jp)
     &             *betamou*fs(ip,jp,jf))*aux1
              fwy(ip)=fwy(ip)+((xk(ip,jf)/sigma)*costet(jp)
     &             *betamou*fs(ip,jp,jf))*aux1
            ENDDO
          ENDDO
        ENDIF
      ENDDO
      RETURN
      END