radia1.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE radia1
!                   *****************
!
     &(lishho)
!
!***********************************************************************
! ARTEMIS   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    COMPUTES THE RADIATION STRESSES AND DRIVING FORCES.
!
!reference  M.W. DINGEMANS, A.C. RADDER AND H.J. DE VRIEND
!+          COMPUTATION OF THE DRIVING FORCES OF WACE-INDUCED
!+          CURRENTS. COASTAL ENGINEERING, 11 (1987) PP 539-563.
!
!history  J-M HERVOUET (LNH)
!+
!+
!+   LINKED TO BIEF 5.0
!
!history  D. AELBRECHT (LNH) ; F. BECQ (LNH)
!+        04/06/1999
!+        V5P1
!+
!
!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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| LISHHO         |<--| SMOOTHING FOR THE WAVE HEIGTH
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_artemis
      USE interface_artemis, ex_radia1 => radia1
!
      USE declarations_special
      IMPLICIT NONE
!
!
      INTEGER I
      INTEGER LISHHO
!
! INTERNAL VARIABLES FOR RADIA1
!
      DOUBLE PRECISION COE , COCO, COSI, SISI
      DOUBLE PRECISION OMEG2
      INTEGER          IRADIA , LISRAD
!
      LOGICAL MAS
!
      INTRINSIC cos, sin
!
!-----------------------------------------------------------------------
!
!     HARD-CODES THE METHOD OF COMPUTATION OF THE RADIATION
!     STRESSES FOLLOWING TESTS : METHOD 2
!
      iradia = 2
!
!-----------------------------------------------------------------------
!
! FOR MEMORY, AND EVEN THOUGH IT IS NOT USED HERE,
! THE FOLLOWING GIVES METHOD 1
!
!=======================================================================
!     RADIATION STRESSES........METHOD 1
!=======================================================================
!
      IF(iradia.EQ.1) THEN
!
      CALL os('X=YZ    ' , x=t3, y=phii, z=phii)
      CALL os('X=X+YZ  ' , x=t3, y=phir, z=phir)
      CALL vector(t2 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , t3 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL vector(t2 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , t2 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL vector(t4 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , t3 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
      CALL vector(t4 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , t4 , t1 , t1 , t1 , t1 , t1 ,
     &            mesh , msk , maskel)
!
      CALL vector(t1 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , t3 , t3 , t3 , t3 , t3 , t3 ,
     &            mesh , msk , maskel)
!
      CALL os('X=Y+Z   ', x=t4, y=t2 , z=t4)
      CALL os('X=Y/Z   ', x=t4, y=t4 , z=t1)
      CALL os('X=Y/Z   ', x=t4, y=t4 , z=t1)
!
      CALL os('X=CY/Z  ', x=t2, y=cg , z=c, c=2.d0  )
      CALL os('X=X+C   ', x=t2, c=-1.d0 )
      omeg2 = omega*omega
      CALL os('X=CX    ' , x=t2 , c=omeg2 )
!
      coe = 1.d0/(8.d0*grav)
!
!=======================================================================
!     RADIATION STRESSES........METHOD 2 (IDENTICAL TO THAT USED IN TOMAWAC)
!=======================================================================
!
      ELSE
!
      CALL os('X=Y     ',x=t3,y=hho)
!
! -------------------------------------------------------------
! SMOOTHES THE WAVE HEIGHT TO ELIMINATE PARASITIC
! OSCILLATIONS
! -------------------------------------------------------------
!
      IF(lishho.GT.0) THEN
        mas = .true.
        CALL filter(t3,mas,t1,t2,am1,'MATMAS          ',
     &              1.d0,sbid,sbid,sbid,sbid,sbid,sbid,
     &              mesh,msk,maskel,lishho)
      ENDIF
!
      CALL os('X=Y     ',x=hho,y=t3)
!
! -------------------------------------------------------------
! COMPUTES STRESSES SXX, SXY AND SYY
! -------------------------------------------------------------
!
      CALL os('X=Y/Z   ' , x=t1 , y=cg , z=c)
      DO i=1,npoin
        coco=cos(inci%R(i))*cos(inci%R(i))
        cosi=cos(inci%R(i))*sin(inci%R(i))
        sisi=sin(inci%R(i))*sin(inci%R(i))
        coe=grav*hho%R(i)*hho%R(i)/8.d0
!
! THE COEFFICIENT 1/8 ABOVE STEMS FROM HHO REPRESENTING THE WAVE
! HEIGHT (ENERGY) IN REGULAR SEAS
!
        sxx%R(i)= sxx%R(i) + (t1%R(i)*(1.d0+coco)-0.5d0)*coe
        sxy%R(i)= sxy%R(i) + (t1%R(i)*cosi)*coe
        syy%R(i)= syy%R(i) + (t1%R(i)*(1.d0+sisi)-0.5d0)*coe
      END DO
      END IF
!
!
!=======================================================================
! SPACIAL GRADIENTS OF RADIATION STRESSES
!=======================================================================
!
!  -----------------------------------------------
!  OPTIONAL SMOOTHING(S) OF THE RADIATION STRESSES
!  -----------------------------------------------
!
      lisrad = 3
!
      CALL os('X=Y     ',x=t3,y=sxx)
      IF(lisrad.GT.0) THEN
        mas = .true.
        CALL filter(t3,mas,t1,t2,am1,'MATMAS          ',
     &              1.d0,sbid,sbid,sbid,sbid,sbid,sbid,
     &              mesh,msk,maskel,lisrad)
      ENDIF
      CALL os('X=Y     ',x=sxx,y=t3)
!
      CALL os('X=Y     ',x=t3,y=sxy)
      IF(lisrad.GT.0) THEN
        mas = .true.
        CALL filter(t3,mas,t1,t2,am1,'MATMAS          ',
     &              1.d0,sbid,sbid,sbid,sbid,sbid,sbid,
     &              mesh,msk,maskel,lisrad)
      ENDIF
      CALL os('X=Y     ',x=sxy,y=t3)
!
      CALL os('X=Y     ',x=t3,y=syy)
      IF(lisrad.GT.0) THEN
        mas = .true.
        CALL filter(t3,mas,t1,t2,am1,'MATMAS          ',
     &              1.d0,sbid,sbid,sbid,sbid,sbid,sbid,
     &              mesh,msk,maskel,lisrad)
      ENDIF
      CALL os('X=Y     ',x=syy,y=t3)

! END OF RADIATION STRESS SMOOTHING(S)
! -------------------------------------------------------
!
!=======================================================================
! DRIVING FORCES FX AND FY FOR WAVE-INDUCED CURRENTS
!=======================================================================
!
      CALL vector(t1 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , t3 , t3 , t3 , t3 , t3 , t3 ,
     &            mesh , msk , maskel )
!
      CALL vector(t2 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , sxx, t4 , t4 , t4 , t4 , t4 ,
     &            mesh , msk , maskel )
      CALL os('X=Y/Z   ',x=t2,y=t2,z=t1)
!
      CALL vector
     & (t3,'=','GRADF          Y',ielm,1.d0,sxy,t4,t4,t4,t4,t4,
     &  mesh , msk , maskel )
      CALL os('X=Y/Z   ',x=t3,y=t3,z=t1)
!     ----------------------------------
!     FORCE FX = - (DSXX/DX + DSXY/DY) / H
!     ----------------------------------
      CALL os('X=Y+Z   ',x=fx,y=t2,z=t3)
      CALL os('X=CY/Z  ',x=fx,y=fx,z=h,c=-1.d0)
!
!     ----------------------------------
!
      CALL vector(t1 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , t3 , t3 , t3 , t3 , t3 , t3 ,
     &            mesh , msk , maskel )
!
      CALL vector
     & (t2,'=','GRADF          X',ielm,1.d0,sxy,t4,t4,t4,t4,t4,
     &  mesh , msk , maskel )
      CALL os('X=Y/Z   ',x=t2,y=t2,z=t1)
!
      CALL vector
     & (t3,'=','GRADF          Y',ielm,1.d0,syy,t4,t4,t4,t4,t4,
     &  mesh , msk , maskel )
      CALL os('X=Y/Z   ',x=t3,y=t3,z=t1)
!
!     ----------------------------------
!     FORCE FY = - (DSXY/DX + DSYY/DY) / H
!     ----------------------------------
      CALL os('X=Y+Z   ',x=fy,y=t2,z=t3)
      CALL os('X=CY/Z  ',x=fy,y=fy,z=h,c=-1.d0)
!
!=======================================================================
!
      RETURN
      END