The TELEMAC-MASCARET system  trunk
caldir.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE caldir
3 ! *****************
4 !
5 !
6 !***********************************************************************
7 ! ARTEMIS V7P3 Aug 2017
8 !***********************************************************************
9 !
10 !brief COMPUTES THE WAVE INCIDENCE
11 !
12 !history C PEYRARD (LNHE)
13 !+
14 !+
15 !+ Creation of DOXYGEN tags for automated documentation and
16 !+ cross-referencing of the FORTRAN sources
17 !
18 !history N.DURAND (HRW)
19 !+ August 2017
20 !+ V7P3
21 !+ PI, DEUPI and PISUR2 now defined in DECLARATIONS_ARTEMIS
22 !
23 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
24 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
25 !
26  USE bief
29 !
31  IMPLICIT NONE
32 !
33  INTEGER I
34 !
35  DOUBLE PRECISION A1, A2 ,ALPHA0, D1, D2, PHI, PHI1, PHI2 ,MODPHI
36  DOUBLE PRECISION TETA01, XU1 ,XU2, XV1, XV2, WT0
37 !
38  INTRINSIC sqrt, atan2, mod, abs, cos, sin, atan
39 !
40 !-----------------------------------------------------------------------
41 !
42  DOUBLE PRECISION, PARAMETER :: ZERO = 1.d-10
43 !
44 !=======================================================================
45 ! SPEEDS AT THE SURFACE (AT T=0 AND T=OMEGA/4)
46 !=======================================================================
47 !
48 ! COMPUTES THE GRADIENTS (PHIR AND PHII)
49 !
50 !
51  CALL vector(u0 , '=' , 'GRADF X' , ielm ,
52  & 1.d0 , phir , sbid, sbid , sbid , sbid , sbid ,
53  & mesh , msk , maskel )
54 !
55  CALL vector(v0 , '=' , 'GRADF Y' , ielm ,
56  & 1.d0 , phir , sbid , sbid , sbid , sbid , sbid ,
57  & mesh , msk , maskel )
58 !
59 ! THE OLD VARIABLE U1 IS STORED IN T3
60 ! BECAUSE IT IS USED TO COMPUTE INCI
61 !
62  CALL vector(t3 , '=' , 'GRADF X' , ielm ,
63  & 1.d0 , phii , sbid , sbid , sbid , sbid , sbid ,
64  & mesh , msk , maskel )
65 !
66 ! THE OLD VARIABLE V1 IS STORED IN T4
67 ! BECAUSE IT IS USED TO COMPUTE INCI
68 !
69  CALL vector(t4 , '=' , 'GRADF Y' , ielm ,
70  & 1.d0 , phii , sbid , sbid , sbid , sbid , sbid ,
71  & mesh , msk , maskel )
72 !
73  CALL vector(t1 , '=' , 'MASBAS ' , ielm ,
74  & 1.d0 , sbid , sbid , sbid , sbid , sbid , sbid ,
75  & mesh , msk , maskel )
76 !
77  CALL os('X=Y/Z ', x=u0, y=u0, z=t1)
78  CALL os('X=Y/Z ', x=v0, y=v0, z=t1)
79  CALL os('X=Y/Z ', x=t3, y=t3, z=t1)
80  CALL os('X=Y/Z ', x=t4, y=t4, z=t1)
81 !
82 !=======================================================================
83 ! COMPUTES WAVE INCIDENCE
84 !=======================================================================
85 !
86 ! U0 (D(PHIR)/DX) : A U1 (D(PHII)/DX): B
87 ! V0 (D(PHIR)/DY) : C V1 (D(PHII)/DY): D
88 ! FROM U= A COS WT + B SIN WT TO : U = A1 COS ( WT - PHI1)
89 ! V= C COS WT + D SIN WT V = A2 COS ( WT - PHI2)
90 !
91  DO i=1,npoin
92  a1 = sqrt( u0%R(i)*u0%R(i) + t3%R(i)*t3%R(i) )
93  phi1 = atan2( t3%R(i),u0%R(i) )
94  a2 = sqrt( v0%R(i)*v0%R(i) + t4%R(i)*t4%R(i) )
95  IF(t4%R(i).EQ.0.d0.AND.v0%R(i).EQ.0.d0) THEN
96  phi2 = 2.d0*atan(1.d0)
97  ELSE
98  phi2 = atan2( t4%R(i),v0%R(i) )
99  ENDIF
100 !
101 ! WRITTEN AS : U = A1 COS ( (WT - PHI1))
102 ! V = A2 COS ( (WT - PHI1) - PHI )
103 ! WHERE PHI IS BETWEEN 0 AND 2*PI
104 !
105  phi = phi2 - phi1
106  IF (phi.LT.0.d0) phi = phi+deupi
107 !
108 ! ESTIMATES THE DIRECTION AND (WT0) WHEN THE ELLIPSE'S MAJOR AXIS
109 ! IS REACHED.
110 ! TREATS INDIVIDUAL CASES (LINEAR POLARISATION)
111 !
112  modphi = mod( phi, pi )
113  IF ( (modphi.LT.zero).OR.((pi-modphi).LT.zero) ) THEN
114  wt0 = phi1
115  IF ( (phi.LT.2d0*zero).OR.((deupi-phi).LT.2d0*zero) )THEN
116  alpha0 = atan2( a2,a1 )
117  ELSE
118 ! (ABS(PHI-PI).LT.2D0*ZERO)
119  alpha0 = deupi - atan2( a2,a1 )
120  ENDIF
121  ELSE
122 ! GENERAL CASE: ELLIPTIC POLARISATION
123 ! TAN(2*(WT0 - PHI1)) = A2**2*SIN(2*PHI)/(A1**2+A2**2*COS(2*PHI))
124  teta01 = atan2( (a2*a2*sin(2*phi)) ,
125  & (a1*a1 + a2*a2*cos(2*phi)) ) / 2.d0
126  xu1 = a1 * cos( teta01)
127  xv1 = a2 * cos( teta01 - phi )
128  xu2 = -a1 * sin( teta01)
129  xv2 = -a2 * sin( teta01 - phi )
130  d1 = xu1*xu1 + xv1*xv1
131  d2 = xu2*xu2 + xv2*xv2
132  IF (d2.GT.d1) THEN
133  teta01 = teta01 + pisur2
134  xu1 = xu2
135  xv1 = xv2
136  ENDIF
137  wt0 = teta01 + phi1
138  alpha0 = atan2( xv1,xu1 )
139  ENDIF
140  inci%R(i) = alpha0
141  t2%R(i) = wt0
142  ENDDO
143 !
144 ! FREE SURFACE IN PHASE WITH ALPHA0
145 ! INCIDENCE IS CONSIDERED POSITIVE WHEN THE FREE SURFACE IS
146 ! POSITIVE.
147 !
148  DO i=1,npoin
149  a1 = -(phii%R(i)*cos(t2%R(i))-phir%R(i)*sin(t2%R(i)))
150  IF (a1.LT.0.d0) THEN
151  IF (inci%R(i).GE.0.d0) THEN
152  inci%R(i) = inci%R(i) - pi
153  ELSE
154  inci%R(i) = inci%R(i) + pi
155  ENDIF
156  ENDIF
157  ENDDO
158 
159 
160 !=======================================================================
161 !
162  RETURN
163  END
type(bief_obj), target maskel
double precision, dimension(:), pointer y
type(bief_obj), target phii
type(bief_obj), target v0
type(bief_obj), pointer t2
type(bief_mesh), target mesh
subroutine caldir
Definition: caldir.f:4
type(bief_obj), target inci
type(bief_obj), pointer t4
type(bief_obj), target phir
double precision, dimension(:), pointer x
type(bief_obj), target sbid
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
type(bief_obj), pointer t3
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
type(bief_obj), pointer t1
type(bief_obj), target u0
Definition: bief.f:3