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