windiss1.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE windiss1
!                   *****************
     &  (fwx, fwy, npoin2, xk, ndire, fs,nf)
!  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,
     &   teta, betam, decal, xkappa, deupi, gravit, roair, roeau,
     &   twold, usold, z0old
!
      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 DTETAR, ZLOGMU, C1, CONST, SIGMA, AUX1, CPHAS
      DOUBLE PRECISION XX, USO, OMEGA, BETAWIN
      INTEGER IP, JF, JP
      dtetar=deupi/dble(ndire)
      DO ip=1,npoin2
        fwx(ip) = 0.d0
        fwy(ip) = 0.d0
      ENDDO
      !
      c1 = deupi * (roair/roeau) * (betam/xkappa**2)
!.....LOOP ON THE DISCRETISED FREQUENCIES
!     """"""""""""""""""""""""""""""""""""""""""""
      DO jf=1,nf

        const=c1*freq(jf)
        sigma=deupi*freq(jf)
        aux1=dfreq(jf)*dtetar
!
!.......LOOP ON THE DISCRETISED DIRECTIONS
!       """"""""""""""""""""""""""""""""""""""""""""
        DO jp=1,ndire
          DO ip=1,npoin2
!.......COMPUTES THE FREQUENCIES (OMEGA AND UETOILE/CPHASE)
            cphas = deupi * freq(jf) / xk(ip,jf)
            omega = gravit * z0old(ip) / cphas**2
            uso  = usold(ip) / cphas + decal
            xx = uso *  cos(teta(jp)-twold(ip))
            zlogmu = log(omega) + xkappa/xx
            IF(zlogmu.LT.0.d0) THEN
              betawin = const*omega*exp(xkappa/xx)*
     &             zlogmu**4*xx**2
              fwx(ip)=fwx(ip)+((xk(ip,jf)/sigma)*sintet(jp)
     &             *betawin*fs(ip,jp,jf))*aux1
              fwy(ip)=fwy(ip)+((xk(ip,jf)/sigma)*costet(jp)
     &             *betawin*fs(ip,jp,jf))*aux1
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END