!+ Function to calculate number denisty of particles.
!
! Method:
!	Straightforward.
!
! 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. 
!
!- ---------------------------------------------------------------------
      FUNCTION NUMBER_PARTICLE(IERR, I_DISTRIBUTION, I_SHAPE_PARTICLE
     &   , N_SIZE, RADIUS, 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_POINT
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/distribution.pcf'
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/math_cnst.ccf'
      INCLUDE 'incl/max_size.acf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     I_DISTRIBUTION
!             NUMBER OF DISTRIBUTION
     &   , I_SHAPE_PARTICLE
!             SHAPE OF PARTICLE
     &   , N_SIZE
!             NUMBER OF PARTICLE SIZES
     &   , N_MODE
!             NUMBER OF MODES
      REAL*8	!, INTENT(OUT)
     &     NUMBER_PARTICLE
!             DIFFERENTIAL POPULATION
      REAL*8	!, INTENT(IN)
     &     RADIUS(0: NPD_SIZE)
!             RADII IN DISTRIBUTION
     &   , NUMBER(0: NPD_SIZE)
!             POPULATION VALUES
     &   , D2_NUMBER(0: NPD_SIZE)
!             2ND DERIVATIVE OF NUMBER
     &   , WEIGHT_MODE(NPD_MODE)
!             WEIGHTS FOR EACH MODE
     &   , N0_LN
!             SCALING FOR LOG NORMAL
     &   , LN_R0_LN(NPD_MODE)
!             LOG OF LOG NORM REF RADIUS
     &   , LN_SIGMA_LN(NPD_MODE)
!             LOG OF LOG NORM STANDARD DEVIATION
     &   , 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_POINT
!             RADIUS FOR EVALUATION
!     PARAMETERS FOR HEYMSFIELD AND PLATT'S DISTRIBUTIONS.
      REAL*8	!, INTENT(IN)
     &     D_HP_MIN
!             MINIMUM RADIUS OF CRYSTAL
     &   , 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
!
!     LOCAL VARIABLES.
      INTEGER
     &     N_KNOT
!             NUMBER OF KNOT POINTS
     &   , I_LOW
!             BISECTION VARIABLE
     &   , I_HIGH
!             BISECTION VARIABLE
     &   , I_SPLIT
!             BISECTION VARIABLE
     &   , I
!             LOOP VARIABLE
      REAL*8
     &     LOG_RADIUS
!             LOG OF RADIUS
     &   , MAX_DIMENSION
!             MAXIMUM DIMENSION OF PARTICLE
      REAL*8
     &     LENGTH
!             LENGTH OF CRYSTAL
     &   , DISCRIMINANT
!             DISCRIMINANT
     &   , CONVERSION
!             CONVERSION FACTOR
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     SPLINE_EVALUATE
!
!
!
!     A CUTOFF IS PLACED AT A PRESCRIBED LIMIT.
      IF (RADIUS_POINT.GT.MAX_RADIUS_SCATTERER) THEN
         NUMBER_PARTICLE=0.0D+00
         RETURN
      ENDIF
!     EVALUATE THE DIFFERENTIAL POPULATION AS DIRECTED 
!     BY I_DISTRIBUTION.
      IF (I_DISTRIBUTION.EQ.IP_EXTERNAL) THEN
         IF (N_SIZE.GT.1) THEN
            N_KNOT=N_SIZE+1
            CALL SPLINE_EVALUATE(IERR, N_KNOT, RADIUS, NUMBER, D2_NUMBER
     &         , RADIUS_POINT, NUMBER_PARTICLE)
            IF (IERR.NE.I_NORMAL) THEN
               IF (IERR.EQ.I_ERR_RANGE) THEN
                  NUMBER_PARTICLE=0.0D+00
                  IERR=I_NORMAL
               ELSE
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
            ENDIF
!
!           THE SPLINE FIT COULD GIVE NEGATIVE RESULTS, AND DATA MAY
!           HAVE STRETCHES OVER WHICH THE POPULATION IS ZERO: THE FIT
!           MUST THEREFORE BE ADJUSTED TO GIVE ZERO IN SUCH CASES.
            IF (NUMBER_PARTICLE.LT.0.0D+00) NUMBER_PARTICLE=0.0D+00
!           CHOP TO FIND RADIUS INTERVAL FOR CHECK.
            I_LOW=0
            I_HIGH=N_KNOT
1           I_SPLIT=(I_LOW+I_HIGH)/2
               IF (RADIUS_POINT.LT.RADIUS(I_SPLIT)) THEN
                  I_HIGH=I_SPLIT
               ELSE
                  I_LOW=I_SPLIT
               ENDIF
               IF ((I_HIGH-I_LOW).GT.1) THEN
                  GOTO 1
               ENDIF
            IF ( (NUMBER(I_LOW).LE.0.0D+00)
     &         .AND.(NUMBER(I_HIGH).LE.0.0D+00) )
     &         THEN
               NUMBER_PARTICLE=0.0D+00
            ENDIF
!
         ELSE
            NUMBER_PARTICLE=NUMBER(N_SIZE)
         ENDIF
      ELSE IF (I_DISTRIBUTION.EQ.IP_LOG_NORMAL) THEN
         LOG_RADIUS=LOG(RADIUS_POINT)
         NUMBER_PARTICLE=0.0D+00
         DO I=1, N_MODE
            NUMBER_PARTICLE=NUMBER_PARTICLE+WEIGHT_MODE(I)
     &         *(1.0D+00/RADIUS_POINT)
     &         *(N0_LN/(LN_SIGMA_LN(I)*SQRT(2.0D+00*PI)))
     &         *EXP(-(LOG_RADIUS-LN_R0_LN(I))**2
     &         /(2.0D+00*LN_SIGMA_LN(I)**2))
         ENDDO
      ELSE IF (I_DISTRIBUTION.EQ.IP_MODIFIED_GAMMA) THEN
         NUMBER_PARTICLE=0.0D+00
         DO I=1, N_MODE
            LOG_RADIUS=LOG(RADIUS_POINT/RM_MG(I))
            NUMBER_PARTICLE=NUMBER_PARTICLE+WEIGHT_MODE(I)
     &         *COEFF_MG(I)*EXP((ALPHA_MG(I)-1.0D+00)*LOG_RADIUS
     &         -EXP(BETA_MG(I)*LOG_RADIUS))
         ENDDO
      ELSE IF (I_DISTRIBUTION.EQ.IP_HEYMSFIELD_PLATT) THEN
!        FIND THE MAXIMUM DIMENSION OF PARTICLE TO APPLY IN THE
!        DISTRIBUTION.
!        THE FORMULAE FOR THE ASPECT RATIO ARE TAKEN FROM AUER AND VEAL
!        J. ATM. SCI. [27] p. 919 (1970).
         DISCRIMINANT=1.004004D+00-9.36D+03*(2.0D+00*RADIUS_POINT
     &      +8.479D-06)
         IF (DISCRIMINANT.GT.0.0D+00) THEN
            LENGTH=2.1368D-04*(1.002D+00-SQRT(DISCRIMINANT))
            IF (LENGTH.GT.2.0D-04) THEN
               LENGTH=4.7185D+06*EXP(2.415D+00*LOG(RADIUS_POINT))
               CONVERSION=1.13952D+07*EXP(1.415D+00*LOG(RADIUS_POINT))
            ELSE
               CONVERSION=2.0D+00/(1.002D+00-4.68D+03*LENGTH)
            ENDIF
         ELSE
            LENGTH=4.7185D+06*EXP(2.415D+00*LOG(RADIUS_POINT))
            CONVERSION=1.13952D+07*EXP(1.415D+00*LOG(RADIUS_POINT))
         ENDIF
         MAX_DIMENSION=MAX(2.0D+00*RADIUS_POINT, LENGTH)
         IF (MAX_DIMENSION.LT.D_HP_MIN) THEN
            NUMBER_PARTICLE=0.0D+00
         ELSE IF (MAX_DIMENSION.LE.D_HP_0) THEN
            NUMBER_PARTICLE=A_HP_1*EXP(B_HP_1*LOG(MAX_DIMENSION))
         ELSE IF (RADIUS_POINT.LT.MAX_RADIUS_SCATTERER) THEN
            NUMBER_PARTICLE=A_HP_2*EXP(B_HP_2*LOG(MAX_DIMENSION))
         ELSE
            NUMBER_PARTICLE=0.0D+00
         ENDIF
!        SINCE THIS DISTRIBUTION IS GIVEN IN TERMS OF THE MAXIMUM 
!        DIMENSION OF THE PARTICLE, WE NOW HAVE THE NUMBER DENSITY IN 
!        A UNIT RANGE OF MAXIMUM DIMENSIONS, AND THIS MUST BE SCALED 
!        FOR A UNIT RANGE OF RADII.
         NUMBER_PARTICLE=NUMBER_PARTICLE*CONVERSION
      ENDIF
!
!
      RETURN
      END
