      SUBROUTINE DMIESS(  RO,      RFR,     RFI,     THETD,     JX,
     2                    QEXT,    QSCAT,   CTBRQS,  ELTRMX,    PI,
     3                    TAU,     CSTHT,   SI2THT,  ACAP,      IT,
     4                    LL,      R,       RE2,     TMAG2,     WVNO  )
C                                                         
C THIS VERSION OF THE FAMOUS DMIESS CODE IS ADAPTED FOR THE VAX
C ITS ACCURACY HAS BEEN CHECKED BY COMPARISON TO CRAY RUNS.
C              -C.P. MCKAY
C
      IMPLICIT REAL  (A-H,O-Z)
C
      COMPLEX    FNAP,      FNBP,      ACAP(LL),     W,
     2              FNA,       FNB,       RF,           RRF,
     3              RRFX,      WM1,       FN1,          FN2,
     4              TC1,       TC2,       WFN(2),       Z(4),
     5              K1,        K2,        K3,           SIN,
     6              COS,       RC,        U(8),         DH1,
     7              DH2,       DH4,       P24H24,       P24H21,
     8              PSTORE,    HSTORE,    DUMMY,        DUMSQ
C
      DIMENSION         W(3,9000)
      DIMENSION     T(5),      TA(4),     TB(2),        TC(2),
     2              TD(2),     TE(2),     PI( 3,IT ),   TAU( 3,IT ),
     3              CSTHT(IT), THETD(IT), SI2THT(IT),   ELTRMX( 4,IT,2 )
C
      EQUIVALENCE   ( WFN(1),TA(1) ),     (  FNA,TB(1) ),
     2              (    FNB,TC(1) ),     ( FNAP,TD(1) ),
     3              (   FNBP,TE(1) )
C
C
C    THIS SUBROUTINE COMPUTES MIE SCATTERING BY A STRATIFIED SPHERE
C    I.E. A PARTICLE CONSISTING OF A SPHERICAL CORE SURROUNDED BY A
C    SPHERICAL SHELL.  THE BASIC CODE USED WAS THAT DESCRIBED IN THE
C    REPORT: " SUBROUTINES FOR COMPUTING THE PARAMETERS OF THE
C    ELECTROMAGNETIC RADIATION SCATTERED BY A SPHERE " J.V. DAVE,
C    I B M SCIENTIFIC CENTER, PALO ALTO , CALIFORNIA.
C    REPORT NO. 320 - 3236 .. MAY 1968 .
C
C    THE MODIFICATIONS FOR STRATIFIED SPHERES ARE DESCRIBED IN
C        TOON AND ACKERMAN, APPL. OPTICS, IN PRESS, 1981
C
C
C    THE PARAMETERS IN THE CALLING STATEMENT ARE DEFINED AS FOLLOWS :
C      RO IS THE OUTER (SHELL) RADIUS;
C      R  IS THE CORE RADIUS;
C      RFR,   RFI  ARE THE REAL AND IMAGINARY PARTS OF THE SHELL INDEX
C                  OF REFRACTION IN THE FORM (RFR - I* RFI);
C      RE2, TMAG2  ARE THE INDEX PARTS FOR THE CORE;
C          ( WE ASSUME SPACE HAS UNIT INDEX. )
C      THETD(J): ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT
C          AND THE SCATTERED RADIATION.  THETD(J) IS< OR= 90.0
C          IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH
C          SUPPLEMENTARY VALUE, SEE COMMENTS BELOW ON ELTRMX;
C      JX: TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE
C          REQUIRED.  JX SHOULD NOT EXCEED IT UNLESS THE DIMENSIONS
C          STATEMENTS ARE APPROPRIATEDLY MODIFIED;
C
C      THE DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN"LIGHT
C          SCATTERING BY SMALL PARTICLES,H.C.VAN DE HULST, JOHN WILEY '
C          SONS, INC., NEW YORK, 1957" .
C      QEXT: EFFIECIENCY FACTOR FOR EXTINCTION,VAN DE HULST,P.14 ' 127.
C      QSCAT: EFFIECINCY FACTOR FOR SCATTERING,V.D. HULST,P.14 ' 127.
C      CTBRQS: AVERAGE(COSINE THETA) * QSCAT,VAN DE HULST,P.128
C      ELTRMX(I,J,K): ELEMENTS OF THE TRANSFORMATION MATRIX F,V.D.HULST
C          ,P.34,45 ' 125. I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1..
C          I = 3: ELEMENT S SUB 21.. I = 4: ELEMENT D SUB 21..
C      ELTRMX(I,J,1) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR
C          THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT
C          OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) ..
C
C      IT: IS THE DIMENSION OF THETD, ELTRMX, CSTHT, PI, TAU, SI2THT,
C          IT MUST CORRESPOND EXACTLY TO THE SECOND DIMENSION OF ELTRMX.
C      LL: IS THE DIMENSION OF ACAP
C          IN THE ORIGINAL PROGRAM THE DIMENSION OF ACAP WAS 7000.
C          FOR CONSERVING SPACE THIS SHOULD BE NOT MUCH HIGHER THAN
C          THE VALUE, N=1.1*(NREAL**2 + NIMAG**2)**.5 * X + 1
C      WVNO: 2*PI / WAVELENGTH
C
C      THIS SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE OF
C      DOWNWARD RECURRENCE RELATIONSHIP.
C
C      TA(1): REAL PART OF WFN(1).  TA(2): IMAGINARY PART OF WFN(1).
C      TA(3): REAL PART OF WFN(2).  TA(4): IMAGINARY PART OF WFN(2).
C      TB(1): REAL PART OF FNA.     TB(2): IMAGINARY PART OF FNA.
C      TC(1): REAL PART OF FNB.     TC(2): IMAGINARY PART OF FNB.
C      TD(1): REAL PART OF FNAP.    TD(2): IMAGINARY PART OF FNAP.
C      TE(1): REAL PART OF FNBP.    TE(2): IMAGINARY PART OF FNBP.
C      FNAP, FNBP  ARE THE PRECEDING VALUES OF FNA, FNB RESPECTIVELY.
C
C   IF THE CORE IS SMALL SCATTERING IS COMPUTED FOR THE SHELL ONLY

c     print*,'debut dmiess ',second(0.)

      IFLAG = 1
      IF ( R/RO .LT. 1.0E-06 )   IFLAG = 2
C
      IF ( JX .LE. IT )   GO TO 20
         WRITE( 6,7 )
         WRITE( 6,6 )
         STOP 30
   20 RF =  CMPLX( RFR,  -RFI )
      RC =  CMPLX( RE2,-TMAG2 )
      X  =  RO * WVNO
      K1 =  RC * WVNO
      K2 =  RF * WVNO
      ZET=0.0
      K3 =  CMPLX( WVNO, ZET )
      Z(1) =  K2 * RO
      Z(2) =  K3 * RO
      Z(3) =  K1 * R
      Z(4) =  K2 * R
      X1   = REAL( Z(1) )
      Y1   = AIMAG( Z(1) )
      X4   = REAL( Z(4) )
      Y4   = AIMAG( Z(4) )
c      X1   = REAL( Z(1) )
c      Y1   = AIMAG( Z(1) )
c      X4   = REAL( Z(4) )
c      Y4   = AIMAG( Z(4) )
      RRF  =  1.0 / RF
      RX   =  1.0 / X
      RRFX =  RRF * RX
      T(1) =  ( X**2 ) * ( RFR**2 + RFI**2 )
      T(1) =  SQRT( T(1) )
      NMX1 =  1.10 * T(1)

C
      IF ( NMX1 .LE. LL-1 )   GO TO 21
         WRITE(6,8)
	 PRINT*,'LL  =',LL
         STOP 32
   21 NMX2 = T(1)
      IF ( NMX1 .GT.  150 )   GO TO 22
         NMX1 = 150
         NMX2 = 135
C
   22 ACAP( NMX1+1 )  =  ( 0.0,0.0 )
      IF ( IFLAG .EQ. 2 )   GO TO 26
         DO 29   N = 1,3
   29    W( N,NMX1+1 )  =  ( 0.0,0.0 )
   26 CONTINUE
      DO 23   N = 1,NMX1
         NN = NMX1 - N + 1
         ACAP(NN) = (NN+1) * RRFX - 1.0 / ( (NN+1) * RRFX + ACAP(NN+1) )
         IF ( IFLAG .EQ. 2 )   GO TO 23
            DO 31   M = 1,3
   31       W( M,NN ) = (NN+1) / Z(M+1)  -
     1                       1.0 / (  (NN+1) / Z(M+1)  +  W( M,NN+1 )  )
   23 CONTINUE
      DO 30   J = 1,JX
      IF ( THETD(J) .LT. 0.0 )  THETD(J) =  ABS( THETD(J) )
      IF ( THETD(J) .GT. 0.0 )  GO TO 24
      CSTHT(J)  = 1.0
      SI2THT(J) = 0.0
      GO TO 30
   24 IF ( THETD(J) .GE. 90.0 )  GO TO 25
      T(1)      =  ( 3.14159265359 * THETD(J) ) / 180.0
      CSTHT(J)  =  COS( T(1) )
      SI2THT(J) =  1.0 - CSTHT(J)**2
      GO TO 30
   25 IF ( THETD(J) .GT. 90.0 )  GO TO 28
      CSTHT(J)  =  0.0
      SI2THT(J) =  1.0
      GO TO 30
   28 WRITE( 6,5 )  THETD(J)
      WRITE( 6,6 )
      STOP 34
   30 CONTINUE
C
      DO 35  J = 1,JX
      PI(1,J)  =  0.0
      PI(2,J)  =  1.0
      TAU(1,J) =  0.0
      TAU(2,J) =  CSTHT(J)
   35 CONTINUE
C
C   INITIALIZATION OF HOMOGENEOUS SPHERE
      T(1)   =  COS(X)
      T(2)   =  SIN(X)
      WM1    =  CMPLX( T(1),-T(2) )
      WFN(1) =  CMPLX( T(2), T(1) )
      WFN(2) =  RX * WFN(1) - WM1
      IF ( IFLAG .EQ. 2 )   GO TO 560
      N = 1
C
C INITIALIZATION PROCEDURE FOR STRATIFIED SPHERE BEGINS HERE
C
      SINX1   =  SIN( X1 )
      SINX4   =  SIN( X4 )
      COSX1   =  COS( X1 )
      COSX4   =  COS( X4 )
      EY1     =  EXP( Y1 )
      E2Y1    =  EY1 * EY1
      EY4     =  EXP( Y4 )
      EY1MY4  =  EXP( Y1 - Y4 )
      EY1PY4  =  EY1 * EY4
      EY1MY4  =  EXP( Y1 - Y4 )
      AA  =  SINX4 * ( EY1PY4 + EY1MY4 )
      BB  =  COSX4 * ( EY1PY4 - EY1MY4 )
      CC  =  SINX1 * ( E2Y1 + 1.0 )
      DD  =  COSX1 * ( E2Y1 - 1.0 )
      DENOM   =  1.0  +  E2Y1 * ( 4.0 * SINX1 * SINX1 - 2.0 + E2Y1 )
      REALP   =  ( AA * CC  +  BB * DD ) / DENOM
      AMAGP   =  ( BB * CC  -  AA * DD ) / DENOM
      DUMMY   =  CMPLX( REALP, AMAGP )
      AA  =  SINX4 * SINX4 - 0.5
      BB  =  COSX4 * SINX4
      P24H24  =  0.5 + CMPLX( AA,BB ) * EY4 * EY4
      AA  =  SINX1 * SINX4  -  COSX1 * COSX4
      BB  =  SINX1 * COSX4  +  COSX1 * SINX4
      CC  =  SINX1 * SINX4  +  COSX1 * COSX4
      DD  = -SINX1 * COSX4  +  COSX1 * SINX4
      P24H21  =  0.5 * CMPLX( AA,BB ) * EY1 * EY4  +
     2           0.5 * CMPLX( CC,DD ) * EY1MY4
      DH4  =  Z(4) / ( 1.0 + ( 0.0,1.0 ) * Z(4) )  -  1.0 / Z(4)
      DH1  =  Z(1) / ( 1.0 + ( 0.0,1.0 ) * Z(1) )  -  1.0 / Z(1)
      DH2  =  Z(2) / ( 1.0 + ( 0.0,1.0 ) * Z(2) )  -  1.0 / Z(2)
      PSTORE  =  ( DH4 + N / Z(4) )  *  ( W(3,N) + N / Z(4) )
      P24H24  =  P24H24 / PSTORE
      HSTORE  =  ( DH1 + N / Z(1) )  *  ( W(3,N) + N / Z(4) )
      P24H21  =  P24H21 / HSTORE
      PSTORE  =  ( ACAP(N) + N / Z(1) )  /  ( W(3,N) + N / Z(4) )
      DUMMY   =  DUMMY * PSTORE
      DUMSQ   =  DUMMY * DUMMY
C
C NOTE:  THE DEFINITIONS OF U(I) IN THIS PROGRAM ARE NOT THE SAME AS
C
C          USUB1 = U(1)                       USUB2 = U(5)
C          USUB3 = U(7)                       USUB4 = DUMSQ
C          USUB5 = U(2)                       USUB6 = U(3)
C          USUB7 = U(6)                       USUB8 = U(4)
C          RATIO OF SPHERICAL BESSEL FTN TO SPHERICAL HENKAL FTN = U(8)
C
      U(1) =  K3 * ACAP(N)  -  K2 * W(1,N)
      U(2) =  K3 * ACAP(N)  -  K2 * DH2
      U(3) =  K2 * ACAP(N)  -  K3 * W(1,N)
      U(4) =  K2 * ACAP(N)  -  K3 * DH2
      U(5) =  K1 *  W(3,N)  -  K2 * W(2,N)
      U(6) =  K2 *  W(3,N)  -  K1 * W(2,N)
      U(7) =  ( 0.0,-1.0 )  *  ( DUMMY * P24H21 - P24H24 )
      U(8) =  TA(3) / WFN(2)
C
      FNA  =  U(8) * ( U(1)*U(5)*U(7)  +  K1*U(1)  -  DUMSQ*K3*U(5) ) /
     2               ( U(2)*U(5)*U(7)  +  K1*U(2)  -  DUMSQ*K3*U(5) )
      FNB  =  U(8) * ( U(3)*U(6)*U(7)  +  K2*U(3)  -  DUMSQ*K2*U(6) ) /
     2               ( U(4)*U(6)*U(7)  +  K2*U(4)  -  DUMSQ*K2*U(6) )
      GO TO 561
  560 TC1  =  ACAP(1) * RRF  +  RX
      TC2  =  ACAP(1) * RF   +  RX
      FNA  =  ( TC1 * TA(3)  -  TA(1) ) / ( TC1 * WFN(2)  -  WFN(1) )
      FNB  =  ( TC2 * TA(3)  -  TA(1) ) / ( TC2 * WFN(2)  -  WFN(1) )
  561 CONTINUE
      FNAP = FNA
      FNBP = FNB
      T(1) = 1.50
C
C    FROM HERE TO THE STATMENT NUMBER 90, ELTRMX(I,J,K) HAS
C    FOLLOWING MEANING:
C    ELTRMX(1,J,K): REAL PART OF THE FIRST COMPLEX AMPLITUDE.
C    ELTRMX(2,J,K): IMAGINARY PART OF THE FIRST COMPLEX AMPLITUDE.
C    ELTRMX(3,J,K): REAL PART OF THE SECOND COMPLEX AMPLITUDE.
C    ELTRMX(4,J,K): IMAGINARY PART OF THE SECOND COMPLEX AMPLITUDE.
C    K = 1 : FOR THETD(J) AND K = 2 : FOR 180.0 - THETD(J)
C    DEFINITION OF THE COMPLEX AMPLITUDE: VAN DE HULST,P.125.
      TB(1) = T(1) * TB(1)
      TB(2) = T(1) * TB(2)
      TC(1) = T(1) * TC(1)
      TC(2) = T(1) * TC(2)
      DO 60 J = 1,JX
          ELTRMX(1,J,1) = TB(1) * PI(2,J) + TC(1) * TAU(2,J)
          ELTRMX(2,J,1) = TB(2) * PI(2,J) + TC(2) * TAU(2,J)
          ELTRMX(3,J,1) = TC(1) * PI(2,J) + TB(1) * TAU(2,J)
          ELTRMX(4,J,1) = TC(2) * PI(2,J) + TB(2) * TAU(2,J)
          ELTRMX(1,J,2) = TB(1) * PI(2,J) - TC(1) * TAU(2,J)
          ELTRMX(2,J,2) = TB(2) * PI(2,J) - TC(2) * TAU(2,J)
          ELTRMX(3,J,2) = TC(1) * PI(2,J) - TB(1) * TAU(2,J)
          ELTRMX(4,J,2) = TC(2) * PI(2,J) - TB(2) * TAU(2,J)
   60 CONTINUE
C
      QEXT   = 2.0 * ( TB(1) + TC(1))
      QSCAT  = ( TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 ) / 0.75
      CTBRQS = 0.0
      N = 2
   65 T(1) = 2*N - 1
      T(2) =   N - 1
      T(3) = 2*N + 1
      DO 70  J = 1,JX
          PI(3,J)  = ( T(1) * PI(2,J) * CSTHT(J) - N * PI(1,J) ) / T(2)
          TAU(3,J) = CSTHT(J) * ( PI(3,J) - PI(1,J) )  -
     1                          T(1) * SI2THT(J) * PI(2,J)  +  TAU(1,J)
   70 CONTINUE
C
C   HERE SET UP HOMOGENEOUS SPHERE
      WM1    =  WFN(1)
      WFN(1) =  WFN(2)
      WFN(2) =  T(1) * RX * WFN(1)  -  WM1
      IF ( IFLAG .EQ. 2 )   GO TO 1000
C
C   HERE SET UP STRATIFIED SPHERE
C
      DH2  =  - N / Z(2)  +  1.0 / ( N / Z(2) - DH2 )
      DH4  =  - N / Z(4)  +  1.0 / ( N / Z(4) - DH4 )
      DH1  =  - N / Z(1)  +  1.0 / ( N / Z(1) - DH1 )
      PSTORE  =  ( DH4 + N / Z(4) )  *  ( W(3,N) + N / Z(4) )
      P24H24  =  P24H24 / PSTORE
      HSTORE  =  ( DH1 + N / Z(1) )  *  ( W(3,N) + N / Z(4) )
      P24H21  =  P24H21 / HSTORE
      PSTORE  =  ( ACAP(N) + N / Z(1) )  /  ( W(3,N) + N / Z(4) )
      DUMMY   =  DUMMY * PSTORE
      DUMSQ   =  DUMMY * DUMMY
C
      U(1) =  K3 * ACAP(N)  -  K2 * W(1,N)
      U(2) =  K3 * ACAP(N)  -  K2 * DH2
      U(3) =  K2 * ACAP(N)  -  K3 * W(1,N)
      U(4) =  K2 * ACAP(N)  -  K3 * DH2
      U(5) =  K1 *  W(3,N)  -  K2 * W(2,N)
      U(6) =  K2 *  W(3,N)  -  K1 * W(2,N)
      U(7) =  ( 0.0,-1.0 )  *  ( DUMMY * P24H21 - P24H24 )
      U(8) =  TA(3) / WFN(2)
C
      FNA  =  U(8) * ( U(1)*U(5)*U(7)  +  K1*U(1)  -  DUMSQ*K3*U(5) ) /
     2               ( U(2)*U(5)*U(7)  +  K1*U(2)  -  DUMSQ*K3*U(5) )
      FNB  =  U(8) * ( U(3)*U(6)*U(7)  +  K2*U(3)  -  DUMSQ*K2*U(6) ) /
     2               ( U(4)*U(6)*U(7)  +  K2*U(4)  -  DUMSQ*K2*U(6) )
C
 1000 CONTINUE
      TC1  =  ACAP(N) * RRF  +  N * RX
      TC2  =  ACAP(N) * RF   +  N * RX
      FN1  =  ( TC1 * TA(3)  -  TA(1) ) / ( TC1 * WFN(2) - WFN(1) )
      FN2  =  ( TC2 * TA(3)  -  TA(1) ) / ( TC2 * WFN(2) - WFN(1) )
      M    =  WVNO * R
      IF ( N .LT. M )   GO TO 1002
      IF ( IFLAG .EQ. 2 )   GO TO 1001
C!!!!!!!!!!!! WARNING MODIF PERSO
C     IF (  ABS(  ( FN1-FNA ) / FN1  )  .LT.  1.0E-09   .AND.
C    1      ABS(  ( FN2-FNB ) / FN2  )  .LT . 1.0E-09  )       IFLAG = 2

      IF (  ABS(  ( FN1-FNA ) / FN1  )  .LT.  1.0E-4   .AND.
     1      ABS(  ( FN2-FNB ) / FN2  )  .LT . 1.0E-4  )       IFLAG = 2
      IF ( IFLAG .EQ. 1 )   GO TO 1002
 1001 FNA  =  FN1
      FNB  =  FN2
 1002 CONTINUE
      T(5)  =  N
      T(4)  =  T(1) / ( T(5) * T(2) )
      T(2)  =  (  T(2) * ( T(5) + 1.0 )  ) / T(5)
C
      CTBRQS  =  CTBRQS  +  T(2) * ( TD(1) * TB(1)  +  TD(2) * TB(2)  +
     1                               TE(1) * TC(1)  +  TE(2) * TC(2) )
     2                   +  T(4) * ( TD(1) * TE(1)  +  TD(2) * TE(2) )
      QEXT    =   QEXT  +  T(3) * ( TB(1) + TC(1) )
      T(4)    =  TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2
      QSCAT   =  QSCAT  +  T(3) * T(4)
      T(2)    =  N * (N+1)
      T(1)    =  T(3) / T(2)
      K = (N/2)*2
      DO 80 J = 1,JX
       ELTRMX(1,J,1) = ELTRMX(1,J,1)+T(1)*(TB(1)*PI(3,J)+TC(1)*TAU(3,J))
       ELTRMX(2,J,1) = ELTRMX(2,J,1)+T(1)*(TB(2)*PI(3,J)+TC(2)*TAU(3,J))
       ELTRMX(3,J,1) = ELTRMX(3,J,1)+T(1)*(TC(1)*PI(3,J)+TB(1)*TAU(3,J))
       ELTRMX(4,J,1) = ELTRMX(4,J,1)+T(1)*(TC(2)*PI(3,J)+TB(2)*TAU(3,J))
      IF ( K .EQ. N )   GO TO 75
       ELTRMX(1,J,2) = ELTRMX(1,J,2)+T(1)*(TB(1)*PI(3,J)-TC(1)*TAU(3,J))
       ELTRMX(2,J,2) = ELTRMX(2,J,2)+T(1)*(TB(2)*PI(3,J)-TC(2)*TAU(3,J))
       ELTRMX(3,J,2) = ELTRMX(3,J,2)+T(1)*(TC(1)*PI(3,J)-TB(1)*TAU(3,J))
       ELTRMX(4,J,2) = ELTRMX(4,J,2)+T(1)*(TC(2)*PI(3,J)-TB(2)*TAU(3,J))
      GO TO 80
   75  ELTRMX(1,J,2) =ELTRMX(1,J,2)+T(1)*(-TB(1)*PI(3,J)+TC(1)*TAU(3,J))
       ELTRMX(2,J,2) =ELTRMX(2,J,2)+T(1)*(-TB(2)*PI(3,J)+TC(2)*TAU(3,J))
       ELTRMX(3,J,2) =ELTRMX(3,J,2)+T(1)*(-TC(1)*PI(3,J)+TB(1)*TAU(3,J))
       ELTRMX(4,J,2) =ELTRMX(4,J,2)+T(1)*(-TC(2)*PI(3,J)+TB(2)*TAU(3,J))
   80 CONTINUE
C
C!!!!!!!!!!!! WARNING MODIF PERSO
C      IF ( T(4) .LT. 1.0E-14 )   GO TO 100
      IF ( T(4) .LT. 1.0E-4 )   GO TO 100
      N = N + 1
      DO 90 J = 1,JX
         PI(1,J)   =   PI(2,J)
         PI(2,J)   =   PI(3,J)
         TAU(1,J)  =  TAU(2,J)
         TAU(2,J)  =  TAU(3,J)
   90 CONTINUE
      FNAP  =  FNA
      FNBP  =  FNB
c     print*,'NMX2 =',nmx2
      IF ( N .LE. NMX2 )   GO TO 65
         WRITE( 6,8 )
         STOP 36
  100 DO 120 J = 1,JX
      DO 120 K = 1,2
         DO  115  I= 1,4
         T(I)  =  ELTRMX(I,J,K)
  115    CONTINUE
         ELTRMX(2,J,K)  =      T(1)**2  +  T(2)**2
         ELTRMX(1,J,K)  =      T(3)**2  +  T(4)**2
         ELTRMX(3,J,K)  =  T(1) * T(3)  +  T(2) * T(4)
         ELTRMX(4,J,K)  =  T(2) * T(3)  -  T(4) * T(1)
  120 CONTINUE
      T(1)    =    2.0 * RX**2
      QEXT    =   QEXT * T(1)
      QSCAT   =  QSCAT * T(1)
      CTBRQS  =  2.0 * CTBRQS * T(1)
C
c     print*,'NMX1= ',nmx1,'  LL=',ll
c     print*,'fin dmiess ',second(0.)
      RETURN
C
    5 FORMAT( 10X,' THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN
     1 90.0 DEGREES. IT IS ', E15.4 )
    6 FORMAT( // 10X, 'PLEASE READ COMMENTS.' // )
    7 FORMAT( // 10X, 'THE VALUE OF THE ARGUMENT JX IS GREATER THAN IT')
    8 FORMAT( // 10X, 'THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST
     1 GET DETAILED OUTPUT AND MODIFY SUBROUTINE' // )
C
      END
