The TELEMAC-MASCARET system  trunk
qnlin1.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE qnlin1
3 ! *****************
4 !
5  &( tstot , tsder , iangnl, nf , ndire , npoin2, f ,
6  & xkmoy , taux1 , taux2 , taux3 , taux4 , taux5 , dfini )
7 !
8 !***********************************************************************
9 ! TOMAWAC V6P3 24/06/2011
10 !***********************************************************************
11 !
12 !brief COMPUTES THE CONTRIBUTION OF THE NON-LINEAR INTERACTIONS
13 !+ SOURCE TERM (FREQUENCY QUADRUPLETS) USING THE DIA METHOD
14 !+ ("DISCRETE INTERACTION APPROXIMATION") PROPOSED BY
15 !+ HASSELMANN AND HASSELMANN (1985).
16 !+
17 !+
18 !+ PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES
19 !+ FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS
20 !+ ARE EVENLY DISTRIBUTED OVER [0;2.PI].
21 !
22 !note THIS SUBROUTINE USES THE OUTPUT FROM 'PRENL1' TO OPTIMISE
23 !+ THE COMPUTATIONS FOR DIA.
24 !
25 !reference HASSELMANN S., HASSELMANN K. ET AL.(1985) :
26 !+ "COMPUTATIONS AND PARAMETERIZATIONS OF THE NONLINEAR
27 !+ ENERGY TRANSFER IN GRAVITY-WAVE SPECTRUM. PART1 :
28 !+ A NEW METHOD FOR EFFICIENT COMPUTATION OF THE EXACT
29 !+ NON-LINEAR TRANSFER INTEGRAL". JPO, VOL 15, PP 1369-1377.
30 !reference HASSELMANN S., HASSELMANN K. ET AL.(1985) :
31 !+ "COMPUTATIONS AND PARAMETERIZATIONS OF THE NONLINEAR
32 !+ ENERGY TRANSFER IN GRAVITY-WAVE SPECTRUM. PART2 :
33 !+ PARAMETERIZATIONS OF THE NONLINEAR ENERGY TRANSFER
34 !+ FOR APPLICATION IN WAVE MODELS". JPO, VOL 15, PP 1378-1391.
35 !
36 !history M. BENOIT
37 !+ 07/06/95
38 !+ V1P0
39 !+ CREATED
40 !
41 !history M. BENOIT
42 !+ 26/06/96
43 !+ V1P2
44 !+ MODIFIED
45 !
46 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
47 !+ 13/07/2010
48 !+ V6P0
49 !+ Translation of French comments within the FORTRAN sources into
50 !+ English comments
51 !
52 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
53 !+ 21/08/2010
54 !+ V6P0
55 !+ Creation of DOXYGEN tags for automated documentation and
56 !+ cross-referencing of the FORTRAN sources
57 !
58 !history G.MATTAROLO (EDF - LNHE)
59 !+ 24/06/2011
60 !+ V6P1
61 !+ Translation of French names of the variables in argument
62 !
63 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
64 !| COEFNL |-->| COEFFICIENTS USED FOR DIA METHOD
65 !| DEPTH |-->| WATER DEPTH
66 !| DFINI |<->| WORK TABLE
67 !| F |-->| DIRECTIONAL SPECTRUM
68 !| F1 |-->| FIRST DISCRETIZED FREQUENCY
69 !| IANGNL |-->| ANGULAR INDICES TABLE
70 !| NF |-->| NUMBER OF FREQUENCIES
71 !| NDIRE |-->| NUMBER OF DIRECTIONS
72 !| NPOIN2 |-->| NUMBER OF POINTS IN 2D MESH
73 !| PROINF |-->| LOGICAL INDICATING INFINITE DEPTH ASSUMPTION
74 !| RAISF |-->| FREQUENTIAL RATIO
75 !| TAUX1 |<->| WORK TABLE
76 !| TAUX2 |<->| WORK TABLE
77 !| TAUX3 |<->| WORK TABLE
78 !| TAUX4 |<->| WORK TABLE
79 !| TAUX5 |<->| WORK TABLE
80 !| TSDER |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
81 !| TSTOT |<->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
82 !| XKMOY |-->| AVERAGE WAVE NUMBER
83 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
84 !
85  USE interface_tomawac, ex_qnlin1 => qnlin1
86  USE declarations_tomawac, ONLY : depth, proinf, f1, raisf, coefnl
87  IMPLICIT NONE
88 !
89 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
90 !
91  INTEGER, INTENT(IN) :: NPOIN2,NDIRE,NF
92  INTEGER, INTENT(IN) :: IANGNL(ndire,8)
93  DOUBLE PRECISION, INTENT(IN) :: F(npoin2,ndire,nf)
94  DOUBLE PRECISION, INTENT(IN) :: XKMOY(npoin2)
95  DOUBLE PRECISION, INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
96  DOUBLE PRECISION, INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
97  DOUBLE PRECISION, INTENT(INOUT) :: TAUX1(npoin2),TAUX2(npoin2)
98  DOUBLE PRECISION, INTENT(INOUT) :: TAUX3(npoin2)
99  DOUBLE PRECISION, INTENT(INOUT) :: TAUX4(npoin2),TAUX5(npoin2)
100  DOUBLE PRECISION, INTENT(INOUT) :: DFINI(npoin2)
101 !
102 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
103 !
104  INTEGER JBP0 , JFP0 , JFP1 , JFM0 , JFM1 , JFP , JFM
105  INTEGER JBP1 , JB , JBM0 , JBM1 , IMAGE , JP
106  INTEGER JPP0 , JPP1 , JPM0 , JPM1 , IP , KAUX , JF
107  INTEGER JFMIN , JFMAX
108  DOUBLE PRECISION COEFP0, COEFP1, COEFM0, COEFM1, COEFJF, XXFAC
109  DOUBLE PRECISION FMOIN , FPLUS , TERM1 , TERM2 , US1PL4, US1ML4
110  DOUBLE PRECISION C1 , C2 , C3 , C4 , C5 , C6
111  DOUBLE PRECISION D1 , D2 , D3 , D4 , D5 , D6
112  DOUBLE PRECISION C1SQ , C2SQ , C3SQ , C4SQ , C5SQ , C6SQ
113  DOUBLE PRECISION C7 , C8 , D7 , D8 , C7SQ , C8SQ
114  DOUBLE PRECISION TERM3 , FDEJF , FREQ
115 !
116 !-----------------------------------------------------------------------
117 !
118 ! RECOVERS THE COEFFICIENTS COMPUTED IN 'PRENL1'
119 !
120  c1 = coefnl( 1)
121  c2 = coefnl( 2)
122  c3 = coefnl( 3)
123  c4 = coefnl( 4)
124  c5 = coefnl( 5)
125  c6 = coefnl( 6)
126  c7 = coefnl( 7)
127  c8 = coefnl( 8)
128  jfp = int(coefnl( 9)+1.d-7)
129  jfm = int(coefnl(10)-1.d-7)
130  us1pl4= coefnl(11)
131  us1ml4= coefnl(12)
132  jfmin = nint(coefnl(13))
133  jfmax = nint(coefnl(14))
134  c1sq = c1**2
135  c2sq = c2**2
136  c3sq = c3**2
137  c4sq = c4**2
138  c5sq = c5**2
139  c6sq = c6**2
140  c7sq = c7**2
141  c8sq = c8**2
142 !
143 ! CORRECTION FACTOR FOR FINITE WATER DEPTH
144 !
145  IF(.NOT.proinf) THEN
146  DO ip=1,npoin2
147  term1 = max(0.75d0*depth(ip)*xkmoy(ip),0.5d0)
148  dfini(ip) = 1.d0+(5.5d0/term1)*(1.d0-0.833d0*term1)
149  & /exp(min(1.25d0*term1,7.d2))
150  ENDDO
151  ENDIF
152 !
153 ! FIRST LOOP ON THE FREQUENCIES
154 !
155  DO jf=jfmin,jfmax
156 !
157 ! COMPUTES THE CONSIDERED FREQUENCY
158 !
159  freq = f1*raisf**(jf-1)
160 !
161 ! GETS THE INDICES OF THE FREQUENCIES EITHER SIDE OF THE
162 ! 'MAX' FREQUENCY: FREQ(JFP0)
163 !
164  jfp0=jf+jfp
165  jfp1=jfp0+1
166 !
167 ! GETS THE INDICES OF THE FREQUENCIES EITHER SIDE OF THE
168 ! 'MIN' FREQUENCY: FREQ(JFM0)
169 !
170  jfm0=jf+jfm-1
171  jfm1=jfm0+1
172 !
173 ! LIMITS THE INDICES TO NF AND TAKES INTO ACCOUNT ANALYTICALLY
174 ! THE SPECTRUM TAIL (DECREASE IN -TAILF).
175 !
176  CALL cqueue( jfp1 , jbp1 , coefp1 )
177  CALL cqueue( jfp0 , jbp0 , coefp0 )
178  CALL cqueue( jf , jb , coefjf )
179  CALL cqueue( jfm1 , jbm1 , coefm1 )
180  CALL cqueue( jfm0 , jbm0 , coefm0 )
181 !
182 ! INTERPOLATION COEFFICIENTS FOR THE MODIFIED SPECTRUM
183 !
184  d1=c1*coefp0*us1pl4
185  d2=c2*coefp0*us1pl4
186  d3=c3*coefp1*us1pl4
187  d4=c4*coefp1*us1pl4
188  d5=c5*coefm0*us1ml4
189  d6=c6*coefm0*us1ml4
190  d7=c7*coefm1*us1ml4
191  d8=c8*coefm1*us1ml4
192 !
193 ! COMPUTES THE MULTIPLICATIVE COEFFICIENT (IN F**11) AND TAKES
194 ! INTO ACCOUNT THE CORRECTION TERM IN FINITE DEPTH
195 !
196  xxfac= 3000.d0*freq**11
197  IF(proinf) THEN
198  DO ip=1,npoin2
199  taux1(ip) = xxfac
200  ENDDO
201  ELSE
202  DO ip=1,npoin2
203  taux1(ip) = dfini(ip)*xxfac
204  ENDDO
205  ENDIF
206 !
207 ! SECOND LOOP ON ANGULAR SYMMETRY
208 !
209  DO image=1,2
210 !
211  kaux=(image-1)*4
212 !
213 ! THIRD LOOP ON THE DIRECTIONS
214 !
215  DO jp=1,ndire
216 !
217  jpp0 = iangnl(jp,kaux+1)
218  jpp1 = iangnl(jp,kaux+2)
219  jpm0 = iangnl(jp,kaux+3)
220  jpm1 = iangnl(jp,kaux+4)
221 !
222  IF (jfm0.LT.1) THEN
223 !
224 !........./-------------------------------------------------------/
225 !........./ AT LEAST ONE OF THE FREQUENCIES IS LOWER THAN FREQ(1) /
226 !........./ THE SPECTRUM F- WITH FREQUENCY (1-XLAMD).FREQ IS ZERO /
227 !........./-------------------------------------------------------/
228 !
229  DO ip=1,npoin2
230  fdejf = f(ip,jp,jb )*coefjf
231  fplus = f(ip,jpp0,jbp0)*d1 + f(ip,jpp1,jbp0)*d2
232  & + f(ip,jpp0,jbp1)*d3 + f(ip,jpp1,jbp1)*d4
233 !
234  term1 = fdejf*fplus
235  term3 = taux1(ip)*fdejf
236 !
237  taux2(ip) = term1*term3
238  taux3(ip) = 2.d0*term1*taux1(ip)
239  taux5(ip) = fdejf*us1pl4*term3
240  ENDDO ! IP
241 !
242  IF (jb.EQ.jf) THEN
243 !
244  DO ip=1,npoin2
245  tstot(ip,jp ,jf )=tstot(ip,jp ,jf )-taux2(ip)*2.d0
246  tsder(ip,jp ,jf )=tsder(ip,jp ,jf )-taux3(ip)*2.d0
247  ENDDO ! IP
248 !
249  IF (jbp0.EQ.jfp0) THEN
250 !
251  DO ip=1,npoin2
252  tstot(ip,jpp0,jfp0)=tstot(ip,jpp0,jfp0)+taux2(ip)*c1
253  tstot(ip,jpp1,jfp0)=tstot(ip,jpp1,jfp0)+taux2(ip)*c2
254  tsder(ip,jpp0,jfp0)=tsder(ip,jpp0,jfp0)
255  & +taux5(ip)*c1sq
256  tsder(ip,jpp1,jfp0)=tsder(ip,jpp1,jfp0)
257  & +taux5(ip)*c2sq
258  ENDDO ! IP
259 !
260  IF (jbp1.EQ.jfp1) THEN
261 !
262  DO ip=1,npoin2
263  tstot(ip,jpp0,jfp1)=tstot(ip,jpp0,jfp1)
264  & +taux2(ip)*c3
265  tstot(ip,jpp1,jfp1)=tstot(ip,jpp1,jfp1)
266  & +taux2(ip)*c4
267  tsder(ip,jpp0,jfp1)=tsder(ip,jpp0,jfp1)
268  & +taux5(ip)*c3sq
269  tsder(ip,jpp1,jfp1)=tsder(ip,jpp1,jfp1)
270  & +taux5(ip)*c4sq
271  ENDDO ! IP
272 !
273  ENDIF
274  ENDIF
275  ENDIF
276 !
277  ELSE
278 !
279 !........./--------------------------------------------------------/
280 !........./ FREQUENCIES F-, F, F+ MAY HAVE ENERGY /
281 !........./--------------------------------------------------------/
282 !
283  DO ip=1,npoin2
284  fdejf = f(ip,jp,jb )*coefjf
285  fplus = f(ip,jpp0,jbp0)*d1 + f(ip,jpp1,jbp0)*d2
286  & + f(ip,jpp0,jbp1)*d3 + f(ip,jpp1,jbp1)*d4
287  fmoin = f(ip,jpm0,jbm0)*d5 + f(ip,jpm1,jbm0)*d6
288  & + f(ip,jpm0,jbm1)*d7 + f(ip,jpm1,jbm1)*d8
289 !
290  term1 = fdejf*(fplus+fmoin)
291  term2 = 2.d0*fplus*fmoin
292  term3 = taux1(ip)*fdejf
293 !
294  taux2(ip) = (term1-term2)*term3
295  taux3(ip) = (2.d0*term1-term2)*taux1(ip)
296  taux5(ip) = (fdejf-2.d0*fmoin)*us1pl4*term3
297  taux4(ip) = (fdejf-2.d0*fplus)*us1ml4*term3
298  ENDDO ! IP
299 !
300  IF (jbm0.EQ.jfm0) THEN
301 !
302  DO ip=1,npoin2
303  tstot(ip,jpm0,jfm0)=tstot(ip,jpm0,jfm0)+taux2(ip)*c5
304  tstot(ip,jpm1,jfm0)=tstot(ip,jpm1,jfm0)+taux2(ip)*c6
305  tsder(ip,jpm0,jfm0)=tsder(ip,jpm0,jfm0)+taux4(ip)*c5sq
306  tsder(ip,jpm1,jfm0)=tsder(ip,jpm1,jfm0)+taux4(ip)*c6sq
307  ENDDO ! IP
308 !
309  IF (jbm1.EQ.jfm1) THEN
310 !
311  DO ip=1,npoin2
312  tstot(ip,jpm0,jfm1)=tstot(ip,jpm0,jfm1)+taux2(ip)*c7
313  tstot(ip,jpm1,jfm1)=tstot(ip,jpm1,jfm1)+taux2(ip)*c8
314  tsder(ip,jpm0,jfm1)=tsder(ip,jpm0,jfm1)
315  & +taux4(ip)*c7sq
316  tsder(ip,jpm1,jfm1)=tsder(ip,jpm1,jfm1)
317  & +taux4(ip)*c8sq
318  ENDDO ! IP
319 !
320  IF (jb.EQ.jf) THEN
321 !
322  DO ip=1,npoin2
323  tstot(ip,jp ,jf )=tstot(ip,jp ,jf )
324  & -taux2(ip)*2.d0
325  tsder(ip,jp ,jf )=tsder(ip,jp ,jf )
326  & -taux3(ip)*2.d0
327  ENDDO ! IP
328 !
329  IF (jbp0.EQ.jfp0) THEN
330 !
331  DO ip=1,npoin2
332  tstot(ip,jpp0,jfp0)=tstot(ip,jpp0,jfp0)
333  & +taux2(ip)*c1
334  tstot(ip,jpp1,jfp0)=tstot(ip,jpp1,jfp0)
335  & +taux2(ip)*c2
336  tsder(ip,jpp0,jfp0)=tsder(ip,jpp0,jfp0)
337  & +taux5(ip)*c1sq
338  tsder(ip,jpp1,jfp0)=tsder(ip,jpp1,jfp0)
339  & +taux5(ip)*c2sq
340  ENDDO ! IP
341 !
342  IF (jbp1.EQ.jfp1) THEN
343 !
344  DO ip=1,npoin2
345  tstot(ip,jpp0,jfp1)=tstot(ip,jpp0,jfp1)
346  & +taux2(ip)*c3
347  tstot(ip,jpp1,jfp1)=tstot(ip,jpp1,jfp1)
348  & +taux2(ip)*c4
349  tsder(ip,jpp0,jfp1)=tsder(ip,jpp0,jfp1)
350  & +taux5(ip)*c3sq
351  tsder(ip,jpp1,jfp1)=tsder(ip,jpp1,jfp1)
352  & +taux5(ip)*c4sq
353  ENDDO ! IP
354 !
355  ENDIF
356  ENDIF
357  ENDIF
358  ENDIF
359  ENDIF
360 !
361  ENDIF
362 !
363  ENDDO ! JP
364 !
365  ENDDO ! IMAGE
366 !
367  ENDDO ! JF
368 !
369 !-----------------------------------------------------------------------
370 !
371  RETURN
372  END
declarations_tomawac::depth
double precision, dimension(:), pointer depth
Definition: declarations_tomawac.f:1010
declarations_tomawac::proinf
logical proinf
Definition: declarations_tomawac.f:788
qnlin1
subroutine qnlin1(TSTOT, TSDER, IANGNL, NF, NDIRE, NPOIN2, F, XKMOY, TAUX1, TAUX2, TAUX3, TAUX4, TAUX5, DFINI)
Definition: qnlin1.f:8
declarations_tomawac::coefnl
double precision, dimension(:), pointer coefnl
Definition: declarations_tomawac.f:1079
declarations_tomawac::raisf
double precision raisf
Definition: declarations_tomawac.f:445
declarations_tomawac
Definition: declarations_tomawac.f:3
cqueue
subroutine cqueue(JFRE, JBIS, COEF1)
Definition: cqueue.f:7
declarations_tomawac::f1
double precision f1
Definition: declarations_tomawac.f:441
interface_tomawac
Definition: interface_tomawac.f:5
Generated on Tue Jan 5 2021 14:36:02 for The TELEMAC-MASCARET system by doxygen 1.8.13