Philippe 07/03/2019: IO bugfix: io_set_mnhversion must be called by all the processes
[MNH-git_open_source-lfs.git] / src / MNH / mode_tmat.f90
1 !MNH_LIC Copyright 2000-2019 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !     algorithme initial créé par Michael Mishchenko (2000)
7 !
8 !     algorithme modifié par Corinne Burlaud (2000) puis Olivier Brunau (2002)
9 !     calcul des paramètres radar en fonction de la forme des gouttes, du taux
10 !     de pluies, de l'élévation du radar, de l'oscillation des gouttes, 
11 !     de la longueur d'onde émise, pour une distribution de taille
12 !     de goutte réelle ou théorique.
13 !     l'essentiel des paramètres modifiables sont demandés à l'utilisateur
14 !     pendant l'execution du programme.
15 !
16 !     Modif par Olivier Caumont (04/2008) pour interfaçage avec diagnostic 
17 !     radar de Méso-NH.
18 !     P. Wautelet 22/01/2019: replace double precision declarations by real(kind(0.0d0)) (to allow compilation by NAG compiler)
19 !
20 !****************************************************************************
21
22 !     This test result was calculated by the code in its
23 !     curent setting.
24 !
25 !     ICHOICE=1  NCHECK=1
26 !     RAT= .963711
27 !     PROLATE SPHEROIDS, A/B=   .5000000
28 !     LAM=  6.283185   MRR= .1500D+01   MRI= .2000D-01
29 !     ACCURACY OF COMPUTATIONS DDELT =  .10D-02
30 !     EQUAL-SURFACE-AREA-SPHERE RADIUS= 10.0000
31 !     thet0= 56.00  thet= 65.00  phi0=114.00  phi=128.00  alpha=145.00
32 !     beta= 52.00
33 !     AMPLITUDE MATRIX
34 !     S11=-.50941D+01 + i* .24402D+02
35 !     S12=-.19425D+01 + i* .19971D+01
36 !     S21=-.11521D+01 + i*-.30977D+01
37 !     S22=-.69323D+01 + i* .24748D+02
38 !     PHASE MATRIX
39 !     650  .3172  -17.9846   10.0498  -12.7580
40 !     -21.1462  631.6322 -127.3059   87.2144
41 !     6    .8322  132.6131  635.2767  -34.7730
42 !     -9.6629  -78.1229   51.4094  643.1738
43 !     time =     .03 min
44 !
45 !   CALCULATION OF THE AMPLITUDE AND PHASE MATRICES FOR                 
46 !   A PARTICLE WITH AN AXIALLY SYMMETRIC SHAPE                   
47 !                                                                      
48 !   This version of the code uses REAL variables,          
49 !   is applicable to spheroids, finite circular cylinders,            
50 !   Chebyshev particles, and generalized Chebyshev particles
51 !   (distorted water drops), and must be used along with the 
52 !   accompanying file ampld.par.f.                
53 !                                                                      
54 !   The code has been developed by Michael Mishchenko at the NASA      
55 !   Goddard Institute for Space Studies, New York.  The development    
56 !   of the code was supported by the NASA Radiation Science Program
57 !   and the NASA FIRE III program.                                                  
58 !   The code may be used without limitations in any not-for-      
59 !   profit scientific research.  The only request is that in any       
60 !   publication using the code the source of the code be acknowledged  
61 !   and relevant references be made.                                   
62 !                                                                       
63 !   The computational method is based on the T-matrix approach         
64 !   [P. C. Waterman, Phys. Rev. D 3, 825 (1971)], also known as        
65 !   the extended boundary condition method.  The method was            
66 !   improved in the following papers:                         
67 !                                                                       
68 !  1.  M. I. Mishchenko and L. D. Travis, T-matrix computations       
69 !      of light scattering by large spheroidal particles,             
70 !      Opt. Commun., vol. 109, 16-21 (1994).                          
71 !                                                                      
72 !   2.  M. I. Mishchenko, L. D. Travis, and A. Macke, Scattering       
73 !       of light by polydisperse, randomly oriented, finite            
74 !       circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).     
75 !                                                                      
76 !   3.  D. J. Wielaard, M. I. Mishchenko, A. Macke, and B. E. Carlson, 
77 !       Improved T-matrix computations for large, nonabsorbing and     
78 !       weakly absorbing nonspherical particles and comparison         
79 !       with geometrical optics approximation, Appl. Opt., vol. 36,    
80 !       4305-4313 (1997).                                             
81 !                                                                      
82 !   A general review of the T-matrix approach can be found in          
83 !                                                                      
84 !   4.  M. I. Mishchenko, L. D. Travis, and D. W. Mackowski,           
85 !       T-matrix computations of light scattering by nonspherical      
86 !       particles:  a review, J. Quant. Spectrosc. Radiat.             
87 !       Transfer, vol. 55, 535-575 (1996).                             
88 !                                                                      
89 !   Additional useful information is contained in the paper
90 !                                                                      
91 !   5.  M. I. Mishchenko and L. D. Travis, Capabilities and            
92 !       limitations of a current FORTRAN implementation of the         
93 !       T-matrix method for randomly oriented, rotationally            
94 !       symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer,   
95 !       vol. 60, 309-324 (1998).                                       
96 !                                                                      
97 !   The definitions and notation used can also be found in
98 !
99 !   6.  M. I. Mishchenko, J. W. Hovenier, and L. D. Travis, Concepts,
100 !       terms, notation.  In "Light Scattering by Nonspherical
101 !       Particles:  Theory, Measurements, and Applications," 
102 !       edited by M. I. Mishchenko, J. W. Hovenier, and L. D. Travis,
103 !      Acedemic Press, San Diego, 1999, pp. 3-27
104 !
105 !   and especially
106 !
107 !   7.  M. I. Mishchenko, Calculation of the amplitude matrix
108 !       for a nonspherical particle in a fixed orientation,
109 !       Appl. Opt. vol. 39, 1026-1031 (2000).
110 !
111 !   Copies of these papers are available upon request from Michael     
112 !   Mishchenko.  Please send your request to crmim@giss.nasa.gov.      
113 !   They are also available in the PDF format at 
114 !   http://www/giss/nasa.gov/~crmim (button "Publications On-Line")
115 !                                                                      
116 !   One of the main advantages of the T-matrix method is that the      
117 !   T-matrix for a given nonspherical particle needs to be computed    
118 !   only once and then can be used any number of times for computing   
119 !   the amplitude and phase matrices for any directions of the incident 
120 !   and scattered beams.  This makes the T-matrix method           
121 !   extremely convenient and efficient in averaging over particle      
122 !   orientations and/or directions of incidence (or scattering).       
123 !                                                                       
124 !   The use of extended precision variables (Ref. 1) can significantly
125 !   increase the maximum convergent equivalent-sphere size parameter 
126 !   and make it well larger than 100 (depending on refractive index    
127 !   and aspect ratio).  The extended-precision code is also available. 
128 !   However, the use of extended precision variables results in a      
129 !   greater consumption of CPU time.                                   
130 !   On IBM RISC workstations, that code is approximately               
131 !   five times slower than this double-precision code.  The            
132 !   CPU time difference between the double-precision and extended-     
133 !   precision codes can be larger on supercomputers.                   
134 !
135 !                                                              
136 !   As described in Ref. 3 above, these changes will affect the        
137 !   performance of the code only for nonabsorbing or weakly            
138 !   absorbing particles (imaginary part of the refractive              
139 !   index smaller than about 0.001).                                   
140
141 !   
142 !      
143 !   INPUT PARAMETERS:                                                  
144 !                                                                     
145 !      AXI - equivalent-sphere radius                                  
146 !      RAT = 1 - particle size is specified in terms of the            
147 !                equal-volume-sphere radius                             
148 !      RAT.ne.1 - particle size is specified in terms of the           
149 !                equal-surface-area-sphere radius                      
150 !      LAM - WAVELENGTH OF INCIDENT LIGHT                                      
151 !      MRR and MRI - real and imaginary parts of the refractive        
152 !                  index                                               
153 !      EPS and NP - specify the shape of the particles.                
154 !             For spheroids NP=-1 and EPS is the ratio of the          
155 !                 horizontal to rotational axes.  EPS is larger than   
156 !                 1 for oblate spheroids and smaller than 1 for        
157 !                 prolate spheroids.                                   
158 !             For cylinders NP=-2 and EPS is the ratio of the          
159 !                 diameter to the length.                              
160 !             For Chebyshev particles NP must be even and positive and 
161 !                 is the degree of the Chebyshev polynomial, while     
162 !                 EPS is the deformation parameter (Ref. 5).                   
163 !             For generalized Chebyshev particles (describing the shape
164 !                 of distorted water drops) NP=-3.  The coefficients
165 !                 of the Chebyshev polynomial expansion of the particle
166 !                 shape (Ref. 7) are specified in subroutine DROP.
167 !      DDELT - accuracy of the computations                            
168 !      NDGS - parameter controlling the number of division points      
169 !             in computing integrals over the particle surface (Ref. 5).       
170 !             For compact particles, the recommended value is 2.       
171 !             For highly aspherical particles larger values (3, 4,...) 
172 !             may be necessary to obtain convergence.                  
173 !             The code does not check convergence over this parameter. 
174 !             Therefore, control comparisons of results obtained with  
175 !             different NDGS-values are recommended.                   
176 !      ALPHA and BETA - Euler angles (in degrees) specifying the orientation 
177 !      of the scattering particle relative to the laboratory reference
178 !         frame (Refs. 6 and 7).
179 !      THET0 - zenith angle of the incident beam in degrees
180 !      THET - zenith angle of the scattered beam in degrees    
181 !      PHI0 - azimuth angle of the incident beam in degrees    
182 !      PHI - azimuth angle of the scattered beam in degrees   
183 !            (Refs. 6 and 7)
184 !      ALPHA, BETA, THET0, THET, PHI0, and PHI are specified at the end of
185 !      the main program before the line                                    
186 !                                                                      
187 !       "CALL AMPL (NMAX,...)"                     
188 !                                                                      
189 !      The part of the main program following the line 
190 !                                                                      
191 !       "COMPUTATION OF THE AMPLITUDE AND PHASE MATRICES"               
192 !                                                                      
193 !      can be repeated any number of times.  At this point the T-matrix 
194 !      for the given scattering particle has already     
195 !      been fully computed and can be repeatedly used in computations  
196 !      for any directions of illumination and scattering and any particle
197 !      orientations.              
198 !                                                                       
199 !   OUTPUT PARAMETERS:                                                 
200 !                                                                      
201 !      Elements of the 2x2 amplitude matrix       
202 !      Elements of the 4x4 phase matrix
203 !                                                                      
204 !   Note that LAM and AXI must be given in the same units of length        
205 !   (e.g., microns).                                                          
206 !                                                                      
207 !   The convergence of the T-matrix method for particles with          
208 !   different sizes, refractive indices, and aspect ratios can be      
209 !   dramatically different.  Usually, large sizes and large aspect     
210 !   ratios cause problems.  The user of this code                      
211 !   should "play" a little with different input parameters in          
212 !   order to get an idea of the range of applicability of this         
213 !   technique.  Sometimes decreasing the aspect ratio                  
214 !   from 3 to 2 can increase the maximum convergent equivalent-        
215 !   sphere size parameter by a factor of several.                      
216 !                                                                      
217 !   The CPU time required rapidly increases with increasing ratio      
218 !   radius/wavelength and/or with increasing particle asphericity.     
219 !   This should be taken into account in planning massive computations.
220 !                                                                      
221 !   Execution can be automatically terminated if dimensions of certain 
222 !   arrays are not big enough or if the convergence procedure decides  
223 !   that the accuracy of double-precision variables is insufficient  
224 !   to obtain a converged T-matrix solution for given particles.       
225 !   In all cases, a message appears explaining                         
226 !   the cause of termination.                                          
227 !                                                                      
228 !   The message                                                        
229 !        "WARNING:  W IS GREATER THAN 1"                               
230 !   means that the single-scattering albedo exceeds the maximum        
231 !   possible value 1.  If W is greater than 1 by more than             
232 !   DDELT, this message can be an indication of numerical              
233 !   instability caused by extreme values of particle parameters.       
234 !                                                                      
235 !   The message "WARNING: NGAUSS=NPNG1" means that convergence over    
236 !   the parameter NG (see Ref. 2) cannot be obtained for the NPNG1     
237 !   value specified in the PARAMETER statement in the file ampld.par.f.
238 !   Often this is not a serious problem, especially for compact
239 !   particles.
240 !                                                                      
241 !   Larger and/or more aspherical particles may require larger
242 !   values of the parameters NPN1, NPN4, and NPNG1 in the file
243 !   ampld.par.f.  It is recommended to keep NPN1=NPN4+25 and
244 !   NPNG1=3*NPN1.  Note that the memory requirement increases
245 !   as the third power of NPN4. If the memory of
246 !   a computer is too small to accomodate the code in its current
247 !   setting, the parameters NPN1, NPN4, and NPNG1 should be
248 !   decreased. However, this will decrease the maximum size parameter
249 !   that can be handled by the code.
250 !
251 !   In some cases any increases of NPN1 will not make the T-matrix     
252 !   computations convergent.  This means that the particle is just     
253 !   too "bad" (extreme size parameter and/or extreme aspect ratio      
254 !   and/or extreme refractive index).                                  
255 !   The main program contains several PRINT statements which are       
256 !   currently commentd out.  If uncommented, these statements will     
257 !   produce numbers which show the convergence rate and can be         
258 !   used to determine whether T-matrix computations for given particle 
259 !   parameters will converge at all, whatever the parameter NPN1 is.   
260 !                                                                       
261 !   Some of the common blocks are used to save memory rather than      
262 !   to transfer data.  Therefore, if a compiler produces a warning     
263 !   message that the lengths of a common block are different in        
264 !   different subroutines, this is not a real problem.                 
265 !                                                                       
266 !   The recommended value of DDELT is 0.001.  For bigger values,       
267 !   false convergence can be obtained.                                 
268 !                                                                       
269 !   In computations for spheres, use EPS=1.000001 instead of EPS=1.    
270 !   EPS=1 can cause overflows in some rare cases.                      
271 !                                                                       
272 !   For some compilers, DACOS must be raplaced by DARCOS and DASIN     
273 !   by DARSIN.                                                         
274 !                                                                       
275 !   I would highly appreciate informing me of any problems encountered 
276 !   with this code.  Please send your message to the following         
277 !   e-mail address:  CRMIM@GISS.NASA.GOV.                              
278 !
279 !   WHILE THE COMPUTER PROGRAM HAS BEEN TESTED FOR A VARIETY OF CASES,
280 !   IT IS NOT INCONCEIVABLE THAT IT CONTAINS UNDETECTED ERRORS. ALSO,
281 !   INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
282 !   VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
283 !   THE AUTHORS AND THEIR ORGANIZATION DISCLAIM ALL LIABILITY FOR
284 !   ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAM. 
285
286
287       SUBROUTINE TMD(Deq,LAM,MRR,MRI,NDGS,OSCIL,ELEV,qbackHH,&
288           qbackVV,kdp,QQEXT,EPS)
289
290
291       USE MODD_TMAT,ONLY : XRT11,XRT12,XRT21,XRT22,&
292                            XIT11,XIT12,XIT21,XIT22,&
293                            XTR1,XTI1,NPN1,NPNG1,NPNG2,NPN2,NPN4,NPN6
294
295       IMPLICIT REAL*8 (A-H,O-Z)
296
297 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1, NPN2=2*NPN1,&
298 !!          NPN4=NPN1, NPN6=NPN4+1)
299       
300
301
302       INTEGER param,oscil
303       REAL*8 SIGBETA,Fbeta
304       REAL*8 PDtotal,PDalpha,PDbeta,Poids
305
306       REAL*8  LAM,MRR,MRI,Deq,X(NPNG2),W(NPNG2),S(NPNG2),SS(NPNG2),&
307              AN(NPN1),R(NPNG2),DR(NPNG2),&
308              DDR(NPNG2),DRR(NPNG2),DRI(NPNG2),ANN(NPN1,NPN1)
309
310       REAL*8 NUMZ,DENZ,NUMTZ,DENTZ,ZDRT,NTotal
311       REAL*8 NUMD,DEND,DELTA,NUMTD,DENTD,DELTAT
312       REAL*8 PAS,KDP
313
314       REAL*8 THET0,THET,Elev,AXI
315
316       REAL*8 MYS11cr,MYS22cr
317
318
319 !!      REAL*8 TR1(NPN2,NPN2),TI1(NPN2,NPN2)
320 !!     REAL*4 RT11(NPN6,NPN4,NPN4),RT12(NPN6,NPN4,NPN4),&
321 !!          RT21(NPN6,NPN4,NPN4),RT22(NPN6,NPN4,NPN4),&
322 ! !         IT11(NPN6,NPN4,NPN4),IT12(NPN6,NPN4,NPN4),&
323 !!          IT21(NPN6,NPN4,NPN4),IT22(NPN6,NPN4,NPN4)
324       COMPLEX*16 S11,S12,S21,S22
325       COMPLEX*16 S11u,S12u,S21u,S22u
326
327       REAL*8 S11carre,S22carre
328       COMPLEX*16 NUMrhoAB,NUMTrhoAB
329       
330 !!      COMMON /CT/ TR1,TI1
331 !!      COMMON /TMAT/ RT11,RT12,RT21,RT22,IT11,IT12,IT21,IT22
332       LOGICAL SORTIE1
333                 
334 !****************************************************************
335
336 ! initialisation des variables
337
338 !c      IF (freq.EQ.1) then
339 !
340 !c     I-Variables depending on the wavelength:
341 !C     Complex refractive index
342 !c         write(8,9887)
343 !c 9887    format(20x,'bande S')         
344 !C     1)  BANDE S
345 !c         MRR=9.006 D0
346 !c         MRI=0.9531 D0
347 !c         LAM=10.7D-2
348 !c         K2=0.9313
349 !C     imag(-K)
350 !c         IMK=0.00688         
351 !c      elseif(freq.EQ.2) then         
352 !c         write(8,9888)
353 !c 9888    format(20x,'bande X')          
354 !C     2) BANDE X
355 !c         MRR=7.813D0
356 !c         MRI=2.427D0
357 !c         LAM=3.20D-2
358 !c         K2=0.9282 
359 !c         IMK=0.0247         
360 !c      endif
361 !      
362 !c      K2=ABS(((MRR+(0,1)*MRI)**2-1.)/((MRR+(0,1)*MRI)**2+2.))**2
363 !C      IMK=AIMAG(((MRR+(0,1)*MRI)**2-1.)/((MRR+(0,1)*MRI)**2+2.))
364 !c      PRINT*,K2,IMK,'>0'
365
366       P=ACOS(-1D0)             !calcul de pi!
367       
368 !****  Lecture du fichier de DSD**************************************       
369 !*********************spectre théorique**********************
370 !*     constante en m-3.mm-1
371 !c      N0=8000D0
372 !      
373 !*     taux de precipitation en mm/h
374       
375 !!c      R_lu(1)=0.1D0
376 !c      R_lu(2)=0.5D0
377 !c      R_lu(3)=1D0
378 !c      R_lu(4)=5D0
379 !c      DO I=1,5   
380 !c         R_lu(I+4)=R_lu(I+3)+5D0
381 !c      ENDDO
382 !      
383 !c      DO I=1,5
384 !c         R_lu(I+9)=R_lu(I+8)+10D0
385 !c      ENDDO
386 !      
387 !c      R_lu(15)=100D0
388 !c      R_lu(16)=125D0
389 !c      R_lu(17)=150D0
390 !      
391 !      
392 !c      DO echant=1,17
393 !c         Rmmh(echant)=R_lu(echant)
394 !c         LAMDAf(echant)=4.1D0*(Rmmh(echant)**(-0.21D0))
395 !c      enddo
396 !      
397 !c      AB=1
398 !      
399 !c      do echant=1,AB
400 !c     Rmmh(echant)=0
401 !*     initial diameter:Deq=0.1 millimètre
402 !*     initial radius:AXI=Deq/2
403 !         
404 !c         AXI=0.05D0
405 !         
406 !c         NDeq(echant,1)=1
407 !c         Deq_lu(echant,1)=Deq
408 !c     print*,Deq_lu(echant,1)
409 !c     calcul de zthéorique
410 !c         Deq_lu2(echant,1)=Deq_lu(echant,1)**6
411 !c         somm(echant)=somm(echant)+
412 !c     &        (NDeq(echant,1)*Deq_lu2(echant,1))
413 !                           
414 !c     recalcul de Rmmh ********************
415 !               
416 !               
417 !c     Vit(echant,I)=3.78*Deq_lu(echant,I)**0.67
418 !c     Volgtte(echant,I)= (4/3)*P*(AXI*1D-3)**3
419 !               
420 !c     Rgtte(echant,I)=Volgtte(echant,I)*Vit(echant,I)*
421 !c     &             (NDeq(echant,I)*1D3)*(deltaD*1D-3)
422 !c     Rmmh(echant)=Rmmh(echant)+Rgtte(echant,I)*1D3*3600
423 !               
424 !               
425 !******************************************
426 !c         AXI=AXI+0.05D0
427 !            
428 !c         sample(echant,1)=echant
429 !c         Ztheorique(echant)=10*lOG10(somm(echant))
430 !c      enddo
431 !      
432 !
433 !*****************************************************************
434 !*****************début du fichier de sortie**********************
435           
436       THET0=90D0-Elev
437       do param=1,2
438
439          IF(param.EQ.1) then
440             THET=THET0
441          else 
442 !c param=2
443             THET=180-THET0
444          endif
445          
446 !********Initialisation ***************
447          ALPHA=0D0
448          BETA=0D0
449          
450          Ntotal=0
451          NUMTZ=0
452          DENTZ=0
453          NUMTD=0
454          DENTD=0
455          MYS11cr=0D0
456          MYS22cr=0D0
457          NUMTrhoAB=0D0
458 !c     ZHHT=0D0
459 !c     ZVVT=0D0
460 !         
461 !c     NDeqtot=0
462 !c???????
463 !         
464 !****************************************
465 !*********debut de la boucle pour chaque diamètre de goutte mesuré
466 !*********32 diamètres par echantillon réel, 64 pour echantillon théorique
467 !         
468 !C     Diametre equivolumetrique en mètre           
469 !c         Deq=Deq_lu(echant,1)*1D-3
470 !            
471 !         AXI=Deq/2
472 !         
473 !c         IF(PARF.EQ."PB70") THEN
474 !c     détermination de la forme de la goutte (Pruppacher & Beard, 1970)
475 !c            if (gtte.EQ.2.AND.Deq.GE.5D-4) then               
476 !c     oblate case
477 !c               EPS=1./(1.03-62.*Deq)       
478 !c            ELSE 
479 !c*     spherical case
480 !c               EPS=1.
481 !c            ENDIF               
482 !c         ELSE
483 !c     détermination de la forme de la goutte (Andsager et al., 1999)
484 !c            if (gtte.EQ.2) then               
485 !c     oblate case
486 !c               if(Deq.ge.1.1E-3.and.Deq.le.4.4E-3) THEN
487 !c                  EPS=1./(1.012-14.4*Deq-1.03E4*Deq**2)
488 !c                  IF(EPS.LT.1.) EPS=1.
489 !c               ELSE
490 !c                  EPS=1.0048+.57*Deq-2.628E4*Deq**2+3.682E6*Deq**3
491 !c     &                 -1.677E8*Deq**4
492 !c                  IF(EPS.LT.0.2) THEN
493 !c                     EPS=5.
494 !c                  ELSE
495 !c                     EPS=1./EPS
496 !c                  ENDIF
497 !c               ENDIF
498 !c           ELSE 
499 !c*     spherical case
500 !c               EPS=1.
501 !c            ENDIF               
502 !c         ENDIF
503 !C         print*,eps
504 !c     WRITE(*,854) 1
505 !c     WRITE(*,855) Deq
506 !c     854  FORMAT('Resultat pour la',I3,' eme goutte')
507 !c     855  FORMAT('Deq=',F11.7)
508 !      
509 !************algorithme NASA*************
510       
511          DDELT=0.001D0 
512 !c     NDGS=2
513 !         
514 !c     NCHECK=1
515 !         
516 !C     IF (NP.EQ.-1) NCHECK=1
517 !C     IF (NP.GT.0.AND.(-1)**NP.EQ.1) NCHECK=1
518 !         
519 !c     WRITE (10,5454) NCHECK
520 !c     5454 FORMAT ('  NCHECK=',I1)
521 !         
522 !!c     IF(ABS(RAT-1D0).GT.1D-8.AND.NP.EQ.-1) CALL SAREA(EPS,RAT)
523 !c     IF(ABS(RAT-1D0).GT.1D-8.AND.NP.GE.0) CALL SURFCH(NP,EPS,RAT)
524 !C     IF(NP.EQ.-3) CALL DROP(RAT)
525 !         
526 !c     8000 FORMAT ('RAT=',F8.6)
527 !         
528 !C     IF(NP.EQ.-1.AND.EPS.GE.1D0) WRITE(10,7000)EPS
529 !C     IF(NP.EQ.-1.AND.EPS.LT.1D0) WRITE(10,7001)EPS
530 !C     IF(NP.GE.0) WRITE(10,7100)NP,EPS
531 !C     IF(NP.EQ.-3)WRITE(10,7160) 
532 !C     WRITE(10,7400) LAM,MRR,MRI
533 !C     WRITE(10,7200) DDELT
534          
535 ! 7000    FORMAT('OBLATE SPHEROIDS, A/B=',F11.7)
536 ! 7001    FORMAT('PROLATE SPHEROIDS, A/B=',F11.7)
537 ! 7100    FORMAT('CHEBYSHEV PARTICLES, T',I1,'(',F5.2,')')
538 ! 7160    FORMAT('GENERALIZED CHEBYSHEV PARTICLES')
539 ! 7200    FORMAT ('ACCURACY OF COMPUTATIONS DDELT = ',D8.2)
540 ! 7400    FORMAT('LAM=',F10.6,3X,'MRR=',D10.4,3X,'MRI=',D10.4)
541          
542          DDELT=0.1D0*DDELT
543          
544 !c     IF (ABS(RAT-1D0).LE.1D-6)WRITE(*,8003)AXI 
545 !c     IF (ABS(RAT-1D0).GT.1D-6) WRITE(*,8004) AXI
546 !         
547 !c     8003 FORMAT('EQUAL-VOLUME-SPHERE RADIUS=',F9.7)
548 !c     8004 FORMAT('EQUAL-SURFACE-AREA-SPHERE RADIUS=',F9.7)
549          
550       A=AXI
551       XEV=2D0*P*A/LAM
552       IXXX=XEV+4.05D0*XEV**0.333333D0
553       INM1=MAX0(4,IXXX)
554       
555 !C     IF (INM1.GE.NPN1)WRITE(10,7333) NPN1
556       IF (INM1.GE.NPN1) STOP
557       
558 ! 7333 FORMAT('CONVERGENCE IS NOT OBTAINED FOR NPN1=',I3,  &
559 !          '.  EXECUTION TERMINATED')
560       
561       QEXT1=0D0
562       QSCA1=0D0
563       
564 !****  Initialisation      
565       NMA=INM1
566       SORTIE1=.FALSE.
567       
568       DO WHILE ((NMA.LE.NPN1).AND.(SORTIE1.EQV..FALSE.))
569          NMAX=NMA
570 !c         MMAX=1
571          NGAUSS=NMAX*NDGS
572          
573 !C     IF (NGAUSS.GT.NPNG1) WRITE(10,7340) NGAUSS
574          IF (NGAUSS.GT.NPNG1) STOP
575          
576 !c 7340    FORMAT('NGAUSS =',I3,' I.E. IS GREATER THAN NPNG1.',
577 !c     &        '  EXECUTION TERMINATED')
578 !c 7334    FORMAT(' NMAX =', I3,'  DSCA=',D8.2,'   DEXT=',D8.2)
579          
580          CALL CONST(NGAUSS,NMAX,X,W,AN,ANN,S,SS)
581          CALL VARY(LAM,MRR,MRI,A,EPS,NGAUSS,X,P,PPI,PIR,PII,R,&
582              DR,DDR,DRR,DRI,NMAX)
583          CALL TMATR0(NGAUSS,X,W,AN,ANN,PPI,PIR,PII,R,DR,&
584              DDR,DRR,DRI,NMAX)
585          QEXT=0D0
586          QSCA=0D0
587          
588          DO N=1,NMAX
589             N1=N+NMAX
590             TR1NN=XTR1(N,N)
591             TI1NN=XTI1(N,N)
592             TR1NN1=XTR1(N1,N1)
593             TI1NN1=XTI1(N1,N1)
594             DN1=FLOAT(2*N+1)
595             QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN+&
596                 TR1NN1*TR1NN1+TI1NN1*TI1NN1)
597             QEXT=QEXT+(TR1NN+TR1NN1)*DN1
598 !            TR1NN=TR1(N,N)
599 !            TI1NN=TI1(N,N)
600 !            TR1NN1=TR1(N1,N1)
601 !            TI1NN1=TI1(N1,N1)
602 !            DN1=FLOAT(2*N+1)
603 !            QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN+&
604 !                TR1NN1*TR1NN1+TI1NN1*TI1NN1)
605 !            QEXT=QEXT+(TR1NN+TR1NN1)*DN1
606          ENDDO
607             
608          DSCA=ABS((QSCA1-QSCA)/QSCA)
609          DEXT=ABS((QEXT1-QEXT)/QEXT)
610          QEXT1=QEXT
611          QSCA1=QSCA
612          
613 !c     PRINT 7334, NMAX,DSCA,DEXT
614          
615          IF(.not.(DSCA.LE.DDELT.AND.DEXT.LE.DDELT)) THEN
616 !C     IF (NMA.EQ.NPN1) WRITE(10,7333) NPN1
617             IF (NMA.EQ.NPN1) STOP      
618          ELSE
619             SORTIE1=.TRUE.
620          ENDIF
621          
622          NMA=NMA+1
623          
624       ENDDO
625       
626       NNNGGG=NGAUSS+1
627 !C      MMAX=NMAX
628       
629 !C     IF (NGAUSS.EQ.NPNG1) WRITE(10,7336) 
630       
631 !****  Initialisation
632       SORTIE1=.FALSE.
633       NGAUS=NNNGGG
634       
635       IF (NGAUSS.NE.NPNG1) THEN
636          DO WHILE ((SORTIE1.EQV..FALSE.).AND.(NGAUS.LE.NPNG1))
637             NGAUSS=NGAUS
638 !c            NGGG=2*NGAUSS
639  7336       FORMAT('WARNING: NGAUSS=NPNG1')
640  7337       FORMAT(' NG=',I3,'  DSCA=',D8.2,'   DEXT=',D8.2)
641             CALL CONST(NGAUSS,NMAX,X,W,AN,ANN,S,SS)
642             CALL VARY(LAM,MRR,MRI,A,EPS,NGAUSS,X,P,PPI,PIR,PII,R,&
643                 DR,DDR,DRR,DRI,NMAX)
644             CALL TMATR0(NGAUSS,X,W,AN,ANN,PPI,PIR,PII,R,DR,&
645                 DDR,DRR,DRI,NMAX)
646             QEXT=0D0
647             QSCA=0D0
648             
649             DO N=1,NMAX
650                N1=N+NMAX
651                TR1NN=XTR1(N,N)
652                TI1NN=XTI1(N,N)
653                TR1NN1=XTR1(N1,N1)
654                TI1NN1=XTI1(N1,N1)
655                DN1=FLOAT(2*N+1)
656                QSCA=QSCA+DN1*(TR1NN*TR1NN+TI1NN*TI1NN+&
657                     TR1NN1*TR1NN1+TI1NN1*TI1NN1)
658                QEXT=QEXT+(TR1NN+TR1NN1)*DN1
659             ENDDO
660             
661             DSCA=ABS((QSCA1-QSCA)/QSCA)
662             DEXT=ABS((QEXT1-QEXT)/QEXT)
663             
664 !c     PRINT 7337, NGGG,DSCA,DEXT
665             
666             QEXT1=QEXT
667             QSCA1=QSCA
668             
669             IF(.NOT.(DSCA.LE.DDELT.AND.DEXT.LE.DDELT)) THEN
670 !C     IF (NGAUS.EQ.NPNG1) WRITE(10,7336)
671             ELSE
672                SORTIE1=.TRUE.
673             ENDIF
674             
675             NGAUS=NGAUS+1
676             
677          ENDDO
678          
679       ENDIF
680       
681       QSCA=0D0
682       QEXT=0D0
683       NNM=NMAX*2
684       
685       DO N=1,NNM
686          QEXT=QEXT+XTR1(N,N)
687       ENDDO
688       
689       DO N2=1,NMAX
690          NN2=N2+NMAX
691          DO N1=1,NMAX
692             NN1=N1+NMAX
693             ZZ1=XTR1(N1,N2)
694             XRT11(1,N1,N2)=ZZ1
695             ZZ2=XTI1(N1,N2)
696             XIT11(1,N1,N2)=ZZ2
697             ZZ3=XTR1(N1,NN2)
698             XRT12(1,N1,N2)=ZZ3
699             ZZ4=XTI1(N1,NN2)
700             XIT12(1,N1,N2)=ZZ4
701             ZZ5=XTR1(NN1,N2)
702             XRT21(1,N1,N2)=ZZ5
703             ZZ6=XTI1(NN1,N2)
704             XIT21(1,N1,N2)=ZZ6
705             ZZ7=XTR1(NN1,NN2)
706             XRT22(1,N1,N2)=ZZ7
707             ZZ8=XTI1(NN1,NN2)
708             XIT22(1,N1,N2)=ZZ8
709             QSCA=QSCA+ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4+&
710                  ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8
711          ENDDO
712       ENDDO
713       
714 !c     PRINT 7800,0,ABS(QEXT),QSCA,NMAX
715       
716       DO  M=1,NMAX
717          CALL TMATR(M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,&
718                DDR,DRR,DRI,NMAX)
719          NM=NMAX-M+1
720          M1=M+1
721          QSC=0D0
722          DO N2=1,NM
723             NN2=N2+M-1
724             N22=N2+NM
725             DO N1=1,NM
726                NN1=N1+M-1
727                N11=N1+NM
728                ZZ1=XTR1(N1,N2)
729                XRT11(M1,NN1,NN2)=ZZ1
730                ZZ2=XTI1(N1,N2)
731                XIT11(M1,NN1,NN2)=ZZ2
732                ZZ3=XTR1(N1,N22)
733                XRT12(M1,NN1,NN2)=ZZ3
734                ZZ4=XTI1(N1,N22)
735                XIT12(M1,NN1,NN2)=ZZ4
736                ZZ5=XTR1(N11,N2)
737                XRT21(M1,NN1,NN2)=ZZ5
738                ZZ6=XTI1(N11,N2)
739                XIT21(M1,NN1,NN2)=ZZ6
740                ZZ7=XTR1(N11,N22)
741                XRT22(M1,NN1,NN2)=ZZ7
742                ZZ8=XTI1(N11,N22)
743                XIT22(M1,NN1,NN2)=ZZ8
744                QSC=QSC+(ZZ1*ZZ1+ZZ2*ZZ2+ZZ3*ZZ3+ZZ4*ZZ4+&
745                   ZZ5*ZZ5+ZZ6*ZZ6+ZZ7*ZZ7+ZZ8*ZZ8)*2D0
746             ENDDO
747          ENDDO
748          
749          NNM=2*NM
750          QXT=0D0
751          
752          DO N=1,NNM
753             QXT=QXT+XTR1(N,N)*2D0
754          ENDDO
755          
756          QSCA=QSCA+QSC
757          QEXT=QEXT+QXT
758          
759 !c     PRINT 7800,M,ABS(QXT),QSC,NMAX
760 ! 7800    FORMAT(' m=',I3,'  qxt=',D12.6,'  qsc=',D12.6,&
761 !            '  nmax=',I3)
762          
763       ENDDO
764       
765 !c      WALB=-QSCA/QEXT
766       
767 !C     IF (WALB.GT.1D0+DDELT)WRITE(10,9111) 
768 ! 9111 FORMAT ('WARNING: W IS GREATER THAN 1')
769 !*****************************************************
770 !C     COMPUTATION OF THE AMPLITUDE AND PHASE MATRICES
771 !***************************************************
772       
773 !*     Initialisation
774       ALPHA=0D0
775       PDalpha=1D0/9D0
776       Poids=0D0
777       PDtotal=0
778       S11=0D0
779       S12=0D0
780       S21=0D0
781       S22=0D0
782       
783       S11carre=0D0
784       S22carre=0D0
785       NUMrhoAB=0D0
786       
787       
788 !c     NDeq(echant,I)=NDeq(echant,I)*1D3
789       
790       if (param.EQ.1) then
791 !c     calcul paramètres de propagation
792          
793 !C     PHI0 restera toujours fixe
794          PHI0=0D0
795          PHI=PHI0
796          
797 !*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
798          if(oscil.EQ.1) then
799 !c     algorithme de l'oscillation des gouttes
800             
801 !*     ecart type de beta
802             IF (Deq.LE.2D-3) THEN
803                SIGBETA=90D0*EXP(-0.95D0*((Deq*1D3)**2))
804             ELSE
805                SIGBETA=2D0
806             ENDIF
807             
808             PAS=SIGBETA/10D0
809 !******************************************
810             BETA=0D0
811             DO WHILE (BETA.LE.(2*SIGBETA))
812                Fbeta=EXP(-(BETA**2)/(2*SIGBETA**2))/&
813                    (SQRT(2*P*SIGBETA**2))
814                Poids=Poids+Fbeta
815                BETA=BETA+PAS
816             ENDDO
817 !*****************************************
818             
819 !***********boucle  pour la variation de alpha
820             
821             DO WHILE (ALPHA.LT.90)
822                BETA=0D0
823                
824 !***********boucle    pour la variation de beta
825                
826                DO WHILE (BETA.LE.(2*SIGBETA))
827                   
828                   Fbeta=EXP(-(BETA**2)/(2*SIGBETA**2))/&
829                   (SQRT(2*P*SIGBETA**2))
830                   PDbeta=Fbeta/Poids 
831                   
832                   
833 !C     AMPLITUDE MATRIX [Eqs. (2)-(4) of Ref. 6]
834                   CALL AMPL(NMAX,LAM,THET0,THET,PHI0,PHI,ALPHA,BETA,&
835                       S11u,S12u,S21u,S22u)     
836                   
837                   S11=S11+S11u*PDalpha*PDbeta
838                   S12=S12+S12u*PDalpha*PDbeta
839                   S21=S21+S21u*PDalpha*PDbeta
840                   S22=S22+S22u*PDalpha*PDbeta
841                   
842                   
843                   BETA=BETA+PAS
844                   
845                   PDtotal=PDbeta*PDalpha+PDtotal
846                   
847                ENDDO
848                
849 !*     end de la boucle de variation de beta
850                
851                ALPHA=ALPHA+10D0
852                
853             ENDDO
854             
855 !*     end de la boucle de variation de alpha
856             
857             S11=S11/PDtotal
858             S12=S12/PDtotal
859             S21=S21/PDtotal
860             S22=S22/PDtotal
861             
862             
863 !*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
864 !c     fin oscillation des gouttes
865 !c            print*,'youve been here'
866
867          elseif(oscil.EQ.2) then
868 !C     AMPLITUDE MATRIX [Eqs. (2)-(4) of Ref. 6]
869             CALL AMPL(NMAX,LAM,THET0,THET,PHI0,PHI,ALPHA,BETA,&
870                 S11,S12,S21,S22)    
871          endif      
872          
873 !c         DIFFER=DBLE(S22-S11)*LAM**2*1D8          
874          
875          
876 !C     Attenuation (H) en dB/km
877 !C     Ah=4343*LAM*AIMAG(S22)
878          
879 !C         Ah=8686*LAM*AIMAG(S22)
880 !C         print*,aimag(s22),Deq
881          QQEXT=8.*LAM*AIMAG(S22)/P/Deq**2
882
883 !C     A reference en db/km
884 !C     Aref=(0.4343*P**2*Deq**3/LAM)*AIMAG(-K)
885 !C     LAM et Deq sont en cm
886             
887 !C         Aref=(0.8686*P**2*(Deq*1D2)**3/(LAM*1D2))*IMK
888 !         
889 !C     Ah-v=4343*LAM*AIMAG(S22-S11)
890             
891 !C         Ahv=8686*LAM*AIMAG(S22-S11)
892             
893 !C     Changement de phase KDP en degres/km
894             
895 !C         KDP=(180/P)*1D3*LAM*DBLE(S22-S11) 
896 !C         KDP=Deq*DBLE(S22-S11)
897          KDP=LAM**2*DBLE(S22-S11)/2./P            
898          
899 !*************************************************
900
901       elseif(param.EQ.2) then
902 !c     calcul paramètres de retrodiffusion
903          
904 !C     PHI0 restera toujours fixe
905          PHI0=0D0
906          PHI=180-PHI0
907          
908 !*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
909          if(oscil.EQ.1) then
910             
911 !c     algorithme de l'oscillation des gouttes
912             
913 !*     ecart type de beta
914             IF (Deq.LE.2D-3) THEN
915                SIGBETA=90D0*EXP(-0.95D0*((Deq*1D3)**2))
916             ELSE
917                SIGBETA=2D0
918             ENDIF
919             
920             PAS=SIGBETA/10D0
921 !******************************************
922             BETA=0D0
923             DO WHILE (BETA.LE.(2*SIGBETA))
924                Fbeta=EXP(-(BETA**2)/(2*SIGBETA**2))/&
925                    (SQRT(2*P*SIGBETA**2))
926                Poids=Poids+Fbeta
927                BETA=BETA+PAS
928             ENDDO
929 !*****************************************
930             
931 !***********boucle  pour la variation de alpha
932                
933             DO WHILE (ALPHA.LT.90)
934                BETA=0D0
935                
936 !***********boucle    pour la variation de beta
937                
938                DO WHILE (BETA.LE.(2*SIGBETA))
939                   
940                   Fbeta=EXP(-(BETA**2)/(2*SIGBETA**2))/&
941                       (SQRT(2*P*SIGBETA**2))
942                   PDbeta=Fbeta/Poids 
943                   
944                   
945 !C     AMPLITUDE MATRIX [Eqs. (2)-(4) of Ref. 6]
946                   CALL AMPL(NMAX,LAM,THET0,THET,PHI0,PHI,ALPHA,BETA,&
947                       S11u,S12u,S21u,S22u)     
948                   
949                   S11=S11+S11u*PDalpha*PDbeta
950                   S12=S12+S12u*PDalpha*PDbeta
951                   S21=S21+S21u*PDalpha*PDbeta
952                   S22=S22+S22u*PDalpha*PDbeta
953                   
954 !c                  print*,'NUMrhoA',NUMrhoAB,pdbeta,pdalpha
955
956                   NUMrhoAB=NUMrhoAB+S11u*CONJG(S22u)*PDalpha*PDbeta
957                   
958                   S11carre=S11carre+((DBLE(S11u))**2+(AIMAG(S11u))**2)*&
959                       PDalpha*PDbeta
960                   
961                   S22carre=S22carre+((DBLE(S22u))**2+(AIMAG(S22u))**2)*&
962                       PDalpha*PDbeta
963                   
964                   BETA=BETA+PAS
965 !c                  print*,'NUMrhoA',NUMrhoAB,pdbeta,pdalpha
966                   
967                   PDtotal=PDbeta*PDalpha+PDtotal
968                   
969                ENDDO
970                
971 !*     end de la boucle de variation de beta
972                
973                ALPHA=ALPHA+10D0
974                
975             ENDDO
976             
977 !*     end de la boucle de variation de alpha
978             
979 !c            print*,'pdtotal',pdtotal
980             
981             S11=S11/PDtotal
982             S12=S12/PDtotal
983             S21=S21/PDtotal
984             S22=S22/PDtotal
985             
986             S11carre=S11carre/PDtotal
987             S22carre=S22carre/PDtotal
988             NUMrhoAB=NUMrhoAB/PDtotal
989             
990 !*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
991 !c     fin oscillation des gouttes
992 !c            print*,'youve been there'
993
994          elseif(oscil.EQ.2) then
995                
996 !C     AMPLITUDE MATRIX [Eqs. (2)-(4) of Ref. 6]
997             CALL AMPL(NMAX,LAM,THET0,THET,PHI0,PHI,ALPHA,BETA,&
998                 S11,S12,S21,S22)    
999             
1000             NUMrhoAB=NUMrhoAB+S11*CONJG(S22)
1001             S11carre=S11carre+((DBLE(S11))**2+(AIMAG(S11))**2)
1002             S22carre=S22carre+((DBLE(S22))**2+(AIMAG(S22))**2)
1003          endif
1004             
1005 !C     ZDR=|Shh|2/|Svv|2=|HH|2/|VV|2
1006 !C     ZDRtotal=Somme(|Shh|2)/Somme(|Svv|2)
1007          NUMZ=(DBLE(S22))**2+(AIMAG(S22))**2
1008          DENZ=(DBLE(S11))**2+(AIMAG(S11))**2
1009 !c     ZDR=NUMZ/DENZ
1010 !*************
1011 !c     ZDR=10*LOG10(ZDR)
1012 !*************
1013 !c     PRINT*,'ZDR (DB) =',ZDR
1014 !               
1015 !C     RHO
1016                
1017                
1018 !c         DENrho=S11carre*S22carre
1019          
1020 !c         RHO=NUMrhoAB/(DENrho**0.5)
1021          
1022 !c         norme=((DBLE(NUMrhoAB))**2+(AIMAG(NUMrhoAB))**2)**0.5
1023          
1024 !C************
1025 !c         RHOhv=((DBLE(RHO))**2+(AIMAG(RHO))**2)**0.5
1026          
1027 !C***********
1028
1029
1030 !C     Calcul de ZH et ZV a partir de la matrice d'amplitude
1031                
1032 !C         ZHH=(LAM**4)*(((DBLE(S22))**2+(AIMAG(S22))**2)*4*P)/(P**5*K2)
1033          
1034 !C         ZVV=(LAM**4)*(((DBLE(S11))**2+(AIMAG(S11))**2)*4*P)/(P**5*K2)
1035           
1036          QBACKHH=16.*ABS(S22)**2/Deq**2
1037          QBACKVV=16.*ABS(S11)**2/Deq**2
1038 !c         print*,Deq,qback
1039
1040 !c     print *, 10*LOG10(ZHH)+180
1041           
1042 !C     DELTA EN RADIANS
1043 !****************
1044          NUMD=AIMAG(S22*CONJG(S11))
1045          DEND=DBLE(S22*CONJG(S11))  
1046          DELTA=NUMD/DEND
1047          
1048 !**************
1049          DELTA=ATAN(DELTA)
1050 !C     DELTA EN DEGRES
1051          DELTA=DELTA*180/P
1052          
1053          
1054          NTotal=Ntotal+1
1055          
1056          NUMTZ=NUMTZ+NUMZ
1057          DENTZ=DENTZ+DENZ
1058          
1059          
1060          NUMTD=NUMTD+NUMD
1061          DENTD=DENTD+DEND
1062          
1063 !c         ZHHT=ZHHT+ZHH
1064 !c         ZVVT=ZVVT+ZVV
1065          
1066          NUMTrhoAB=NUMTrhoAB+NUMrhoAB
1067          
1068          MYS11cr=MYS11cr+S11carre
1069          
1070          MYS22cr=MYS22cr+S22carre
1071 !***********************************************
1072             
1073             
1074 !C     PHASE MATRIX [Eqs. (13)-(29) of Ref. 6]
1075 !            
1076 !c      Z11=0.5D0*(S11*CONJG(S11)+S12*CONJG(S12)
1077 !c     &     +S21*CONJG(S21)+S22*CONJG(S22))
1078 !c      Z12=0.5D0*(S11*CONJG(S11)-S12*CONJG(S12)
1079 !c     &     +S21*CONJG(S21)-S22*CONJG(S22))
1080 !c      Z13=-S11*CONJG(S12)-S22*CONJG(S21)
1081 !c      Z14=(0D0,1D0)*(S11*CONJG(S12)-S22*CONJG(S21))
1082 !c      Z21=0.5D0*(S11*CONJG(S11)+S12*CONJG(S12)
1083 !c     &     -S21*CONJG(S21)-S22*CONJG(S22))
1084 !c      Z22=0.5D0*(S11*CONJG(S11)-S12*CONJG(S12)
1085 !c     &     -S21*CONJG(S21)+S22*CONJG(S22))
1086 !c      Z23=-S11*CONJG(S12)+S22*CONJG(S21)
1087 !c      Z24=(0D0,1D0)*(S11*CONJG(S12)+S22*CONJG(S21))
1088 !c      Z31=-S11*CONJG(S21)-S22*CONJG(S12)
1089 !c      Z32=-S11*CONJG(S21)+S22*CONJG(S12)
1090 !c      Z33=S11*CONJG(S22)+S12*CONJG(S21)
1091 !c      Z34=(0D0,-1D0)*(S11*CONJG(S22)+S21*CONJG(S12))
1092 !c      Z41=(0D0,1D0)*(S21*CONJG(S11)+S22*CONJG(S12))
1093 !c      Z42=(0D0,1D0)*(S21*CONJG(S11)-S22*CONJG(S12))
1094 !c      Z43=(0D0,-1D0)*(S22*CONJG(S11)-S12*CONJG(S21))
1095 !c      Z44=S22*CONJG(S11)-S12*CONJG(S21)
1096       
1097 !C     WRITE (10,5000) 
1098 !C     WRITE (10,5001) Z11,Z12,Z13,Z14
1099 !C     WRITE (10,5001) Z21,Z22,Z23,Z24
1100 !C     WRITE (10,5001) Z31,Z32,Z33,Z34
1101 !C     WRITE (10,5001) Z41,Z42,Z43,Z44
1102             
1103 ! 5000 FORMAT ('PHASE MATRIX')
1104 ! 5001 FORMAT (4F11.7)
1105             
1106 !****mclock gives the time of execution (CPU)in centieme of second
1107 !c            ITIME=MCLOCK()
1108 !c            TIME=FLOAT(ITIME)/6000D0
1109 !C     WRITE(10,1001)TIME
1110 !c 1001       FORMAT (' time =',F8.2,' min')
1111             
1112 !*************************************************************
1113             
1114 !***   enddo de la boucle des gouttes
1115 !                  
1116 !c     ZHHT=ZHHT/NTotal
1117 !c     ZVVT=ZVVT/NTotal
1118 !         
1119 !C     ZH  et ZV en mm6.m-3
1120 !         
1121 !c      ZHHT=ZHHT*1D18
1122 !c      ZVVT=ZVVT*1D18
1123 !         
1124 !C     ZH et ZV en dBZ
1125 !         
1126 !c      ZHHT=10*LOG10(ZHHT)
1127 !c      ZVVT=10*LOG10(ZVVT)
1128       
1129       ZDRT= NUMTZ/DENTZ
1130       ZDRT=10*LOG10(ZDRT)
1131       
1132       DELTAT= NUMTD/DENTD
1133       DELTAT=ATAN(DELTAT)
1134       DELTAT=DELTAT*180/P
1135       
1136       NUMTrhoAB=NUMTrhoAB/NTOTAL
1137       
1138       MYS11cr=MYS11cr/NTOTAL
1139       MYS22cr=MYS22cr/NTOTAL
1140 !c      DENTrhoAB=MYS11cr*MYS22cr
1141 !c      normTOT=((DBLE(NUMTrhoAB))**2+(AIMAG(NUMTrhoAB))**2)**0.5
1142       
1143       endif
1144       
1145 !c      RHOT=NUMTrhoAB/((DENTrhoAB)**0.5)
1146       
1147 !c      RHOhvT=((DBLE(RHOT))**2+(AIMAG(RHOT))**2)**0.5
1148       
1149 !c     ecriture du fichier
1150                
1151 !c     ENDDO
1152       
1153 !c     CLOSE(8)
1154       
1155       ENDDO
1156       
1157       RETURN
1158       END
1159
1160 !************fin du main program***********************************
1161
1162 !C********************************************************************
1163 !C   CALCULATION OF THE AMPLITUDE MATRIX                             *
1164 !C******************************************************************** 
1165       SUBROUTINE AMPL(NMAX,DLAM,TL,TL1,PL,PL1,ALPHA,BETA,&
1166                       VV,VH,HV,HH)  
1167
1168       USE MODD_TMAT,ONLY :XRT11,XRT12,XRT21,XRT22,XIT11,XIT12,&
1169                           XIT21,XIT22,NPN1,NPN4,NPN6
1170
1171 !c      INCLUDE 'ampld.par.f'
1172 !!      Parameter (NPN1=200,NPN4=NPN1, NPN6=NPN4+1)
1173       IMPLICIT REAL*8 (A-B,D-H,O-Z), COMPLEX*16 (C)
1174       REAL*8 AL(3,2),AL1(3,2),AP(2,3),AP1(2,3),B(3,3),&
1175             R(2,2),R1(2,2),C(3,2),CA,CB,CT,CP,CTP,CPP,CT1,CP1,&
1176             CTP1,CPP1
1177 !C       REAL*8 ZDR,NUM,DEN
1178       REAL*8 DV1(NPN6),DV2(NPN6),DV01(NPN6),DV02(NPN6)
1179 !!      REAL*4 TR11(NPN6,NPN4,NPN4),TR12(NPN6,NPN4,NPN4),&
1180 !!          TR21(NPN6,NPN4,NPN4),TR22(NPN6,NPN4,NPN4),&
1181 !!          TI11(NPN6,NPN4,NPN4),TI12(NPN6,NPN4,NPN4),&
1182 !!          TI21(NPN6,NPN4,NPN4),TI22(NPN6,NPN4,NPN4)
1183       COMPLEX*16 CAL(NPN4,NPN4),VV,VH,HV,HH
1184 !!      COMMON /TMAT/ TR11,TR12,TR21,TR22,TI11,TI12,TI21,TI22
1185
1186       IF (ALPHA.LT.0D0.OR.ALPHA.GT.360D0.OR.&
1187           BETA.LT.0D0.OR.BETA.GT.180D0.OR.&
1188           TL.LT.0D0.OR.TL.GT.180D0.OR.&
1189           TL1.LT.0D0.OR.TL1.GT.180D0.OR.&
1190           PL.LT.0D0.OR.PL.GT.360D0.OR.&
1191           PL1.LT.0D0.OR.PL1.GT.360D0) THEN 
1192 !C      WRITE (10,2000)
1193          STOP
1194       ENDIF  
1195
1196 ! 2000 FORMAT ('AN ANGULAR PARAMETER IS OUTSIDE ITS',&
1197 !             ' ALLOWABLE RANGE')
1198       PIN=ACOS(-1D0)
1199
1200       PIN2=PIN*0.5D0
1201       PI=PIN/180D0
1202       ALPH=ALPHA*PI
1203       BET=BETA*PI
1204       THETL=TL*PI
1205       PHIL=PL*PI
1206       THETL1=TL1*PI
1207       PHIL1=PL1*PI
1208
1209       EPSILON=1D-7
1210       IF (THETL.LT.PIN2) THETL=THETL+EPSILON
1211       IF (THETL.GT.PIN2) THETL=THETL-EPSILON
1212       IF (THETL1.LT.PIN2) THETL1=THETL1+EPSILON
1213       IF (THETL1.GT.PIN2) THETL1=THETL1-EPSILON
1214       IF (PHIL.LT.PIN) PHIL=PHIL+EPSILON
1215       IF (PHIL.GT.PIN) PHIL=PHIL-EPSILON
1216       IF (PHIL1.LT.PIN) PHIL1=PHIL1+EPSILON
1217       IF (PHIL1.GT.PIN) PHIL1=PHIL1-EPSILON
1218       
1219 !C_____________COMPUTE THETP, PHIP, THETP1, AND PHIP1, EQS. (8), (19), AND (20)
1220
1221       CB=COS(BET)
1222       SB=SIN(BET)
1223       CT=COS(THETL)
1224       ST=SIN(THETL)
1225       CP=COS(PHIL-ALPH)
1226       SP=SIN(PHIL-ALPH)
1227       CTP=CT*CB+ST*SB*CP
1228       THETP=ACOS(CTP)
1229       CPP=CB*ST*CP-SB*CT
1230       SPP=ST*SP
1231       PHIP=ATAN(SPP/CPP)
1232
1233       IF (PHIP.GT.0D0.AND.SP.LT.0D0) PHIP=PHIP+PIN
1234       IF (PHIP.LT.0D0.AND.SP.GT.0D0) PHIP=PHIP+PIN
1235       IF (PHIP.LT.0D0) PHIP=PHIP+2D0*PIN
1236
1237       CT1=COS(THETL1)
1238       ST1=SIN(THETL1)
1239       CP1=COS(PHIL1-ALPH)
1240       SP1=SIN(PHIL1-ALPH)
1241       CTP1=CT1*CB+ST1*SB*CP1
1242       THETP1=ACOS(CTP1)
1243       CPP1=CB*ST1*CP1-SB*CT1
1244       SPP1=ST1*SP1
1245       PHIP1=ATAN(SPP1/CPP1)
1246       IF (PHIP1.GT.0D0.AND.SP1.LT.0D0) PHIP1=PHIP1+PIN
1247       IF (PHIP1.LT.0D0.AND.SP1.GT.0D0) PHIP1=PHIP1+PIN
1248       IF (PHIP1.LT.0D0) PHIP1=PHIP1+2D0*PIN
1249
1250 !C____________COMPUTE MATRIX BETA, EQ. (21)
1251
1252       CA=COS(ALPH)
1253       SA=SIN(ALPH)
1254       B(1,1)=CA*CB
1255       B(1,2)=SA*CB
1256       B(1,3)=-SB
1257       B(2,1)=-SA
1258       B(2,2)=CA
1259       B(2,3)=0D0
1260       B(3,1)=CA*SB
1261       B(3,2)=SA*SB
1262       B(3,3)=CB
1263
1264 !C____________COMPUTE MATRICES AL AND AL1, EQ. (14) 
1265
1266       CP=COS(PHIL)
1267       SP=SIN(PHIL)
1268       CP1=COS(PHIL1)
1269       SP1=SIN(PHIL1)
1270       AL(1,1)=CT*CP
1271       AL(1,2)=-SP
1272       AL(2,1)=CT*SP
1273       AL(2,2)=CP
1274       AL(3,1)=-ST
1275       AL(3,2)=0D0
1276       AL1(1,1)=CT1*CP1
1277       AL1(1,2)=-SP1
1278       AL1(2,1)=CT1*SP1
1279       AL1(2,2)=CP1
1280       AL1(3,1)=-ST1
1281       AL1(3,2)=0D0
1282
1283 !C____________COMPUTE MATRICES AP^(-1) AND AP1^(-1), EQ. (15) 
1284
1285       CT=CTP
1286       ST=SIN(THETP) 
1287       CP=COS(PHIP)
1288       SP=SIN(PHIP)
1289       CT1=CTP1
1290       ST1=SIN(THETP1)
1291       CP1=COS(PHIP1)
1292       SP1=SIN(PHIP1)
1293       AP(1,1)=CT*CP
1294       AP(1,2)=CT*SP
1295       AP(1,3)=-ST  
1296       AP(2,1)=-SP
1297       AP(2,2)=CP 
1298       AP(2,3)=0D0
1299       AP1(1,1)=CT1*CP1
1300       AP1(1,2)=CT1*SP1
1301       AP1(1,3)=-ST1   
1302       AP1(2,1)=-SP1
1303       AP1(2,2)=CP1 
1304       AP1(2,3)=0D0
1305
1306 !C____________COMPUTE MATRICES R AND R^(-1), EQ. (13)
1307       DO I=1,3
1308           DO J=1,2
1309           X=0D0
1310               DO K=1,3
1311               X=X+B(I,K)*AL(K,J)
1312               ENDDO
1313           C(I,J)=X
1314           ENDDO
1315       ENDDO
1316
1317       DO I=1,2
1318           DO J=1,2
1319           X=0D0
1320               DO K=1,3
1321               X=X+AP(I,K)*C(K,J)
1322               ENDDO
1323           R(I,J)=X
1324           ENDDO
1325       ENDDO
1326
1327       DO I=1,3
1328           DO J=1,2
1329           X=0D0
1330               DO K=1,3
1331               X=X+B(I,K)*AL1(K,J)
1332               ENDDO
1333           C(I,J)=X
1334           ENDDO
1335       ENDDO
1336
1337       DO I=1,2
1338           DO J=1,2
1339           X=0D0
1340               DO K=1,3
1341               X=X+AP1(I,K)*C(K,J)
1342               ENDDO
1343           R1(I,J)=X
1344           ENDDO
1345       ENDDO
1346
1347       D=1D0/(R1(1,1)*R1(2,2)-R1(1,2)*R1(2,1))
1348       X=R1(1,1)
1349       R1(1,1)=R1(2,2)*D
1350       R1(1,2)=-R1(1,2)*D
1351       R1(2,1)=-R1(2,1)*D
1352       R1(2,2)=X*D
1353
1354       CI=(0D0,1D0)
1355
1356       DO NN=1,NMAX
1357           DO N=1,NMAX
1358           CN=CI**(NN-N-1)
1359           DNN=FLOAT((2*N+1)*(2*NN+1)) 
1360           DNN=DNN/FLOAT( N*NN*(N+1)*(NN+1) ) 
1361           RN=SQRT(DNN)
1362           CAL(N,NN)=CN*RN
1363           ENDDO
1364       ENDDO
1365     
1366       DCTH0=CTP
1367       DCTH=CTP1 
1368       PH=PHIP1-PHIP
1369       VV=(0D0,0D0)
1370       VH=(0D0,0D0)
1371       HV=(0D0,0D0)
1372       HH=(0D0,0D0)
1373
1374       DO M=0,NMAX
1375          M1=M+1
1376          NMIN=MAX(M,1)
1377          
1378          CALL VIGAMPL(DCTH, NMAX, M, DV1, DV2)
1379          CALL VIGAMPL(DCTH0, NMAX, M, DV01, DV02)
1380     
1381          FC=2D0*COS(M*PH)
1382          FS=2D0*SIN(M*PH)
1383          
1384          DO NN=NMIN,NMAX
1385             DV1NN=M*DV01(NN)
1386             DV2NN=DV02(NN)
1387             DO N=NMIN,NMAX
1388                DV1N=M*DV1(N)
1389                DV2N=DV2(N)
1390                
1391                CT11=CMPLX(XRT11(M1,N,NN),XIT11(M1,N,NN))
1392                CT22=CMPLX(XRT22(M1,N,NN),XIT22(M1,N,NN))
1393                
1394                IF (M.EQ.0) THEN
1395                   
1396                   CN=CAL(N,NN)*DV2N*DV2NN
1397                   
1398                   VV=VV+CN*CT22  
1399                   HH=HH+CN*CT11
1400                   
1401                ELSE
1402                   
1403                   CT12=CMPLX(XRT12(M1,N,NN),XIT12(M1,N,NN))
1404                   CT21=CMPLX(XRT21(M1,N,NN),XIT21(M1,N,NN))
1405                   
1406                   CN1=CAL(N,NN)*FC
1407                   CN2=CAL(N,NN)*FS
1408                   
1409                   D11=DV1N*DV1NN
1410                   D12=DV1N*DV2NN
1411                   D21=DV2N*DV1NN
1412                   D22=DV2N*DV2NN
1413                   
1414                   VV=VV+(CT11*D11+CT21*D21+CT12*D12+CT22*D22)*CN1   
1415                   
1416                   VH=VH+(CT11*D12+CT21*D22+CT12*D11+CT22*D21)*CN2
1417                   
1418                   HV=HV-(CT11*D21+CT21*D11+CT12*D22+CT22*D12)*CN2   
1419                   
1420                   HH=HH+(CT11*D22+CT21*D12+CT12*D21+CT22*D11)*CN1      
1421                ENDIF
1422             ENDDO
1423          ENDDO
1424       ENDDO
1425
1426       DK=2D0*PIN/DLAM
1427       VV=VV/DK
1428       VH=VH/DK
1429       HV=HV/DK
1430       HH=HH/DK
1431       CVV=VV*R(1,1)+VH*R(2,1)
1432       CVH=VV*R(1,2)+VH*R(2,2)
1433       CHV=HV*R(1,1)+HH*R(2,1)
1434       CHH=HV*R(1,2)+HH*R(2,2)
1435       VV=R1(1,1)*CVV+R1(1,2)*CHV
1436       VH=R1(1,1)*CVH+R1(1,2)*CHH
1437       HV=R1(2,1)*CVV+R1(2,2)*CHV
1438       HH=R1(2,1)*CVH+R1(2,2)*CHH
1439       
1440 !C      WRITE (10,1005) TL,TL1,PL,PL1,ALPHA,BETA 
1441 !C      WRITE (10,1006)
1442 !C      WRITE (10,1101) VV
1443 !C      WRITE (10,1102) VH
1444 !C      WRITE (10,1103) HV
1445 !C      WRITE (10,1104) HH
1446 ! 1101 FORMAT ('S11=',D11.5,' + i*',D11.5)
1447 ! 1102 FORMAT ('S12=',D11.5,' + i*',D11.5)
1448 ! 1103 FORMAT ('S21=',D11.5,' + i*',D11.5)
1449 ! 1104 FORMAT ('S22=',D11.5,' + i*',D11.5)
1450 !
1451 !C     ZDR=|Shh|2/|Svv|2=|HH|2/|VV|2
1452 !C     ZDRtotal=Somme(|Shh|2)/Somme(|Svv|2)
1453 !C        NUM=(REAL(HH))**2+(AIMAG(HH))**2
1454 !C        DEN=(REAL(VV))**2+(AIMAG(VV))**2
1455 !C        ZDR=NUM/DEN
1456 !C       PRINT*,'ZDR=',ZDR
1457 !C        NUMT=NUMT+NUM
1458 !C        DENT=DENT+DEN
1459 !C        ZDRT=NUMT/DENT
1460         
1461
1462
1463
1464 ! 1005 FORMAT ('thet0=',F6.2,'  thet=',F6.2,'  phi0=',F6.2,&
1465 !             '  phi=',F6.2,'  alpha=',F6.2,'  beta=',F6.2)
1466 ! 1006 FORMAT ('AMPLITUDE MATRIX')
1467       RETURN
1468       END
1469       
1470 !C*****************************************************************
1471 !C     Calculation of the functions                               *
1472 !C     DV1(N)=dvig(0,m,n,arccos x)/sin(arccos x)                  *
1473 !C     and                                                        *
1474 !C     DV2(N)=[d/d(arccos x)] dvig(0,m,n,arccos x)                *
1475 !C     1.LE.N.LE.NMAX                                             *
1476 !C     0.LE.X.LE.1                                                *
1477 !C*****************************************************************
1478       SUBROUTINE VIGAMPL(X, NMAX, M, DV1, DV2)
1479
1480       USE MODD_TMAT,   ONLY: NPN1,NPN4,NPN6
1481 !!      Parameter (NPN1=200,NPN4=NPN1, NPN6=NPN4+1)
1482       IMPLICIT REAL*8 (A-H,O-Z)
1483       REAL*8 DV1(NPN6), DV2(NPN6)
1484
1485       DO N=1,NMAX
1486       DV1(N)=0D0
1487       DV2(N)=0D0
1488       ENDDO
1489
1490       DX=ABS(X)
1491
1492       IF (ABS(1D0-DX).GT.1D-10) THEN
1493       A=1D0
1494       QS=SQRT(1D0-X*X)
1495       QS1=1D0/QS
1496       DSI=QS1
1497            IF (M.EQ.0) THEN
1498            D1=1D0
1499            D2=X  
1500                DO N=1,NMAX  
1501                QN=FLOAT(N)
1502                QN1=FLOAT(N+1)
1503                QN2=FLOAT(2*N+1)
1504                D3=(QN2*X*D2-QN*D1)/QN1 
1505                DER=QS1*(QN1*QN/QN2)*(-D1+D3)
1506                DV1(N)=D2*DSI
1507                DV2(N)=DER
1508                D1=D2
1509                D2=D3
1510                ENDDO
1511            ELSE
1512            QMM=FLOAT(M*M)
1513                DO I=1,M
1514                I2=I*2
1515                A=A*SQRT(FLOAT(I2-1)/FLOAT(I2))*QS
1516                ENDDO
1517            D1=0D0
1518            D2=A 
1519                DO N=M,NMAX
1520                QN=FLOAT(N)
1521                QN2=FLOAT(2*N+1)
1522                QN1=FLOAT(N+1)
1523                QNM=SQRT(QN*QN-QMM)
1524                QNM1=SQRT(QN1*QN1-QMM)
1525                D3=(QN2*X*D2-QNM*D1)/QNM1
1526                DER=QS1*(-QN1*QNM*D1+QN*QNM1*D3)/QN2
1527                DV1(N)=D2*DSI
1528                DV2(N)=DER
1529                D1=D2
1530                D2=D3
1531                ENDDO
1532            ENDIF
1533       ELSE
1534           IF (M.EQ.1) THEN
1535
1536               DO N=1,NMAX
1537               DN=FLOAT(N*(N+1))
1538               DN=0.5D0*SQRT(DN)
1539                   IF (X.LT.0D0) DN=DN*(-1)**(N+1)
1540               DV1(N)=DN
1541                   IF (X.LT.0D0) DN=-DN
1542               DV2(N)=DN
1543               ENDDO
1544            ENDIF
1545       ENDIF
1546
1547       RETURN
1548       END 
1549
1550 !C**********************************************************************
1551
1552       SUBROUTINE CONST(NGAUSS,NMAX,X,W,AN,ANN,S,SS)
1553
1554       USE MODD_TMAT,   ONLY: NPN1,NPNG1,NPNG2
1555       IMPLICIT REAL*8 (A-H,O-Z)
1556 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1)
1557       REAL*8 X(NPNG2),W(NPNG2),&
1558              S(NPNG2),SS(NPNG2),&
1559              AN(NPN1),ANN(NPN1,NPN1),DD(NPN1)
1560  
1561       DO N=1,NMAX
1562          NN=N*(N+1)
1563          AN(N)=FLOAT(NN)
1564          D=SQRT(FLOAT(2*N+1)/FLOAT(NN))
1565          DD(N)=D
1566          DO N1=1,N
1567             DDD=D*DD(N1)*0.5D0
1568             ANN(N,N1)=DDD
1569             ANN(N1,N)=DDD
1570          ENDDO
1571       ENDDO
1572       
1573       NG=2*NGAUSS
1574       
1575       CALL GAUSS(NG,0,0,X,W)
1576       
1577       DO I=1,NGAUSS
1578          Y=X(I)
1579          Y=1D0/(1D0-Y*Y)
1580          SS(I)=Y
1581          SS(NG-I+1)=Y
1582          Y=SQRT(Y)
1583          S(I)=Y
1584          S(NG-I+1)=Y
1585       ENDDO
1586       
1587       RETURN
1588       END
1589  
1590 !C**********************************************************************
1591  
1592       SUBROUTINE VARY(LAM,MRR,MRI,A,EPS,NGAUSS,X,P,PPI,PIR,PII,&
1593                      R,DR,DDR,DRR,DRI,NMAX)
1594        USE MODD_TMAT,   ONLY: NPN1,NPNG1,NPNG2
1595               
1596 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1)
1597       IMPLICIT REAL*8 (A-H,O-Z)
1598       REAL*8  X(NPNG2),R(NPNG2),DR(NPNG2),MRR,MRI,LAM,&
1599              Z(NPNG2),ZR(NPNG2),ZI(NPNG2),&
1600              DDR(NPNG2),DRR(NPNG2),DRI(NPNG2)
1601              
1602 !!      COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1603
1604       NG=NGAUSS*2
1605
1606 !c      IF (NP.GT.0) CALL RSP2(X,NG,A,EPS,NP,R,DR)
1607 !C      IF (NP.EQ.-1) 
1608       CALL RSP1(X,NG,NGAUSS,A,EPS,R,DR)
1609 !c      IF (NP.EQ.-3) CALL RSP4(X,NG,A,R,DR)
1610 !*****PI=2*PI/lambda
1611       PI=P*2D0/LAM 
1612 !*****PPI=(2*PI/lamda)**2  
1613       PPI=PI*PI
1614       PIR=PPI*MRR
1615       PII=PPI*MRI
1616       V=1D0/(MRR*MRR+MRI*MRI)
1617       PRR=MRR*V
1618       PRI=-MRI*V
1619       TA=0D0
1620
1621
1622       DO I=1,NG
1623          VV=SQRT(R(I))
1624          V=VV*PI
1625          TA=MAX(TA,V)
1626          VV=1D0/V
1627          DDR(I)=VV
1628          DRR(I)=PRR*VV
1629          DRI(I)=PRI*VV
1630          V1=V*MRR
1631          V2=V*MRI
1632          Z(I)=V
1633          ZR(I)=V1
1634          ZI(I)=V2
1635       ENDDO
1636
1637
1638 !C      IF (NMAX.GT.NPN1) WRITE (10,9000) NMAX,NPN1 
1639       IF (NMAX.GT.NPN1) STOP
1640
1641  9000 FORMAT(' NMAX = ',I2,', i.e., greater than ',I3)
1642       TB=TA*SQRT(MRR*MRR+MRI*MRI)
1643       TB=MAX(TB,FLOAT(NMAX))
1644       NNMAX1=1.2D0*SQRT(MAX(TA,FLOAT(NMAX)))+3D0
1645       NNMAX2=(TB+4D0*(TB**0.33333D0)+1.2D0*SQRT(TB))
1646       NNMAX2=NNMAX2-NMAX+5
1647       CALL BESS(Z,ZR,ZI,NG,NMAX,NNMAX1,NNMAX2)
1648
1649       RETURN
1650       END
1651  
1652 !C**********************************************************************
1653  
1654       SUBROUTINE RSP1(X,NG,NGAUSS,REV,EPS,R,DR)
1655       IMPLICIT REAL*8 (A-H,O-Z)
1656       REAL*8 X(NG),R(NG),DR(NG)
1657
1658       A=REV*EPS**(1D0/3D0)
1659       AA=A*A
1660       EE=EPS*EPS
1661       EE1=EE-1D0
1662
1663       DO I=1,NGAUSS
1664          C=X(I)
1665          CC=C*C
1666          SS=1D0-CC
1667          S=SQRT(SS)
1668          RR=1D0/(SS+EE*CC)
1669          R(I)=AA*RR
1670          R(NG-I+1)=R(I)
1671          DR(I)=RR*C*S*EE1
1672          DR(NG-I+1)=-DR(I)
1673       ENDDO
1674
1675       RETURN
1676       END
1677  
1678
1679 !C*********************************************************************
1680  
1681       SUBROUTINE BESS (X,XR,XI,NG,NMAX,NNMAX1,NNMAX2)
1682
1683       USE MODD_TMAT,ONLY : XJ,XY,XJR,XJI,XDJ,XDY,XDJR,XDJI,NPN1,NPNG1,NPNG2
1684 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1)
1685       IMPLICIT REAL*8 (A-H,O-Z)
1686       REAL*8 X(NG),XR(NG),XI(NG),&
1687 !!             J(NPNG2,NPN1),Y(NPNG2,NPN1),JR(NPNG2,NPN1),&
1688 !!             JI(NPNG2,NPN1),DJ(NPNG2,NPN1),DY(NPNG2,NPN1),&
1689 !!             DJR(NPNG2,NPN1),DJI(NPNG2,NPN1),&
1690              AJ(NPN1),AY(NPN1),AJR(NPN1),AJI(NPN1),&
1691              ADJ(NPN1),ADY(NPN1),ADJR(NPN1),&
1692              ADJI(NPN1)
1693 !!      COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1694  
1695       DO I=1,NG
1696       XX=X(I)
1697
1698       CALL RJB(XX,AJ,ADJ,NMAX,NNMAX1)
1699       CALL RYB(XX,AY,ADY,NMAX)
1700
1701       YR=XR(I)
1702       YI=XI(I)
1703
1704       CALL CJB(YR,YI,AJR,AJI,ADJR,ADJI,NMAX,NNMAX2)
1705
1706           DO N=1,NMAX
1707           XJ(I,N)=AJ(N)
1708           XY(I,N)=AY(N)
1709           XJR(I,N)=AJR(N)
1710           XJI(I,N)=AJI(N)
1711           XDJ(I,N)=ADJ(N)
1712           XDY(I,N)=ADY(N)
1713           XDJR(I,N)=ADJR(N)
1714           XDJI(I,N)=ADJI(N)
1715           ENDDO
1716       ENDDO
1717
1718       RETURN
1719       END
1720  
1721 !C**********************************************************************
1722  
1723       SUBROUTINE RJB(X,Y,U,NMAX,NNMAX)
1724       IMPLICIT REAL*8 (A-H,O-Z)
1725       REAL*8 Y(NMAX),U(NMAX),Z(800)
1726
1727       L=NMAX+NNMAX
1728       XX=1D0/X
1729       Z(L)=1D0/(FLOAT(2*L+1)*XX)
1730       L1=L-1
1731
1732       DO I=1,L1
1733       I1=L-I
1734       Z(I1)=1D0/(FLOAT(2*I1+1)*XX-Z(I1+1))
1735       ENDDO
1736
1737       Z0=1D0/(XX-Z(1))
1738       Y0=Z0*COS(X)*XX
1739       Y1=Y0*Z(1)
1740       U(1)=Y0-Y1*XX
1741       Y(1)=Y1
1742
1743       DO I=2,NMAX
1744       YI1=Y(I-1)
1745       YI=YI1*Z(I)
1746       U(I)=YI1-FLOAT(I)*YI*XX
1747       Y(I)=YI
1748       ENDDO
1749
1750       RETURN
1751       END
1752  
1753 !C**********************************************************************
1754  
1755       SUBROUTINE RYB(X,Y,V,NMAX)
1756       IMPLICIT REAL*8 (A-H,O-Z)
1757       REAL*8 Y(NMAX),V(NMAX)
1758
1759       C=COS(X)
1760       S=SIN(X)
1761       X1=1D0/X
1762       X2=X1*X1
1763       X3=X2*X1
1764       Y1=-C*X2-S*X1
1765       Y(1)=Y1
1766       Y(2)=(-3D0*X3+X1)*C-3D0*X2*S
1767       NMAX1=NMAX-1
1768
1769       DO I=2,NMAX1
1770       Y(I+1)=FLOAT(2*I+1)*X1*Y(I)-Y(I-1)
1771       ENDDO
1772
1773       V(1)=-X1*(C+Y1)
1774       DO  I=2,NMAX
1775       V(I)=Y(I-1)-FLOAT(I)*X1*Y(I)
1776       ENDDO
1777
1778       RETURN
1779       END
1780  
1781 !C**********************************************************************
1782 !C                                                                     *
1783 !C   CALCULATION OF SPHERICAL BESSEL FUNCTIONS OF THE FIRST KIND       *
1784 !C   J=JR+I*JI OF COMPLEX ARGUMENT X=XR+I*XI OF ORDERS FROM 1 TO NMAX  *
1785 !C   BY USING BACKWARD RECURSION. PARAMETR NNMAX DETERMINES NUMERICAL  *
1786 !C   ACCURACY. U=UR+I*UI - FUNCTION (1/X)(D/DX)(X*J(X))                *
1787 !C                                                                     *
1788 !C**********************************************************************
1789  
1790       SUBROUTINE CJB (XR,XI,YR,YI,UR,UI,NMAX,NNMAX)
1791
1792       USE MODD_TMAT,ONLY:NPN1
1793 !!      Parameter (NPN1=200)
1794       IMPLICIT REAL*8 (A-H,O-Z)
1795       REAL*8 YR(NMAX),YI(NMAX),UR(NMAX),UI(NMAX)
1796       REAL*8 CYR(NPN1),CYI(NPN1),CZR(1200),CZI(1200)
1797
1798       L=NMAX+NNMAX
1799       XRXI=1D0/(XR*XR+XI*XI)
1800       CXXR=XR*XRXI
1801       CXXI=-XI*XRXI 
1802       QF=1D0/FLOAT(2*L+1)
1803       CZR(L)=XR*QF
1804       CZI(L)=XI*QF
1805       L1=L-1
1806
1807       DO I=1,L1
1808       I1=L-I
1809       QF=FLOAT(2*I1+1)
1810       AR=QF*CXXR-CZR(I1+1)
1811       AI=QF*CXXI-CZI(I1+1)
1812       ARI=1D0/(AR*AR+AI*AI)
1813       CZR(I1)=AR*ARI
1814       CZI(I1)=-AI*ARI
1815       ENDDO   
1816
1817       AR=CXXR-CZR(1)
1818       AI=CXXI-CZI(1)
1819       ARI=1D0/(AR*AR+AI*AI)
1820       CZ0R=AR*ARI
1821       CZ0I=-AI*ARI
1822       CR=COS(XR)*COSH(XI)
1823       CI=-SIN(XR)*SINH(XI)
1824       AR=CZ0R*CR-CZ0I*CI
1825       AI=CZ0I*CR+CZ0R*CI
1826       CY0R=AR*CXXR-AI*CXXI
1827       CY0I=AI*CXXR+AR*CXXI
1828       CY1R=CY0R*CZR(1)-CY0I*CZI(1)
1829       CY1I=CY0I*CZR(1)+CY0R*CZI(1)
1830       AR=CY1R*CXXR-CY1I*CXXI
1831       AI=CY1I*CXXR+CY1R*CXXI
1832       CU1R=CY0R-AR
1833       CU1I=CY0I-AI
1834       CYR(1)=CY1R
1835       CYI(1)=CY1I
1836 !c      CUR(1)=CU1R
1837 !c      CUI(1)=CU1I
1838       YR(1)=CY1R
1839       YI(1)=CY1I
1840       UR(1)=CU1R
1841       UI(1)=CU1I
1842
1843       DO I=2,NMAX
1844       QI=FLOAT(I)
1845       CYI1R=CYR(I-1)
1846       CYI1I=CYI(I-1)
1847       CYIR=CYI1R*CZR(I)-CYI1I*CZI(I)
1848       CYII=CYI1I*CZR(I)+CYI1R*CZI(I)
1849       AR=CYIR*CXXR-CYII*CXXI
1850       AI=CYII*CXXR+CYIR*CXXI
1851       CUIR=CYI1R-QI*AR
1852       CUII=CYI1I-QI*AI
1853       CYR(I)=CYIR
1854       CYI(I)=CYII
1855 !c      CUR(I)=CUIR
1856 !c      CUI(I)=CUII
1857       YR(I)=CYIR
1858       YI(I)=CYII
1859       UR(I)=CUIR
1860       UI(I)=CUII
1861       ENDDO   
1862
1863       RETURN
1864       END
1865  
1866 !C**********************************************************************
1867  
1868       SUBROUTINE TMATR0(NGAUSS,X,W,AN,ANN,PPI,PIR,PII,R,DR,DDR,&
1869                        DRR,DRI,NMAX)
1870
1871       USE MODD_TMAT,ONLY : XJ,XY,XJR,XJI,XDJ,XDY,XDJR,XDJI,&
1872                            XQR,XQI,XRGQR,XRGQI,&
1873                            NPN1,NPNG1,NPNG2,NPN2,NPN4,NPN6
1874 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1, NPN2=2*NPN1,&   
1875 !!          NPN4=NPN1,NPN6=NPN4+1)
1876       IMPLICIT REAL*8 (A-H,O-Z)
1877       REAL*8  X(NPNG2),W(NPNG2),AN(NPN1),&
1878              R(NPNG2),DR(NPNG2),SIG(NPN2),&
1879 !!             J(NPNG2,NPN1),Y(NPNG2,NPN1),&
1880 !!             JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1),&
1881 !!             DY(NPNG2,NPN1),DJR(NPNG2,NPN1),&
1882 !!             DJI(NPNG2,NPN1),
1883              DDR(NPNG2),DRR(NPNG2),&
1884              D1(NPNG2,NPN1),D2(NPNG2,NPN1),&
1885              DRI(NPNG2),RR(NPNG2),&
1886              DV1(NPN1),DV2(NPN1)
1887         REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: R11,R12,R21,R22
1888         REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: I11,I12,I21,I22
1889         REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: RG11,RG12,RG21,RG22
1890         REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: IG11,IG12,IG21,IG22
1891
1892      REAL*8  ANN(NPN1,NPN1),&
1893 !!             QR(NPN2,NPN2),QI(NPN2,NPN2),&
1894 !!             RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),&
1895              TQR(NPN2,NPN2),TQI(NPN2,NPN2),&
1896              TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2)
1897              
1898 !!      COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
1899 !!      COMMON /CTT/ QR,QI,RGQR,RGQI
1900
1901       MM1=1
1902       NNMAX=NMAX+NMAX
1903 !c      NG=2*NGAUSS
1904 !c      NGSS=NG
1905 !c      FACTOR=1D0
1906
1907 !c      IF (NCHECK.EQ.1) THEN
1908       NGSS=NGAUSS
1909       FACTOR=2D0
1910 !c      ENDIF
1911
1912       SI=1D0
1913
1914
1915
1916 ! Allocation des tableaux :
1917 ALLOCATE(R11(NPN1,NPN1))
1918 ALLOCATE(R12(NPN1,NPN1))
1919 ALLOCATE(R21(NPN1,NPN1))
1920 ALLOCATE(R22(NPN1,NPN1))
1921 ALLOCATE(I11(NPN1,NPN1))
1922 ALLOCATE(I12(NPN1,NPN1))
1923 ALLOCATE(I21(NPN1,NPN1))
1924 ALLOCATE(I22(NPN1,NPN1))
1925 ALLOCATE(RG11(NPN1,NPN1))
1926 ALLOCATE(RG12(NPN1,NPN1))
1927 ALLOCATE(RG21(NPN1,NPN1))
1928 ALLOCATE(RG22(NPN1,NPN1))
1929 ALLOCATE(IG11(NPN1,NPN1))
1930 ALLOCATE(IG12(NPN1,NPN1))
1931 ALLOCATE(IG21(NPN1,NPN1))
1932 ALLOCATE(IG22(NPN1,NPN1))
1933
1934
1935       DO N=1,NNMAX
1936       SI=-SI
1937       SIG(N)=SI
1938       ENDDO
1939
1940    20 DO I=1,NGAUSS
1941       I1=NGAUSS+I
1942       I2=NGAUSS-I+1
1943       CALL VIG ( X(I1), NMAX, 0, DV1, DV2)
1944           DO N=1,NMAX
1945           SI=SIG(N)
1946           DD1=DV1(N)
1947           DD2=DV2(N)
1948           D1(I1,N)=DD1
1949           D2(I1,N)=DD2
1950           D1(I2,N)=DD1*SI
1951           D2(I2,N)=-DD2*SI
1952           ENDDO
1953       ENDDO
1954
1955    30 DO I=1,NGSS
1956            RR(I)=W(I)*R(I)
1957       ENDDO
1958  
1959       DO N1=MM1,NMAX
1960       AN1=AN(N1)
1961           DO N2=MM1,NMAX
1962           AN2=AN(N2)
1963           AR12=0D0
1964           AR21=0D0
1965           AI12=0D0
1966           AI21=0D0
1967           GR12=0D0
1968           GR21=0D0
1969           GI12=0D0
1970           GI21=0D0
1971               IF (.NOT.(SIG(N1+N2).LT.0D0)) THEN
1972                   DO I=1,NGSS
1973                   D1N1=D1(I,N1)
1974                   D2N1=D2(I,N1)
1975                   D1N2=D1(I,N2)
1976                   D2N2=D2(I,N2)
1977                   A12=D1N1*D2N2
1978                   A21=D2N1*D1N2
1979                   A22=D2N1*D2N2
1980                   AA1=A12+A21
1981
1982                   QJ1=XJ(I,N1)
1983                   QY1=XY(I,N1)
1984                   QJR2=XJR(I,N2)
1985                   QJI2=XJI(I,N2)
1986                   QDJR2=XDJR(I,N2)
1987                   QDJI2=XDJI(I,N2)
1988                   QDJ1=XDJ(I,N1)
1989                   QDY1=XDY(I,N1)
1990  
1991                   C1R=QJR2*QJ1
1992                   C1I=QJI2*QJ1
1993                   B1R=C1R-QJI2*QY1
1994                   B1I=C1I+QJR2*QY1
1995  
1996                   C2R=QJR2*QDJ1
1997                   C2I=QJI2*QDJ1
1998                   B2R=C2R-QJI2*QDY1
1999                   B2I=C2I+QJR2*QDY1
2000
2001                   DDRI=DDR(I)
2002                   C3R=DDRI*C1R
2003                   C3I=DDRI*C1I
2004                   B3R=DDRI*B1R
2005                   B3I=DDRI*B1I
2006
2007                   C4R=QDJR2*QJ1
2008                   C4I=QDJI2*QJ1
2009                   B4R=C4R-QDJI2*QY1
2010                   B4I=C4I+QDJR2*QY1
2011   
2012                   DRRI=DRR(I)
2013                   DRII=DRI(I)
2014                   C5R=C1R*DRRI-C1I*DRII
2015                   C5I=C1I*DRRI+C1R*DRII
2016                   B5R=B1R*DRRI-B1I*DRII
2017                   B5I=B1I*DRRI+B1R*DRII
2018  
2019                   URI=DR(I)
2020                   RRI=RR(I)
2021  
2022                   F1=RRI*A22
2023                   F2=RRI*URI*AN1*A12
2024                   AR12=AR12+F1*B2R+F2*B3R
2025                   AI12=AI12+F1*B2I+F2*B3I
2026                   GR12=GR12+F1*C2R+F2*C3R
2027                   GI12=GI12+F1*C2I+F2*C3I
2028                     
2029                   F2=RRI*URI*AN2*A21
2030                   AR21=AR21+F1*B4R+F2*B5R
2031                   AI21=AI21+F1*B4I+F2*B5I
2032                   GR21=GR21+F1*C4R+F2*C5R
2033                   GI21=GI21+F1*C4I+F2*C5I
2034                   ENDDO
2035               ENDIF
2036           AN12=ANN(N1,N2)*FACTOR
2037           R12(N1,N2)=AR12*AN12
2038           R21(N1,N2)=AR21*AN12
2039           I12(N1,N2)=AI12*AN12
2040           I21(N1,N2)=AI21*AN12
2041           RG12(N1,N2)=GR12*AN12
2042           RG21(N1,N2)=GR21*AN12
2043           IG12(N1,N2)=GI12*AN12
2044           IG21(N1,N2)=GI21*AN12
2045           ENDDO
2046       ENDDO
2047
2048       TPIR=PIR
2049       TPII=PII
2050       TPPI=PPI
2051  
2052       NM=NMAX
2053
2054       DO N1=MM1,NMAX
2055       K1=N1-MM1+1
2056       KK1=K1+NM
2057           DO N2=MM1,NMAX
2058           K2=N2-MM1+1
2059           KK2=K2+NM
2060  
2061           TAR12= I12(N1,N2)
2062           TAI12=-R12(N1,N2)
2063           TGR12= IG12(N1,N2)
2064           TGI12=-RG12(N1,N2)
2065  
2066           TAR21=-I21(N1,N2)
2067           TAI21= R21(N1,N2)
2068           TGR21=-IG21(N1,N2)
2069           TGI21= RG21(N1,N2)
2070  
2071           TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12
2072           TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12
2073           TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12
2074           TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12
2075  
2076           TQR(K1,KK2)=0D0
2077           TQI(K1,KK2)=0D0
2078           TRGQR(K1,KK2)=0D0
2079           TRGQI(K1,KK2)=0D0
2080  
2081           TQR(KK1,K2)=0D0
2082           TQI(KK1,K2)=0D0
2083           TRGQR(KK1,K2)=0D0
2084           TRGQI(KK1,K2)=0D0
2085
2086           TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21
2087           TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21
2088           TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21
2089           TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21
2090           ENDDO
2091       ENDDO
2092  
2093       NNMAX=2*NM
2094
2095       DO N1=1,NNMAX
2096           DO N2=1,NNMAX
2097           XQR(N1,N2)=TQR(N1,N2)
2098           XQI(N1,N2)=TQI(N1,N2)
2099           XRGQR(N1,N2)=TRGQR(N1,N2)
2100           XRGQI(N1,N2)=TRGQI(N1,N2)
2101           ENDDO
2102       ENDDO
2103 DEALLOCATE(R12)
2104 DEALLOCATE(R21)
2105 DEALLOCATE(R22)
2106 DEALLOCATE(I11)
2107 DEALLOCATE(I12)
2108 DEALLOCATE(I21)
2109 DEALLOCATE(I22)
2110 DEALLOCATE(RG11)
2111 DEALLOCATE(RG12)
2112 DEALLOCATE(RG21)
2113 DEALLOCATE(RG22)
2114 DEALLOCATE(IG11)
2115 DEALLOCATE(IG12)
2116 DEALLOCATE(IG21)
2117 DEALLOCATE(IG22)
2118
2119       CALL TT(NMAX)
2120
2121       RETURN
2122       END
2123  
2124 !C**********************************************************************
2125  
2126       SUBROUTINE TMATR(M,NGAUSS,X,W,AN,ANN,S,SS,PPI,PIR,PII,R,DR,DDR,&
2127                       DRR,DRI,NMAX)
2128
2129       USE MODD_TMAT,ONLY : XJ,XY,XJR,XJI,XDJ,XDY,XDJR,XDJI,&
2130                            XQR,XQI,XRGQR,XRGQI,&
2131                            NPN1,NPNG1,NPNG2,NPN2,NPN4,NPN6
2132
2133 !c      INCLUDE 'ampld.par.f'
2134 !!      Parameter (NPN1=200, NPNG1=600, NPNG2=2*NPNG1, NPN2=2*NPN1, & 
2135 !!           NPN4=NPN1, NPN6=NPN4+1)
2136       IMPLICIT REAL*8 (A-H,O-Z)
2137       REAL*8  X(NPNG2),W(NPNG2),AN(NPN1),S(NPNG2),SS(NPNG2),&
2138              R(NPNG2),DR(NPNG2),SIG(NPN2),&
2139 !!             J(NPNG2,NPN1),Y(NPNG2,NPN1),&
2140 !!             JR(NPNG2,NPN1),JI(NPNG2,NPN1),DJ(NPNG2,NPN1),&
2141 !!             DY(NPNG2,NPN1),DJR(NPNG2,NPN1),&
2142 !!             DJI(NPNG2,NPN1),DDR(NPNG2),
2143              DRR(NPNG2),&
2144              D1(NPNG2,NPN1),D2(NPNG2,NPN1),&
2145              DRI(NPNG2),DS(NPNG2),DSS(NPNG2),RR(NPNG2),&
2146              DV1(NPN1),DV2(NPN1)
2147       REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: R11,R12,R21,R22
2148       REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: I11,I12,I21,I22
2149       REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: RG11,RG12,RG21,RG22
2150       REAL(KIND=8),DIMENSION(:,:),ALLOCATABLE :: IG11,IG12,IG21,IG22
2151
2152       REAL*8  ANN(NPN1,NPN1),&
2153 !!             QR(NPN2,NPN2),QI(NPN2,NPN2),&
2154 !!             RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),&
2155              TQR(NPN2,NPN2),TQI(NPN2,NPN2),&
2156              TRGQR(NPN2,NPN2),TRGQI(NPN2,NPN2)
2157 !!      COMMON /CBESS/ J,Y,JR,JI,DJ,DY,DJR,DJI
2158 !!      COMMON /CTT/ QR,QI,RGQR,RGQI
2159
2160       MM1=M
2161       QM=FLOAT(M)
2162       QMM=QM*QM
2163 !c      NG=2*NGAUSS
2164 !c      NGSS=NG
2165 !c      FACTOR=1D0
2166
2167 !c      IF (NCHECK.EQ.1) THEN
2168          NGSS=NGAUSS
2169          FACTOR=2D0
2170 !c      ENDIF
2171
2172       SI=1D0
2173       NM=NMAX+NMAX
2174
2175 ! Allocation des tableaux :
2176 ALLOCATE(R11(NPN1,NPN1))
2177 ALLOCATE(R12(NPN1,NPN1))
2178 ALLOCATE(R21(NPN1,NPN1))
2179 ALLOCATE(R22(NPN1,NPN1))
2180 ALLOCATE(I11(NPN1,NPN1))
2181 ALLOCATE(I12(NPN1,NPN1))
2182 ALLOCATE(I21(NPN1,NPN1))
2183 ALLOCATE(I22(NPN1,NPN1))
2184 ALLOCATE(RG11(NPN1,NPN1))
2185 ALLOCATE(RG12(NPN1,NPN1))
2186 ALLOCATE(RG21(NPN1,NPN1))
2187 ALLOCATE(RG22(NPN1,NPN1))
2188 ALLOCATE(IG11(NPN1,NPN1))
2189 ALLOCATE(IG12(NPN1,NPN1))
2190 ALLOCATE(IG21(NPN1,NPN1))
2191 ALLOCATE(IG22(NPN1,NPN1))
2192
2193
2194
2195       DO N=1,NM
2196          SI=-SI
2197          SIG(N)=SI
2198       ENDDO
2199
2200    20 DO I=1,NGAUSS
2201          I1=NGAUSS+I
2202          I2=NGAUSS-I+1
2203          CALL VIG(X(I1),NMAX,M,DV1,DV2)
2204          DO N=1,NMAX
2205             SI=SIG(N)
2206             DD1=DV1(N)
2207             DD2=DV2(N)
2208             D1(I1,N)=DD1
2209             D2(I1,N)=DD2
2210             D1(I2,N)=DD1*SI
2211             D2(I2,N)=-DD2*SI
2212          ENDDO
2213       ENDDO
2214
2215  30   DO I=1,NGSS
2216          WR=W(I)*R(I)
2217          DS(I)=S(I)*QM*WR
2218          DSS(I)=SS(I)*QMM
2219          RR(I)=WR
2220       ENDDO
2221  
2222       DO N1=MM1,NMAX
2223          AN1=AN(N1)
2224          DO N2=MM1,NMAX
2225             AN2=AN(N2)
2226             AR11=0D0
2227             AR12=0D0
2228             AR21=0D0
2229             AR22=0D0
2230             AI11=0D0
2231             AI12=0D0
2232             AI21=0D0
2233             AI22=0D0
2234             GR11=0D0
2235             GR12=0D0
2236             GR21=0D0
2237             GR22=0D0
2238             GI11=0D0
2239             GI12=0D0
2240             GI21=0D0
2241             GI22=0D0
2242             SI=SIG(N1+N2)
2243             DO  I=1,NGSS
2244                D1N1=D1(I,N1)
2245                D2N1=D2(I,N1)
2246                D1N2=D1(I,N2)
2247                D2N2=D2(I,N2)
2248                A11=D1N1*D1N2
2249                A12=D1N1*D2N2
2250                A21=D2N1*D1N2
2251                A22=D2N1*D2N2
2252                AA1=A12+A21
2253                AA2=A11*DSS(I)+A22
2254                QJ1=XJ(I,N1)
2255                QY1=XY(I,N1)
2256                QJR2=XJR(I,N2)
2257                QJI2=XJI(I,N2)
2258                QDJR2=XDJR(I,N2)
2259                QDJI2=XDJI(I,N2)
2260                QDJ1=XDJ(I,N1)
2261                QDY1=XDY(I,N1)
2262
2263                C1R=QJR2*QJ1
2264                C1I=QJI2*QJ1
2265                B1R=C1R-QJI2*QY1
2266                B1I=C1I+QJR2*QY1
2267                
2268                C2R=QJR2*QDJ1
2269                C2I=QJI2*QDJ1
2270                B2R=C2R-QJI2*QDY1
2271                B2I=C2I+QJR2*QDY1
2272                
2273                DDRI=DDR(I)
2274                C3R=DDRI*C1R
2275                C3I=DDRI*C1I
2276                B3R=DDRI*B1R
2277                B3I=DDRI*B1I
2278  
2279                C4R=QDJR2*QJ1
2280                C4I=QDJI2*QJ1
2281                B4R=C4R-QDJI2*QY1
2282                B4I=C4I+QDJR2*QY1
2283                
2284                DRRI=DRR(I)
2285                DRII=DRI(I)
2286                C5R=C1R*DRRI-C1I*DRII
2287                C5I=C1I*DRRI+C1R*DRII
2288                B5R=B1R*DRRI-B1I*DRII
2289                B5I=B1I*DRRI+B1R*DRII
2290                
2291                C6R=QDJR2*QDJ1
2292                C6I=QDJI2*QDJ1
2293                B6R=C6R-QDJI2*QDY1
2294                B6I=C6I+QDJR2*QDY1
2295                
2296                C7R=C4R*DDRI
2297                C7I=C4I*DDRI
2298                B7R=B4R*DDRI
2299                B7I=B4I*DDRI
2300                
2301                C8R=C2R*DRRI-C2I*DRII
2302                C8I=C2I*DRRI+C2R*DRII
2303                B8R=B2R*DRRI-B2I*DRII
2304                B8I=B2I*DRRI+B2R*DRII
2305                
2306                URI=DR(I)
2307                DSI=DS(I)
2308 !c               DSSI=DSS(I)
2309                RRI=RR(I)
2310                IF (SI.GT.0D0) THEN    
2311                   F1=RRI*AA2
2312                   F2=RRI*URI*AN1*A12
2313                   AR12=AR12+F1*B2R+F2*B3R
2314                   AI12=AI12+F1*B2I+F2*B3I
2315                   GR12=GR12+F1*C2R+F2*C3R
2316                   GI12=GI12+F1*C2I+F2*C3I
2317                   
2318                   F2=RRI*URI*AN2*A21
2319                   AR21=AR21+F1*B4R+F2*B5R
2320                   AI21=AI21+F1*B4I+F2*B5I
2321                   GR21=GR21+F1*C4R+F2*C5R
2322                   GI21=GI21+F1*C4I+F2*C5I
2323 !c                  IF (NCHECK.NE.1) THEN
2324 !c                     E2=DSI*URI*A11
2325 !c                     E3=E2*AN2
2326 !c                     E2=E2*AN1
2327 !c                     AR22=AR22+E1*B6R+E2*B7R+E3*B8R
2328 !c                     AI22=AI22+E1*B6I+E2*B7I+E3*B8I
2329 !c                     GR22=GR22+E1*C6R+E2*C7R+E3*C8R
2330 !c                     GI22=GI22+E1*C6I+E2*C7I+E3*C8I
2331 !c                  ENDIF
2332                ELSE
2333 !c                  print*,'aa1',aa1
2334                   E1=DSI*AA1
2335                   AR11=AR11+E1*B1R
2336                   AI11=AI11+E1*B1I
2337                   GR11=GR11+E1*C1R
2338                   GI11=GI11+E1*C1I
2339 !c                  IF (NCHECK.EQ.1) THEN
2340                      E2=DSI*URI*A11
2341                      E3=E2*AN2
2342                      E2=E2*AN1
2343                      AR22=AR22+E1*B6R+E2*B7R+E3*B8R
2344                      AI22=AI22+E1*B6I+E2*B7I+E3*B8I
2345                      GR22=GR22+E1*C6R+E2*C7R+E3*C8R
2346                      GI22=GI22+E1*C6I+E2*C7I+E3*C8I
2347 !c                  ENDIF
2348                ENDIF
2349             ENDDO
2350             
2351             AN12=ANN(N1,N2)*FACTOR
2352             R11(N1,N2)=AR11*AN12
2353             R12(N1,N2)=AR12*AN12
2354             R21(N1,N2)=AR21*AN12
2355             R22(N1,N2)=AR22*AN12
2356             I11(N1,N2)=AI11*AN12
2357             I12(N1,N2)=AI12*AN12
2358             I21(N1,N2)=AI21*AN12
2359             I22(N1,N2)=AI22*AN12
2360             RG11(N1,N2)=GR11*AN12
2361             RG12(N1,N2)=GR12*AN12
2362             RG21(N1,N2)=GR21*AN12
2363             RG22(N1,N2)=GR22*AN12
2364             IG11(N1,N2)=GI11*AN12
2365             IG12(N1,N2)=GI12*AN12
2366             IG21(N1,N2)=GI21*AN12
2367             IG22(N1,N2)=GI22*AN12
2368             
2369          ENDDO
2370       ENDDO
2371       
2372       TPIR=PIR
2373       TPII=PII
2374       TPPI=PPI
2375       NM=NMAX-MM1+1
2376
2377       DO N1=MM1,NMAX
2378          K1=N1-MM1+1
2379          KK1=K1+NM
2380          DO N2=MM1,NMAX
2381             K2=N2-MM1+1
2382             KK2=K2+NM
2383             
2384             TAR11=-R11(N1,N2)
2385             TAI11=-I11(N1,N2)
2386             TGR11=-RG11(N1,N2)
2387             TGI11=-IG11(N1,N2)
2388             
2389             TAR12= I12(N1,N2)
2390             TAI12=-R12(N1,N2)
2391             TGR12= IG12(N1,N2)
2392             TGI12=-RG12(N1,N2)
2393             
2394             TAR21=-I21(N1,N2)
2395             TAI21= R21(N1,N2)
2396             TGR21=-IG21(N1,N2)
2397             TGI21= RG21(N1,N2)
2398             
2399             TAR22=-R22(N1,N2)
2400             TAI22=-I22(N1,N2)
2401             TGR22=-RG22(N1,N2)
2402             TGI22=-IG22(N1,N2)
2403             
2404             TQR(K1,K2)=TPIR*TAR21-TPII*TAI21+TPPI*TAR12
2405             TQI(K1,K2)=TPIR*TAI21+TPII*TAR21+TPPI*TAI12
2406             TRGQR(K1,K2)=TPIR*TGR21-TPII*TGI21+TPPI*TGR12
2407             TRGQI(K1,K2)=TPIR*TGI21+TPII*TGR21+TPPI*TGI12
2408             
2409             TQR(K1,KK2)=TPIR*TAR11-TPII*TAI11+TPPI*TAR22
2410             TQI(K1,KK2)=TPIR*TAI11+TPII*TAR11+TPPI*TAI22
2411             TRGQR(K1,KK2)=TPIR*TGR11-TPII*TGI11+TPPI*TGR22
2412             TRGQI(K1,KK2)=TPIR*TGI11+TPII*TGR11+TPPI*TGI22
2413             
2414             TQR(KK1,K2)=TPIR*TAR22-TPII*TAI22+TPPI*TAR11
2415             TQI(KK1,K2)=TPIR*TAI22+TPII*TAR22+TPPI*TAI11
2416             TRGQR(KK1,K2)=TPIR*TGR22-TPII*TGI22+TPPI*TGR11
2417             TRGQI(KK1,K2)=TPIR*TGI22+TPII*TGR22+TPPI*TGI11
2418             
2419             TQR(KK1,KK2)=TPIR*TAR12-TPII*TAI12+TPPI*TAR21
2420             TQI(KK1,KK2)=TPIR*TAI12+TPII*TAR12+TPPI*TAI21
2421             TRGQR(KK1,KK2)=TPIR*TGR12-TPII*TGI12+TPPI*TGR21
2422             TRGQI(KK1,KK2)=TPIR*TGI12+TPII*TGR12+TPPI*TGI21
2423          ENDDO
2424       ENDDO
2425       
2426       NNMAX=2*NM
2427       
2428       DO N1=1,NNMAX
2429          DO N2=1,NNMAX
2430             XQR(N1,N2)=TQR(N1,N2)
2431             XQI(N1,N2)=TQI(N1,N2)
2432             XRGQR(N1,N2)=TRGQR(N1,N2)
2433             XRGQI(N1,N2)=TRGQI(N1,N2)
2434          ENDDO
2435       ENDDO
2436 DEALLOCATE(R12)
2437 DEALLOCATE(R21)
2438 DEALLOCATE(R22)
2439 DEALLOCATE(I11)
2440 DEALLOCATE(I12)
2441 DEALLOCATE(I21)
2442 DEALLOCATE(I22)
2443 DEALLOCATE(RG11)
2444 DEALLOCATE(RG12)
2445 DEALLOCATE(RG21)
2446 DEALLOCATE(RG22)
2447 DEALLOCATE(IG11)
2448 DEALLOCATE(IG12)
2449 DEALLOCATE(IG21)
2450 DEALLOCATE(IG22)
2451
2452       
2453       CALL TT(NM)
2454
2455  
2456       RETURN
2457       END
2458  
2459 !C*****************************************************************
2460  
2461       SUBROUTINE VIG(X,NMAX,M,DV1,DV2)
2462
2463       USE MODD_TMAT, ONLY:NPN1
2464 !!      Parameter (NPN1=200)
2465       IMPLICIT REAL*8 (A-H,O-Z)
2466       REAL*8 DV1(NPN1),DV2(NPN1)
2467  
2468       A=1D0
2469       QS=SQRT(1D0-X*X)
2470       QS1=1D0/QS
2471       DO N=1,NMAX
2472       DV1(N)=0D0
2473       DV2(N)=0D0
2474       ENDDO   
2475
2476       IF (M.NE.0) THEN
2477       QMM=FLOAT(M*M)
2478           DO I=1,M
2479           I2=I*2
2480           A=A*SQRT(FLOAT(I2-1)/FLOAT(I2))*QS
2481           ENDDO   
2482       D1=0D0
2483       D2=A
2484           DO N=M,NMAX
2485           QN=FLOAT(N)
2486           QN2=FLOAT(2*N+1)
2487           QN1=FLOAT(N+1)
2488           QNM=SQRT(QN*QN-QMM)
2489           QNM1=SQRT(QN1*QN1-QMM)
2490           D3=(QN2*X*D2-QNM*D1)/QNM1
2491           DER=QS1*(-QN1*QNM*D1+QN*QNM1*D3)/QN2
2492           DV1(N)=D2
2493           DV2(N)=DER
2494           D1=D2
2495           D2=D3
2496           ENDDO   
2497       ELSE
2498       D1=1D0
2499       D2=X  
2500           DO N=1,NMAX  
2501           QN=FLOAT(N)
2502           QN1=FLOAT(N+1)
2503           QN2=FLOAT(2*N+1)
2504           D3=(QN2*X*D2-QN*D1)/QN1 
2505           DER=QS1*(QN1*QN/QN2)*(-D1+D3)
2506           DV1(N)=D2
2507           DV2(N)=DER
2508           D1=D2
2509           D2=D3
2510           ENDDO
2511       ENDIF
2512
2513       RETURN
2514       END 
2515  
2516 !C**********************************************************************
2517 !C                                                                     *
2518 !C   CALCULATION OF THE MATRIX    T = - RG(Q) * (Q**(-1))              *
2519 !C                                                                     *
2520 !C   INPUT INFORTMATION IS IN COMMON /CTT/                             *
2521 !C   OUTPUT INFORMATION IS IN COMMON /CT/                              *
2522 !C                                                                     *
2523 !C**********************************************************************
2524  
2525       SUBROUTINE TT(NMAX)
2526       USE MODD_TMAT,ONLY : XTR1,XTI1,XQR,XQI,XRGQR,XRGQI,NPN1,NPN2
2527
2528
2529 !!      Parameter (NPN1=200, NPN2=2*NPN1)
2530       IMPLICIT REAL*8 (A-H,O-Z)
2531       REAL*8 F(NPN2,NPN2),&
2532 !!            QR(NPN2,NPN2),QI(NPN2,NPN2),&
2533 !!            RGQR(NPN2,NPN2),RGQI(NPN2,NPN2),&
2534             A(NPN2,NPN2),C(NPN2,NPN2),D(NPN2,NPN2),E(NPN2,NPN2)
2535 !!      REAL*8 TR1(NPN2,NPN2),TI1(NPN2,NPN2)
2536 !c      INTEGER IPVT(NPN2)
2537 !!      COMMON /CT/ TR1,TI1
2538 !!      COMMON /CTT/ QR,QI,RGQR,RGQI
2539
2540        NDIM=NPN2
2541        NNMAX=2*NMAX
2542  
2543 !C        Gaussian elimination
2544
2545           DO N1=1,NNMAX
2546               DO N2=1,NNMAX
2547               F(N1,N2)=XQI(N1,N2)
2548               ENDDO
2549           ENDDO
2550
2551 !c          IF (NCHECK.EQ.1) THEN
2552           CALL INV1(NMAX,F,A)
2553 !c          ELSE
2554 !c          CALL INVERT(NDIM,NNMAX,F,A,COND,IPVT,WORK,B) 
2555 !c          ENDIF
2556  
2557        CALL PROD(XQR,A,C,NDIM,NNMAX)
2558        CALL PROD(C,XQR,D,NDIM,NNMAX)
2559
2560           DO N1=1,NNMAX
2561               DO N2=1,NNMAX
2562               C(N1,N2)=D(N1,N2)+XQI(N1,N2)
2563               ENDDO
2564           ENDDO
2565
2566 !c          IF (NCHECK.EQ.1) THEN
2567           CALL INV1(NMAX,C,XQI)
2568 !c          ELSE
2569 !c          CALL INVERT(NDIM,NNMAX,C,QI,COND,IPVT,WORK,B) 
2570 !c          ENDIF
2571
2572        CALL PROD(A,XQR,D,NDIM,NNMAX)
2573        CALL PROD(D,XQI,XQR,NDIM,NNMAX)
2574  
2575        CALL PROD(XRGQR,XQR,A,NDIM,NNMAX)
2576        CALL PROD(XRGQI,XQI,C,NDIM,NNMAX)
2577        CALL PROD(XRGQR,XQI,D,NDIM,NNMAX)
2578        CALL PROD(XRGQI,XQR,E,NDIM,NNMAX)
2579
2580           DO N1=1,NNMAX
2581               DO N2=1,NNMAX
2582               XTR1(N1,N2)=-A(N1,N2)-C(N1,N2)
2583               XTI1(N1,N2)= D(N1,N2)-E(N1,N2)
2584               ENDDO
2585           ENDDO
2586            RETURN
2587            END
2588  
2589 !C********************************************************************
2590 !C     Calcul of produit de 2 matrices de dimension NDIM*N           *
2591 !C********************************************************************
2592       SUBROUTINE PROD(A,B,C,NDIM,N)
2593
2594       REAL*8 A(NDIM,N),B(NDIM,N),C(NDIM,N),cij
2595
2596       DO I=1,N
2597           DO J=1,N
2598           CIJ=0d0
2599               DO K=1,N
2600               CIJ=CIJ+A(I,K)*B(K,J)
2601               ENDDO
2602           C(I,J)=CIJ
2603           ENDDO
2604       ENDDO
2605
2606       RETURN
2607       END
2608  
2609 !C**********************************************************************
2610  
2611       SUBROUTINE INV1(NMAX,F,A)
2612
2613       USE MODD_TMAT,ONLY : NPN1,NPN2
2614       IMPLICIT REAL*8 (A-H,O-Z)
2615 !!      Parameter (NPN1=200, NPN2=2*NPN1)
2616       REAL*8  A(NPN2,NPN2),F(NPN2,NPN2),B(NPN1),&
2617              WORK(NPN1),Q1(NPN1,NPN1),Q2(NPN1,NPN1),&
2618              P1(NPN1,NPN1),P2(NPN1,NPN1)
2619       INTEGER IPVT(NPN1),IND1(NPN1),IND2(NPN1)
2620
2621       NDIM=NPN1
2622       NN1=(FLOAT(NMAX)-0.1D0)*0.5D0+1D0 
2623       NN2=NMAX-NN1
2624       DO I=1,NMAX
2625       IND1(I)=2*I-1
2626           IF(I.GT.NN1) IND1(I)=NMAX+2*(I-NN1)
2627       IND2(I)=2*I
2628           IF(I.GT.NN2) IND2(I)=NMAX+2*(I-NN2)-1
2629       ENDDO
2630
2631       NNMAX=2*NMAX
2632
2633       DO I=1,NMAX
2634       I1=IND1(I)
2635       I2=IND2(I)
2636           DO J=1,NMAX
2637           J1=IND1(J)
2638           J2=IND2(J)
2639           Q1(J,I)=F(J1,I1)
2640           Q2(J,I)=F(J2,I2)
2641           ENDDO
2642       ENDDO
2643
2644       CALL INVERT(NDIM,NMAX,Q1,P1,COND,IPVT,WORK,B)
2645       CALL INVERT(NDIM,NMAX,Q2,P2,COND,IPVT,WORK,B)
2646
2647       DO I=1,NNMAX
2648           DO J=1,NNMAX
2649           A(J,I)=0D0
2650           ENDDO
2651       ENDDO
2652
2653       DO I=1,NMAX
2654       I1=IND1(I)
2655       I2=IND2(I)
2656           DO J=1,NMAX
2657           J1=IND1(J)
2658           J2=IND2(J)
2659           A(J1,I1)=P1(J,I)
2660           A(J2,I2)=P2(J,I)
2661           ENDDO
2662       ENDDO
2663
2664       RETURN
2665       END
2666  
2667 !C*********************************************************************
2668  
2669       SUBROUTINE INVERT(NDIM,N,A,X,COND,IPVT,WORK,B)
2670       IMPLICIT REAL*8 (A-H,O-Z)
2671       REAL*8 A(NDIM,N),X(NDIM,N),WORK(N),B(N)
2672       INTEGER IPVT(N)
2673
2674       CALL DECOMP (NDIM,N,A,COND,IPVT,WORK)
2675
2676 !C     IF (COND+1D0.EQ.COND) WRITE(10,5) COND
2677
2678 !C     IF (COND+1D0.EQ.COND) STOP
2679
2680 !C    5 FORMAT('THE MATRIX IS SINGULAR FOR THE GIVEN NUMERICAL ACCURACY',&
2681 !C           'COND = ',D12.6)
2682
2683       DO I=1,N
2684           DO J=1,N
2685           B(J)=0D0
2686               IF (J.EQ.I) B(J)=1D0
2687           ENDDO
2688       CALL SOLVE (NDIM,N,A,B,IPVT)
2689            DO J=1,N
2690            X(J,I)=B(J)
2691            ENDDO
2692       ENDDO
2693    
2694       RETURN
2695       END
2696  
2697 !C********************************************************************
2698  
2699       SUBROUTINE DECOMP (NDIM,N,A,COND,IPVT,WORK)
2700       IMPLICIT REAL*8 (A-H,O-Z)
2701       REAL*8 A(NDIM,N),COND,WORK(N)
2702       INTEGER IPVT(N)
2703
2704       IPVT(N)=1
2705 !*****IF 1
2706       IF(N.NE.1) THEN
2707       NM1=N-1
2708       ANORM=0D0
2709           DO J=1,N
2710           T=0D0
2711                DO I=1,N
2712                T=T+ABS(A(I,J))
2713                ENDDO
2714                IF (T.GT.ANORM) ANORM=T
2715           ENDDO
2716
2717           DO K=1,NM1
2718           KP1=K+1
2719           M=K
2720               DO I=KP1,N
2721                    IF (ABS(A(I,K)).GT.ABS(A(M,K))) M=I
2722               ENDDO
2723           IPVT(K)=M
2724               IF (M.NE.K) IPVT(N)=-IPVT(N)
2725           T=A(M,K)
2726           A(M,K)=A(K,K)
2727           A(K,K)=T
2728               IF (T.NE.0d0) THEN
2729                   DO I=KP1,N
2730                   A(I,K)=-A(I,K)/T
2731                   ENDDO
2732                   DO J=KP1,N
2733                   T=A(M,J)
2734                   A(M,J)=A(K,J)
2735                   A(K,J)=T
2736                       IF (T.NE.0D0) THEN
2737                           DO I=KP1,N
2738                           A(I,J)=A(I,J)+A(I,K)*T
2739                           ENDDO
2740                       ENDIF
2741                  ENDDO
2742              ENDIF
2743          ENDDO
2744
2745       DO K=1,N
2746       T=0D0
2747           IF (K.NE.1) THEN
2748           KM1=K-1
2749               DO I=1,KM1
2750               T=T+A(I,K)*WORK(I)
2751               ENDDO
2752           ELSE
2753           EK=1D0
2754               IF (T.LT.0D0) EK=-1D0
2755               IF (A(K,K).EQ.0D0) THEN 
2756 !               COND=1D52
2757               RETURN
2758               ELSE
2759               WORK(K)=-(EK+T)/A(K,K)
2760               ENDIF
2761           ENDIF
2762       ENDDO
2763
2764       DO KB=1,NM1
2765       K=N-KB
2766       T=0D0
2767       KP1=K+1
2768           DO I=KP1,N
2769           T=T+A(I,K)*WORK(K)
2770           ENDDO
2771       WORK(K)=T
2772       M=IPVT(K)
2773           IF (M.EQ.K) THEN
2774           T=WORK(M)
2775           WORK(M)=WORK(K)
2776           WORK(K)=T
2777           ENDIF
2778       ENDDO
2779       YNORM=0D0
2780
2781       DO I=1,N
2782       YNORM=YNORM+ABS(WORK(I))
2783       ENDDO
2784   
2785       CALL SOLVE (NDIM,N,A,WORK,IPVT)
2786       ZNORM=0D0
2787
2788       DO I=1,N
2789       ZNORM=ZNORM+ABS(WORK(I))
2790       ENDDO
2791
2792       COND=ANORM*ZNORM/YNORM
2793       IF (COND.LT.1d0) COND=1D0
2794       RETURN
2795
2796       ELSE
2797 !****ELSE DU IF 1
2798       COND=1D0
2799           IF (A(1,1).EQ.0D0) THEN
2800           COND=1D52
2801           ENDIF
2802       ENDIF
2803 !****END DU IF 1   
2804
2805       RETURN
2806       END
2807  
2808 !C**********************************************************************
2809  
2810       SUBROUTINE SOLVE (NDIM,N,A,B,IPVT)
2811       IMPLICIT REAL*8 (A-H,O-Z)
2812       REAL*8 A(NDIM,N),B(N)
2813       INTEGER IPVT(N)
2814
2815       IF (N.NE.1) THEN
2816       NM1=N-1
2817           DO K=1,NM1
2818           KP1=K+1
2819           M=IPVT(K)
2820           T=B(M)
2821           B(M)=B(K)
2822           B(K)=T
2823               DO I=KP1,N
2824               B(I)=B(I)+A(I,K)*T
2825               ENDDO
2826           ENDDO
2827           DO KB=1,NM1
2828           KM1=N-KB
2829           K=KM1+1
2830           B(K)=B(K)/A(K,K)
2831           T=-B(K)
2832               DO I=1,KM1
2833               B(I)=B(I)+A(I,K)*T
2834               ENDDO
2835           ENDDO
2836       ENDIF
2837
2838       B(1)= B(1)/A(1,1)
2839
2840       RETURN
2841       END
2842  
2843 !C**********************************************************************
2844 !C    CALCULATION OF POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE         *
2845 !C    FORMULA. IF IND1 = 0 - ON INTERVAL (-1,1), IF IND1 = 1 - ON      *
2846 !C    INTERVAL  (0,1). IF  IND2 = 1 RESULTS ARE PRINTED.               *
2847 !C    N - NUMBER OF POINTS                                             *
2848 !C    Z - DIVISION POINTS                                              *
2849 !C    W - WEIGHTS                                                      *
2850 !C**********************************************************************
2851  
2852       SUBROUTINE GAUSS(N,IND1,IND2,Z,W)
2853       IMPLICIT REAL*8 (A-H,P-Z)
2854       REAL*8 Z(N),W(N)
2855       DATA A,B,C /1D0,2D0,3D0/
2856     
2857       IND=MOD(N,2)
2858       K=N/2+IND
2859       F=FLOAT(N)
2860
2861 !*****DO 1
2862       DO I=1,K
2863       M=N+1-I
2864           IF(I.EQ.1) X=A-B/((F+A)*F)
2865           IF(I.EQ.2) X=(Z(N)-A)*4D0+Z(N)
2866           IF(I.EQ.3) X=(Z(N-1)-Z(N))*1.6D0+Z(N-1)
2867           IF(I.GT.3) X=(Z(M+1)-Z(M+2))*C+Z(M+3)
2868           IF(I.EQ.K.AND.IND.EQ.1) X=0D0
2869       NITER=0
2870       CHECK=1D-16
2871       PB=1D17
2872
2873           DO WHILE(ABS(PB).GT.CHECK*ABS(X))
2874           PB=1D0
2875           NITER=NITER+1
2876               IF (NITER.GT.100) CHECK=CHECK*10D0
2877               PC=X
2878               DJ=A
2879               DO J=2,N
2880               DJ=DJ+A
2881               PA=PB
2882               PB=PC
2883               PC=X*PB+(X*PB-PA)*(DJ-A)/DJ
2884               ENDDO
2885           PA=A/((PB-X*PC)*F)
2886           PB=PA*PC*(A-X*X)
2887           X=X-PB
2888           ENDDO  
2889
2890       Z(M)=X
2891       W(M)=PA*PA*(A-X*X)
2892           IF(IND1.EQ.0) W(M)=B*W(M) 
2893           IF(.NOT.(I.EQ.K.AND.IND.EQ.1))THEN
2894           Z(I)=-Z(M)
2895           W(I)=W(M)
2896           ENDIF
2897
2898       ENDDO
2899 !****END DU DO 1
2900 !!      IF(IND2.EQ.1) THEN
2901 !!      
2902 !! 1100 FORMAT(' ***  POINTS AND WEIGHTS OF GAUSSIAN QUADRATURE FORMULA',&
2903 !!      ' OF ',I4,'-TH ORDER')
2904 !!          DO I=1,K
2905 !!c          ZZ=-Z(I)
2906 !!C          WRITE(10,1200) I,ZZ,I,W(I)
2907 !!          ENDDO
2908 !! 1200 FORMAT(' ',4X,'X(',I4,') = ',F17.14,5X,'W(',I4,') = ',F17.14)
2909 !!      ELSE
2910 !!C     PRINT 1300,N
2911 !! 1300 FORMAT(' GAUSSIAN QUADRATURE FORMULA OF ',I4,'-TH ORDER IS USED')
2912 !!      ENDIF
2913
2914       IF(IND1.NE.0) THEN
2915           DO  I=1,N
2916           Z(I)=(A+Z(I))/B
2917           ENDDO 
2918       ENDIF
2919
2920       RETURN
2921       END
2922
2923
2924
2925
2926
2927
2928
2929