5 &( tstot , tsder , iangnl, coefnl, nf , ndire ,
6 & npoin2, f , xkmoy , taux1 , dfini , xcoef )
69 INTEGER,
INTENT(IN) :: NPOIN2, NDIRE , NF
70 INTEGER,
INTENT(IN) :: IANGNL(ndire,16)
71 DOUBLE PRECISION,
INTENT(IN) :: XCOEF
72 DOUBLE PRECISION,
INTENT(IN) :: F(npoin2,ndire,nf), COEFNL(32)
73 DOUBLE PRECISION,
INTENT(IN) :: XKMOY(npoin2)
74 DOUBLE PRECISION,
INTENT(INOUT) :: TAUX1(npoin2), DFINI(npoin2)
75 DOUBLE PRECISION,
INTENT(INOUT) :: TSTOT(npoin2,ndire,nf)
76 DOUBLE PRECISION,
INTENT(INOUT) :: TSDER(npoin2,ndire,nf)
80 INTEGER IP , JP , JF , JFMIN , JFMAX
81 INTEGER JFP1 , JFM2 , JFP3 , JFM4 ,
82 & jf1_0 , jf1_1 , jf2_0 , jf2_1 , jf3_0 , jf3_1 , jf4_0 , jf4_1 ,
83 & jb1_0 , jb1_1 , jb2_0 , jb2_1 , jb3_0 , jb3_1 , jb4_0 , jb4_1 ,
84 & jp1p_0, jp1p_1, jp1m_0, jp1m_1, jp2p_0, jp2p_1, jp2m_0, jp2m_1,
85 & jp3p_0, jp3p_1, jp3m_0, jp3m_1, jp4p_0, jp4p_1, jp4m_0, jp4m_1
86 DOUBLE PRECISION US1PM4, US1MM4, US1PL4, US1ML4
87 DOUBLE PRECISION FREQ , XXFAC , TERM1 , W ,
88 & c01 , c02 , c03 , c04 , c05 , c06 , c07 , c08 ,
89 & c09 , c10 , c11 , c12 , c13 , c14 , c15 , c16 ,
90 & d01 , d02 , d03 , d04 , d05 , d06 , d07 , d08 ,
91 & d09 , d10 , d11 , d12 , d13 , d14 , d15 , d16 ,
92 & c01sq , c02sq , c03sq , c04sq , c05sq , c06sq , c07sq , c08sq ,
93 & c09sq , c10sq , c11sq , c12sq , c13sq , c14sq , c15sq , c16sq ,
94 & cf1_0 , cf1_1 , cf2_0 , cf2_1 , cf3_0 , cf3_1 , cf4_0 , cf4_1 ,
95 & qnl_a , qnl_b , qnl_c , qnl_d
97 DOUBLE PRECISION :: TMP
98 DOUBLE PRECISION,
DIMENSION(:),
POINTER ::
99 & f1plus, f1moin, f2plus, f2moin, f3plus, f3moin, f4plus, f4moin,
100 & s_1p2m, s_1m2p, s_3p4m, s_3m4p, p_1p2m, p_1m2p, p_3p4m, p_3m4p
101 DOUBLE PRECISION,
DIMENSION(:),
POINTER ::
102 & q_apb , q_apc , q_bpd , q_cpd
107 LOGICAL,
SAVE :: DEJA = .false.
108 DOUBLE PRECISION,
ALLOCATABLE,
TARGET,
DIMENSION(:,:),
SAVE ::
113 ALLOCATE(trav_mdia(npoin2,22),stat=ierr)
114 CALL check_allocate(ierr,
'TRAV_MDIA')
117 f1plus => trav_mdia(:,1)
118 f1moin => trav_mdia(:,2)
119 f2plus => trav_mdia(:,3)
120 f2moin => trav_mdia(:,4)
121 f3plus => trav_mdia(:,5)
122 f3moin => trav_mdia(:,6)
123 f4plus => trav_mdia(:,7)
124 f4moin => trav_mdia(:,8)
125 s_1p2m => trav_mdia(:,9)
126 s_1m2p => trav_mdia(:,10)
127 s_3p4m => trav_mdia(:,11)
128 s_3m4p => trav_mdia(:,12)
129 p_1p2m => trav_mdia(:,13)
130 p_1m2p => trav_mdia(:,14)
131 p_3p4m => trav_mdia(:,15)
132 p_3m4p => trav_mdia(:,16)
133 q_apb => trav_mdia(:,17)
134 q_apc => trav_mdia(:,18)
135 q_bpd => trav_mdia(:,19)
136 q_cpd => trav_mdia(:,20)
160 jfp1 = int(coefnl( 9)+1.d-7)
161 jfm2 = int(coefnl(10)-1.d-7)
162 jfp3 = int(coefnl(25)+1.d-7)
163 jfm4 = int(coefnl(26)-1.d-7)
170 jfmin = min(nint(coefnl(13)),nint(coefnl(29)))
171 jfmax = max(nint(coefnl(14)),nint(coefnl(30)))
195 term1 = max(0.75d0*
depth(ip)*xkmoy(ip),0.5d0)
196 dfini(ip) = 1.d0+(5.5d0/term1)*(1.d0-0.833d0*term1)
197 & *exp(-1.25d0*term1)
207 freq = f1*raisf**(jf-1)
240 CALL cqueue( jf1_1 , jb1_1 , cf1_1 )
241 CALL cqueue( jf1_0 , jb1_0 , cf1_0 )
242 CALL cqueue( jf2_1 , jb2_1 , cf2_1 )
243 CALL cqueue( jf2_0 , jb2_0 , cf2_0 )
244 CALL cqueue( jf3_1 , jb3_1 , cf3_1 )
245 CALL cqueue( jf3_0 , jb3_0 , cf3_0 )
246 CALL cqueue( jf4_1 , jb4_1 , cf4_1 )
247 CALL cqueue( jf4_0 , jb4_0 , cf4_0 )
271 xxfac= 0.5d0 * xcoef /
gravit**4 * freq**11
278 taux1(ip) = dfini(ip)*xxfac
285 jp1p_0 = iangnl(jp, 1)
286 jp1p_1 = iangnl(jp, 2)
287 jp2m_0 = iangnl(jp, 3)
288 jp2m_1 = iangnl(jp, 4)
289 jp1m_0 = iangnl(jp, 5)
290 jp1m_1 = iangnl(jp, 6)
291 jp2p_0 = iangnl(jp, 7)
292 jp2p_1 = iangnl(jp, 8)
293 jp3p_0 = iangnl(jp, 9)
294 jp3p_1 = iangnl(jp,10)
295 jp4m_0 = iangnl(jp,11)
296 jp4m_1 = iangnl(jp,12)
297 jp3m_0 = iangnl(jp,13)
298 jp3m_1 = iangnl(jp,14)
299 jp4p_0 = iangnl(jp,15)
300 jp4p_1 = iangnl(jp,16)
308 f1plus(ip) = f(ip,jp1p_0,jb1_0)*d01 + f(ip,jp1p_1,jb1_0)*d02
309 & + f(ip,jp1p_0,jb1_1)*d03 + f(ip,jp1p_1,jb1_1)*d04
312 f1moin(ip) = f(ip,jp1m_0,jb1_0)*d01 + f(ip,jp1m_1,jb1_0)*d02
313 & + f(ip,jp1m_0,jb1_1)*d03 + f(ip,jp1m_1,jb1_1)*d04
316 f2plus(ip) = f(ip,jp2p_0,jb2_0)*d05 + f(ip,jp2p_1,jb2_0)*d06
317 & + f(ip,jp2p_0,jb2_1)*d07 + f(ip,jp2p_1,jb2_1)*d08
320 f2moin(ip) = f(ip,jp2m_0,jb2_0)*d05 + f(ip,jp2m_1,jb2_0)*d06
321 & + f(ip,jp2m_0,jb2_1)*d07 + f(ip,jp2m_1,jb2_1)*d08
324 f3plus(ip) = f(ip,jp3p_0,jb3_0)*d09 + f(ip,jp3p_1,jb3_0)*d10
325 & + f(ip,jp3p_0,jb3_1)*d11 + f(ip,jp3p_1,jb3_1)*d12
328 f3moin(ip) = f(ip,jp3m_0,jb3_0)*d09 + f(ip,jp3m_1,jb3_0)*d10
329 & + f(ip,jp3m_0,jb3_1)*d11 + f(ip,jp3m_1,jb3_1)*d12
332 f4plus(ip) = f(ip,jp4p_0,jb4_0)*d13 + f(ip,jp4p_1,jb4_0)*d14
333 & + f(ip,jp4p_0,jb4_1)*d15 + f(ip,jp4p_1,jb4_1)*d16
336 f4moin(ip) = f(ip,jp4m_0,jb4_0)*d13 + f(ip,jp4m_1,jb4_0)*d14
337 & + f(ip,jp4m_0,jb4_1)*d15 + f(ip,jp4m_1,jb4_1)*d16
341 s_1p2m(ip) = f1plus(ip) + f2moin(ip)
342 s_1m2p(ip) = f1moin(ip) + f2plus(ip)
343 s_3p4m(ip) = f3plus(ip) + f4moin(ip)
344 s_3m4p(ip) = f3moin(ip) + f4plus(ip)
345 p_1p2m(ip) = f1plus(ip) * f2moin(ip)
346 p_1m2p(ip) = f1moin(ip) * f2plus(ip)
347 p_3p4m(ip) = f3plus(ip) * f4moin(ip)
348 p_3m4p(ip) = f3moin(ip) * f4plus(ip)
353 qnl_a=w*( p_1p2m(ip)*s_3p4m(ip) - p_3p4m(ip)*s_1p2m(ip) )
354 qnl_b=w*( p_1p2m(ip)*s_3m4p(ip) - p_3m4p(ip)*s_1p2m(ip) )
355 qnl_c=w*( p_1m2p(ip)*s_3p4m(ip) - p_3p4m(ip)*s_1m2p(ip) )
356 qnl_d=w*( p_1m2p(ip)*s_3m4p(ip) - p_3m4p(ip)*s_1m2p(ip) )
358 q_apb(ip) = qnl_a + qnl_b
359 q_apc(ip) = qnl_a + qnl_c
360 q_bpd(ip) = qnl_b + qnl_d
361 q_cpd(ip) = qnl_c + qnl_d
366 s_1p2m(ip) = s_1p2m(ip) + s_1m2p(ip)
367 s_3p4m(ip) = s_3p4m(ip) + s_3m4p(ip)
368 p_1p2m(ip) = p_1p2m(ip) + p_1m2p(ip)
369 p_3p4m(ip) = p_3p4m(ip) + p_3m4p(ip)
374 IF (jb4_0.EQ.jf4_0)
THEN 377 tstot(ip,jp4m_0,jf4_0)=tstot(ip,jp4m_0,jf4_0)
379 tstot(ip,jp4m_1,jf4_0)=tstot(ip,jp4m_1,jf4_0)
382 tmp = (-f3plus(ip)*s_1p2m(ip) + p_1p2m(ip))*us1ml4*w
383 tsder(ip,jp4m_0,jf4_0)=tsder(ip,jp4m_0,jf4_0)
385 tsder(ip,jp4m_1,jf4_0)=tsder(ip,jp4m_1,jf4_0)
390 tstot(ip,jp4p_0,jf4_0)=tstot(ip,jp4p_0,jf4_0)
392 tstot(ip,jp4p_1,jf4_0)=tstot(ip,jp4p_1,jf4_0)
395 tmp = (-f3moin(ip)*s_1p2m(ip)+ p_1p2m(ip))*us1ml4*w
397 tsder(ip,jp4p_0,jf4_0)=tsder(ip,jp4p_0,jf4_0)
399 tsder(ip,jp4p_1,jf4_0)=tsder(ip,jp4p_1,jf4_0)
404 IF (jb4_1.EQ.jf4_1)
THEN 407 tstot(ip,jp4m_0,jf4_1)=tstot(ip,jp4m_0,jf4_1)
409 tstot(ip,jp4m_1,jf4_1)=tstot(ip,jp4m_1,jf4_1)
412 tmp = (-f3plus(ip)*s_1p2m(ip) + p_1p2m(ip))*us1ml4*w
413 tsder(ip,jp4m_0,jf4_1)=tsder(ip,jp4m_0,jf4_1)
415 tsder(ip,jp4m_1,jf4_1)=tsder(ip,jp4m_1,jf4_1)
420 tstot(ip,jp4p_0,jf4_1)=tstot(ip,jp4p_0,jf4_1)
422 tstot(ip,jp4p_1,jf4_1)=tstot(ip,jp4p_1,jf4_1)
425 tmp = (-f3moin(ip)*s_1p2m(ip) + p_1p2m(ip))*us1ml4*w
426 tsder(ip,jp4p_0,jf4_1)=tsder(ip,jp4p_0,jf4_1)
428 tsder(ip,jp4p_1,jf4_1)=tsder(ip,jp4p_1,jf4_1)
433 IF (jb2_0.EQ.jf2_0)
THEN 436 tstot(ip,jp2m_0,jf2_0)=tstot(ip,jp2m_0,jf2_0)
438 tstot(ip,jp2m_1,jf2_0)=tstot(ip,jp2m_1,jf2_0)
441 tmp = ( f1plus(ip)*s_3p4m(ip) - p_3p4m(ip))*us1mm4*w
442 tsder(ip,jp2m_0,jf2_0)=tsder(ip,jp2m_0,jf2_0)
444 tsder(ip,jp2m_1,jf2_0)=tsder(ip,jp2m_1,jf2_0)
449 tstot(ip,jp2p_0,jf2_0)=tstot(ip,jp2p_0,jf2_0)
451 tstot(ip,jp2p_1,jf2_0)=tstot(ip,jp2p_1,jf2_0)
454 tmp = ( f1moin(ip)*s_3p4m(ip) - p_3p4m(ip))*us1mm4*w
455 tsder(ip,jp2p_0,jf2_0)=tsder(ip,jp2p_0,jf2_0)
457 tsder(ip,jp2p_1,jf2_0)=tsder(ip,jp2p_1,jf2_0)
462 IF (jb2_1.EQ.jf2_1)
THEN 465 tstot(ip,jp2m_0,jf2_1)=tstot(ip,jp2m_0,jf2_1)
467 tstot(ip,jp2m_1,jf2_1)=tstot(ip,jp2m_1,jf2_1)
470 tmp = ( f1plus(ip)*s_3p4m(ip) - p_3p4m(ip))*us1mm4*w
471 tsder(ip,jp2m_0,jf2_1)=tsder(ip,jp2m_0,jf2_1)
473 tsder(ip,jp2m_1,jf2_1)=tsder(ip,jp2m_1,jf2_1)
478 tstot(ip,jp2p_0,jf2_1)=tstot(ip,jp2p_0,jf2_1)
480 tstot(ip,jp2p_1,jf2_1)=tstot(ip,jp2p_1,jf2_1)
483 tmp = ( f1moin(ip)*s_3p4m(ip) - p_3p4m(ip))*us1mm4*w
484 tsder(ip,jp2p_0,jf2_1)=tsder(ip,jp2p_0,jf2_1)
486 tsder(ip,jp2p_1,jf2_1)=tsder(ip,jp2p_1,jf2_1)
491 IF (jb1_0.EQ.jf1_0)
THEN 494 tstot(ip,jp1p_0,jf1_0)=tstot(ip,jp1p_0,jf1_0)
496 tstot(ip,jp1p_1,jf1_0)=tstot(ip,jp1p_1,jf1_0)
499 tmp = ( f2moin(ip)*s_3p4m(ip) - p_3p4m(ip))*us1pm4*w
500 tsder(ip,jp1p_0,jf1_0)=tsder(ip,jp1p_0,jf1_0)
502 tsder(ip,jp1p_1,jf1_0)=tsder(ip,jp1p_1,jf1_0)
507 tstot(ip,jp1m_0,jf1_0)=tstot(ip,jp1m_0,jf1_0)
509 tstot(ip,jp1m_1,jf1_0)=tstot(ip,jp1m_1,jf1_0)
512 tmp = ( f2plus(ip)*s_3p4m(ip) - p_3p4m(ip))*us1pm4*w
513 tsder(ip,jp1m_0,jf1_0)=tsder(ip,jp1m_0,jf1_0)
515 tsder(ip,jp1m_1,jf1_0)=tsder(ip,jp1m_1,jf1_0)
520 IF (jb1_1.EQ.jf1_1)
THEN 523 tstot(ip,jp1p_0,jf1_1)=tstot(ip,jp1p_0,jf1_1)
525 tstot(ip,jp1p_1,jf1_1)=tstot(ip,jp1p_1,jf1_1)
528 tmp = ( f2moin(ip)*s_3p4m(ip) - p_3p4m(ip))*us1pm4*w
529 tsder(ip,jp1p_0,jf1_1)=tsder(ip,jp1p_0,jf1_1)
531 tsder(ip,jp1p_1,jf1_1)=tsder(ip,jp1p_1,jf1_1)
536 tstot(ip,jp1m_0,jf1_1)=tstot(ip,jp1m_0,jf1_1)
538 tstot(ip,jp1m_1,jf1_1)=tstot(ip,jp1m_1,jf1_1)
541 tmp = ( f2plus(ip)*s_3p4m(ip) - p_3p4m(ip))*us1pm4*w
542 tsder(ip,jp1m_0,jf1_1)=tsder(ip,jp1m_0,jf1_1)
544 tsder(ip,jp1m_1,jf1_1)=tsder(ip,jp1m_1,jf1_1)
549 IF (jb3_0.EQ.jf3_0)
THEN 552 tstot(ip,jp3p_0,jf3_0)=tstot(ip,jp3p_0,jf3_0)
554 tstot(ip,jp3p_1,jf3_0)=tstot(ip,jp3p_1,jf3_0)
557 tmp = (-f4moin(ip)*s_1p2m(ip) + p_1p2m(ip))*us1pl4*w
558 tsder(ip,jp3p_0,jf3_0)=tsder(ip,jp3p_0,jf3_0)
560 tsder(ip,jp3p_1,jf3_0)=tsder(ip,jp3p_1,jf3_0)
565 tstot(ip,jp3m_0,jf3_0)=tstot(ip,jp3m_0,jf3_0)
567 tstot(ip,jp3m_1,jf3_0)=tstot(ip,jp3m_1,jf3_0)
570 tmp = (-f4plus(ip)*s_1p2m(ip) + p_1p2m(ip) )*us1pl4*w
571 tsder(ip,jp3m_0,jf3_0)=tsder(ip,jp3m_0,jf3_0)
573 tsder(ip,jp3m_1,jf3_0)=tsder(ip,jp3m_1,jf3_0)
578 IF (jb3_1.EQ.jf3_1)
THEN 581 tstot(ip,jp3p_0,jf3_1)=tstot(ip,jp3p_0,jf3_1)
583 tstot(ip,jp3p_1,jf3_1)=tstot(ip,jp3p_1,jf3_1)
586 tmp = (-f4moin(ip)*s_1p2m(ip) + p_1p2m(ip))*us1pl4*w
587 tsder(ip,jp3p_0,jf3_1)=tsder(ip,jp3p_0,jf3_1)
589 tsder(ip,jp3p_1,jf3_1)=tsder(ip,jp3p_1,jf3_1)
594 tstot(ip,jp3m_0,jf3_1)=tstot(ip,jp3m_0,jf3_1)
596 tstot(ip,jp3m_1,jf3_1)=tstot(ip,jp3m_1,jf3_1)
599 tmp = (-f4plus(ip)*s_1p2m(ip) + p_1p2m(ip))*us1pl4*w
600 tsder(ip,jp3m_0,jf3_1)=tsder(ip,jp3m_0,jf3_1)
602 tsder(ip,jp3m_1,jf3_1)=tsder(ip,jp3m_1,jf3_1)
subroutine qnlin2(TSTOT, TSDER, IANGNL, COEFNL, NF, NDIRE, NPOIN2, F, XKMOY, TAUX1, DFINI, XCOEF)
double precision, dimension(:), pointer depth
subroutine cqueue(JFRE, JBIS, COEF1)