stress.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE stress
!                   *****************
!
     &( tauwav, tstot , f     , usnew , twnew , z0new ,
     &  npoin2, ndire , nf    , xtauw , ytauw , tauhf )
!
!***********************************************************************
! TOMAWAC   V6P1                                   28/06/2011
!***********************************************************************
!
!brief    COMPUTES THE WAVE STRESSES FOR ALL THE NODES
!+                IN THE MESH.
!
!note     THIS SUBROUTINE COMPUTES TAUHF DIRECTLY FROM USTAR AND
!+          ALFA. (THIS WAS DONE VIA A CALL TO 'TAUWHF' PREVIOUSLY).
!
!history  M. BENOIT (EDF/DER/LNH)
!+        03/05/95
!+        V1P0
!+
!
!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  G.MATTAROLO (EDF - LNHE)
!+        28/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| COSTET         |-->| COSINE OF TETA ANGLE
!| DFREQ          |-->| FREQUENCY STEPS BETWEEN DISCRETIZED FREQUENCIES
!| F              |-->| DIRECTIONAL SPECTRUM
!| FREQ           |-->| DISCRETIZED FREQUENCIES
!| NF             |-->| NUMBER OF FREQUENCIES
!| NDIRE          |-->| NUMBER OF DIRECTIONS
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| SINTET         |-->| SINE OF TETA ANGLE
!| TAUHF          |<->| WORK TABLE
!| TAUWAV         |<--| STRESS DUE TO THE WAVES
!| TETA           |-->| DISCRETIZED DIRECTIONS
!| TSTOT          |-->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
!| TWNEW          |-->| WIND DIRECTION
!| USNEW          |-->| FRICTION VELOCITY
!| XKAPPA         |-->| VON KARMAN CONSTANT
!| XTAUW          |<->| WORK TABLE
!| YTAUW          |<->| WORK TABLE
!| Z0NEW          |-->| SURFACE ROUGHNESS LENGTH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_tomawac, ONLY : deupi, gravit,  roair,  roeau,
     &                           freq, dfreq,   teta, sintet, costet,
     &                          decal,xkappa,  betam
!
      USE interface_tomawac, ex_stress => stress
      IMPLICIT NONE
!
!.....VARIABLES IN ARGUMENT
!     """"""""""""""""""""
      INTEGER, INTENT(IN)    ::  NPOIN2, NDIRE , NF
      DOUBLE PRECISION, INTENT(IN)    :: USNEW(npoin2), TWNEW(npoin2)
      DOUBLE PRECISION, INTENT(IN)    :: Z0NEW(npoin2)
      DOUBLE PRECISION, INTENT(IN)    :: TSTOT(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(IN)    :: F(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: TAUWAV(npoin2), TAUHF(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: XTAUW(npoin2) , YTAUW(npoin2)
!.....VARIABLES FROM MODULE TOMAWAC
!     """"""""""""""""""""
!     DECAL           SHIFT GROWING CURVE DUE TO WIND
!     XKAPPA          VON KARMAN CONSTANT
!     BETAM           WIND GENERATION COEFFICIENT
!
!.....LOCAL VARIABLES
!     """""""""""""""""
      INTEGER  IP    , JP    , JF    , JTOT  , J
      DOUBLE PRECISION DTETAR, FRMAX , COEF1 , COEF2 , TTAUHF
      DOUBLE PRECISION USTAR , ALFA  , COSTMP, C1    , C2    , DIREC
      DOUBLE PRECISION Y     , OMEGA , CM    , ZX    , ZARG  , ZMU
      DOUBLE PRECISION ZLOG  , ZBETA , UST   , Z0    , UST2  , ALF2
      DOUBLE PRECISION CONST1, OMEGAM, X0    , YC    , DELY  , AUX
!
!
      dtetar= deupi/dble(ndire)
      frmax = freq(nf)
      coef1 = dtetar*deupi**4*frmax**5/gravit**2
      coef2 = deupi*roeau/roair*dtetar
!
      DO ip=1,npoin2
        xtauw(ip)=0.d0
        ytauw(ip)=0.d0
      ENDDO ! IP
!
!.....INTEGRATES THE SOURCE TERM OVER FREQUENCIES AND DIRECTIONS
!     """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
      DO jf=1,nf
        aux=coef2*freq(jf)*dfreq(jf)
        DO jp=1,ndire
          c1=aux*sintet(jp)
          c2=aux*costet(jp)
          DO ip=1,npoin2
            xtauw(ip)=xtauw(ip)+tstot(ip,jp,jf)*c1
            ytauw(ip)=ytauw(ip)+tstot(ip,jp,jf)*c2
          ENDDO ! IP
        ENDDO ! JP
      ENDDO ! JF
!
!.....COMPUTES THE PARAMETERISED HIGH FREQUENCY PART
!     """""""""""""""""""""""""""""""""""""""""""""""""""
      DO ip=1,npoin2
        tauhf(ip)=0.d0
      ENDDO ! IP
!
      DO jp=1,ndire
        direc=teta(jp)
        DO ip=1,npoin2
          costmp=max(cos(direc-twnew(ip)),0.d0)
          tauhf(ip)=tauhf(ip)+f(ip,jp,nf)*costmp**3
        ENDDO ! IP
      ENDDO ! JP
!
      jtot  = 50
      const1= betam/xkappa**2
      omegam= deupi*frmax
      x0    = 0.05d0
!
      DO ip=1,npoin2
        ustar=usnew(ip)
        alfa =z0new(ip)*gravit/ustar**2
!
!----------------------------------------------OLD SUBROUTINE TAUWHF
!C      CALL TAUWHF
!C   *( TAUHF , USTAR , ALFA  , BETAM , XKAPPA , DECAL , FRMAX , GRAVIT)
!----------------------------------------------
!
!MB.....LIMITING FACTORS TO REPRODUCE WAM4 (SUBROUTINE TAUHF OF PREPROC)
!.......(THIS LIMITATION IS NOT JUSTIFIED A PRIORI. IT IS AN ARTEFACT
!.......OF TAUHF BEING DISCRETISED ON A GRID)
        ust2  = min(ustar,5.d0)
        alf2  = min(alfa,0.11d0)
!MB.....END OF LIMITATION
!
        ust   = max(ust2,0.000001d0)
        z0    = alf2*ust**2/gravit
!
        yc    = max(omegam,x0*gravit/ust)*sqrt(z0/gravit)
        dely  = max((1.d0-yc)/float(jtot),0.d0)
        ttauhf = 0.d0
        DO j=1,jtot
          y     = yc+dble(j-1)*dely
          omega = y*sqrt(gravit/z0)
          cm    = gravit/omega
          zx    = ust/cm +decal
          zarg  = min(xkappa/zx,20.d0)
          zmu   = min(gravit*z0/cm**2*exp(zarg),1.d0)
          zlog  = min(log(zmu),0.d0)
          zbeta = const1*zmu*zlog**4
          ttauhf= ttauhf+zbeta/y*dely
        ENDDO ! J
!----------------------------------------------ANCIENNE SUB TAUWHF
!
        tauhf(ip) = ttauhf*coef1*ustar**2*tauhf(ip)
      ENDDO ! IP
!
!.....TAKES THE PARAMETERISED HIGH FREQUENCY PART INTO ACCOUNT
!     """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
!
      DO ip=1,npoin2
        xtauw(ip) = xtauw(ip) + tauhf(ip)*sin(twnew(ip))
        ytauw(ip) = ytauw(ip) + tauhf(ip)*cos(twnew(ip))
        tauwav(ip)= sqrt(xtauw(ip)**2+ytauw(ip)**2)
      ENDDO ! IP
!
      RETURN
      END