The TELEMAC-MASCARET system  trunk
integ_gaia.f
Go to the documentation of this file.
1 ! *********************
2  SUBROUTINE integ_gaia
3 ! *********************
4 !
5  &(a , b , iein , npoin)
6 !
7 !***********************************************************************
8 ! GAIA
9 !***********************************************************************
10 !
26 !
29 !
30 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
35 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36 !
38  IMPLICIT NONE
39 !
40 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
41 !
42  INTEGER, INTENT(IN) :: NPOIN
43  DOUBLE PRECISION, INTENT(IN) :: A(npoin)
44  DOUBLE PRECISION, INTENT(INOUT) :: IEIN(npoin), B(npoin)
45 !
46 !!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
47 !
48  INTEGER, PARAMETER :: LNA = 15, lnb = 42
49  SAVE einj3,aein,bein
50 !
51  DOUBLE PRECISION EINJ3(lnb,lna)
52  DOUBLE PRECISION AEIN(lna),BEIN(lnb)
53 !
54  DOUBLE PRECISION A1, A2, B1, B2, INT1, INT2
55 !
56  INTEGER I, IT, JT, INF, SUP, MIL
57 !
58 ! **********************************************************************
59 ! TABLES
60 ! **********************************************************************
61 !
62 ! RATIO QSS/QSC
63 ! COMPUTED USING EINSTEIN INTEGRALS
64 ! AND GENERATED BY MAPLE V
65 !
66 ! DATA FOR A= .0
67 !
68  einj3(:,1) = (/
69  & .0000d+00, .4806d+00, .9944d+00, .1535d+01, .2386d+01, .3889d+01,
70  & .7128d+01, .9708d+01, .1422d+02, .1798d+02, .2055d+02, .2383d+02,
71  & .2813d+02, .3399d+02, .4243d+02, .5549d+02, .7811d+02, .1257d+03,
72  & .1628d+03, .2269d+03, .2797d+03, .3155d+03, .3609d+03, .4202d+03,
73  & .5005d+03, .6152d+03, .7913d+03, .1093d+04, .1721d+04, .2206d+04,
74  & .3037d+04, .3716d+04, .4176d+04, .4756d+04, .5511d+04, .6531d+04,
75  & .7982d+04, .1020d+05, .1398d+05, .4627d+05, .9803d+05, .5538d+06/
76  & )
77 !
78 ! DATA FOR A= .1
79 !
80  einj3(:,2) = (/
81  & .0000d+00, .4298d+00, .8772d+00, .1338d+01, .2048d+01, .3271d+01,
82  & .5813d+01, .7778d+01, .1113d+02, .1386d+02, .1569d+02, .1801d+02,
83  & .2101d+02, .2505d+02, .3077d+02, .3944d+02, .5407d+02, .8379d+02,
84  & .1063d+03, .1442d+03, .1747d+03, .1951d+03, .2207d+03, .2536d+03,
85  & .2977d+03, .3595d+03, .4526d+03, .6082d+03, .9202d+03, .1154d+04,
86  & .1545d+04, .1857d+04, .2066d+04, .2326d+04, .2660d+04, .3105d+04,
87  & .3727d+04, .4659d+04, .6207d+04, .1844d+05, .3649d+05, .1758d+06/
88  & )
89 !
90 ! DATA FOR A= .2
91 !
92  einj3(:,3) = (/
93  & .0000d+00, .3881d+00, .7823d+00, .1180d+01, .1782d+01, .2792d+01,
94  & .4821d+01, .6344d+01, .8874d+01, .1089d+02, .1223d+02, .1390d+02,
95  & .1604d+02, .1888d+02, .2282d+02, .2869d+02, .3835d+02, .5729d+02,
96  & .7124d+02, .9408d+02, .1121d+03, .1239d+03, .1386d+03, .1573d+03,
97  & .1820d+03, .2160d+03, .2663d+03, .3481d+03, .5065d+03, .6219d+03,
98  & .8095d+03, .9561d+03, .1052d+04, .1171d+04, .1322d+04, .1521d+04,
99  & .1793d+04, .2193d+04, .2841d+04, .7579d+04, .1400d+05, .5756d+05/
100  & )
101 !
102 ! DATA FOR A= .4
103 !
104  einj3(:,4) = (/
105  & .0000d+00, .3241d+00, .6390d+00, .9459d+00, .1393d+01, .2110d+01,
106  & .3460d+01, .4418d+01, .5934d+01, .7092d+01, .7840d+01, .8753d+01,
107  & .9895d+01, .1137d+02, .1335d+02, .1619d+02, .2063d+02, .2877d+02,
108  & .3441d+02, .4320d+02, .4981d+02, .5405d+02, .5919d+02, .6557d+02,
109  & .7375d+02, .8468d+02, .1001d+03, .1241d+03, .1674d+03, .1970d+03,
110  & .2429d+03, .2771d+03, .2989d+03, .3253d+03, .3580d+03, .3997d+03,
111  & .4551d+03, .5332d+03, .6534d+03, .1408d+04, .2273d+04, .6813d+04/
112  & )
113 !
114 ! DATA FOR A= .6
115 !
116  einj3(:,5) = (/
117  & .0000d+00, .2773d+00, .5367d+00, .7818d+00, .1128d+01, .1659d+01,
118  & .2603d+01, .3238d+01, .4199d+01, .4903d+01, .5346d+01, .5877d+01,
119  & .6525d+01, .7340d+01, .8403d+01, .9866d+01, .1205d+02, .1579d+02,
120  & .1823d+02, .2186d+02, .2446d+02, .2609d+02, .2802d+02, .3036d+02,
121  & .3327d+02, .3703d+02, .4215d+02, .4969d+02, .6239d+02, .7059d+02,
122  & .8262d+02, .9120d+02, .9653d+02, .1028d+03, .1104d+03, .1199d+03,
123  & .1320d+03, .1484d+03, .1724d+03, .3023d+03, .4277d+03, .9411d+03/
124  & )
125 !
126 ! DATA FOR A= .8
127 !
128  einj3(:,6) = (/
129  & .0000d+00, .2419d+00, .4607d+00, .6618d+00, .9376d+00, .1346d+01,
130  & .2033d+01, .2474d+01, .3114d+01, .3564d+01, .3841d+01, .4166d+01,
131  & .4555d+01, .5031d+01, .5635d+01, .6435d+01, .7576d+01, .9407d+01,
132  & .1054d+02, .1213d+02, .1323d+02, .1390d+02, .1467d+02, .1559d+02,
133  & .1670d+02, .1808d+02, .1990d+02, .2244d+02, .2644d+02, .2887d+02,
134  & .3226d+02, .3457d+02, .3597d+02, .3758d+02, .3948d+02, .4177d+02,
135  & .4462d+02, .4832d+02, .5346d+02, .7757d+02, .9715d+02, .1601d+03/
136  & )
137 !
138 ! DATA FOR A= 1.0
139 !
140  einj3(:,7) = (/
141  & .0000d+00, .2142d+00, .4023d+00, .5711d+00, .7968d+00, .1120d+01,
142  & .1638d+01, .1956d+01, .2400d+01, .2702d+01, .2884d+01, .3093d+01,
143  & .3338d+01, .3631d+01, .3992d+01, .4455d+01, .5086d+01, .6039d+01,
144  & .6595d+01, .7344d+01, .7837d+01, .8129d+01, .8461d+01, .8844d+01,
145  & .9296d+01, .9844d+01, .1053d+02, .1145d+02, .1281d+02, .1358d+02,
146  & .1461d+02, .1528d+02, .1567d+02, .1611d+02, .1662d+02, .1722d+02,
147  & .1794d+02, .1884d+02, .2003d+02, .2487d+02, .2817d+02, .3657d+02/
148  & )
149 !
150 ! DATA FOR A= 1.2
151 !
152  einj3(:,8) = (/
153  & .0000d+00, .1919d+00, .3563d+00, .5006d+00, .6892d+00, .9512d+00,
154  & .1354d+01, .1592d+01, .1912d+01, .2122d+01, .2246d+01, .2387d+01,
155  & .2549d+01, .2738d+01, .2965d+01, .3247d+01, .3615d+01, .4141d+01,
156  & .4432d+01, .4807d+01, .5044d+01, .5180d+01, .5332d+01, .5504d+01,
157  & .5701d+01, .5933d+01, .6215d+01, .6572d+01, .7066d+01, .7332d+01,
158  & .7668d+01, .7877d+01, .7997d+01, .8129d+01, .8277d+01, .8447d+01,
159  & .8644d+01, .8880d+01, .9178d+01, .1024d+02, .1084d+02, .1207d+02/
160  & )
161 !
162 ! DATA FOR A= 1.5
163 !
164  einj3(:,9) = (/
165  & .0000d+00, .1658d+00, .3032d+00, .4205d+00, .5694d+00, .7685d+00,
166  & .1058d+01, .1220d+01, .1428d+01, .1560d+01, .1635d+01, .1718d+01,
167  & .1812d+01, .1918d+01, .2041d+01, .2188d+01, .2369d+01, .2609d+01,
168  & .2732d+01, .2881d+01, .2970d+01, .3020d+01, .3074d+01, .3133d+01,
169  & .3198d+01, .3271d+01, .3356d+01, .3456d+01, .3583d+01, .3645d+01,
170  & .3719d+01, .3763d+01, .3786d+01, .3812d+01, .3839d+01, .3869d+01,
171  & .3903d+01, .3941d+01, .3986d+01, .4117d+01, .4174d+01, .4259d+01/
172  & )
173 !
174 ! DATA FOR A= 2.0
175 !
176  einj3(:,10) = (/
177  & .0000d+00, .1349d+00, .2418d+00, .3298d+00, .4372d+00, .5737d+00,
178  & .7582d+00, .8542d+00, .9708d+00, .1040d+01, .1078d+01, .1120d+01,
179  & .1164d+01, .1213d+01, .1268d+01, .1329d+01, .1399d+01, .1483d+01,
180  & .1522d+01, .1566d+01, .1590d+01, .1603d+01, .1616d+01, .1630d+01,
181  & .1645d+01, .1661d+01, .1678d+01, .1697d+01, .1718d+01, .1727d+01,
182  & .1737d+01, .1742d+01, .1745d+01, .1748d+01, .1751d+01, .1754d+01,
183  & .1757d+01, .1760d+01, .1764d+01, .1772d+01, .1774d+01, .1777d+01/
184  & )
185 !
186 ! DATA FOR A= 2.5
187 !
188  einj3(:,11) = (/
189  & .0000d+00, .1135d+00, .2004d+00, .2699d+00, .3524d+00, .4531d+00,
190  & .5821d+00, .6459d+00, .7199d+00, .7621d+00, .7849d+00, .8089d+00,
191  & .8344d+00, .8615d+00, .8907d+00, .9222d+00, .9568d+00, .9955d+00,
192  & .1012d+01, .1031d+01, .1040d+01, .1045d+01, .1050d+01, .1055d+01,
193  & .1060d+01, .1066d+01, .1072d+01, .1077d+01, .1084d+01, .1086d+01,
194  & .1089d+01, .1090d+01, .1091d+01, .1091d+01, .1092d+01, .1093d+01,
195  & .1093d+01, .1094d+01, .1095d+01, .1096d+01, .1097d+01, .1097d+01/
196  & )
197 !
198 ! DATA FOR A= 3.0
199 !
200  einj3(:,12) = (/
201  & .0000d+00, .9783d-01, .1708d+00, .2278d+00, .2939d+00, .3724d+00,
202  & .4689d+00, .5148d+00, .5664d+00, .5950d+00, .6101d+00, .6259d+00,
203  & .6424d+00, .6596d+00, .6777d+00, .6969d+00, .7173d+00, .7392d+00,
204  & .7484d+00, .7580d+00, .7630d+00, .7655d+00, .7681d+00, .7707d+00,
205  & .7733d+00, .7759d+00, .7786d+00, .7813d+00, .7841d+00, .7853d+00,
206  & .7864d+00, .7870d+00, .7872d+00, .7875d+00, .7878d+00, .7881d+00,
207  & .7884d+00, .7887d+00, .7890d+00, .7896d+00, .7897d+00, .7898d+00/
208  & )
209 !
210 ! DATA FOR A= 4.0
211 !
212  einj3(:,13) = (/
213  & .0000d+00, .7657d-01, .1314d+00, .1730d+00, .2196d+00, .2727d+00,
214  & .3345d+00, .3624d+00, .3927d+00, .4088d+00, .4172d+00, .4258d+00,
215  & .4346d+00, .4437d+00, .4530d+00, .4626d+00, .4725d+00, .4828d+00,
216  & .4870d+00, .4913d+00, .4935d+00, .4946d+00, .4957d+00, .4968d+00,
217  & .4979d+00, .4990d+00, .5001d+00, .5012d+00, .5023d+00, .5028d+00,
218  & .5033d+00, .5035d+00, .5036d+00, .5037d+00, .5038d+00, .5039d+00,
219  & .5041d+00, .5042d+00, .5043d+00, .5045d+00, .5046d+00, .5046d+00/
220  & )
221 !
222 ! DATA FOR A= 7.0
223 !
224  einj3(:,14) = (/
225  & .0000d+00, .4618d-01, .7721d-01, .9958d-01, .1235d+00, .1493d+00,
226  & .1771d+00, .1890d+00, .2013d+00, .2076d+00, .2108d+00, .2141d+00,
227  & .2174d+00, .2207d+00, .2240d+00, .2274d+00, .2308d+00, .2343d+00,
228  & .2357d+00, .2371d+00, .2378d+00, .2381d+00, .2385d+00, .2388d+00,
229  & .2392d+00, .2396d+00, .2399d+00, .2403d+00, .2406d+00, .2408d+00,
230  & .2409d+00, .2410d+00, .2410d+00, .2410d+00, .2411d+00, .2411d+00,
231  & .2412d+00, .2412d+00, .2412d+00, .2413d+00, .2413d+00, .2413d+00/
232  & )
233 !
234 ! DATA FOR A=15.0
235 !
236  einj3(:,15) = (/
237  & .0000d+00, .2236d-01, .3656d-01, .4639d-01, .5654d-01, .6702d-01,
238  & .7786d-01, .8230d-01, .8681d-01, .8909d-01, .9023d-01, .9138d-01,
239  & .9254d-01, .9369d-01, .9486d-01, .9602d-01, .9720d-01, .9837d-01,
240  & .9884d-01, .9932d-01, .9955d-01, .9967d-01, .9979d-01, .9991d-01,
241  & .1000d+00, .1001d+00, .1003d+00, .1004d+00, .1005d+00, .1005d+00,
242  & .1006d+00, .1006d+00, .1006d+00, .1006d+00, .1007d+00, .1007d+00,
243  & .1007d+00, .1007d+00, .1007d+00, .1007d+00, .1007d+00, .1007d+00/
244  & )
245 !
246 ! VALUES OF TABLE BEIN
247 !
248  bein = (/
249  & 1.d0 ,.75d0 ,.6d0 ,.5d0 ,.4d0 ,.3d0 ,
250  & .2d0 ,.16d0 ,.12d0 ,.1d0 ,.09d0 ,.08d0 ,
251  & .07d0 ,.06d0 ,.05d0 ,.04d0 ,.03d0 ,.02d0 ,
252  & .016d0 ,.012d0 ,.01d0 ,.009d0 ,.008d0 ,.007d0 ,
253  & .006d0 ,.005d0 ,.004d0 ,.003d0 ,.002d0 ,.0016d0 ,
254  & .0012d0,.001d0 ,.0009d0,.0008d0,.0007d0 ,.0006d0 ,
255  & .0005d0,.0004d0,.0003d0,.0001d0,.00005d0,.00001d0 /
256  & )
257 !
258 ! VALUES OF TABLE AEIN
259 !
260  aein = (/ 0.d0,.1d0,.2d0,.4d0,.6d0,.8d0,1.d0,1.2d0,1.5d0,
261  & 2.d0 ,2.5d0,3.d0,4.d0,7.d0,15.d0 /)
262 !
263  DO i=1,npoin
264 !
265 !
266  IF (b(i).LT.bein(lnb)) b(i)=bein(lnb)
267 ! TREATS B AND A
268 ! =======================
269 !
270 ! CASE WHERE SUSPENSION IS NEGLIGIBLE
271 !
272 !
273  IF (a(i).GT.aein(lna).OR.(b(i).GE.1.d0)) THEN
274  iein(i) = 0.d0
275 !
276 ! VALUES OUTSIDE OF THE TABLE
277 !
278  ELSEIF ((b(i).GT.1.d0).OR.(b(i).LT.bein(lnb)).OR.
279  & (a(i).LT.aein(1))) THEN
280 !
281  100 FORMAT(/,1x,'TBIJFR : B EST SUPERIEUR A 1',/,
282  & 3x,'B = ',g10.4,/,
283  & 3x,'NOEUD = ',i5,/,
284  & 3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
285  & ,' DANS VOTRE CAS',/)
286  101 FORMAT(/,1x,'TBIJFR : B EST INFERIEUR A ',g10.4,
287  & ' LIMITE DE LA ZONE TABULEE',/,
288  & 3x,'B = ',g10.4,/,
289  & 3x,'NOEUD = ',i5,/,
290  & 3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
291  & ,' DANS VOTRE CAS',/)
292  102 FORMAT(/,1x,'TBIJFR : A EST INFERIEUR A ',g10.4,
293  & ' LIMITE DE LA ZONE TABULEE',/,
294  & 3x,'A = ',g10.4,/,
295  & 3x,'NOEUD = ',i5,/,
296  & 3x,'VERIFIER LA PERTINENCE DE LA FORMULE DE BIJKER'
297  & ,' DANS VOTRE CAS',/)
298 !
299  IF (b(i).GT.1.d0) THEN
300  WRITE(lu,100) b(i),i
301  ELSEIF (b(i).LT.bein(lnb)) THEN
302  WRITE(lu,101) bein(lnb),b(i),i
303  ELSEIF (a(i).LT.aein(1)) THEN
304  WRITE(lu,102) aein(1),a(i),i
305  ENDIF
306  CALL plante(0)
307  stop
308 !
309 ! VALUES IN THE TABLE
310 !
311  ELSE
312 !
313 ! LOOKS FOR THE BOUNDS FOR B
314 ! ===============================
315 !
316  inf = 1
317  sup = lnb
318 !----START OF THE DO WHILE LOOP----
319  51 mil = (inf + sup ) / 2
320  IF (bein(mil).GT.b(i)) THEN
321  inf = mil
322  ELSE
323  sup = mil
324  ENDIF
325  IF ((sup - inf).NE.1) GOTO 51
326 !----END OF THE DO WHILE LOOP----
327  jt=sup
328  b1 = bein(inf)
329  b2 = bein(sup)
330 !
331 ! LOOKS FOR THE BOUNDS FOR A
332 ! ===============================
333 !
334  inf = 1
335  sup = lna
336 !----START OF THE DO WHILE LOOP----
337  DO
338  mil = (inf + sup ) / 2
339  IF (aein(mil).LT.a(i)) THEN
340  inf = mil
341  ELSE
342  sup = mil
343  ENDIF
344  IF ((sup - inf).EQ.1) EXIT
345  ENDDO
346 !----END OF THE DO WHILE LOOP----
347  it=sup
348  a1 = aein(inf)
349  a2 = aein(sup)
350 !
351 ! COMPUTES THE IEIN INTEGRAL
352 ! ==========================
353 !
354 ! INTERPOLATION WITH CONSTANT A1
355  int1 = (einj3(jt-1,it-1)-einj3(jt,it-1))/(b1-b2)*
356  & (b(i)-b1) + einj3(jt-1,it-1)
357 ! INTERPOLATION WITH CONSTANT A2
358  int2 = (einj3(jt-1,it) -einj3(jt,it) )/(b1-b2)*
359  & (b(i)-b1) + einj3(jt-1,it)
360 ! RE-INTERPOLATION BETWEEN THE TWO PRECEDING VALUES
361  iein(i) = (int1-int2)/(a1-a2)*(a(i)-a1) + int1
362 !
363  ENDIF
364 !
365  ENDDO !I
366 !
367  RETURN
368  END SUBROUTINE integ_gaia
subroutine integ_gaia(A, B, IEIN, NPOIN)
Definition: integ_gaia.f:7