!+ Subroutine to integrate scattering parameters over a distribution.
!
! Method:
!	The extinction, the scattering and the asymmetry are 
!	initialized to 0. Initial estimates of these quantities are
!	made and the range of integration is extended until further
!	extension does not sensibly alter these estimates. This is
!	done by dividing the range of integration into panels and
!	adding new panels as required. The resolution within a panel 
!	is then increased until the integrals converge.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE SCATTER_INTEGRAL(IERR, RADIUS_0
     &   , I_DISTRIBUTION, I_SHAPE_PARTICLE, I_SCATTER_ALGORITHM
     &   , N_SIZE, RADIUS_DIST, NUMBER, ASPECT, D2_NUMBER, D2_ASPECT
     &   , N_MODE, WEIGHT_MODE
     &   , N0_LN, LN_R0_LN, LN_SIGMA_LN
     &   , COEFF_MG, ALPHA_MG, RM_MG, BETA_MG
     &   , D_HP_MIN, A_HP_1, B_HP_1, A_HP_2, B_HP_2, D_HP_0
     &   , WAVELENGTH, N_RELATIVE
     &   , L_PHASE, N_ANGLE, MU_ANGLE
     &   , EXTINCTION, SCATTERING, ASYMMETRY
     &   , S1_AVE, S2_AVE)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/parm_integ.acf'
      INCLUDE 'incl/prec_integral.tcf'
      INCLUDE 'incl/scatter_algorithm.pcf'
      INCLUDE 'incl/shape_particle.pcf'
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/math_cnst.ccf'
!
!     DUMMY ARGUMENTS
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     I_DISTRIBUTION
!             TYPE OF SIZE DISTRIBUTION
     &   , I_SHAPE_PARTICLE
!             SHAPE OF PARTICLE
     &   , I_SCATTER_ALGORITHM
!             SCATTERING ALGORITHM
     &   , N_SIZE
!             NUMBER OF SIZES
     &   , N_MODE
!             NUMBER OF MODES
     &   , N_ANGLE
!             NUMBER OF SCATTERING ANGLES
      LOGICAL
     &     L_PHASE
!             CALCULATE PHASE FUNCTION
      REAL*8	!, INTENT(IN)
     &     RADIUS_DIST(0: NPD_SIZE)
!             RADII OF SIZE DISTRIBUTION
     &   , NUMBER(0: NPD_SIZE)
!             NUMBER DENSITY OF PARTICLES
     &   , ASPECT(0: NPD_SIZE)
!             ASPECT RATIO OF PARTICLES
     &   , D2_NUMBER(0: NPD_SIZE)
!             SECOND DERIVATIVE OF NUMBER
     &   , D2_ASPECT(0: NPD_SIZE)
!             SECOND DERIVATIVE OF ASPECT
     &   , WEIGHT_MODE(NPD_MODE)
!             WEIGHTING FACTORS FOR EACH MODE
     &   , N0_LN
!             FACTOR FOR LOG-NORMAL DIST.
     &   , LN_R0_LN(NPD_MODE)
!             LOG OF LOG-NORM REF. RADIUS
     &   , LN_SIGMA_LN(NPD_MODE)
!             LOG OF LOG-NORM STD. DEV.
     &   , COEFF_MG(NPD_MODE)
!             COEFFICIENT FOR MODIFIED GAMMA
     &   , ALPHA_MG(NPD_MODE)
!             SHAPE PARAMETER FOR MODIFIED GAMMA
     &   , RM_MG(NPD_MODE)
!             CHARACTERISTIC RADIUS FOR MODIFIED GAMMA
     &   , BETA_MG(NPD_MODE)
!             PARAMETER FOR MODIFIED GAMMA
     &   , RADIUS_0
!             TYPICAL RADIUS
     &   , WAVELENGTH
!             WAVELENGTH OF LIGHT
     &   , MU_ANGLE(NPD_SCATT_ANGLE)
!             COSINES OF SCATTERING ANGLES
!     PARAMETERS FOR HEYMSFIELD AND PLATT'S DISTRIBUTIONS.
      REAL*8	!, INTENT(IN)
     &     D_HP_MIN
!             MINIMUM MAXIMUM DIMENSION
     &   , A_HP_1
!             1ST INTERCEPT
     &   , B_HP_1
!             1ST SLOPE
     &   , A_HP_2
!             2ND INTERCEPT
     &   , B_HP_2
!             2ND SLOPE
     &   , D_HP_0
!             TRANSITIONAL DIMENSION
      COMPLEX*16	!, INTENT(IN)
     &     N_RELATIVE
!             RELATIVE REFGRACTIVE INDEX
      REAL*8	!, INTENT(OUT)
     &     EXTINCTION
!             MEAN EXTINCTION
     &   , SCATTERING
!             MEAN SCATTERING
     &   , ASYMMETRY
!             MEAN ASYMMETRY
      COMPLEX*16	!, INTENT(OUT)
     &     S1_AVE(NPD_SCATT_ANGLE)
!             MEAN VALUE OF S1
     &   , S2_AVE(NPD_SCATT_ANGLE)
!             MEAN VALUE OF S1
!
!
!     LOCAL VARIABLES.
      INTEGER
     &     N_INTEGRAL
!             NUMBER OF INTEGRALS TO EVALUATE
     &   , N_PANEL
!             NUMBER OF PANELS
     &   , I_PANEL
!             INDEX OF PANEL
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , N_POINT
!             NUMBER OF POINTS IN INTERVAL
     &   , I_REFINEMENT
!             INDEX OF REFINEMENT
      LOGICAL
     &     L_ADD_UPPER
!             LOGICAL TO ADD AN UPPER PANEL
     &   , L_ADD_LOWER
!             LOGICAL TO ADD A LOWER PANEL
     &   , L_REFINE
!             LOGICAL TO PERFORM REFINEMENT
      REAL*8
     &     RADIUS_LOW
!             LOWER LIMIT OF INTEGRALS
     &   , RADIUS_HIGH
!             UPPER LIMIT OF INTEGRALS
     &   , RADIUS_PANEL_LOW
!             LOWER LIMIT OF PANEL
     &   , RADIUS_PANEL_HIGH
!             UPPER LIMIT OF PANEL
     &   , INTEGRAL_ESTIMATE(NPD_INTEGRAL)
!             ESTIMATES OF INTEGRALS
     &   , INTEGRAL_PANEL_TEMP(NPD_INTEGRAL)
!             ESTIMATES OF PANEL INTEGRALS
     &   , INTEGRAL_PANEL(NPD_INTEGRAL, NPD_PANEL, NPD_REFINEMENT)
!             INTEGRALS OVER PANELS
     &   , RADIUS(NPD_POINT)
!             RADII OF POINTS OF INTEGRATION
     &   , Y(NPD_POINT, NPD_INTEGRAL)
!             POINT VALUES OF INTEGRAND
     &   , RADIUS_COARSE(NPD_PANEL, NPD_PANEL_POINT)
!             STARTING VALUES FOR ROMBERG STEP
     &   , Y_COARSE(NPD_INTEGRAL, NPD_PANEL, NPD_PANEL_POINT)
!             STARTING VALUES FOR ROMBERG STEP
     &   , REFINEMENT(NPD_INTEGRAL)
!             VALUE OF REFINEMENT
     &   , NUMBER_POINT
!             POINT VALUE OF NUMBER DENSITY
     &   , ASPECT_POINT
!             POINT VALUE OF ASPECT RATIO
      REAL*8
     &     ASYMMETRY_POINT
!             POINT VALUE OF ASYMMETRY
     &   , Q_EXT
!             EXTINCTION EFFICIENCY
     &   , Q_SCATTER
!             SCATTERING EFFICIENCY
     &   , XS_EXT
!             EXTINCTION CROSS-SECTION
     &   , XS_SCATTER
!             SCATTERING CROSS-SECTION
     &   , SIZE
!             SIZE PARAMETER
      COMPLEX*16
     &     S1(NPD_SCATT_ANGLE)
!             ELEMENT OF SCATTERING MATRIX
     &   , S2(NPD_SCATT_ANGLE)
!             ELEMENT OF SCATTERING MATRIX
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     MIE_SCATTER, TAKANO_LIOU_SCATTER
!     FUNCTIONS CALLED:
      REAL*8
     &     NUMBER_PARTICLE
!             FUNCTION FOR NUMBER DENSITY
     &   , ASPECT_PARTICLE
!             ASPECT RATIO OF PARTICLES
     &   , TRAPEZOID
!             INTEGRATION ROUTINE
      EXTERNAL
     &     NUMBER_PARTICLE, ASPECT_PARTICLE, TRAPEZOID
!
!
!
      EXTINCTION=0.0D+00
      SCATTERING=0.0D+00
      ASYMMETRY=0.0D+00
!
!     USE A VARIABLE FOR THE NUMBER OF INTEGRALS TO FACILITATE
!     SUBSEQUENT MODIFICATION.
      N_INTEGRAL=3
!
!     PRELIMINARY INTEGRATION
      DO J=1, N_INTEGRAL
         INTEGRAL_ESTIMATE(J)=0.0D+00
      ENDDO
      RADIUS_LOW=RADIUS_0
      RADIUS_HIGH=RADIUS_0
      N_PANEL=0
      L_ADD_UPPER=.TRUE.
      L_ADD_LOWER=.TRUE.
!
!     INCLUDE EXTRA PANELS IN THE RANGE OF INTEGRATION AT EITHER LIMIT
!     UNTIL THE INCLUSION OF FURTHER PANELS MAKES LITTLE CHANGE IN THE
!     ESTIMATE OF THE INTEGRAL.
!
 1    IF (L_ADD_UPPER) THEN
            RADIUS_PANEL_LOW=RADIUS_HIGH
            RADIUS_PANEL_HIGH=P_PANEL_RATIO*RADIUS_HIGH
            DO I=1, NPD_PANEL_POINT
               RADIUS(I)=RADIUS_PANEL_LOW
     &            *EXP((DBLE(I-1)/DBLE(NPD_PANEL_POINT-1))
     &            *LOG(RADIUS_PANEL_HIGH/RADIUS_PANEL_LOW))
               NUMBER_POINT=NUMBER_PARTICLE(IERR, I_DISTRIBUTION
     &            , I_SHAPE_PARTICLE
     &            , N_SIZE, RADIUS_DIST, NUMBER, D2_NUMBER
     &            , N_MODE, WEIGHT_MODE
     &            , N0_LN, LN_R0_LN, LN_SIGMA_LN
     &            , COEFF_MG, ALPHA_MG, RM_MG, BETA_MG
     &            , D_HP_MIN, A_HP_1, B_HP_1, A_HP_2, B_HP_2, D_HP_0
     &            , RADIUS(I))
               IF (IERR.NE.I_NORMAL) THEN
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE) THEN
                  IF (NUMBER_POINT.GT.0.0D+00) THEN
!                    THE SCATTERING CODE IS NOT CALLED IF THERE ARE
!                    NO PARTICLES AT THIS SIZE.
                     SIZE=2.0D+00*PI*RADIUS(I)/WAVELENGTH
                     CALL MIE_SCATTER(IERR, SIZE, N_RELATIVE
     &                  , Q_SCATTER, Q_EXT, ASYMMETRY_POINT
     &                  , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
!                     print*,Q_SCATTER,Q_EXT,ASYMMETRY_POINT
                     if (IERR.EQ.1)
     &                  print *, i, radius(i),wavelength, number_point

                  ELSE
                     Q_SCATTER=0.0D+00
                     Q_EXT=0.0D+00
                     ASYMMETRY_POINT=0.0
                  ENDIF
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     RETURN
                  ENDIF
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_EXT
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_SCATTER
               ELSE IF (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL)
     &            THEN
                  ASPECT_POINT=ASPECT_PARTICLE(IERR, I_DISTRIBUTION
     &               , I_SHAPE_PARTICLE
     &               , N_SIZE, RADIUS_DIST, ASPECT, D2_ASPECT
     &               , RADIUS(I))
                  IF (IERR.NE.I_NORMAL) RETURN
                  CALL TAKANO_LIOU_SCATTER(WAVELENGTH, RADIUS(I)
     &               , ASPECT_POINT, N_RELATIVE
     &               , XS_SCATTER, XS_EXT, ASYMMETRY_POINT
     &               )
                  Y(I, 1)=NUMBER_POINT*XS_EXT
                  Y(I, 2)=NUMBER_POINT*XS_SCATTER
               ENDIF
               Y(I, 3)=ASYMMETRY_POINT*Y(I, 2)
            ENDDO
!           THE PANEL IS INCLUDED IF ANY ONE OF THE INTEGRALS IS OF
!           SUFFICIENT SIZE.
            L_ADD_UPPER=.FALSE.
            DO J=1, N_INTEGRAL
               INTEGRAL_PANEL_TEMP(J)=TRAPEZOID(NPD_PANEL_POINT
     &            , RADIUS, Y(1, J))
               L_ADD_UPPER=L_ADD_UPPER.OR.(ABS(INTEGRAL_PANEL_TEMP(J))
     &            .GT.TOL_PANEL*ABS(INTEGRAL_ESTIMATE(J)))
            ENDDO
            IF (L_ADD_UPPER) THEN
               N_PANEL=N_PANEL+1
               RADIUS_HIGH=RADIUS_HIGH*P_PANEL_RATIO
               DO I=1, NPD_PANEL_POINT
                  RADIUS_COARSE(N_PANEL, I)=RADIUS(I)
               ENDDO
               DO J=1, N_INTEGRAL
                  INTEGRAL_ESTIMATE(J)=INTEGRAL_ESTIMATE(J)
     &               +INTEGRAL_PANEL_TEMP(J)
                  INTEGRAL_PANEL(J, N_PANEL, 1)=INTEGRAL_PANEL_TEMP(J)
                  DO I=1, NPD_PANEL_POINT
                     Y_COARSE(J, N_PANEL, I)=Y(I, J)
                  ENDDO
               ENDDO
            ENDIF
         ENDIF
!
         IF (L_ADD_LOWER) THEN
            RADIUS_PANEL_HIGH=RADIUS_LOW
            RADIUS_PANEL_LOW=RADIUS_PANEL_HIGH/P_PANEL_RATIO
            DO I=1, NPD_PANEL_POINT
               RADIUS(I)=RADIUS_PANEL_LOW
     &            *EXP((DBLE(I-1)/DBLE(NPD_PANEL_POINT-1))
     &            *LOG(RADIUS_PANEL_HIGH/RADIUS_PANEL_LOW))
               NUMBER_POINT=NUMBER_PARTICLE(IERR, I_DISTRIBUTION
     &            , I_SHAPE_PARTICLE
     &            , N_SIZE, RADIUS_DIST, NUMBER, D2_NUMBER
     &            , N_MODE, WEIGHT_MODE
     &            , N0_LN, LN_R0_LN, LN_SIGMA_LN
     &            , COEFF_MG, ALPHA_MG, RM_MG, BETA_MG
     &            , D_HP_MIN, A_HP_1, B_HP_1, A_HP_2, B_HP_2, D_HP_0
     &            , RADIUS(I))
               IF (IERR.NE.I_NORMAL) THEN
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE) THEN
                  IF (NUMBER_POINT.GT.0.0D+00) THEN
!                    THE SCATTERING CODE IS NOT CALLED IF THERE ARE
!                    NO PARTICLES AT THIS SIZE.
                     SIZE=2.0D+00*PI*RADIUS(I)/WAVELENGTH
                     CALL MIE_SCATTER(IERR, SIZE, N_RELATIVE
     &                  , Q_SCATTER, Q_EXT, ASYMMETRY_POINT
     &                  , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
!                     print*,Q_SCATTER,Q_EXT,ASYMMETRY_POINT
                     if (IERR.EQ.1)
     &                  print *, i, radius(i),wavelength, number_point
                  ELSE
                     Q_SCATTER=0.0D+00
                     Q_EXT=0.0D+00
                     ASYMMETRY_POINT=0.0
                  ENDIF
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     RETURN
                  ENDIF
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_EXT
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_SCATTER
               ELSE IF (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL)
     &            THEN
                  ASPECT_POINT=ASPECT_PARTICLE(IERR, I_DISTRIBUTION
     &               , I_SHAPE_PARTICLE
     &               , N_SIZE, RADIUS_DIST, ASPECT, D2_ASPECT
     &               , RADIUS(I))
                  IF (IERR.NE.I_NORMAL) RETURN
                  CALL TAKANO_LIOU_SCATTER(WAVELENGTH, RADIUS(I)
     &               , ASPECT_POINT, N_RELATIVE
     &               , XS_SCATTER, XS_EXT, ASYMMETRY_POINT
     &               )
                  Y(I, 1)=NUMBER_POINT*XS_EXT
                  Y(I, 2)=NUMBER_POINT*XS_SCATTER
               ENDIF
               Y(I, 3)=ASYMMETRY_POINT*Y(I, 2)
            ENDDO
            DO J=1, N_INTEGRAL
               INTEGRAL_PANEL_TEMP(J)=TRAPEZOID(NPD_PANEL_POINT
     &            , RADIUS, Y(1, J))
            ENDDO
!           THE PANEL IS INCLUDED IF ANY ONE OF THE INTEGRALS IS OF
!           SUFFICIENT SIZE.
            L_ADD_LOWER=.FALSE.
            DO J=1, N_INTEGRAL
               INTEGRAL_PANEL_TEMP(J)=TRAPEZOID(NPD_PANEL_POINT
     &            , RADIUS, Y(1, J))
               L_ADD_LOWER=L_ADD_LOWER.OR.(ABS(INTEGRAL_PANEL_TEMP(J))
     &            .GT.TOL_PANEL*ABS(INTEGRAL_ESTIMATE(J)))
            ENDDO
            IF (L_ADD_LOWER) THEN
               N_PANEL=N_PANEL+1
               RADIUS_LOW=RADIUS_LOW/P_PANEL_RATIO
               DO I=1, NPD_PANEL_POINT
                  RADIUS_COARSE(N_PANEL, I)=RADIUS(I)
               ENDDO
               DO J=1, N_INTEGRAL
                  INTEGRAL_ESTIMATE(J)=INTEGRAL_ESTIMATE(J)
     &               +INTEGRAL_PANEL_TEMP(J)
                  INTEGRAL_PANEL(J, N_PANEL, 1)=INTEGRAL_PANEL_TEMP(J)
                  DO I=1, NPD_PANEL_POINT
                     Y_COARSE(J, N_PANEL, I)=Y(I, J)
                  ENDDO
               ENDDO
            ENDIF
         ENDIF
!
      IF (L_ADD_LOWER.OR.L_ADD_UPPER) GOTO 1
!
!
!     NOW USE ROMBERG EXTRAPOLATION TO DETERMINE THE EXACT VALUE
!     OF THE INTEGRAL.
      DO I_PANEL=1, N_PANEL
!
         I_REFINEMENT=1
         N_POINT=NPD_PANEL_POINT
         DO I=1, N_POINT
            RADIUS(I)=RADIUS_COARSE(I_PANEL, I)
            DO J=1, N_INTEGRAL
               Y(I, J)=Y_COARSE(J, I_PANEL, I)
            ENDDO
         ENDDO
!
2           I_REFINEMENT=I_REFINEMENT+1
!           ADD NEW POINTS BETWEEN OLD ONES.
            DO I=N_POINT, 1, -1
               RADIUS(2*I-1)=RADIUS(I)
               DO J=1, N_INTEGRAL
                  Y(2*I-1, J)=Y(I, J)
               ENDDO
            ENDDO
            N_POINT=2*N_POINT-1
            DO I=2, N_POINT-1, 2
               RADIUS(I)=SQRT(RADIUS(I+1)*RADIUS(I-1))
               NUMBER_POINT=NUMBER_PARTICLE(IERR, I_DISTRIBUTION
     &            , I_SHAPE_PARTICLE
     &            , N_SIZE, RADIUS_DIST, NUMBER, D2_NUMBER
     &            , N_MODE, WEIGHT_MODE
     &            , N0_LN, LN_R0_LN, LN_SIGMA_LN
     &            , COEFF_MG, ALPHA_MG, RM_MG, BETA_MG
     &            , D_HP_MIN, A_HP_1, B_HP_1, A_HP_2, B_HP_2, D_HP_0
     &            , RADIUS(I))
               IF (IERR.NE.I_NORMAL) THEN
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE) THEN
                  IF (NUMBER_POINT.GT.0.0D+00) THEN
!                    THE SCATTERING CODE IS NOT CALLED IF THERE ARE
!                    NO PARTICLES AT THIS SIZE.
                     SIZE=2.0D+00*PI*RADIUS(I)/WAVELENGTH
                     CALL MIE_SCATTER(IERR, SIZE, N_RELATIVE
     &                  , Q_SCATTER, Q_EXT, ASYMMETRY_POINT
     &                  , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
                
                     if (IERR.EQ.1)
     &                  print *, i, radius(i),wavelength, number_point
                  ELSE
                     Q_SCATTER=0.0D+00
                     Q_EXT=0.0D+00
                     ASYMMETRY_POINT=0.0
                  ENDIF
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     RETURN
                  ENDIF
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_EXT
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)*RADIUS(I)*Q_SCATTER
               ELSE IF (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL)
     &            THEN
                  ASPECT_POINT=ASPECT_PARTICLE(IERR, I_DISTRIBUTION
     &               , I_SHAPE_PARTICLE
     &               , N_SIZE, RADIUS_DIST, ASPECT, D2_ASPECT
     &               , RADIUS(I))
                  IF (IERR.NE.I_NORMAL) RETURN
                  CALL TAKANO_LIOU_SCATTER(WAVELENGTH, RADIUS(I)
     &               , ASPECT_POINT, N_RELATIVE
     &               , XS_SCATTER, XS_EXT, ASYMMETRY_POINT
     &               )
                  Y(I, 1)=NUMBER_POINT*XS_EXT
                  Y(I, 2)=NUMBER_POINT*XS_SCATTER
               ENDIF
               Y(I, 3)=ASYMMETRY_POINT*Y(I, 2)
               print*,i,radius(i),Y(i,1),y(i,2),y(i,3)
            ENDDO
            L_REFINE=.FALSE.
            DO J=1, N_INTEGRAL
!              IN RELATIVE TERMS THERE IS LITTLE TO BE GAINED BY USING
!              THE EARLIER VALUE OF THE INTEGRAL FOR THE SUM
!              OVER ALTERNATE POINTS.
               REFINEMENT(J)=TRAPEZOID(N_POINT, RADIUS, Y(1, J))
     &            -INTEGRAL_PANEL(J, I_PANEL, I_REFINEMENT-1)
               INTEGRAL_PANEL(J, I_PANEL, I_REFINEMENT)
     &            =INTEGRAL_PANEL(J, I_PANEL, I_REFINEMENT-1)
     &            +REFINEMENT(J)/0.75D+00
               L_REFINE=L_REFINE.OR.(ABS(REFINEMENT(J))
     &            .GT.TOL_REFINEMENT*ABS(INTEGRAL_ESTIMATE(J)))
            ENDDO
            IF (L_REFINE) THEN
               IF ( (2*N_POINT.LT.NPD_POINT).AND.
     &              (I_REFINEMENT.LT.NPD_REFINEMENT) ) THEN
                  GOTO 2
               ELSE
                  WRITE(IU_ERR, '(/A, 1X, 1PE10.3, 1X, A)')
     &               '+++ WARNING: SCATTERING INTEGRALS AT WAVELENGTH'
     &               ,WAVELENGTH, 'HAVE NOT CONVERGED.'
                  WRITE(IU_ERR, '(A)') 'LAST VALUES WILL BE USED.'
               ENDIF
            ENDIF
!
         EXTINCTION=EXTINCTION
     &      +INTEGRAL_PANEL(1, I_PANEL, I_REFINEMENT)
         SCATTERING=SCATTERING
     &      +INTEGRAL_PANEL(2, I_PANEL, I_REFINEMENT)
         ASYMMETRY=ASYMMETRY
     &      +INTEGRAL_PANEL(3, I_PANEL, I_REFINEMENT)
      ENDDO
!
      ASYMMETRY=ASYMMETRY/SCATTERING
!
      IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
         EXTINCTION=PI*EXTINCTION
         SCATTERING=PI*SCATTERING
      ENDIF
                                
      print*,EXTINCTION,SCATTERING,ASYMMETRY
!
     
      RETURN
      END
