The TELEMAC-MASCARET system  trunk
diffrac.f
Go to the documentation of this file.
1 ! ******************
2  SUBROUTINE diffrac
3 ! ******************
4 !
5  &( cx , cy , ct , xk , cg , npoin2, ndire ,
6  & iff , nf , f , rx , ry , rxx , ryy , neigb )
7 !
8 !***********************************************************************
9 ! TOMAWAC V7P0 25/06/2012
10 !***********************************************************************
11 !
12 !brief COMPUTES DIFFRACTION.
13 !+
14 !+ COMPUTES THE DIFFRACTION TERM AND THE
15 !+ DIFFRACTION-CORRECTED GROUP VELOCITY
16 !
17 !history E. KRIEZI (LNH)
18 !+ 04/12/2006
19 !+ V5P5
20 !+
21 !
22 !history G.MATTAROLO (EDF - LNHE)
23 !+ 23/06/2012
24 !+ V6P2
25 !+ Modification for V6P2
26 !
27 !history J-M HERVOUET (EDF R&D, LNHE)
28 !+ 10/12/2012
29 !+ V6P3
30 !+ 4 subroutines GRAD-... inlined and removed from the tomawac library
31 !+ and then optimised.
32 !
33 !history J-M HERVOUET (EDF R&D, LNHE)
34 !+ 21/03/2013
35 !+ V6P3
36 !+ Now velocities modified, not fully computed.
37 !+ A preliminary call of conwac is requested.
38 !
39 !history J-M HERVOUET (EDF - LNHE)
40 !+ 08/01/2014
41 !+ V7P0
42 !+ CALL PARCOM suppressed by using new argument ASSPAR in VECTOR
43 !
44 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45 !| CG |-->| DISCRETIZED GROUP VELOCITY
46 !| COSTET |-->| COSINE OF TETA ANGLE
47 !| CX |<->| ADVECTION FIELD ALONG X(OR PHI)
48 !| CY |<->| ADVECTION FIELD ALONG Y(OR LAMBDA)
49 !| CT |<->| ADVECTION FIELD ALONG TETA
50 !| DFREQ |-->| FREQUENCY STEPS BETWEEN DISCRETIZED FREQUENCIES
51 !| DIFFRA |-->| IF >0 DIFFRACTION IS CONSIDERED
52 !| IF = 1 MSE FORMULATION
53 !| IF = 2 RMSE FORMULATION
54 !| F |<->| VARIANCE DENSITY DIRECTIONAL SPECTRUM
55 !| FLTDIF |-->| IF TRUE, LOCAL AMPLITUDES ARE FILTERED
56 !| FREQ |-->| DISCRETIZED FREQUENCIES
57 !| IFF |-->| FREQUENCY INDEX
58 !| MAXNSP |-->| CONSTANT FOR MESHFREE TECHNIQUE
59 !| NB_CLOSE |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
60 !| NBOR |-->| GLOBAL NUMBER OF BOUNDARY POINTS
61 !| NEIGB |-->| NEIGHBOUR POINTS FOR MESHFREE METHOD
62 !| NF |-->| NUMBER OF FREQUENCIES
63 !| NDIRE |-->| NUMBER OF DIRECTIONS
64 !| NPOIN2 |-->| NUMBER OF POINTS IN 2D MESH
65 !| NPTFR |-->| NUMBER OF BOUNDARY POINTS
66 !| PROINF |-->| LOGICAL INDICATING INFINITE DEPTH ASSUMPTION
67 !| RX |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
68 !| RXX |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
69 !| RY |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
70 !| RYY |-->| ARRAY USED IN THE MESHFREE TECHNIQUE
71 !| SINTET |-->| SINE OF TETA ANGLE
72 !| SPHE |-->| LOGICAL INDICATING SPHERICAL COORD ASSUMPTION
73 !| XK |-->| DISCRETIZED WAVE NUMBER
74 !| XKONPT |<--| ARRAY USED FOR COMPUTING DIFFRACTION PARAMETER
75 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
76 !
77  USE bief
78  USE declarations_tomawac, ONLY : deupi,st1,st0,ielm2,sa,mesh,
79  & freq , costet, sintet, dfreq, sphe, proinf, nbor, nptfr,
80  & sccg,sdelta,sxkonpt,sddx,sddy,
81  & sqrdelta,sqrccg,frdk,frda,scda,
82  & l_delta,deja_diffrac,
83  & optder, fltdif, diffra, nb_close, f2difm, ampli, div , maxnsp
84 !
85  USE declarations_special
86  USE interface_tomawac, ex_diffrac => diffrac
87  IMPLICIT NONE
88 !
89 !
90 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
91 !
92  INTEGER, INTENT(IN) :: NF,NDIRE,NPOIN2,IFF
93  INTEGER, INTENT(IN) :: NEIGB(npoin2,maxnsp)
94  DOUBLE PRECISION, INTENT(INOUT) :: CX(npoin2,ndire)
95  DOUBLE PRECISION, INTENT(INOUT) :: CY(npoin2,ndire)
96  DOUBLE PRECISION, INTENT(INOUT) :: CT(npoin2,ndire)
97  DOUBLE PRECISION, INTENT(IN) :: CG(npoin2,nf),XK(npoin2,nf)
98  DOUBLE PRECISION, INTENT(IN) :: F(npoin2,ndire,nf)
99  DOUBLE PRECISION, INTENT(IN) :: RX(maxnsp,npoin2)
100  DOUBLE PRECISION, INTENT(IN) :: RY(maxnsp,npoin2)
101  DOUBLE PRECISION, INTENT(IN) :: RXX(maxnsp,npoin2)
102  DOUBLE PRECISION, INTENT(IN) :: RYY(maxnsp,npoin2)
103 
104 !
105 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
106 !
107  INTEGER I, IP, IPOIN
108  DOUBLE PRECISION CDELTA,DELTAN
109 !
110 !***********************************************************************
111 !
112 ! TODO: THERE ARE ENOUGH ARRAYS ELSEWHERE, THIS IS USELESS...
113 !
114  IF(.NOT.deja_diffrac)THEN
115  ALLOCATE(sqrdelta(npoin2))
116  ALLOCATE(sqrccg(npoin2))
117  ALLOCATE(frdk(npoin2,2))
118  ALLOCATE(frda(npoin2,2))
119  ALLOCATE(scda(npoin2,3))
120  ALLOCATE(l_delta(npoin2))
121  deja_diffrac=.true.
122  ENDIF
123 !
124 !-----------------------------------------------------------------------
125 ! INFINITE WATER DEPTH ...
126 !-----------------------------------------------------------------------
127 !
128  IF(proinf) THEN
129 !
130  WRITE(lu,*) ''
131  WRITE(lu,*) '***************************************'
132  WRITE(lu,*) ' ATTENTION : DIFFRACTION IS NOT TAKEN '
133  WRITE(lu,*) ' INTO ACCOUNT IN THE CASE OF INFINITE '
134  WRITE(lu,*) ' WATER DEPTH '
135  WRITE(lu,*) '***************************************'
136  CALL plante(1)
137  stop
138 !
139 !-----------------------------------------------------------------------
140 ! FINITE DEPTH
141 !-----------------------------------------------------------------------
142 !
143  ELSE
144 !
145 ! ------------------------------------------------------------------
146 ! ... CARTESIAN COORDINATES
147 ! ------------------------------------------------------------------
148 !
149  IF(.NOT.sphe) THEN
150 !
151 ! DIFFRACTION IS TAKEN INTO ACCOUNT
152 !
153 ! CCG VECTOR COMPUTATION
154 !
155  DO ipoin=1,npoin2
156  sccg%R(ipoin) = cg(ipoin,iff)*deupi*freq(iff)/xk(ipoin,iff)
157  sxkonpt%R(ipoin)=1.d0/(xk(ipoin,iff)**2)
158  sqrccg(ipoin)=sqrt(abs(sccg%R(ipoin)))
159  ENDDO
160 !
161 ! INVERSE OF INTEGRALS OF TEST FUNCTIONS
162 !
163  CALL vector(st0,'=','MASBAS ',ielm2,1.d0,
164  & st1,st1,st1,st1,st1,st1,mesh,.false.,st1,
165  & asspar=.true.)
166  CALL os('X=1/Y ',x=st0,y=st0)
167 !
168 ! LOOP OVER THE DIRECTIONS
169 !
170  DO ip = 1,ndire
171 !
172 ! COMPUTATION OF LOCAL AMPLITUDES OF DIRECTIONAL SPECTRA
173 !
174  DO ipoin = 1,npoin2
175  ampli(ipoin) = sqrt(2.d0*f(ipoin,ip,iff)*dfreq(iff)
176  & *deupi/ndire)
177  IF(diffra.EQ.2)THEN
178  ampli(ipoin)=ampli(ipoin)*xk(ipoin,iff)*sqrccg(ipoin)
179  ENDIF
180  ENDDO
181 !
182 ! Filtering the local amplitudes of directional spectra
183 !
184  IF(fltdif) CALL filt_sa
185 !
186  CALL vector(st1,'=','GRADF X',ielm2,1.d0,sa,
187  & st0,st0,st0,st0,st0,mesh,.false.,st0,
188  & asspar=.true.)
189  DO i=1,npoin2
190  frda(i,1)=st1%R(i)*st0%R(i)
191  ENDDO
192  CALL vector(st1,'=','GRADF Y',ielm2,1.d0,sa,
193  & st0,st0,st0,st0,st0,mesh,.false.,st0,
194  & asspar=.true.)
195  DO i=1,npoin2
196  frda(i,2)=st1%R(i)*st0%R(i)
197  ENDDO
198 !
199 ! DIFFRA=1 - Mean Slope Equation model
200 ! DIFFRA=2 - Revised Mean Slope Equation model
201 !
202  IF(diffra.EQ.1) THEN
203 !
204  CALL vector(st1,'=','GRADF X',ielm2,1.d0,sccg,
205  & st0,st0,st0,st0,st0,mesh,.false.,st0,
206  & asspar=.true.)
207  DO i=1,npoin2
208  frdk(i,1)=st1%R(i)*st0%R(i)
209  ENDDO
210  CALL vector(st1,'=','GRADF Y',ielm2,1.d0,sccg,
211  & st0,st0,st0,st0,st0,mesh,.false.,st0,
212  & asspar=.true.)
213  DO i=1,npoin2
214  frdk(i,2)=st1%R(i)*st0%R(i)
215  ENDDO
216 !
217  ELSE
218 !
219  CALL vector(st1,'=','GRADF X',ielm2,1.d0,sxkonpt,
220  & st0,st0,st0,st0,st0,mesh,.false.,st0,
221  & asspar=.true.)
222  DO i=1,npoin2
223  frdk(i,1)=st1%R(i)*st0%R(i)
224  ENDDO
225  CALL vector(st1,'=','GRADF Y',ielm2,1.d0,sxkonpt,
226  & st0,st0,st0,st0,st0,mesh,.false.,st0,
227  & asspar=.true.)
228  DO i=1,npoin2
229  frdk(i,2)=st1%R(i)*st0%R(i)
230  ENDDO
231 !
232  ENDIF
233 !
234 !
235 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
236 ! CALCUL DE D2A/DX2 + D2A/DY2 AVEC LA METHODE FREEMESH
237 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238 !
239  IF(optder.EQ.1) THEN
240 !
241  DO ipoin = 1,npoin2
242 ! calculate first and second derivative of A ( FFD=A)
243 ! TODO: note JMH: pour first il y a .FALSE. !!!
244 ! et seul SCDA(IPOIN,3) est utilise
245  CALL rpi_intr(neigb,nb_close,
246  & rx(1,ipoin),ry(1,ipoin),
247  & rxx(1,ipoin),ryy(1,ipoin),
248  & npoin2,ipoin,maxnsp,ampli,
249  & frda(ipoin,1),frda(ipoin,2),
250  & scda(ipoin,1),scda(ipoin,2),
251  & scda(ipoin,3),.false.,.true.)
252  ENDDO
253 !
254  ELSEIF(optder.EQ.2) THEN
255 !
256 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
257 ! CALCUL DE D2A/DX2 + D2A/DY2 AVEC LA METHODE BETE
258 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259 !
260 ! FRDA not very convenient, copy required here...
261 !
262  DO i=1,npoin2
263  sddx%R(i)=frda(i,1)
264  sddy%R(i)=frda(i,2)
265  ENDDO
266  CALL vector(st1,'=','GRADF X',ielm2,1.d0,sddx,
267  & st0,st0,st0,st0,st0,mesh,.false.,st0,
268  & asspar=.true.)
269  DO i=1,npoin2
270  scda(i,3)=st1%R(i)*st0%R(i)
271  ENDDO
272  CALL vector(st1,'=','GRADF Y',ielm2,1.d0,sddy,
273  & st0,st0,st0,st0,st0,mesh,.false.,st0,
274  & asspar=.true.)
275  DO i=1,npoin2
276  scda(i,3)=scda(i,3)+st1%R(i)*st0%R(i)
277  ENDDO
278 !
279  ELSE
280  WRITE(lu,*) 'OPTDER=',optder,' NOT TREATED'
281  CALL plante(1)
282  stop
283  ENDIF
284 !
285 ! DIFFRA=1 - Mean Slope Equation model
286 ! DIFFRA=2 - Revised Mean Slope Equation model
287 !
288  IF(diffra.EQ.1)THEN
289  DO ipoin = 1,npoin2
290  div(ipoin)=sccg%R(ipoin)*scda(ipoin,3)
291  & + frdk(ipoin,1)*frda(ipoin,1)
292  & + frdk(ipoin,2)*frda(ipoin,2)
293  ENDDO
294  ELSE
295  DO ipoin = 1,npoin2
296  div(ipoin)=sxkonpt%R(ipoin)*scda(ipoin,3)
297  & + frdk(ipoin,1)*frda(ipoin,1)
298  & + frdk(ipoin,2)*frda(ipoin,2)
299  ENDDO
300  ENDIF
301 !
302 ! Calculating Delta=div/A
303 !
304  DO ipoin = 1,npoin2
305  l_delta(ipoin)=.true.
306  IF(f(ipoin,ip,iff).LE.f2difm) THEN
307  sdelta%R(ipoin) = 0.d0
308  l_delta(ipoin)=.false.
309  sqrdelta(ipoin) =1.d0
310  ELSE
311 ! DIFFRA=1 - Mean Slope Equation model
312 ! DIFFRA=2 - Revised Mean Slope Equation model
313  IF(diffra.EQ.1) THEN
314  sdelta%R(ipoin)=div(ipoin)*sxkonpt%R(ipoin)/
315  & (sccg%R(ipoin)*ampli(ipoin))
316  ELSE
317  sdelta%R(ipoin)=(div(ipoin)/ampli(ipoin))
318  ENDIF
319 !
320  IF(sdelta%R(ipoin).LE.-1.d0) THEN
321 ! TODO: JMH: discutable !!!!!!
322  sqrdelta(ipoin) =1.d0
323  l_delta(ipoin)=.false.
324  sdelta%R(ipoin)= 0.d0
325  ELSE
326  sqrdelta(ipoin) = sqrt(1.d0+sdelta%R(ipoin))
327  l_delta(ipoin)=.true.
328  ENDIF
329 ! TODO: JMH: discutable !!!!!
330  IF(sqrdelta(ipoin).LE.f2difm) THEN
331  sqrdelta(ipoin) =1.d0
332  l_delta(ipoin)=.false.
333  sdelta%R(ipoin)= 0.d0
334  ENDIF
335  ENDIF
336  ENDDO
337 !
338  DO i = 1,nptfr
339  ipoin = nbor(i)
340  l_delta(ipoin)=.false.
341 ! TODO: JMH: discutable !!!
342  sdelta%R(ipoin)= 0.d0
343  ENDDO
344 !
345 ! DELTA GRADIENT COMPUTATION
346 !
347  CALL vector(st1,'=','GRADF X',ielm2,1.d0,sdelta,
348  & st0,st0,st0,st0,st0,mesh,.false.,st0,
349  & asspar=.true.)
350  CALL ov('X=YZ ', x=sddx%R, y=st1%R, z=st0%R, dim1=npoin2)
351  CALL vector(st1,'=','GRADF Y',ielm2,1.d0,sdelta,
352  & st0,st0,st0,st0,st0,mesh,.false.,st0,
353  & asspar=.true.)
354  CALL ov('X=YZ ', x=sddy%R, y=st1%R, z=st0%R, dim1=npoin2)
355 !
356 ! calculation of CG_n =CG(1+delta)^0.5
357 ! and of modified transfer rates Cx,Cy,Ctheta
358 !
359  DO ipoin=1,npoin2
360  IF(l_delta(ipoin)) THEN
361  deltan = -sintet(ip)*sddy%R(ipoin)+costet(ip)*sddx%R(ipoin)
362  cdelta = cg(ipoin,iff)/sqrdelta(ipoin)/2.d0
363  ct(ipoin,ip)=ct(ipoin,ip)*sqrdelta(ipoin)-cdelta*deltan
364  cx(ipoin,ip)=cx(ipoin,ip)*sqrdelta(ipoin)
365  cy(ipoin,ip)=cy(ipoin,ip)*sqrdelta(ipoin)
366  ENDIF
367  ENDDO
368 !
369  ENDDO ! IP
370 !
371 ! ----------------------------------------------------------------
372 ! ... AND SPHERICAL COORDINATES
373 ! ----------------------------------------------------------------
374 !
375  ELSE
376 !
377  WRITE(lu,*) ''
378  WRITE(lu,*) '***************************************'
379  WRITE(lu,*) ' ATTENTION : THE PRESENT VERSION OF '
380  WRITE(lu,*) ' TOMAWAC CANNOT SIMULATE DIFFRACTION '
381  WRITE(lu,*) ' WHEN SPHERICAL COORDINATES ARE SET '
382  WRITE(lu,*) '***************************************'
383  CALL plante(1)
384  stop
385 !
386 ! ENDIF (Finite depth)
387  ENDIF
388 !ENDIF (Cartesian coordinates)
389  ENDIF
390 !
391 !-----------------------------------------------------------------------
392 !
393  RETURN
394  END
ov
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7
declarations_special
Definition: declarations_special.F:3
filt_sa
subroutine filt_sa
Definition: filt_sa.f:4
rpi_intr
subroutine rpi_intr(NEIGB, NB_CLOSE, RX, RY, RXX, RYY, NPOIN2, I, MAXNSP, FFD, FIRDIV1, FIRDIV2, SECDIV1, SECDIV2, SECDIV3, FRSTDIV, SCNDDIV)
Definition: rpi_intr.f:9
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_tomawac
Definition: declarations_tomawac.f:3
declarations_tomawac::st0
type(bief_obj), target st0
Definition: declarations_tomawac.f:95
declarations_tomawac::deupi
double precision deupi
Definition: declarations_tomawac.f:1298
declarations_tomawac::sa
type(bief_obj), target sa
Definition: declarations_tomawac.f:215
diffrac
subroutine diffrac(CX, CY, CT, XK, CG, NPOIN2, NDIRE, IFF, NF, F, RX, RY, RXX, RYY, NEIGB)
Definition: diffrac.f:8
declarations_tomawac::st1
type(bief_obj), target st1
Definition: declarations_tomawac.f:95
vector
subroutine vector(VEC, OP, FORMUL, IELM1, XMUL, F, G, H, U, V, W, MESH, MSK, MASKEL, LEGO, ASSPAR)
Definition: vector.f:7
os
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
declarations_tomawac::ielm2
integer ielm2
Definition: declarations_tomawac.f:874
interface_tomawac
Definition: interface_tomawac.f:5
declarations_tomawac::mesh
type(bief_mesh), target mesh
Definition: declarations_tomawac.f:211
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:35:57 for The TELEMAC-MASCARET system by doxygen 1.8.13