!+ Subroutine to calculate properties of a size distribution.
!
! Method:
!       The total number density, the effective radius and the volume 
!       fraction 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 SIZE_INTEGRAL(IERR, RADIUS_0
     &   , I_DISTRIBUTION, I_SHAPE_PARTICLE
     &   , 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
     &   , NUMBER_TOTAL, RADIUS_EFF, VOLUME_FRACTION)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/parm_integ.acf'
      INCLUDE 'incl/shape_particle.pcf'
      INCLUDE 'incl/prec_integral.tcf'
      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 BEING CONSIDERED
     &   , N_SIZE
!             NUMBER OF SIZES
     &   , N_MODE
!             NUMBER OF MODES
      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)
!             WEIGHTS 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)
!             PARAMETER 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
!     PARAMETERS FOR HEYMSFIELD AND PLATT'S DISTRIBUTIONS.
      REAL*8	!, INTENT(IN)
     &     D_HP_MIN
!             LEAST VALE OF 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
      REAL*8	!, INTENT(OUT)
     &     NUMBER_TOTAL
!             TOTAL NUMBER OF PARTICLES/UNIT VOL.
     &   , RADIUS_EFF
!             EFFECTIVE RADIUS
     &   , VOLUME_FRACTION
!             VOLUME FRACTION
!
!     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
     &   , SINE_PI_3
!             SINE(PI/3)
      REAL*8
     &     SECOND_MOMENT
!             SECOND MOMENT OF RADIUS
     &   , THIRD_MOMENT
!             THIRD MOMENT OF RADIUS
     &   , ONE_HALF_MOMENT
!             ONE AND A HALFTH MOMENT OF RADIUS
!
!     FUNCTIONS CALLED:
      REAL*8
     &     NUMBER_PARTICLE
!             FUNCTION FOR NUMBER DENSITY
     &   , ASPECT_PARTICLE
!             ASPECT RATIO OF PARTICLE
     &   , TRAPEZOID
!             INTEGRATION ROUTINE
      EXTERNAL
     &     NUMBER_PARTICLE, ASPECT_PARTICLE, TRAPEZOID
!
!
!
      SINE_PI_3=5.0D-01*SQRT(3.0D+00)
      NUMBER_TOTAL=0.0D+00
      SECOND_MOMENT=0.0D+00
      THIRD_MOMENT=0.0D+00
      IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
         ONE_HALF_MOMENT=0.0D+00
      ENDIF
!     USE A VARIABLE FOR THE NUMBER OF INTEGRALS TO FACILITATE
!     SUBSEQUENT MODIFICATION.
      IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
         N_INTEGRAL=3
      ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
         N_INTEGRAL=4
      ENDIF
!
!     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_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
                  Y(I, 2)=Y(I, 1)*RADIUS(I)
               ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_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
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
     &               *(SINE_PI_3+2.0D+00*ASPECT_POINT)
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)**3*ASPECT_POINT
                  Y(I, 4)=NUMBER_POINT*RADIUS(I)**3*EXP(1.5D+00
     &               *LOG(SINE_PI_3+2.0D+00*ASPECT_POINT))
               ENDIF
               Y(I, 3)=NUMBER_POINT
            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_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
                  Y(I, 2)=Y(I, 1)*RADIUS(I)
               ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_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
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
     &               *(SINE_PI_3+2.0D+00*ASPECT_POINT)
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)**3*ASPECT_POINT
                  Y(I, 4)=NUMBER_POINT*RADIUS(I)**3*EXP(1.5D+00
     &               *LOG(SINE_PI_3+2.0D+00*ASPECT_POINT))
               ENDIF
               Y(I, 3)=NUMBER_POINT
            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_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
                  Y(I, 2)=Y(I, 1)*RADIUS(I)
               ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_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
                  Y(I, 1)=NUMBER_POINT*RADIUS(I)*RADIUS(I)
     &               *(SINE_PI_3+2.0D+00*ASPECT_POINT)
                  Y(I, 2)=NUMBER_POINT*RADIUS(I)**3*ASPECT_POINT
                  Y(I, 4)=NUMBER_POINT*RADIUS(I)**3*EXP(1.5D+00
     &               *LOG(SINE_PI_3+2.0D+00*ASPECT_POINT))
               ENDIF
               Y(I, 3)=NUMBER_POINT
            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)/7.5D-01
               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_STDOUT, '(/A)')
     &               '+++ WARNING: INTEGRAL HAS NOT CONVERGED.'
                  WRITE(IU_STDOUT, '(A)') 'LAST VALUE WILL BE USED.'
               ENDIF
            ENDIF
!
         SECOND_MOMENT=SECOND_MOMENT
     &      +INTEGRAL_PANEL(1, I_PANEL, I_REFINEMENT)
         THIRD_MOMENT=THIRD_MOMENT
     &      +INTEGRAL_PANEL(2, I_PANEL, I_REFINEMENT)
         NUMBER_TOTAL=NUMBER_TOTAL
     &      +INTEGRAL_PANEL(3, I_PANEL, I_REFINEMENT)
         IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
            ONE_HALF_MOMENT=ONE_HALF_MOMENT
     &         +INTEGRAL_PANEL(4, I_PANEL, I_REFINEMENT)
         ENDIF
      ENDDO
!
      IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
         RADIUS_EFF=THIRD_MOMENT/SECOND_MOMENT
         VOLUME_FRACTION=4.0D+00*PI*THIRD_MOMENT/3.0D+00
      ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
         RADIUS_EFF=SQRT(3.0D+00/(2.0D+00*PI))
     &      *ONE_HALF_MOMENT/SECOND_MOMENT
         VOLUME_FRACTION=6.0D+00*SINE_PI_3*THIRD_MOMENT
      ENDIF
!
!
      RETURN
      END
