!+ Subroutine to obtain information about size and shape 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.
!
!- ---------------------------------------------------------------------
      SUBROUTINE PARTICLE_SIZE(IERR, I_DISTRIBUTION, I_SHAPE_PARTICLE
     &   , N_SIZE, RADIUS, 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
     &   , RADIUS_0
     &   )
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/distribution.pcf'
      INCLUDE 'incl/distribution.ncf'
      INCLUDE 'incl/shape_particle.pcf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/def_sct_io.icf'
!
!     DUMMY VARIABLES.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
     &   , I_DISTRIBUTION
!             NUMBER OF DISTRIBUTION
     &   , I_SHAPE_PARTICLE
!             SHAPE OF PARTICLE
     &   , N_SIZE
!             NUMBER OF SIZES IN DISTRIB.
     &   , N_MODE
!             NUMBER OF MODES IN DISTRIBUTION
      REAL*8	!, INTENT(OUT)
     &     WEIGHT_MODE(NPD_MODE)
!             WEIGHTS FOR EACH MODE
      REAL*8	!, INTENT(OUT)
     &     RADIUS(0: NPD_SIZE)
!             RADII OF SIZE DISTRIBUTION
     &   , NUMBER(0: NPD_SIZE)
!             POPULATIONS IN DISTRIBUTION
     &   , ASPECT(0: NPD_SIZE)
!             ASPECT RATIO OF PARTICLE
     &   , D2_NUMBER(0: NPD_SIZE)
!             2ND DERIVATIVE FOR SPLINE
     &   , D2_ASPECT(0: NPD_SIZE)
!             2ND DERIVATIVE FOR SPLINE
     &   , N0_LN
!             SCALING FACTOR FOR LOG NORM.
     &   , LN_R0_LN(NPD_MODE)
!             REFERENCE RADIUS FOR LOG NORM.
     &   , LN_SIGMA_LN(NPD_MODE)
!             STAND. DEV. OF LOG NORM.
     &   , COEFF_MG(NPD_MODE)
!             PARAMETER FOR MODIFIED GAMMA
     &   , ALPHA_MG(NPD_MODE)
!             "
     &   , RM_MG(NPD_MODE)
!             "
     &   , BETA_MG(NPD_MODE)
!             "
     &   , RADIUS_0
!             TYPICAL RADIUS FOR INTEGRAL
!     PARAMETERS FOR HEYMSFIELD AND PLATT'S DISTRIBUTIONS.
      REAL*8	!, INTENT(OUT)
     &     D_HP_MIN
!             LEAST 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
!
!
!     LOCAL VARIABLES.
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
      EXTERNAL
     &     LOCK_CODE
!
      CHARACTER
     &     LINE*80
!             LINE OF INPUT FILE
      INTEGER
     &     IOS
!             I/O ERROR FLAG
     &   , I
!             LOOP VARIABLE
     &   , N_KNOT
!             NUMBER OF KNOT POINTS
     &   , I_MAX
!             INDEX OF MODE OF NUMBER
      LOGICAL
     &     L_DATA_REGION
!             TRUE IN REGIONS OF DATA
      REAL*8
     &     NUMBER_MAX
!             PEAK VALUE OF NUMBER
     &   , NUMBER_TOT
!             TOTAL NUMBER DENSITY
     &   , WORK(0: NPD_SIZE)
!             WORKING ARRAY
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     OPEN_FILE_IN, SPLINE_FIT
!
!     FUNCTIONS CALLED:
      REAL*8
     &     GAMMA_FNC
!             GAMMA-FUNCTION
      EXTERNAL
     &     GAMMA_FNC
!
!
      INCLUDE 'incl/distrib_name.dcf'
!
!
!
      WRITE (IU_USER, '(/A)') 'SPECIFY THE SIZE DISTRIBUTION FOR '
     &   //'AVERAGING.'
      WRITE(IU_USER, '(/A/)') 'TYPE ONE OF THE FOLLOWING:'
      DO I=1, NPD_DISTRIBUTION
         WRITE(IU_USER, '(6X, I5, A2, 1X, A)')
     &      I, '.)',NAME_DISTRIBUTION(I)
      ENDDO
      WRITE(IU_USER, '(/)')
1     READ(IU_STDIN, *, IOSTAT=IOS) I_DISTRIBUTION
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(/A)') '+++ UNRECOGNIZED RESPONSE: '
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RE-TYPE.'
            GOTO 1
         ENDIF
      ENDIF
      IF ( (I_DISTRIBUTION.LE.0)
     &   .OR.(I_DISTRIBUTION.GT.NPD_DISTRIBUTION) ) THEN
         WRITE(IU_ERR, '(A)') '+++ ILLEGAL RESPONSE: '
     &      //'THIS NUMBER IS OUTSIDE THE SUPPORTED RANGE.'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RE-TYPE.'
            GOTO 1
         ENDIF
      ENDIF
!
      IF (I_DISTRIBUTION.EQ.IP_EXTERNAL) THEN
!        INPUT READ FROM FILE.
         CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &      , 'ENTER NAME OF FILE OF DISTRIBUTION.')
         IF (IERR.NE.I_NORMAL) RETURN
!        INITIALIZE DISTRIBUTION AND READ IN IT, DETERMINING THE MODAL
!        RADIUS FOR LATER SCALING.
         L_DATA_REGION=.FALSE.
         I_SHAPE_PARTICLE=IP_SHAPE_NULL
         N_SIZE=0
         RADIUS(N_SIZE)=0.0D+00
         NUMBER(N_SIZE)=0.0D+00
         ASPECT(N_SIZE)=1.0D+00
!        FIND RADIUS WHERE NUMBER IS MAXIMAL.
         NUMBER_MAX=0.0D+00
         I_MAX=1
2        READ(IU_FILE_IN, '(A)', END=4) LINE
         IF (L_DATA_REGION) THEN
            IF (LINE(1:4).NE.'*END') THEN
               BACKSPACE(IU_FILE_IN)
               N_SIZE=N_SIZE+1
               IF (N_SIZE.GT.NPD_SIZE) THEN
                  WRITE(IU_ERR, '(/A)') '*** ERROR: THERE ARE '
     &               //'TOO MANY ENTRIES IN THE SIZE DISTRIBUTION.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
                  READ(IU_FILE_IN, *) RADIUS(N_SIZE), NUMBER(N_SIZE)
               ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
                  READ(IU_FILE_IN, *) RADIUS(N_SIZE), NUMBER(N_SIZE)
     &               , ASPECT(N_SIZE)
               ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_NULL) THEN
                  WRITE(IU_ERR, '(A)') 
     &               '*** ERROR: SHAPE OF PARTICLE IS NOT SET.'
                  IERR=I_ERR_IO
                  RETURN
               ENDIF
!              ELIMINATE VERY CLOSE POINTS WHICH WILL GIVE INDETERMINATE
!              DIVISION IN THE SPLINE FITTING PROGRAM.
               IF (ABS(RADIUS(N_SIZE)-RADIUS(N_SIZE-1)).LT.TOL_TEST)
     &            N_SIZE=N_SIZE-1
               IF (NUMBER(N_SIZE).GT.NUMBER_MAX) THEN
                  I_MAX=N_SIZE
                  NUMBER_MAX=NUMBER(N_SIZE)
               ENDIF
            ELSE
               L_DATA_REGION=.FALSE.
            ENDIF
         ELSE
            IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
               L_DATA_REGION=.TRUE.
            ELSE IF (LINE(1:7).EQ.'*SPHERE') THEN
               I_SHAPE_PARTICLE=IP_SHAPE_SPHERE
            ELSE IF (LINE(1:7).EQ.'*HEXCYL') THEN
               I_SHAPE_PARTICLE=IP_SHAPE_HEXCYL
            ENDIF
         ENDIF

         GOTO 2
!
4        CLOSE(IU_FILE_IN)
!
!        N_KNOT IS THE NUMBER OF SPLINE KNOT POINTS. N_SIZE IS THE
!        NUMBER OF POINTS IN THE DISTRIBUTION AS READ IN, EXCLUDING
!        THE ZERO POINT.
         N_KNOT=N_SIZE+1
         RADIUS_0=RADIUS(I_MAX)
!
!        ESTABLISH A SPLINE FIT TO THE DATA.
         IF (N_SIZE.GT.1) THEN
            CALL SPLINE_FIT(N_KNOT, RADIUS, NUMBER, D2_NUMBER, WORK)
            CALL SPLINE_FIT(N_KNOT, RADIUS, ASPECT, D2_ASPECT, WORK)
         ENDIF
!
      ELSE IF (I_DISTRIBUTION.EQ.IP_LOG_NORMAL) THEN
!
!        LOG_NORMAL DISTRIBUTION.
!
!        SPHERES ARE ASSUMED HERE.
         I_SHAPE_PARTICLE=IP_SHAPE_SPHERE
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER OVERALL NUMBER DENSITY'
5        READ(IU_STDIN, *, IOSTAT=IOS) N0_LN
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 5
            ENDIF
         ENDIF
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER NUMBER OF MODES'
10       READ(IU_STDIN, *, IOSTAT=IOS) N_MODE
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 10
            ENDIF
         ENDIF
!
         RADIUS_0=0.0D+00
         DO I=1, N_MODE
            WRITE(IU_USER, '(/A, I3)')
     &         'ENTER VALUES OF WEIGHT, R0 AND SIGMA FOR MODE ', I
11           READ(IU_STDIN, *, IOSTAT=IOS) WEIGHT_MODE(I), LN_R0_LN(I)
     &         , LN_SIGMA_LN(I)
            IF (IOS.NE.0) THEN
               WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
               IF (LOCK_CODE(.TRUE.)) THEN
                  INCLUDE 'batch_error.finc'
               ELSE
                  WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
                  GOTO 11
               ENDIF
            ENDIF
            RADIUS_0=RADIUS_0+WEIGHT_MODE(I)*LN_R0_LN(I)
!           SINCE LOGS OF R0 AND SIGMA ARE NEEDED WE FORM THESE.
            LN_R0_LN(I)=LOG(LN_R0_LN(I))
            LN_SIGMA_LN(I)=LOG(LN_SIGMA_LN(I))
         ENDDO
!
      ELSE IF (I_DISTRIBUTION.EQ.IP_MODIFIED_GAMMA) THEN
!
!        SPHERES ARE ASSUMED HERE.
         I_SHAPE_PARTICLE=IP_SHAPE_SPHERE
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER TOTAL NUMBER DENSITY.'
6        READ(IU_STDIN, *, IOSTAT=IOS)
     &      NUMBER_TOT 
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 6
            ENDIF
         ENDIF
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER NUMBER OF MODES.'
12       READ(IU_STDIN, *, IOSTAT=IOS) N_MODE
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 12
            ENDIF
         ENDIF
!
         RADIUS_0=0.0D+00
         DO I=1, N_MODE
            WRITE(IU_USER, '(/A, /A, I3)')
     &         'ENTER VALUES OF WEIGHTING, SHAPE PARAMETER, '
     &         //'CHARACTERISTIC RADIUS AND SECOND POWER'
     &         , 'FOR MODE ', I
13           READ(IU_STDIN, *, IOSTAT=IOS)
     &         WEIGHT_MODE(I), ALPHA_MG(I), RM_MG(I), BETA_MG(I)
            IF (IOS.NE.0) THEN
               WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
               IF (LOCK_CODE(.TRUE.)) THEN
                  INCLUDE 'batch_error.finc'
               ELSE
                  WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
                  GOTO 13
               ENDIF
            ENDIF
!           THE TYPICAL RADIUS IS SET TO THE MODE OF THE DISTRIBUTION.
            RADIUS_0=RADIUS_0+WEIGHT_MODE(I)*RM_MG(I)
!           PRECALCULATE THE COEFFICIENT FOR FUTURE USE.
            COEFF_MG(I)=NUMBER_TOT*BETA_MG(I)
     &         /(RM_MG(I)*GAMMA_FNC(IERR, (ALPHA_MG(I)/BETA_MG(I))))
         ENDDO
!
!
      ELSE IF (I_DISTRIBUTION.EQ.IP_HEYMSFIELD_PLATT) THEN
!
!        HEXAGONAL COLUMNS ARE ASSUMED HERE.
         I_SHAPE_PARTICLE=IP_SHAPE_HEXCYL
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER THE SMALLEST VALUE OF THE  MAMIMUM '
     &      //'DIMENSION OF THE CRYSTAL.'
7        READ(IU_STDIN, *, IOSTAT=IOS) D_HP_MIN
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 7
            ENDIF
         ENDIF
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER THE NUMBER DENSITY AT 100 MICRONS '
     &      //'AND THE 1ST SLOPE.'
8        READ(IU_STDIN, *, IOSTAT=IOS) A_HP_1, B_HP_1
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 8
            ENDIF
         ENDIF
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER THE NUMBER DENSITY AT 1000 MICRONS '
     &      //'AND THE 2ND SLOPE.'
9        READ(IU_STDIN, *, IOSTAT=IOS) A_HP_2, B_HP_2
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ ERRONEOUS INPUT:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE_ENTER.'
               GOTO 9
            ENDIF
         ENDIF
!        FORM THE ACTUAL INTERCEPT PARAMETERS FROM THE NORMALIZED
!        NUMBER DENSITIES AND DETERMINE THE TRANSITIONAL DIMENSION.
         A_HP_1=A_HP_1/(1.0D-04**B_HP_1)
         A_HP_2=A_HP_2/(1.0D-03**B_HP_2)
         IF (ABS(B_HP_1-B_HP_2).GT.TOL_TEST) THEN
            D_HP_0=EXP(LOG(A_HP_2/A_HP_1)/(B_HP_1-B_HP_2))
         ELSE
            D_HP_0=5.0D-04
         ENDIF
         IF ( (D_HP_0.LE.1.0D-04).OR.(D_HP_0.GE.1.0D-03) ) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR: THE TRANSITIONAL DIMENSION IS OUT OF RANGE.'
            IERR=I_ERR_FATAL
            RETURN
         ENDIF
!        FOR WANT OF A MORE SUITABLE VALUE OF THE TYPICAL RADIUS WE TAKE
!        100 MICRONS.
         RADIUS_0=1.0D-04

      ENDIF
!
!
      RETURN
      END
