The TELEMAC-MASCARET system  trunk
prenl3.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE prenl3
3 ! *****************
4 !
5 !
6 !***********************************************************************
7 ! TOMAWAC V6P1 24/06/2011
8 !***********************************************************************
9 !
10 !brief PREPARES THE COMPUTATION FOR THE NON-LINEAR INTERACTION
11 !+ SOURCE TERM BETWEEN QUADRUPLETS USING THE GQM METHOD
12 !+ ("GAUSSIAN QUADRATURE METHOD") PROPOSED BY LAVRENOV
13 !+ (2001)
14 !
15 !+ PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES
16 !+ FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS
17 !+ ARE EVENLY DISTRIBUTED OVER [0;2.PI].
18 !
19 !+note THIS SUBROUTINE IS TO BE USED IN CONJONCT
20 !+ SUBROUTINE QNLIN3, WHICH IT OPTIMISES.
21 !
22 !reference LAVRENOV, I.V. (2001):
23 !+ "EFFECT OF WIND WAVE PARAMETER FLUCTUATION ON THE NONLINEAR
24 !+ SPECTRUM EVOLUTION". J. PHYS. OCEANOGR. 31, 861-873.
25 !
26 !history E. GAGNAIRE-RENOU
27 !+ 04/2011
28 !+ V6P1
29 !+ CREATED
30 !
31 !history G.MATTAROLO (EDF - LNHE)
32 !+ 24/06/2011
33 !+ V6P1
34 !+ Translation of French names of the variables in argument
35 !
36 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
37 !| DIMBUF |-->| VARIABLE FOR SPECTRUM INTERPOLATION
38 !| ELIM |<--| FRACTION OF ELIMINATED CONFIGURATIONS
39 !| F_COEF |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
40 !| F_POIN |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
41 !| F_PROJ |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
42 !| FREQ |-->| DISCRETIZED FREQUENCIES
43 !| IDCONF |<--| WORK TABLE
44 !| IQ_OM1 |<--| NUMBER OF INTEGRATION POINT ON OMEGA1
45 !| K_IF1 |<--| WORK TABLE
46 !| K_IF2 |<--| WORK TABLE
47 !| K_IF3 |<--| WORK TABLE
48 !| K_1M |<--| WORK TABLE
49 !| K_1M2M |<--| WORK TABLE
50 !| K_1M2P |<--| WORK TABLE
51 !| K_1M3M |<--| WORK TABLE
52 !| K_1M3P |<--| WORK TABLE
53 !| K_1P |<--| WORK TABLE
54 !| K_1P2M |<--| WORK TABLE
55 !| K_1P2P |<--| WORK TABLE
56 !| K_1P3M |<--| WORK TABLE
57 !| K_1P3P |<--| WORK TABLE
58 !| LBUF |-->| VARIABLE FOR SPECTRUM INTERPOLATION
59 !| NCONF |<--| NUMBER OF RETAINED CONFIGURATIONS
60 !| NCONFM |-->| MAXIMUM NUMBER OF CONFIGURATIONS
61 !| NF |-->| NUMBER OF FREQUENCIES
62 !| NF1 |-->| NUMBER OF INTEGRATION POINT ON OMEGA1
63 !| NQ_OM2 |-->| NUMBER OF INTEGRATION POINT ON OMEGA2
64 !| NQ_TE1 |-->| SETTING FOR INTEGRATION ON THETA1
65 !| NT |-->| NUMBER OF DIRECTIONS
66 !| NT1 |-->| NUMBER OF INTEGRATION POINT ON THETA1
67 !| RAISF |-->| FREQUENTIAL RATIO
68 !| SEUIL1 |-->| THRESHOLD1 FOR CONFIGURATIONS ELIMINATION
69 !| SEUIL2 |-->| THRESHOLD2 FOR CONFIGURATIONS ELIMINATION
70 !| T_POIN |<--| WORK TABLE FOR SPECTRUM INTERPOLATION
71 !| TAILF |-->| SPECTRUM QUEUE FACTOR
72 !| TB_FAC |<--| WORK TABLE
73 !| TB_SCA |<--| SCALE COEFFICIENT
74 !| TB_TMP |<--| WORK TABLE
75 !| TB_TPM |<--| WORK TABLE
76 !| TB_V14 |<--| WORK TABLE for F1
77 !| TB_V24 |<--| WORK TABLE
78 !| TB_V34 |<--| WORK TABLE
79 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
80 
81 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
82 !
83 !Missing the description of the variables in argument
84 !
85 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
86 !
88 !
90  USE interface_tomawac, ex_prenl3 => prenl3
91  IMPLICIT NONE
92 !
93 !
94 !
95 !.....LOCAL VARIABLES
96 ! """""""""""""""""
97  INTEGER JF , JT , JF1 , JT1 , NF1P1 , IAUX
98  INTEGER IQ_TE1, IQ_OM2
99  DOUBLE PRECISION EPSI_A, AUX , CCC , DENO , AAA , DP2SG
100  DOUBLE PRECISION V1 , V1_4 , DV1 , DTETAR
101  DOUBLE PRECISION V2 , V2_4 , V3 , V3_4
102  DOUBLE PRECISION W2 , W2_M , W2_1 , W_MIL , W_RAD
103  DOUBLE PRECISION RK0 , XK0 , YK0 , RK1 , XK1 , YK1
104  DOUBLE PRECISION RK2 , XK2P , YK2P , XK2M , YK2M
105  DOUBLE PRECISION RK3 , XK3P , YK3P , XK3M , YK3M
106  DOUBLE PRECISION D01P , C_D01P, S_D01P, D0AP , C_D0AP, S_D0AP
107  DOUBLE PRECISION GA2P , C_GA2P, S_GA2P, GA3P , C_GA3P, S_GA3P
108  DOUBLE PRECISION, ALLOCATABLE :: F1SF(:)
109 !
110 !.....Variables related to the Gaussian quadratures
111  DOUBLE PRECISION W_CHE_TE1, W_CHE_OM2, C_LEG_OM2
112  DOUBLE PRECISION, ALLOCATABLE :: X_CHE_TE1(:),X_CHE_OM2(:),
113  & x_leg_om2(:),w_leg_om2(:)
114 !
115 !.....Variables related to the configuration selection
116  DOUBLE PRECISION TEST1 , TEST2
117  DOUBLE PRECISION, ALLOCATABLE :: MAXCLA(:)
118  INTEGER NT
119 !
120 !.....External functions (USE INTERFACE)
121 ! DOUBLE PRECISION COUPLE
122 ! EXTERNAL COUPLE
123 !
124  nt=ndire
125  dtetar=deupi/dble(nt)
126 !
127 !=======================================================================
128 ! INITIALISATION OF AUXILIAIRY TABLES FOR SPECTRUM INTERPOLATION
129 !=======================================================================
130  DO jf=1,dimbuf
131  f_poin(jf)=0
132  t_poin(jf)=0
133  f_coef(jf)=0.d0
134  f_proj(jf)=0.d0
135  tb_sca(jf)=0.0d0
136  ENDDO
137  DO jf=1,lbuf
138  f_poin(jf)=1
139  f_coef(jf)=0.0d0
140  f_proj(jf)=0.0d0
141  ENDDO
142  DO jf=1,nf
143  iaux=lbuf+jf
144  f_poin(iaux)=jf
145  f_coef(iaux)=1.0d0
146  f_proj(iaux)=1.0d0
147  ENDDO
148  aux=1.d0/raisf**tailf
149  DO jf=1,lbuf
150  iaux=lbuf+nf+jf
151  f_poin(iaux)=nf
152  f_coef(iaux)=aux**jf
153  f_proj(iaux)=0.0d0
154  ENDDO
155 !
156  DO jt=lbuf,1,-1
157  t_poin(jt)=nt-mod(lbuf-jt,nt)
158  ENDDO
159  DO jt=1,nt
160  t_poin(lbuf+jt)=jt
161  ENDDO
162  DO jt=1,lbuf
163  t_poin(lbuf+nt+jt)=mod(jt-1,nt)+1
164  ENDDO
165 !======================================================================
166 !
167 !=======================================================================
168 ! COMPUTES SCALE COEFFICIENTS FOR THE COUPLING COEFFICIENT
169 !=======================================================================
170  dp2sg=deupi*deupi/gravit
171  DO jf=1,lbuf
172  aux=freq(1)/raisf**(lbuf-jf+1)
173  tb_sca(jf)=(dp2sg*aux**2)**6/(deupi**3*aux)
174  ENDDO
175  DO jf=1,nf
176  tb_sca(lbuf+jf)=(dp2sg*freq(jf)**2)**6/(deupi**3*freq(jf))
177  ENDDO
178  DO jf=1,lbuf
179  iaux=lbuf+nf+jf
180  aux=freq(nf)*raisf**jf
181  tb_sca(iaux)=(dp2sg*aux**2)**6/(deupi**3*aux)
182  ENDDO
183 !=======================================================================
184 !
185 !=======================================================================
186 ! COMPUTES VALUES FOR GAUSSIAN QUADRATURES
187 !=======================================================================
188  ALLOCATE(x_che_te1(1:nq_te1),x_che_om2(1:nq_om2))
189  ALLOCATE(x_leg_om2(1:nq_om2),w_leg_om2(1:nq_om2))
190 !
191 !.....Abscissa and weight (constant) for Gauss-Chebyshev
192  DO iq_te1=1,nq_te1
193  x_che_te1(iq_te1)=cos(pi*(dble(iq_te1)-0.5d0)/dble(nq_te1))
194  ENDDO
195  w_che_te1=pi/dble(nq_te1)
196  DO iq_om2=1,nq_om2
197  x_che_om2(iq_om2)=cos(pi*(dble(iq_om2)-0.5d0)/dble(nq_om2))
198  ENDDO
199  w_che_om2=pi/dble(nq_om2)
200 !
201 !.....Abscissa et weight for Gauss-Legendre
202  CALL gauleg
203  &( w_leg_om2 , x_leg_om2 , nq_om2 )
204  DO iq_om2=1,nq_om2
205  x_leg_om2(iq_om2)=0.25d0*(1.d0+x_leg_om2(iq_om2))**2
206  ENDDO
207 !=======================================================================
208 !
209 !
210 !=======================================================================
211 ! COMPUTES VALUES OF RATIO F1/F AS FUNCTION OF THE IQ_OM1 INDICATOR
212 !=======================================================================
213  nf1p1=nf1+1
214  ALLOCATE(f1sf(1:nf1p1))
215  CALL f1f1f1 ( f1sf , nf1 , iq_om1)
216 !=======================================================================
217 !
218 ! ==================================================
219 ! STARTS LOOP 1 OVER THE RATIOS F1/F0
220 ! ==================================================
221  DO jf1=1,nf1
222 ! ---------Computes and stores v1=f1/f0 and v1**4
223  v1=(f1sf(jf1+1)+f1sf(jf1))/2.d0
224  k_if1(jf1)=nint(dble(lbuf)+log(v1)/log(raisf))
225  v1_4=v1**4
226  tb_v14(jf1)=v1_4
227 ! ---------Computes and stores dv1=df1/f0
228  dv1=f1sf(jf1+1)-f1sf(jf1)
229 ! ---------Computes the A parameter
230  aaa=((1.d0+v1)**4-4.d0*(1.d0+v1_4))/(8.d0*v1**2)
231 !
232 !
233 !
234 ! =================================================
235 ! STARTS LOOP 2 OVER THE DELTA_1+ VALUES
236 ! =================================================
237  DO jt1=1,nt1
238 !
239 !......Computes the Delta1+ values (=Theta_1-Theta_0) between 0 and Pi.
240  IF (jt1.LE.nq_te1) THEN
241 ! ---------First interval : X from -1 to A
242  iq_te1=jt1
243  c_d01p=(-1.d0+aaa)/2.d0+(1.d0+aaa)/2.d0*x_che_te1(iq_te1)
244  ccc=dv1*sqrt((aaa-c_d01p)/(1.d0-c_d01p))*w_che_te1
245  ELSE
246 ! ---------Second interval : X from A to 1
247  iq_te1=jt1-nq_te1
248  c_d01p=( 1.d0+aaa)/2.d0+(1.d0-aaa)/2.d0*x_che_te1(iq_te1)
249  ccc=dv1*sqrt((c_d01p-aaa)/(1.d0+c_d01p))*w_che_te1
250  ENDIF
251  s_d01p=sqrt(1.d0-c_d01p*c_d01p)
252  d01p =acos(c_d01p)
253  k_1p(jt1,jf1)=lbuf+nint(d01p/dtetar)
254  k_1m(jt1,jf1)=lbuf-nint(d01p/dtetar)
255 !
256 ! ---------Computes Epsilon_a
257  epsi_a=2.d0*sqrt(1.d0+v1_4+2.d0*v1*v1*c_d01p)/(1.d0+v1)**2
258 ! ---------Computes Delta_A+ and its cosinus
259  c_d0ap=(1.d0-v1_4+0.25d0*epsi_a**2*(1.d0+v1)**4)
260  & /(epsi_a*(1.d0+v1)**2)
261  s_d0ap=sqrt(1.0d0-c_d0ap*c_d0ap)
262  d0ap = acos(c_d0ap)
263 !
264 !.......Integration over OMEGA2 depending on EPS_A
265  IF (epsi_a.LT.1.d0) THEN
266 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267 !........Case of a single singularity (in OMEGA2-)
268 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269  w2_m=0.5d0*(1.d0-epsi_a/2.d0)
270  w2_1=0.5d0
271 !
272  w_rad=w2_1-w2_m
273  c_leg_om2=sqrt(w_rad)
274 !
275 ! ----------------------------------------------------
276 !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
277 !........Case of a single singularity (in OMEGA2-)
278 !........Integration over OMEGA2 via GAUSS-LEGENDRE quadrature
279 ! ----------------------------------------------------
280  DO iq_om2=1,nq_om2
281 ! ---------Computes W2, V2, and V3
282  w2=w2_m+w_rad*x_leg_om2(iq_om2)
283  v2=w2*(1.d0+v1)
284  v2_4=v2**4
285  tb_v24(iq_om2,jt1,jf1)=v2_4
286  k_if2(iq_om2,jt1,jf1) = nint(dble(lbuf)
287  & + log(v2)/log(raisf))
288  v3=1.d0+v1-v2
289  v3_4=v3**4
290  tb_v34(iq_om2,jt1,jf1)=v3_4
291  k_if3(iq_om2,jt1,jf1) = nint(dble(lbuf)
292  & + log(v3)/log(raisf))
293 ! ---------Computes Gamma_2+ et Gamma_3+ angles
294  c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
295  c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
296  s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
297  ga2p =acos(c_ga2p)
298  c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a
299  & /(1.d0-w2)**2
300  c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
301  s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
302  ga3p =acos(c_ga3p)
303 ! Shifting of the direction indexes - Config. +Delta1 (SIG=1)
304  k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar
305  & +dble(lbuf))
306  k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar
307  & +dble(lbuf))
308  k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar
309  & +dble(lbuf))
310  k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar
311  & +dble(lbuf))
312 ! Shifting of the direction indexes - Config. -Delta1 (SIG=-1)
313  k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar
314  & +dble(lbuf))
315  k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar
316  & +dble(lbuf))
317  k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar
318  & +dble(lbuf))
319  k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar
320  & +dble(lbuf))
321 !
322 !.........Computes the coupling coefficients (only for Delta_1+ )
323  rk0=1.d0
324  rk1=v1*v1
325  rk2=v2*v2
326  rk3=(1.d0+v1-v2)**2
327  xk0 = rk0
328  yk0 = 0.0d0
329  xk1 = rk1*c_d01p
330  yk1 = rk1*s_d01p
331  xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
332  yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
333  xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
334  yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
335  xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
336  yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
337  xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
338  yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
339  tb_tpm(iq_om2,jt1,jf1)=couple
340  &( xk0 , yk0 , xk1 , yk1 , xk2p , yk2p , xk3m , yk3m )
341  tb_tmp(iq_om2,jt1,jf1)=couple
342  &( xk0 , yk0 , xk1 , yk1 , xk2m , yk2m , xk3p , yk3p )
343 !
344 !.........Computes the multiplicative coefficient for QNL4
345  deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2)
346  & *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a))
347  & *((w2-0.5d0)**2+0.25d0*(1.d0-epsi_a)) )
348  tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2))
349  & /(1.d0+v1)**5 * w_leg_om2(iq_om2)*c_leg_om2* ccc
350  ENDDO
351 ! -----------------------------------------------
352 !........END OF THE LOOP 3 OVER OMEGA_2 (CASE Epsilon_A < 1)
353 ! -----------------------------------------------
354 !
355  ELSE
356 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357 !........STARTS LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
358 !........Case of two singularities (in OMEGA2- and OMEGA2_1)
359 !........Integration over OMEGA2 via GAUSS-CHEBYSCHEV quadrature
360 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - -
361  w2_m=0.5d0*(1.d0-epsi_a/2.d0)
362  w2_1=0.5d0*(1.d0-sqrt(epsi_a-1.d0))
363 !
364  w_mil=(w2_m+w2_1)/2.d0
365  w_rad=(w2_1-w2_m)/2.d0
366 !
367  DO iq_om2=1,nq_om2
368 ! ---------Computes W2, V2, and V3
369  w2=w_mil+w_rad*x_che_om2(iq_om2)
370  v2=w2*(1.d0+v1)
371  v2_4=v2**4
372  tb_v24(iq_om2,jt1,jf1)=v2_4
373  k_if2(iq_om2,jt1,jf1)=nint(dble(lbuf)
374  & +log(v2)/log(raisf))
375  v3=1.d0+v1-v2
376  v3_4=v3**4
377  tb_v34(iq_om2,jt1,jf1)=v3_4
378  k_if3(iq_om2,jt1,jf1)=nint(dble(lbuf)
379  & +log(v3)/log(raisf))
380 ! ---------Computes Gamma_2+ et Gamma_3+ angles
381  c_ga2p=(epsi_a**2/4.d0+w2**4-(1.d0-w2)**4)/(epsi_a*w2*w2)
382  c_ga2p=max(min(c_ga2p,1.d0),-1.d0)
383  s_ga2p=sqrt(1.d0-c_ga2p*c_ga2p)
384  ga2p =acos(c_ga2p)
385  c_ga3p=(epsi_a**2/4.d0-w2**4+(1.d0-w2)**4)/epsi_a
386  & /(1.d0-w2)**2
387  c_ga3p=max(min(c_ga3p,1.d0),-1.d0)
388  s_ga3p=sqrt(1.d0-c_ga3p*c_ga3p)
389  ga3p =acos(c_ga3p)
390 ! Shifts the direction indexes - Config. +Delta1 (SIG=1)
391  k_1p2p(iq_om2,jt1,jf1)=nint(( d0ap+ga2p)/dtetar
392  & +dble(lbuf))
393  k_1p3m(iq_om2,jt1,jf1)=nint(( d0ap-ga3p)/dtetar
394  & +dble(lbuf))
395  k_1p2m(iq_om2,jt1,jf1)=nint(( d0ap-ga2p)/dtetar
396  & +dble(lbuf))
397  k_1p3p(iq_om2,jt1,jf1)=nint(( d0ap+ga3p)/dtetar
398  & +dble(lbuf))
399 ! Shifts the direction indexes - Config. -Delta1 (SIG=-1)
400  k_1m2p(iq_om2,jt1,jf1)=nint((-d0ap+ga2p)/dtetar
401  & +dble(lbuf))
402  k_1m3m(iq_om2,jt1,jf1)=nint((-d0ap-ga3p)/dtetar
403  & +dble(lbuf))
404  k_1m2m(iq_om2,jt1,jf1)=nint((-d0ap-ga2p)/dtetar
405  & +dble(lbuf))
406  k_1m3p(iq_om2,jt1,jf1)=nint((-d0ap+ga3p)/dtetar
407  & +dble(lbuf))
408 !
409 !.........Computes the coupling coefficients (only for Delta_1+ )
410  rk0=1.d0
411  rk1=v1*v1
412  rk2=v2*v2
413  rk3=(1.d0+v1-v2)**2
414  xk0 = rk0
415  yk0 = 0.0d0
416  xk1 = rk1*c_d01p
417  yk1 = rk1*s_d01p
418  xk2p = rk2*(c_d0ap*c_ga2p-s_d0ap*s_ga2p)
419  yk2p = rk2*(s_d0ap*c_ga2p+c_d0ap*s_ga2p)
420  xk2m = rk2*(c_d0ap*c_ga2p+s_d0ap*s_ga2p)
421  yk2m = rk2*(s_d0ap*c_ga2p-c_d0ap*s_ga2p)
422  xk3p = rk3*(c_d0ap*c_ga3p-s_d0ap*s_ga3p)
423  yk3p = rk3*(s_d0ap*c_ga3p+c_d0ap*s_ga3p)
424  xk3m = rk3*(c_d0ap*c_ga3p+s_d0ap*s_ga3p)
425  yk3m = rk3*(s_d0ap*c_ga3p-c_d0ap*s_ga3p)
426  tb_tpm(iq_om2,jt1,jf1)=couple
427  &( xk0 , yk0 , xk1 , yk1 , xk2p , yk2p , xk3m , yk3m )
428  tb_tmp(iq_om2,jt1,jf1)=couple
429  &( xk0 , yk0 , xk1 , yk1 , xk2m , yk2m , xk3p , yk3p )
430 !
431 !.........Computes the multiplicative coefficient for QNL4
432  deno=2.d0*sqrt( (0.5d0*(1.d0+epsi_a/2.d0)-w2)
433  & *((w2-0.5d0)**2+0.25d0*(1.d0+epsi_a))
434  & *(0.5d0*(1.d0+sqrt(epsi_a-1.d0))-w2) )
435  tb_fac(iq_om2,jt1,jf1)=1.d0/(deno*v1*w2*(1.d0-w2))
436  & /(1.d0+v1)**5 * w_che_om2* ccc
437 !
438  ENDDO
439 ! -----------------------------------------------
440 !........END OF LOOP 3 OVER OMEGA_2 (CASE Epsilon_A > 1)
441 ! -----------------------------------------------
442 !
443  ENDIF
444 !
445 !
446  ENDDO
447 ! =================================================
448 ! END OF LOOP 2 OVER THE DELTA_1+ VALUES
449 ! =================================================
450 !
451  ENDDO
452 ! ==================================================
453 ! END OF LOOP 1 OVER THE F1/F0 RATIOS
454 ! ==================================================
455  DEALLOCATE(f1sf)
456  DEALLOCATE(x_che_te1)
457  DEALLOCATE(x_che_om2)
458  DEALLOCATE(x_leg_om2)
459  DEALLOCATE(w_leg_om2)
460 !
461 !
462 !
463 !
464 ! ===========================================================
465 ! POST-PROCESSING TO ELIMINATE PART OF THE CONFIGURATIONS
466 ! ===========================================================
467 !
468 !.....It looks, for every value of the ratio V1, for the maximum value
469 !.....of FACTOR*COUPLING : it is stored in the local table NAXCLA(.)
470 ! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
471  ALLOCATE(maxcla(1:nf1))
472  DO jf1=1,nf1
473  aux=0.0d0
474  DO jt1=1,nt1
475  DO iq_om2=1,nq_om2
476  aaa=tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1)
477  IF (aaa.GT.aux) aux=aaa
478  ccc=tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1)
479  IF (ccc.GT.aux) aux=ccc
480  ENDDO
481  ENDDO
482  maxcla(jf1)=aux
483  ENDDO
484 !
485 !.....It looks for the max V1 value
486 ! """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""
487  aux=0.0d0
488  DO jf1=1,nf1
489  IF (maxcla(jf1).GT.aux) aux=maxcla(jf1)
490  ENDDO
491  test1=seuil1*aux
492 !
493 !.....Set to zero the coupling coefficients not used
494 ! """""""""""""""""""""""""""""""""""""""""""""""""""""
495  nconf=0
496  DO jf1=1,nf1
497  test2 =seuil2*maxcla(jf1)
498  DO jt1=1,nt1
499  DO iq_om2=1,nq_om2
500  aaa=tb_fac(iq_om2,jt1,jf1)*tb_tpm(iq_om2,jt1,jf1)
501  ccc=tb_fac(iq_om2,jt1,jf1)*tb_tmp(iq_om2,jt1,jf1)
502  IF ((aaa.GT.test1.OR.aaa.GT.test2).OR.
503  & (ccc.GT.test1.OR.ccc.GT.test2)) THEN
504  nconf=nconf+1
505  idconf(nconf,1)=jf1
506  idconf(nconf,2)=jt1
507  idconf(nconf,3)=iq_om2
508  ENDIF
509  ENDDO
510  ENDDO
511  ENDDO
512  DEALLOCATE(maxcla)
513 !
514 !.....It counts the fraction of the eliminated configurations
515 ! """"""""""""""""""""""""""""""""""""""""""""""""""""""
516  elim=(1.d0-dble(nconf)/dble(nconfm))*100.d0
517  WRITE(lu,*) 'PART DE CONFIGURATIONS ANNULEES : ',elim,' %'
518 !
519  RETURN
520  END
integer, dimension(:,:), allocatable k_1m
double precision, dimension(:,:,:), allocatable tb_fac
double precision, dimension(:,:,:), allocatable tb_tpm
double precision, dimension(:), allocatable tb_v14
integer, dimension(:,:,:), allocatable k_1m2p
integer, dimension(dimbuf) f_poin
double precision, dimension(:), pointer freq
integer, dimension(:,:), allocatable idconf
integer, dimension(:,:,:), allocatable k_if2
integer, dimension(:,:,:), allocatable k_1m3p
double precision, dimension(:,:,:), allocatable tb_v34
integer, dimension(:,:,:), allocatable k_1p3p
integer, dimension(:,:,:), allocatable k_1m3m
double precision, dimension(:,:,:), allocatable tb_tmp
double precision, dimension(dimbuf) f_proj
double precision function couple(XK1, YK1, XK2, YK2, XK3, YK3, XK4, YK4)
Definition: couple.f:7
integer, dimension(:,:,:), allocatable k_1p2p
subroutine f1f1f1(F1SF, NF1, IQ_OM1)
Definition: f1f1f1.f:7
integer, dimension(:,:,:), allocatable k_1m2m
double precision, dimension(dimbuf) f_coef
integer, dimension(:), allocatable k_if1
integer, dimension(:,:,:), allocatable k_1p2m
integer, parameter dimbuf
integer, dimension(:,:), allocatable k_1p
double precision, dimension(dimbuf) tb_sca
subroutine prenl3
Definition: prenl3.f:4
integer, dimension(:,:,:), allocatable k_if3
subroutine gauleg(W_LEG, X_LEG, NPOIN)
Definition: gauleg.f:7
integer, dimension(:,:,:), allocatable k_1p3m
integer, dimension(dimbuf) t_poin
double precision, dimension(:,:,:), allocatable tb_v24
integer, parameter lbuf