The TELEMAC-MASCARET system  trunk
qnlin3.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE qnlin3
3 ! *****************
4 !
5  &( t1tot,t1der, f2, n1poin2, n1plan, n1f )
6 !
7 !***********************************************************************
8 ! TOMAWAC V6P1 24/06/2011
9 !***********************************************************************
10 !
11 !brief COMPUTES THE CONTRIBUTION OF THE NON-LINEAR INTERACTIONS
12 !+ SOURCE TERM BETWEEN QUADRUPLETS USING THE GQM METHOD
13 !+ ("GAUSSIAN QUADRATURE METHOD") PROPOSED BY LAVRENOV
14 !+ (2001)
15 !+
16 !+ PROCEDURE SPECIFIC TO THE CASE WHERE THE FREQUENCIES
17 !+ FOLLOW A GEOMETRICAL PROGRESSION AND THE DIRECTIONS
18 !+ ARE EVENLY DISTRIBUTED OVER [0;2.PI].
19 !
20 !note THIS SUBROUTINE USES THE OUTPUT FROM 'PRENL3' TO OPTIMISE
21 !+ THE COMPUTATIONS FOR DIA.
22 !
23 !reference LAVRENOV, I.V. (2001):
24 !+ "EFFECT OF WIND WAVE PARAMETER FLUCTUATION ON THE NONLINEAR
25 !+ SPECTRUM EVOLUTION". J. PHYS. OCEANOGR. 31, 861-873.
26 !
27 !history E. GAGNAIRE-RENOU
28 !+ 04/2011
29 !+ V6P1
30 !+ CREATED
31 !
32 !history G.MATTAROLO (EDF - LNHE)
33 !+ 24/06/2011
34 !+ V6P1
35 !+ Translation of French names of the variables in argument
36 !
37 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38 !| DIMBUF |-->| VARIABLE FOR SPECTRUM INTERPOLATION
39 !| F |-->| DIRECTIONAL SPECTRUM
40 !| F_COEF |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
41 !| F_POIN |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
42 !| F_PROJ |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
43 !| FSEUIL |-->| WORK TABLE
44 !| K_IF1 |-->| WORK TABLE
45 !| K_IF2 |-->| WORK TABLE
46 !| K_IF3 |-->| WORK TABLE
47 !| K_1M |-->| WORK TABLE
48 !| K_1M2M |-->| WORK TABLE
49 !| K_1M2P |-->| WORK TABLE
50 !| K_1M3M |-->| WORK TABLE
51 !| K_1M3P |-->| WORK TABLE
52 !| K_1P |-->| WORK TABLE
53 !| K_1P2M |-->| WORK TABLE
54 !| K_1P2P |-->| WORK TABLE
55 !| K_1P3M |-->| WORK TABLE
56 !| K_1P3P |-->| WORK TABLE
57 !| IDCONF |-->| WORK TABLE
58 !| LBUF |-->| VARIABLE FOR SPECTRUM INTERPOLATION
59 !| NB_NOD |-->| NUMBER OF POINTS IN 2D MESH
60 !| NCONF |-->| NUMBER OF RETAINED CONFIGURATIONS
61 !| NCONFM |-->| MAXIMUM NUMBER OF CONFIGURATIONS
62 !| NF |-->| NUMBER OF FREQUENCIES
63 !| NF1 |-->| NUMBER OF INTEGRATION POINT ON OMEGA1
64 !| NQ_OM2 |-->| NUMBER OF INTEGRATION POINT ON OMEGA2
65 !| NT |-->| NUMBER OF DIRECTIONS
66 !| NT1 |-->| NUMBER OF INTEGRATION POINT ON TETA1
67 !| RAISF |-->| FREQUENTIAL RATIO
68 !| SEUIL |-->| THRESHOLD0 FOR CONFIGURATIONS ELIMINATION (GQM)
69 !| T_POIN |-->| WORK TABLE FOR SPECTRUM INTERPOLATION
70 !| TB_FAC |-->| WORK TABLE
71 !| TB_SCA |-->| SCALE COEFFICIENT
72 !| TB_TMP |-->| WORK TABLE
73 !| TB_TPM |-->| WORK TABLE
74 !| TB_V14 |-->| WORK TABLE
75 !| TB_V24 |-->| WORK TABLE
76 !| TB_V34 |-->| WORK TABLE
77 !| TSDER |<->| DERIVED PART OF THE SOURCE TERM CONTRIBUTION
78 !| TSTOT |<->| TOTAL PART OF THE SOURCE TERM CONTRIBUTION
79 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
80 !
81  USE interface_tomawac, ex_qnlin3 => qnlin3
83 !
84  IMPLICIT NONE
85 !
86 !.....VARIABLES IN ARGUMENT
87 ! """"""""""""""""""""
88 ! A WAY TO CHANGE DIMENSION OF ARRAY
89 ! TSTOT,TSDER AND F HAVE ONLY ONE DIMENSION
90 
91  INTEGER, INTENT(IN) :: N1POIN2,N1PLAN,N1F
92  DOUBLE PRECISION, INTENT(INOUT) :: T1TOT( n1poin2,n1plan,n1f)
93  DOUBLE PRECISION, INTENT(INOUT) :: T1DER( n1poin2,n1plan,n1f)
94  DOUBLE PRECISION, INTENT(INOUT) :: F2(n1poin2,n1plan,n1f)
95 !
96 !.....LOCAL VARIABLES
97 ! """""""""""""""""
98  INTEGER IP , JF , JT , JF1 , JT1 , IQ_OM2,
99  & jfm0 , jfm1 , jfm2 , jfm3 , ixf1 , ixf2 ,
100  & ixf3 , jfmin , jfmax , iconf
101  INTEGER KT1P , KT1M , JT1P , JT1M , KT1P2P, KT1P2M,
102  & kt1p3p, kt1p3m, kt1m2p, kt1m2m, kt1m3p, kt1m3m,
103  & jt1p2p, jt1p2m, jt1p3p, jt1p3m, jt1m2p, jt1m2m,
104  & jt1m3p, jt1m3m
105  DOUBLE PRECISION V1_4 , V2_4 , V3_4 , Q_2P3M, Q_2M3P, FACTOR,
106  & t_2p3m, t_2m3p, s_2p3m, s_2m3p, scal_t, t2p3m ,
107  & t2m3p , sp0 , sp1p , sp1m , sp1p2p, sp1p2m,
108  & sp1p3p, sp1p3m, sp1m2p, sp1m2m, sp1m3p, sp1m3m,
109  & cf0 , cp0 , cf1 , cp1 , cf2 , cp2 ,
110  & cf3 , cp3 , q2pd0 , q2pd1 , q2pd2p, q2pd3m,
111  & q2md0 , q2md1 , q2md2m, q2md3p,
112  & aux00 , aux01 , aux02 , aux03 , aux04 , aux05 ,
113  & aux06 , aux07 , aux08 , aux09 , aux10
114 !
115 !=======================================================================
116 ! COMPUTES THE GENERALIZED MIN AND MAX FREQUENCIES : INSTEAD OF GOING
117 ! FROM 1 TO NF IN FREQ(JF) FOR THE MAIN FREQUENCY, IT GOES FROM JFMIN
118 ! TO JFMAX
119 ! JFMIN IS GIVEN BY Fmin=FREQ(1) /Gamma_min
120 ! JFMAX IS GIVEN BY Fmax=FREQ(NF)*Gamma_max
121 ! TESTS HAVE SHOWN THAT IT CAN BE ASSUMED Gamma_min=1. (JFMIN=1) AND
122 ! Gamma_max=1.3 (JFMAX>NF) TO OBTAIN IMPROVED RESULTS
123 !=======================================================================
124  jfmin= 1-int(log(1.0d0)/log(raisf))
125  jfmax=nf+int(log(1.3d0)/log(raisf))
126 !
127 !=======================================================================
128 ! COMPUTES THE SPECTRUM THRESHOLD VALUES (BELOW WHICH QNL4 IS NOT
129 ! CALCULATED). THE THRESHOLD IS SET WITHIN 0 AND 1.
130 !=======================================================================
131  DO ip=1,npoin2
132  aux00=0.0d0
133  DO jf=1,nf
134  DO jt=1,ndire
135  IF (f2(ip,jt,jf).GT.aux00) aux00=f2(ip,jt,jf)
136  ENDDO
137  ENDDO
138  tra37(ip)=aux00*seuil
139  ENDDO
140 !=======================================================================
141 !
142 !
143 !
144 !
145 ! ==================================================
146 ! STARTS LOOP 1 OVER THE SELECTED CONFIGURATIONS
147 ! ==================================================
148  DO iconf=1,nconf
149 ! ---------selected configuration characteristics
150  jf1 =idconf(iconf,1)
151  jt1 =idconf(iconf,2)
152  iq_om2=idconf(iconf,3)
153 !
154 ! ---------Recovers V1**4=(f1/f0)**4
155  v1_4 =tb_v14(jf1)
156 ! ---------Recovers the shift of the frequency index on f1
157  ixf1 =k_if1(jf1)
158 ! ---------Recovers the direction indexes for Delat1
159  kt1p =k_1p(jt1,jf1)
160  kt1m =k_1m(jt1,jf1)
161 ! ---------Recovers V2**4=(f2/f0)**4 and V3**4=(f3/f0)**4
162  v2_4 =tb_v24(iq_om2,jt1,jf1)
163  v3_4 =tb_v34(iq_om2,jt1,jf1)
164 ! ---------Recovers the frequency indexes shift on f2 and f3
165  ixf2 =k_if2(iq_om2,jt1,jf1)
166  ixf3 =k_if3(iq_om2,jt1,jf1)
167 ! ---------Recovers the direction indexes shift
168  kt1p2p=k_1p2p(iq_om2,jt1,jf1)
169  kt1p2m=k_1p2m(iq_om2,jt1,jf1)
170  kt1p3p=k_1p3p(iq_om2,jt1,jf1)
171  kt1p3m=k_1p3m(iq_om2,jt1,jf1)
172  kt1m2p=k_1m2p(iq_om2,jt1,jf1)
173  kt1m2m=k_1m2m(iq_om2,jt1,jf1)
174  kt1m3p=k_1m3p(iq_om2,jt1,jf1)
175  kt1m3m=k_1m3m(iq_om2,jt1,jf1)
176 ! ---------Recovers the coupling coefficients
177  t2p3m =tb_tpm(iq_om2,jt1,jf1)
178  t2m3p =tb_tmp(iq_om2,jt1,jf1)
179 ! ---------Recovers the multiplicative factor of QNL4
180  factor=tb_fac(iq_om2,jt1,jf1)
181 !
182 ! = = = = = = = = = = = = = = = = = = = = = = = = =
183 ! STARTS LOOP 2 OVER THE SPECTRUM FREQUENCIES
184 ! = = = = = = = = = = = = = = = = = = = = = = = = =
185  DO jf=jfmin,jfmax
186 !
187 !.........Recovers the coefficient for the coupling factor
188 !.........Computes the coupling coefficients for the case +Delta1 (SIG=1)
189  scal_t=tb_sca(lbuf+jf)*factor
190  t_2p3m=t2p3m*scal_t
191  t_2m3p=t2m3p*scal_t
192 !
193 !.........Frequency indexes and coefficients
194  jfm0=f_poin(jf+lbuf)
195  cf0 =f_coef(jf+lbuf)
196  cp0 =f_proj(jf+lbuf)
197  jfm1=f_poin(jf+ixf1)
198  cf1 =f_coef(jf+ixf1)
199  cp1 =f_proj(jf+ixf1)
200  jfm2=f_poin(jf+ixf2)
201  cf2 =f_coef(jf+ixf2)
202  cp2 =f_proj(jf+ixf2)
203  jfm3=f_poin(jf+ixf3)
204  cf3 =f_coef(jf+ixf3)
205  cp3 =f_proj(jf+ixf3)
206 !
207 ! -------------------------------------------------
208 ! STARTS LOOP 3 OVER THE SPECTRUM DIRECTIONS
209 ! -------------------------------------------------
210  DO jt=1,ndire
211 !
212 !...........Direction indexes
213 ! direct config (+delta1) (sig =1)
214  jt1p =t_poin(jt+kt1p)
215  jt1p2p=t_poin(jt+kt1p2p)
216  jt1p2m=t_poin(jt+kt1p2m)
217  jt1p3p=t_poin(jt+kt1p3p)
218  jt1p3m=t_poin(jt+kt1p3m)
219 ! image config (-delta1)
220  jt1m =t_poin(jt+kt1m)
221  jt1m2p=t_poin(jt+kt1m2p)
222  jt1m2m=t_poin(jt+kt1m2m)
223  jt1m3p=t_poin(jt+kt1m3p)
224  jt1m3m=t_poin(jt+kt1m3m)
225 !
226 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227 ! STARTS LOOP 4 OVER THE MESH NODES
228 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229  DO ip=1,npoin2
230 !
231  sp0=f2(ip,jt,jfm0)*cf0
232 !
233  IF (sp0.GT.tra37(ip)) THEN
234 !
235 ! Config. +Delta1 (SIG=1)
236 ! =======================
237 !...............Computes the spectrum values in 1, 2, 3
238  sp1p =f2(ip,jt1p ,jfm1)*cf1
239  sp1p2p=f2(ip,jt1p2p,jfm2)*cf2
240  sp1p3m=f2(ip,jt1p3m,jfm3)*cf3
241  sp1p2m=f2(ip,jt1p2m,jfm2)*cf2
242  sp1p3p=f2(ip,jt1p3p,jfm3)*cf3
243 !
244 !...............Computes auxiliary products and variables
245  aux01=sp0*v1_4+sp1p
246  aux02=sp0*sp1p
247  aux03=sp1p2p*sp1p3m
248  aux04=sp1p2p*v3_4+sp1p3m*v2_4
249  aux05=sp1p2m*sp1p3p
250  aux06=sp1p2m*v3_4+sp1p3p*v2_4
251  aux07=aux02*v3_4
252  aux08=aux02*v2_4
253 !
254 !...............Computes the components of the transfer term
255  s_2p3m=aux03*aux01-aux02*aux04
256  s_2m3p=aux05*aux01-aux02*aux06
257  q_2p3m=t_2p3m*s_2p3m
258  q_2m3p=t_2m3p*s_2m3p
259  aux00 =q_2p3m+q_2m3p
260 !
261 !...............Computes the components of the derived terms (dQ/dF)
262  q2pd0 =t_2p3m*(aux03*v1_4 - sp1p*aux04)*cf0
263  q2pd1 =t_2p3m*(aux03 - sp0 *aux04)*cf1
264  q2pd2p=t_2p3m*(aux01*sp1p3m - aux07 )*cf2
265  q2pd3m=t_2p3m*(aux01*sp1p2p - aux08 )*cf3
266  q2md0 =t_2m3p*(aux05*v1_4 - sp1p*aux06)*cf0
267  q2md1 =t_2m3p*(aux03 - sp0 *aux06)*cf1
268  q2md2m=t_2m3p*(aux01*sp1p3p - aux07 )*cf2
269  q2md3p=t_2m3p*(aux01*sp1p2m - aux08 )*cf3
270  aux09=q2pd0+q2md0
271  aux10=q2pd1+q2md1
272 !
273 !...............Sum of Qnl4 term in the table TSTOT
274  t1tot(ip,jt ,jfm0)=t1tot(ip,jt ,jfm0)+aux00 *cp0
275  t1tot(ip,jt1p ,jfm1)=t1tot(ip,jt1p ,jfm1)+aux00 *cp1
276  t1tot(ip,jt1p2p,jfm2)=t1tot(ip,jt1p2p,jfm2)-q_2p3m*cp2
277  t1tot(ip,jt1p2m,jfm2)=t1tot(ip,jt1p2m,jfm2)-q_2m3p*cp2
278  t1tot(ip,jt1p3m,jfm3)=t1tot(ip,jt1p3m,jfm3)-q_2p3m*cp3
279  t1tot(ip,jt1p3p,jfm3)=t1tot(ip,jt1p3p,jfm3)-q_2m3p*cp3
280 !
281 !...............Sum of the term dQnl4/dF in the table T1DER
282  t1der(ip,jt ,jfm0)=t1der(ip,jt ,jfm0)+aux09 *cp0
283  t1der(ip,jt1p ,jfm1)=t1der(ip,jt1p ,jfm1)+aux10 *cp1
284  t1der(ip,jt1p2p,jfm2)=t1der(ip,jt1p2p,jfm2)-q2pd2p*cp2
285  t1der(ip,jt1p2m,jfm2)=t1der(ip,jt1p2m,jfm2)-q2md2m*cp2
286  t1der(ip,jt1p3m,jfm3)=t1der(ip,jt1p3m,jfm3)-q2pd3m*cp3
287  t1der(ip,jt1p3p,jfm3)=t1der(ip,jt1p3p,jfm3)-q2md3p*cp3
288 !
289 ! Config. -Delta1 (SIG=-1)
290 ! ========================
291 !...............Computes the spectrum values in 1, 2, 3
292  sp1m =f2(ip,jt1m ,jfm1)*cf1
293  sp1m2p=f2(ip,jt1m2p,jfm2)*cf2
294  sp1m3m=f2(ip,jt1m3m,jfm3)*cf3
295  sp1m2m=f2(ip,jt1m2m,jfm2)*cf2
296  sp1m3p=f2(ip,jt1m3p,jfm3)*cf3
297 !
298 !...............Computes auxiliary products and variables
299  aux01=sp0*v1_4+sp1m
300  aux02=sp0*sp1m
301  aux03=sp1m2p*sp1m3m
302  aux04=sp1m2p*v3_4+sp1m3m*v2_4
303  aux05=sp1m2m*sp1m3p
304  aux06=sp1m2m*v3_4+sp1m3p*v2_4
305  aux07=aux02*v3_4
306  aux08=aux02*v2_4
307 !
308 !...............Computes the transfer term components
309  s_2p3m=aux03*aux01-aux02*aux04
310  s_2m3p=aux05*aux01-aux02*aux06
311  q_2p3m=t_2m3p*s_2p3m
312  q_2m3p=t_2p3m*s_2m3p
313  aux00 =q_2p3m+q_2m3p
314 !
315 !...............Computes the derived terms components (dQ/dF)
316  q2pd0 =t_2p3m*(aux03*v1_4 - sp1m*aux04)*cf0
317  q2pd1 =t_2p3m*(aux03 - sp0 *aux04)*cf1
318  q2pd2p=t_2p3m*(aux01*sp1m3m - aux07 )*cf2
319  q2pd3m=t_2p3m*(aux01*sp1m2p - aux08 )*cf3
320  q2md0 =t_2m3p*(aux05*v1_4 - sp1m*aux06)*cf0
321  q2md1 =t_2m3p*(aux03 - sp0 *aux06)*cf1
322  q2md2m=t_2m3p*(aux01*sp1m3p - aux07 )*cf2
323  q2md3p=t_2m3p*(aux01*sp1m2m - aux08 )*cf3
324  aux09=q2pd0+q2md0
325  aux10=q2pd1+q2md1
326 !
327 !...............Sum of Qnl4 term in the table T1TOT
328  t1tot(ip,jt ,jfm0)=t1tot(ip,jt ,jfm0)+aux00 *cp0
329  t1tot(ip,jt1m ,jfm1)=t1tot(ip,jt1m ,jfm1)+aux00 *cp1
330  t1tot(ip,jt1m2p,jfm2)=t1tot(ip,jt1m2p,jfm2)-q_2p3m*cp2
331  t1tot(ip,jt1m2m,jfm2)=t1tot(ip,jt1m2m,jfm2)-q_2m3p*cp2
332  t1tot(ip,jt1m3m,jfm3)=t1tot(ip,jt1m3m,jfm3)-q_2p3m*cp3
333  t1tot(ip,jt1m3p,jfm3)=t1tot(ip,jt1m3p,jfm3)-q_2m3p*cp3
334 !
335 !...............Sum of the term dQnl4/dF in the table T1DER
336  t1der(ip,jt ,jfm0)=t1der(ip,jt ,jfm0)+aux09 *cp0
337  t1der(ip,jt1m ,jfm1)=t1der(ip,jt1m ,jfm1)+aux10 *cp1
338  t1der(ip,jt1m2p,jfm2)=t1der(ip,jt1m2p,jfm2)-q2pd2p*cp2
339  t1der(ip,jt1m2m,jfm2)=t1der(ip,jt1m2m,jfm2)-q2md2m*cp2
340  t1der(ip,jt1m3m,jfm3)=t1der(ip,jt1m3m,jfm3)-q2pd3m*cp3
341  t1der(ip,jt1m3p,jfm3)=t1der(ip,jt1m3p,jfm3)-q2md3p*cp3
342 !
343  ENDIF
344 !
345  ENDDO
346 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
347 ! END OF LOOP 4 OVER THE MESH NODES
348 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
349 !
350  ENDDO
351 ! -------------------------------------------------
352 ! END OF LOOP 3 OVER THE SPECTRUM DIRECTIONS
353 ! -------------------------------------------------
354 !
355  ENDDO
356 ! = = = = = = = = = = = = = = = = = = = = = = = = =
357 ! END OF LOOP 2 OVER THE SPECTRUM FREQUENCIES
358 ! = = = = = = = = = = = = = = = = = = = = = = = = =
359 !
360  ENDDO
361 ! ==================================================
362 ! END OF LOOP 1 OVER THE SELECTED CONFIGURATIONS
363 ! ==================================================
364 !
365  RETURN
366  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
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
integer, dimension(:,:,:), allocatable k_1p2p
subroutine qnlin3(T1TOT, T1DER, F2, N1POIN2, N1PLAN, N1F)
Definition: qnlin3.f:7
integer, dimension(:,:,:), allocatable k_1m2m
double precision, dimension(dimbuf) f_coef
integer, dimension(:), allocatable k_if1
integer, dimension(:,:,:), allocatable k_1p2m
integer, dimension(:,:), allocatable k_1p
double precision, dimension(:), pointer tra37
double precision, dimension(dimbuf) tb_sca
integer, dimension(:,:,:), allocatable k_if3
integer, dimension(:,:,:), allocatable k_1p3m
integer, dimension(dimbuf) t_poin
double precision, dimension(:,:,:), allocatable tb_v24
integer, parameter lbuf