caldir.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE caldir
!                   *****************
!
!
!***********************************************************************
! ARTEMIS   V7P3                                     Aug 2017
!***********************************************************************
!
!brief    COMPUTES THE WAVE INCIDENCE
!
!history  C PEYRARD (LNHE)
!+
!+
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  N.DURAND (HRW)
!+        August 2017
!+        V7P3
!+   PI, DEUPI and PISUR2 now defined in DECLARATIONS_ARTEMIS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_telemac
      USE declarations_artemis
!
      USE declarations_special
      IMPLICIT NONE
!
      INTEGER I
!
      DOUBLE PRECISION A1, A2 ,ALPHA0, D1, D2, PHI, PHI1, PHI2 ,MODPHI
      DOUBLE PRECISION TETA01, XU1 ,XU2, XV1, XV2, WT0
!
      INTRINSIC sqrt, atan2, mod, abs, cos, sin, atan
!
!-----------------------------------------------------------------------
!
      DOUBLE PRECISION, PARAMETER :: ZERO = 1.d-10
!
!=======================================================================
! SPEEDS AT THE SURFACE (AT T=0 AND T=OMEGA/4)
!=======================================================================
!
! COMPUTES THE GRADIENTS (PHIR AND PHII)
!
!
      CALL vector(u0 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , phir , sbid, sbid , sbid , sbid , sbid ,
     &            mesh , msk , maskel )
!
      CALL vector(v0 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , phir , sbid , sbid , sbid , sbid , sbid ,
     &            mesh , msk , maskel )
!
!     THE OLD VARIABLE U1 IS STORED IN T3
!     BECAUSE IT IS USED TO COMPUTE INCI
!
      CALL vector(t3 , '=' , 'GRADF          X' , ielm ,
     &            1.d0 , phii , sbid , sbid , sbid , sbid , sbid ,
     &            mesh , msk , maskel )
!
!     THE OLD VARIABLE V1 IS STORED IN T4
!     BECAUSE IT IS USED TO COMPUTE INCI
!
      CALL vector(t4 , '=' , 'GRADF          Y' , ielm ,
     &            1.d0 , phii , sbid , sbid , sbid , sbid , sbid ,
     &            mesh , msk , maskel )
!
      CALL vector(t1 , '=' , 'MASBAS          ' , ielm ,
     &            1.d0 , sbid , sbid , sbid , sbid , sbid , sbid ,
     &            mesh , msk , maskel )
!
      CALL os('X=Y/Z   ', x=u0, y=u0, z=t1)
      CALL os('X=Y/Z   ', x=v0, y=v0, z=t1)
      CALL os('X=Y/Z   ', x=t3, y=t3, z=t1)
      CALL os('X=Y/Z   ', x=t4, y=t4, z=t1)
!
!=======================================================================
! COMPUTES WAVE INCIDENCE
!=======================================================================
!
!        U0 (D(PHIR)/DX) : A      U1 (D(PHII)/DX): B
!        V0 (D(PHIR)/DY) : C      V1 (D(PHII)/DY): D
! FROM U= A COS WT + B SIN WT  TO : U = A1 COS ( WT - PHI1)
!      V= C COS WT + D SIN WT       V = A2 COS ( WT - PHI2)
!
      DO i=1,npoin
        a1 = sqrt( u0%R(i)*u0%R(i) + t3%R(i)*t3%R(i) )
        phi1 = atan2( t3%R(i),u0%R(i) )
        a2 = sqrt( v0%R(i)*v0%R(i) + t4%R(i)*t4%R(i) )
        IF(t4%R(i).EQ.0.d0.AND.v0%R(i).EQ.0.d0) THEN
          phi2 = 2.d0*atan(1.d0)
        ELSE
          phi2 = atan2( t4%R(i),v0%R(i) )
        ENDIF
!
! WRITTEN AS : U = A1 COS ( (WT - PHI1))
!              V = A2 COS ( (WT - PHI1) - PHI )
! WHERE PHI IS BETWEEN 0 AND 2*PI
!
        phi = phi2 - phi1
        IF (phi.LT.0.d0)   phi = phi+deupi
!
! ESTIMATES THE DIRECTION AND (WT0) WHEN THE ELLIPSE'S MAJOR AXIS
! IS REACHED.
! TREATS INDIVIDUAL CASES (LINEAR POLARISATION)
!
        modphi = mod( phi, pi )
        IF ( (modphi.LT.zero).OR.((pi-modphi).LT.zero) ) THEN
          wt0 = phi1
          IF ( (phi.LT.2d0*zero).OR.((deupi-phi).LT.2d0*zero) )THEN
            alpha0 = atan2( a2,a1 )
          ELSE
!                  (ABS(PHI-PI).LT.2D0*ZERO)
            alpha0 = deupi - atan2( a2,a1 )
          ENDIF
        ELSE
! GENERAL CASE: ELLIPTIC POLARISATION
!        TAN(2*(WT0 - PHI1)) = A2**2*SIN(2*PHI)/(A1**2+A2**2*COS(2*PHI))
          teta01 = atan2( (a2*a2*sin(2*phi)) ,
     &                    (a1*a1 + a2*a2*cos(2*phi)) ) / 2.d0
          xu1 = a1 * cos( teta01)
          xv1 = a2 * cos( teta01 - phi )
          xu2 = -a1 * sin( teta01)
          xv2 = -a2 * sin( teta01 - phi )
          d1 = xu1*xu1 + xv1*xv1
          d2 = xu2*xu2 + xv2*xv2
          IF (d2.GT.d1) THEN
            teta01 = teta01 + pisur2
            xu1    = xu2
            xv1    = xv2
          ENDIF
          wt0    = teta01 + phi1
          alpha0 = atan2( xv1,xu1 )
        ENDIF
        inci%R(i)  = alpha0
        t2%R(i) = wt0
      ENDDO
!
! FREE SURFACE IN PHASE WITH ALPHA0
! INCIDENCE IS CONSIDERED POSITIVE WHEN THE FREE SURFACE IS
! POSITIVE.
!
      DO i=1,npoin
        a1 = -(phii%R(i)*cos(t2%R(i))-phir%R(i)*sin(t2%R(i)))
        IF (a1.LT.0.d0) THEN
          IF (inci%R(i).GE.0.d0) THEN
            inci%R(i) = inci%R(i) - pi
          ELSE
            inci%R(i) = inci%R(i) + pi
          ENDIF
        ENDIF
      ENDDO


!=======================================================================
!
      RETURN
      END