radiat.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE radiat
!                   *****************
!
     &     ( fx1    , fy1, xk1   ,fs, cg1,
     &  cgsuc1,dsxxdx, dsxydx, dsxydy, dsyydy)
!
!***********************************************************************
! TOMAWAC   V7P0
!***********************************************************************
!
!brief    COMPUTES THE RADIATION STRESSES AND DRIVING FORCES
!+                FOR THE GENERATION OF WAVE-INDUCED CURRENTS.
!+
!+           (SEE NOTES FOR METHODOLOGY)
!code
!+  THE RESULT OF THIS COMPUTATION IS GIVEN AS :
!+       FI = - 1/D D( SIJ )/D( XJ )    UNIT : M/S**2
!
!note     COMPUTATION ACCORDING TO THE "THEORICAL" FORMULATION, WITH
!+          COMPUTATION OF THE TERMS IN THE TENSOR OF THE RADIATION
!+          STRESSES, AND THEN THEIR GRADIENTS IN SPACE.
!
!history  M. BENOIT (EDF/DER/LNH)
!+        13/12/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)
!+        27/06/2011
!+        V6P1
!+   Translation of French names of the variables in argument
!
!history  J-M HERVOUET (EDF - LNHE)
!+        08/01/2014
!+        V7P0
!+   CALL PARCOM suppressed by using new argument ASSPAR in VECTOR
!
!history  C. VILLARET (HR-WALLINGFORD)
!+        15/09/2014
!+        V7P0
!+   Cancellation of radiation stresses below a given depth hmin.
!
!history T. Fouquet (made by C Raoul)
!        01/02/2018
!        V7P3
!        Forces calculated in spherical coordinates
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| CG1            |-->| DISCRETIZED GROUP VELOCITY
!| CGSUC1         |<--| WORK TABLE
!| DEPTH1         |-->| WATER DEPTH
!| DSXXDX         |<->| WORK TABLE
!| DSXYDX         |<->| WORK TABLE
!| DSXYDY         |<->| WORK TABLE
!| DSYYDY         |<->| WORK TABLE
!| FS             |-->| DIRECTIONAL SPECTRUM
!| FX             |<--| DRIVING FORCE ALONG X
!| FY             |<--| DRIVING FORCE ALONG Y
!| NPOIN2         |-->| NUMBER OF POINTS IN 2D MESH
!| SXX            |<--| RADIATION STRESS ALONG XX
!| SXY            |<--| RADIATION STRESS ALONG XY
!| SYY            |<--| RADIATION STRESS ALONG YY
!| XK1            |-->| DISCRETIZED WAVE NUMBER
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_tomawac
!
      USE interface_tomawac, ex_radiat => radiat
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      DOUBLE PRECISION, INTENT(IN)    :: FS(npoin2,ndire,nf)
      DOUBLE PRECISION, INTENT(IN)    :: CG1(npoin2,nf),XK1(npoin2,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: CGSUC1(npoin2,nf)
      DOUBLE PRECISION, INTENT(INOUT) :: FX1(npoin2),FY1(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: DSXXDX(npoin2),DSXYDX(npoin2)
      DOUBLE PRECISION, INTENT(INOUT) :: DSXYDY(npoin2),DSYYDY(npoin2)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER JP,JF,IP
      DOUBLE PRECISION COEF,COEF2,COCO,SISI,SICO,OMEGA,DTETAR
!
!     CV: MINIMUM WATER DEPTH FOR TIDAL FLATS TREATMENT
!     NEW KEYWORD IS NEEDED
!
      DOUBLE PRECISION HMIN
!
!-----------------------------------------------------------------------
!
      hmin=0.1d0
!
      dtetar=deupi/ndire
!
      DO ip=1,npoin2
        sxx(ip) = 0.d0
        sxy(ip) = 0.d0
        syy(ip) = 0.d0
      ENDDO
!
!     COMPUTES THE WORKING ARRAY N = CG/C
!
      DO jf = 1,nf
        omega=deupi*freq(jf)
        DO ip=1,npoin2
          cgsuc1(ip,jf)=cg1(ip,jf)*xk1(ip,jf)/omega
        ENDDO
      ENDDO
!
!     COMPUTES THE RADIATION STRESSES INTEGRATED OVER THE SPECTRUM
!     SUMS UP THE DISCRETISED PART OF THE SPECTRUM
!
      DO jp=1,ndire
        coco=costet(jp)**2
        sisi=sintet(jp)**2
        sico=sintet(jp)*costet(jp)
        DO jf=1,nf
          coef=gravit*dfreq(jf)*dtetar
          DO ip=1,npoin2
            coef2=coef*fs(ip,jp,jf)
            sxx(ip)=sxx(ip)+(cgsuc1(ip,jf)*(1.d0+sisi)-0.5d0)*coef2
            sxy(ip)=sxy(ip)+(cgsuc1(ip,jf)*sico             )*coef2
            syy(ip)=syy(ip)+(cgsuc1(ip,jf)*(1.d0+coco)-0.5d0)*coef2
          ENDDO
        ENDDO
      ENDDO
!
!     COMPUTES THE GRADIENTS IN SPACE OF THE RADIATION STRESSES
!
!
!     INVERSE OF INTEGRALS OF TEST FUNCTIONS
!
      CALL vector(st0,'=','MASBAS          ',ielm2,1.d0,
     &            st1,st1,st1,st1,st1,st1,mesh,.false.,st1,
     &            asspar=.true.)
      CALL os('X=1/Y   ',x=st0,y=st0)
!
!     DERIVATIVE IN X
!
      CALL ov('X=Y     ',x=t4,y=sxx,dim1=npoin2)

      CALL vector
     & (st1,'=','GRADF          X',ielm2,1.d0,st4,
     &  st3,st3,st3,st3,st3,mesh,.false.,st3,asspar=.true.)
!
      CALL ov('X=Y     ',x=t4,y=sxy,dim1=npoin2)
      CALL vector
     & (st2,'=','GRADF          X',ielm2,1.d0,st4,
     &  st3,st3,st3,st3,st3,mesh,.false.,st3,asspar=.true.)
!
      CALL ov('X=YZ    ',x=dsxxdx, y=t1, z=t0, dim1=npoin2)
      CALL ov('X=YZ    ',x=dsxydx, y=t2, z=t0, dim1=npoin2)
!
!     DERIVATIVE IN Y
!
      CALL ov('X=Y     ',x=t4, y=syy, dim1=npoin2)
      CALL vector
     & (st1,'=','GRADF          Y',ielm2,1.d0,st4,
     &  st3,st3,st3,st3,st3,mesh,.false.,st3,asspar=.true.)
!
      CALL ov('X=Y     ',x=t4, y=sxy, dim1=npoin2)
      CALL vector
     & (st2,'=','GRADF          Y',ielm2,1.d0,st4,
     &  st3,st3,st3,st3,st3,mesh,.false.,st3,asspar=.true.)
!
      CALL ov('X=YZ    ',x=dsyydy, y=t1, z=t0, dim1=npoin2)
      CALL ov('X=YZ    ',x=dsxydy, y=t2, z=t0, dim1=npoin2)
!
!     COMPUTES THE DRIVING FORCES FOR WAVE-INDUCED CURRENTS

      IF (.NOT.sphe) THEN
!                               +-----------------------------+
!.............................. ! CARTESIAN COORDINATE SYSTEM !
!                               +-----------------------------+
        DO ip=1,npoin2
          IF(depth(ip).GE.hmin) THEN
            fx1(ip)= - (dsxxdx(ip)+dsxydy(ip))/depth(ip)
            fy1(ip)= - (dsxydx(ip)+dsyydy(ip))/depth(ip)
          ELSE
            fx1(ip)=0.d0
            fy1(ip)=0.d0
          ENDIF
        ENDDO
!
      ELSE
!                               +-----------------------------+
!.............................. ! SPHERICAL COORDINATE SYSTEM !
!                               +-----------------------------+
        DO ip=1,npoin2
          tra31(ip)=1.0d0/cosf(ip)
        ENDDO
        CALL ov('X=CYZ    ',x=dsxxdx, y=dsxxdx, z=tra31, c=sr,
     &          dim1=npoin2)
        CALL ov('X=CYZ    ',x=dsxydx, y=dsxydx, z=tra31, c=sr,
     &          dim1=npoin2)
        CALL ov('X=CY     ',x=dsyydy, y=dsyydy, c=sr, dim1=npoin2)
        CALL ov('X=CY     ',x=dsxydy, y=dsxydy, c=sr, dim1=npoin2)

        DO ip=1,npoin2
          IF(depth(ip).GE.hmin) THEN
            fx1(ip)= - (dsxxdx(ip)+dsxydy(ip))*gradeg/depth(ip)
            fy1(ip)= - (dsxydx(ip)+dsyydy(ip))*gradeg/depth(ip)
          ELSE
            fx1(ip)=0.d0
            fy1(ip)=0.d0
          ENDIF
        ENDDO
      ENDIF


!

!
!-----------------------------------------------------------------------
!
      RETURN
      END