!+ Program to perform scattering calculations.
!
! Method:
!	A file of wavelengths for the calculation and a file of
!	refractive indices are read in. A distribution and a 
!	scattering algorithm are supplied. The program calculates
!	single scattering properties at each of the wavelengths given.
!
! 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.
!
!- ---------------------------------------------------------------------
      PROGRAM SCATTER
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/def_sct_io.icf'
      INCLUDE 'incl/scatter_algorithm.pcf'
      INCLUDE 'incl/scatter_algorithm.ncf'
      INCLUDE 'incl/shape_particle.pcf'
      INCLUDE 'incl/aerosol_component.pcf'
      INCLUDE 'incl/aerosol_name.ncf'
      INCLUDE 'incl/scatter_pp.pcf'
      INCLUDE 'incl/scatter_pp.ncf'
      INCLUDE 'incl/distribution.pcf'
      INCLUDE 'incl/distribution.ncf'
      INCLUDE 'incl/math_cnst.ccf'
      INCLUDE 'incl/file_type.pcf'
      INCLUDE 'incl/error.pcf'
!
      CHARACTER
     &     FILE_OUT*80
!             NAME OF OUTPUT FILE
     &   , FILE_PHASE*80
!             FILE CONTAINING PHASE FUNCTION
     &   , LINE*80
!             LINE OF INPUT
     &   , NAME_COMPONENT*20
!             NAME OF SCATTERING COMPONENT
     &   , CHAR_YN
!             CHARACTER RESPONSE VARIABLE
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
      EXTERNAL
     &     LOCK_CODE
!
      INTEGER
     &     IERR
!             ERROR FLAG
     &   , IOS
!             I/O ERROR FLAG
     &   , N_WAVELENGTH
!             NUMBER OF WAVELENGTHS
     &   , N_REFRACT
!             NUMBER OF REFRACTIVE INDEX PAIRS
     &   , N_REFRACT_WATER
!             NUMBER OF REF. INDEX PAIRS FOR WATER
     &   , N_HUMIDITIES
!             NUMBER OF HUMIDITIES
     &   , ISIZE
!
     &   , N_SIZE
!             NUMBER OF ELEMENTS OF SIZE DIST.
     &   , N_ANGLE
!             NUMBER OF SCATTERING ANGLES
     &   , I_DISTRIBUTION
!             DISTRIBUTION TYPE
     &   , I_SHAPE_PARTICLE
!             SHAPE OF PARTICLE
     &   , LENGTH_OUT
!             LENGTH OF NAME OF OUTPUT FILE
     &   , I_SCATTER_TYPE
!             INDEX OF SCATTERING TYPE
     &   , ICHEM_TYPE
!             CHEMICAL TYPE OF AEROSOL
     &   , I_SCATTER_ALGORITHM
!             INDEXING NUMBER OF ALGORITHM
     &   , I_COMPONENT
!             INDEX OF SCATTERING COMPONENT
     &   , I_OUTPUT_TYPE
!             TYPE OF OUTPUT FILE
     &   , I
!             LOOP INDEX
     &   , J
!             LOOP INDEX
     &   , K
!             LOOP INDEX
      LOGICAL
     &     L_DATA_REGION
!             FLAG FOR INPUT OF DATA
     &   , L_PHASE
!             CALCULATE PHASE FUNCTION
     &   , L_EXIST
!             FILE EXISTENCE FLAG
     &   , L_HUMID
!             COMPUTE FOR VARIED HUMIDITIES
     &   , L_INTERP_SPLINE
!             USE CUBIC SPLINES TO INTERPOLATE THE REFRACTIVE INDEX
     &   ,L_ELLIE
!             Compute phase angle for external distribution with N_SIZE ne 1
      REAL*8
     &     SIZE
!             SIZE PARAMETER
     &   , WAVELENGTH(NPD_WAVELENGTH_SCAT)
!             WAVELENGTHS FOR SCATTERING PARAMETERS
     &   , WAVELENGTH_REFRACT(NPD_REFRACT)
!             WAVELENGTHS OF REFRACTIVE INDICES
     &   , WAVELENGTH_REFRACT_WATER(NPD_REFRACT)
!             WAVELENGTHS OF REFRACTIVE INDICES OF WATER
     &   , RE_REFRACT(NPD_REFRACT)
!             REAL PART OF REFRACTIVE INDEX
     &   , RE_REFRACT_WATER(NPD_REFRACT)
!             REAL PART OF WATER'S REF. INDEX
     &   , IM_REFRACT(NPD_REFRACT)
!             IMAGINARY PART OF REFRACTIVE INDEX
     &   , IM_REFRACT_WATER(NPD_REFRACT)
!             IMAG. PART OF WATER'S REF. INDEX
     &   , D2_RE_REFRACT(NPD_REFRACT)
!             SECOND DERIV. OF RE(REFRACT)
     &   , D2_RE_REFRACT_WATER(NPD_REFRACT)
!             SECOND DERIV. OF RE(REFRACT) WATER
     &   , D2_IM_REFRACT(NPD_REFRACT)
!             SECOND DERIV. OF IM(REFRACT)
     &   , D2_IM_REFRACT_WATER(NPD_REFRACT)
!             SECOND DERIV. OF IM(REFRACT) WATER
     &   , HUMIDITY(NPD_HUMIDITIES)
!             HUMIDITIES
     &   , NUMBER_TOTAL
!             NUMBER DENSITY OF PARTICLES
     &   , RADIUS_EFF
!             EFFECTIVE RADIUS
     &   , RADIUS_TYPICAL
!             TYPICAL RADIUS
     &   , NUMBER_0
!             NUMBER AT PEAK RADIUS
     &   , NUMBER_EFF
!             NUMBER AT EFFECTIVE RADIUS
     &   , VOLUME_FRACTION
!             VOLUME FRACTION OF SPHERES
     &   , NUMBER_TOTAL_DRY
!             NUMBER DENSITY OF DRY PARTICLES
     &   , RADIUS_EFF_DRY
!             DRY EFFECTIVE RADIUS
     &   , VOLUME_FRACTION_DRY
!             VOLUME FRACTION OF DRY SPHERES
     &   , GROWTH_FACTOR
!             HUMIDITY GROWTH FACTOR.
     &   , DILUTION_FACTOR
!             HUMIDITY DILUTION FACTOR.
!
!     PARAMETERS FOR SIZE DISTRIBUTION READ FROM FILE
      REAL*8
     &     RADIUS(0: NPD_SIZE)
!             RADII OF SIZE DISTRIBUTION
     &   , NUMBER(0: NPD_SIZE)
!             POPULATION SIZE OF DISTRIBUTION
     &   , ASPECT(0: NPD_SIZE)
!             ASPECT RATIOS OF CRYSTALS
     &   , D2_NUMBER(0: NPD_SIZE)
!             SECOND DERIV. OF NUMBER
     &   , D2_ASPECT(0: NPD_SIZE)
!             SECOND DERIV. OF ASPECT
!
!     PARAMETERS FOR ANALYTIC SIZE DISTRIBUTIONS
      INTEGER
     &     N_MODE
!             NUMBER OF MODES
      REAL*8
     &     WEIGHT_MODE(NPD_MODE)
!             WEIGHTING FACTORS FOR MODES
      REAL*8
     &     N0_LN
!             NUMBER DENSITY OF LOG-NORMAL DISTRIBUTION
     &   , LN_R0_LN(NPD_MODE)
!             LOG OF LOG-NORMAL REFERENCE RADIUS
     &   , LN_SIGMA_LN(NPD_MODE)
!             LOG OF LOG-NORMAL STANDARD DEVIATION
     &   , LN_R0_LN_DRY
!             DRY VALUE OF LN_R0_LN
     &   , LN_SIGMA_LN_DRY
!             DRY VALUE OF LN_SIGMA_LN
     &   , COEFF_MG(NPD_MODE)
!             COEFFICIENT OF MODIFIED GAMMA DISTRIBUTION
     &   , ALPHA_MG(NPD_MODE)
!             PARAMETER OF MODIFIED GAMMA DISTRIBUTION
     &   , RM_MG(NPD_MODE)
!             MODE RADIUS OF MODIFIED GAMMA DISTRIBUTION
     &   , BETA_MG(NPD_MODE)
!             PARAMETER OF MODIFIED GAMMA DISTRIBUTION
     &   , RADIUS_0
!             TYPICAL RADIUS OF DISTRIBUTION
!     PARAMETERS FOR HEYMSFIELD AND PLATT'S DISTRIBUTIONS
      REAL*8	!, INTENT(IN)
     &     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
      REAL*8
     &     Q_SCATTER
!             SCATTERING EFFICIENCY
     &   , Q_EXT
!             EXTINCTION EFFICIENCY
     &   , ANGLE(NPD_SCATT_ANGLE)
!             SCATTERING ANGLES
     &   , MU_ANGLE(NPD_SCATT_ANGLE)
!             COSINES OF SCATTERING ANGLES
     &   , XS_SCATTER
!             SCATTERING CROSS-SECTION
     &   , XS_EXT
!             EXTINCTION CROSS-SECTION
      COMPLEX*16
     &     S1(NPD_SCATT_ANGLE)
!             MATRIX ELEMENT
     &   , S2(NPD_SCATT_ANGLE)
!             MATRIX ELEMENT
      REAL*8
     &     ABSORPTION
!             AVERAGED ABSORPTION
     &   , EXTINCTION
!             AVERAGED EXTINCTION
     &   , SCATTERING
!             AVERAGED SCATTERING
     &   , ASYMMETRY
!             AVERAGED ASYMMETRY PARAMETER
      REAL*8
     &     WORK(NPD_REFRACT)
!             WORKING ARRAY
      COMPLEX*16
     &     N_RELATIVE
!             RELATIVE INDEX OF REFRACTION
     &   , N_RELATIVE_WATER
!             REFRACTIVE INDEX OF WATER
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     OPEN_FILE_IN, SPLINE_FIT, PARTICLE_SIZE, SIZE_INTEGRAL
     &   , GROW_PARTICLES, MIE_SCATTER, TAKANO_LIOU_SCATTER
     &   , SCATTER_INTEGRAL, ADT_INTEGRAL 
!
!     FUNCTIONS CALLED:
      COMPLEX*16
     &     REFRACTIVE_INDEX
!             FUNCTION TO FIND REFRACTIVE INDEX
      REAL*8
     &     NUMBER_PARTICLE
!             FUNCTION FOR NUMBER DENSITY
      EXTERNAL
     &     REFRACTIVE_INDEX, NUMBER_PARTICLE
!
!
      INCLUDE 'incl/distrib_name.dcf'
      INCLUDE 'incl/scatter_name.dcf'
      INCLUDE 'incl/aerosol_name.dcf'
      INCLUDE 'incl/scatter_algorithm.dcf'
      DATA IERR/I_NORMAL/
!
!
!
      WRITE(IU_USER, '(/A)') 'This program generates a file of single '
     &   //'scattering coefficients'
      WRITE(IU_USER, '(A)') 'at a range of frequencies held in a file'
     &   //' prepared by the user.'
      WRITE(IU_USER, '(A)') 'a table of refractive indices must also '
     &   //'be provided.'
      WRITE(IU_USER, '(/)')
!
!     ASK WHETHER HUMIDITY EFFECTS ARE REQUIRED AND SET THE TYPE OF
!     OUTPUT FILE.
!
      WRITE(IU_USER, '(/A)') 'Is growth with humidity required?'
10    READ(IU_STDIN, '(A)') CHAR_YN
      IF ( (CHAR_YN.EQ.'Y').OR.(CHAR_YN.EQ.'y') ) THEN
         L_HUMID = .TRUE.
         I_OUTPUT_TYPE=IT_FILE_MIE_HUMID
         WRITE(IU_USER, '(A)') 'A list of the humidities must also '
     &      //'be supplied.'
      ELSE IF ( (CHAR_YN.EQ.'N').OR.(CHAR_YN.EQ.'n') ) THEN
         L_HUMID = .FALSE.
         I_OUTPUT_TYPE=IT_FILE_MIE_DRY
         N_HUMIDITIES = 1
      ELSE
         WRITE(IU_ERR, '(/A)') '+++ UNRECOGNISED RESPONSE:'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A/)') 'PLEASE RE-ENTER.'
            GOTO 10
         ENDIF
      ENDIF
!
      print*,N_HUMIDITIES

      WRITE(IU_USER, '(/)')
!
!     SPECIFICATION OF WAVELENGTHS.
      CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &   , 'SPECIFY NAME OF FILE OF WAVELENGTHS.')
      IF (IERR.NE.I_NORMAL) STOP
      L_DATA_REGION=.FALSE.
      N_WAVELENGTH=0
20    READ(IU_FILE_IN, '(A)', END=21) LINE
         IF (L_DATA_REGION) THEN
            IF (LINE(1:4).NE.'*END') THEN
               BACKSPACE(IU_FILE_IN)
               N_WAVELENGTH=N_WAVELENGTH+1
               IF (N_WAVELENGTH.GT.NPD_WAVELENGTH_SCAT) THEN
                  WRITE(IU_ERR, '(/A)') '*** ERROR: THERE ARE '
     &               //'TOO MANY WAVELENGTHS IN THE INPUT FILE.'
                  IERR=I_ERR_FATAL
                  STOP
               ENDIF
               READ(IU_FILE_IN, *) WAVELENGTH(N_WAVELENGTH)
            ELSE
               L_DATA_REGION=.FALSE.
            ENDIF
         ELSE
            IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
               L_DATA_REGION=.TRUE.
            ENDIF
         ENDIF
!
      GOTO 20
!
21    CLOSE(IU_FILE_IN)
!
!     SPECIFICATION OF REFRACTIVE INDICES.
      CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &   , 'SPECIFY NAME OF FILE OF REFRACTIVE INDICES.')
      IF (IERR.NE.I_NORMAL) STOP
      L_DATA_REGION=.FALSE.
      N_REFRACT=0
30    READ(IU_FILE_IN, '(A)', END=31) LINE
         IF (L_DATA_REGION) THEN
            IF (LINE(1:4).NE.'*END') THEN
               BACKSPACE(IU_FILE_IN)
               N_REFRACT=N_REFRACT+1
               IF (N_REFRACT.GT.NPD_REFRACT) THEN
                  WRITE(IU_ERR, '(/A)') '*** ERROR: THERE ARE '
     &               //'TOO MANY WAVELENGTHS IN THE INPUT FILE.'
                  IERR=I_ERR_FATAL
                  STOP
               ENDIF
               READ(IU_FILE_IN, *) WAVELENGTH_REFRACT(N_REFRACT)
     &            , RE_REFRACT(N_REFRACT), IM_REFRACT(N_REFRACT)
            ELSE
               L_DATA_REGION=.FALSE.
            ENDIF
         ELSE
            IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
               L_DATA_REGION=.TRUE.
            ENDIF
         ENDIF
!
      GOTO 30
!
31    CLOSE(IU_FILE_IN)
!
!     ASK WHETHER CUBIC SPLINES ARE TO BE USED FOR INTERPOLATING THE
!     REFRACTIVE INDEX OF THE PARTICLE TO THE REQUIRED WAVELENGTHS.
!
!     NOTE THAT THE REFRACTIVE INDEX OF WATER, WHEN USED FOR
!     AEROSOLS THAT GROW WITH HUMIDITY, WILL ALWAYS USE CUBIC SPLINES.
!
      WRITE(IU_USER, '(/A)')
     &   'Which way should the refractive index be interpolated?'
      WRITE(IU_USER, '(/A)')
     &   'Answer C for cubic splines, L for linear interpolation.'
40    READ(IU_STDIN, '(A)') CHAR_YN
      IF ( (CHAR_YN.EQ.'C').OR.(CHAR_YN.EQ.'c') ) THEN
         L_INTERP_SPLINE = .TRUE.
      ELSE IF ( (CHAR_YN.EQ.'L').OR.(CHAR_YN.EQ.'l') ) THEN
         L_INTERP_SPLINE = .FALSE.
      ELSE
         WRITE(IU_ERR, '(/A/)')
     &      '+++ UNRECOGNISED RESPONSE: PLEASE RE-ENTER.'
         GOTO 40
      ENDIF
!
!     SPECIFICATION OF REFRACTIVE INDICES FOR WATER.
!     THIS IS ONLY REQUIRED IF HUMIDITY EFFECTS ARE BEING TAKEN
!     INTO ACCOUNT.
!
      IF (L_HUMID) THEN
         CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &      , 'SPECIFY FILENAME OF WATER REFRACTIVE INDICES.')
         print*,IERR,I_NORMAL
         IF (IERR.NE.I_NORMAL) STOP
         L_DATA_REGION=.FALSE.
         N_REFRACT_WATER=0
50       READ(IU_FILE_IN, '(A)', END=51) LINE
            IF (L_DATA_REGION) THEN
               IF (LINE(1:4).NE.'*END') THEN
                  BACKSPACE(IU_FILE_IN)
                  N_REFRACT_WATER=N_REFRACT_WATER+1
!                  print*,N_REFRACT_WATER
                  IF (N_REFRACT_WATER.GT.NPD_REFRACT) THEN
                     WRITE(IU_ERR, '(/A)') '*** ERROR: FILE OF '
     &                  //'REFRACTIVE INDICES FOR WATER IS TOO LONG.'
                     IERR=I_ERR_FATAL
                     STOP
                  ENDIF
                  READ(IU_FILE_IN, *)
     &               WAVELENGTH_REFRACT_WATER(N_REFRACT_WATER)
     &               , RE_REFRACT_WATER(N_REFRACT_WATER)
     &               , IM_REFRACT_WATER(N_REFRACT_WATER)
               ELSE
                  L_DATA_REGION=.FALSE.
               ENDIF
            ELSE
               IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
                  L_DATA_REGION=.TRUE.
               ENDIF
            ENDIF
!
         GOTO 50
!
51       CLOSE(IU_FILE_IN)
!
!        SPECIFICATION OF HUMIDITIES.
!
         CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &      , 'SPECIFY NAME OF FILE OF RELATIVE HUMIDITIES .')
         IF (IERR.NE.I_NORMAL) STOP
c         print*,L_data_region
         L_DATA_REGION=.FALSE.
         N_HUMIDITIES=0
c         print*,N_HUMIDITIES
60       READ(IU_FILE_IN, '(A)', END=61) LINE
         print*,L_data_region
!         L_data_region=.TRUE.
         IF (L_DATA_REGION) THEN
            IF (LINE(1:4).NE.'*END') THEN
               BACKSPACE(IU_FILE_IN)
               N_HUMIDITIES=N_HUMIDITIES+1
c               PRINT*,'N_HUMIDITIES=',N_HUMIDITIES
               IF (N_HUMIDITIES.GT.NPD_HUMIDITIES) THEN
                  WRITE(IU_ERR, '(/A)') 
     &               '*** ERROR: THERE ARE TOO MANY HUMIDITIES: '
     &               //'INCREASE NPD_HUMIDITIES AND RECOMPILE.'
                  IERR=I_ERR_FATAL
                  STOP
               ENDIF
               READ(IU_FILE_IN, *) HUMIDITY(N_HUMIDITIES)
c               PRINT*,HUMIDITY(N_HUMIDITIES)
            ELSE
               L_DATA_REGION=.FALSE.
            ENDIF
         ELSE
            IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
               L_DATA_REGION=.TRUE.
            ENDIF
         ENDIF
!
         GOTO 60
!
61       CLOSE(IU_FILE_IN)
!
!        FIND THE SECOND DERIVATIVES FOR THE SPLINE FITS TO THE
!        REFRACTIVE INDEX OF WATER.
!
         CALL SPLINE_FIT(N_REFRACT_WATER, WAVELENGTH_REFRACT_WATER
     &      , RE_REFRACT_WATER, D2_RE_REFRACT_WATER, WORK)
         CALL SPLINE_FIT(N_REFRACT_WATER, WAVELENGTH_REFRACT_WATER
     &      , IM_REFRACT_WATER, D2_IM_REFRACT_WATER, WORK)
!
!        END OF BLOCK OF CODE FOR HUMIDITY EFFECTS ONLY.
!
      ENDIF
!
      IF (L_INTERP_SPLINE) THEN
!
!        FIND THE SECOND DERIVATIVES FOR THE SPLINE FITS TO THE 
!        REFRACTIVE INDEX.
         CALL SPLINE_FIT(N_REFRACT, WAVELENGTH_REFRACT
     &      , RE_REFRACT, D2_RE_REFRACT, WORK)
         CALL SPLINE_FIT(N_REFRACT, WAVELENGTH_REFRACT
     &      , IM_REFRACT, D2_IM_REFRACT, WORK)
!
      ELSE
!
!        IF CUBIC SPLINES WERE NOT REQUESTED, ZERO THE SECOND
!        DERIVATIVES TO OBTAIN LINEAR INTERPOLATION INSTEAD.
         DO J = 1, NPD_REFRACT
            D2_RE_REFRACT(J) = 0.0D+00
            D2_IM_REFRACT(J) = 0.0D+00
         ENDDO
      ENDIF
!
!     SET THE SIZE DISTRIBUTION.
!     IF HUMIDITY VARIATION IS BEING USED, THIS IS THE SIZE DISTRIBUTION
!     OF THE DRY AEROSOL, BEFORE HUMIDITY GROWTH TAKES PLACE.
!
      CALL 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
     &   )
      IF (IERR.NE.I_NORMAL) THEN
         IERR=I_ERR_FATAL
         STOP
      ENDIF
!
!     CHECK THAT THE TYPE OF DISTRIBUTION IS CONSISTENT WITH THE
!     ASSUMPTIONS MADE IN THE PROGRAM.
!
      IF ( L_HUMID .AND. I_DISTRIBUTION .NE. IP_LOG_NORMAL ) THEN
         WRITE(IU_USER,'(/A)') 'HUMIDITY EFFECTS ARE ONLY AVAILABLE '
     &      //'WITH LOGNORMAL SIZE DISTRIBUTIONS AT PRESENT.'
         STOP
      ENDIF
!
!     SELECT THE SCATTERING ALGORITHM.
      WRITE(IU_USER, '(/A)')
     &   'SELECT THE ALGORITHM FOR THE SCATTERING CALCULATION.'
      DO I=1, NPD_ALGORITHM_SCATTER
         WRITE(IU_USER, '(3X, I3, A3, A40)')
     &      I, '.) ', NAME_SCATTER_ALGORITHM(I)
      ENDDO
      WRITE(IU_USER, '(/)')
70    READ(IU_STDIN, *) I_SCATTER_ALGORITHM
!
!     ADT IS A PERMITTED OPTION ONLY WITH A GAMMA DISTRIBUTION.
      IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_ADT) THEN
         IF (I_DISTRIBUTION.NE.IP_MODIFIED_GAMMA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '+++ ADT IS PERMITTED ONLY WITH A MODIFIED '
     &         //'GAMMA-DISTRIBUTION.'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A/)') 'PLEASE RE-ENTER.'
               GOTO 70
            ENDIF
         ENDIF
      ENDIF
      
!     print*,I_SCATTER_ALGORITHM
!     print*,I_SHAPE_PARTICLE
!      print*,IP_SHAPE_SPHERE

      print*,N_SIZE
!
!     DETERMINE THE CONSISTENCY OF THE SCATTERING ALGORITHM WITH THE
!     SHAPE OF THE PARTICLE.
      IF ( ( (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE).AND.
     &       (I_SHAPE_PARTICLE.NE.IP_SHAPE_SPHERE) ).OR.
     &     ( (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_ADT).AND.
     &       (I_SHAPE_PARTICLE.NE.IP_SHAPE_SPHERE) ).OR.
     &     ( (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL).AND.
     &       (I_SHAPE_PARTICLE.NE.IP_SHAPE_HEXCYL) ) ) THEN
         WRITE(IU_ERR, '(/A)') 
     &      '*** ERROR: THE SCATTERING ALGORITHM SELECTED '
     &      //'IS NOT APPLICABLE TO THIS SHAPE.'
         STOP
      ENDIF
!
!     DETERMINE WHETHER THE PHASE FUNCTION IS TO BE CALCULATED. THIS IS
!     ONLY PERMITTED FOR AN EXTERNAL DISTRIBUTION WITH ONE SIZE AND
!     ONLY FOR CERTAIN ALGORITHMS.

!  ********* Ellie's extra algorithm to allow pahse function for N_size.ge.1


      WRITE(IU_USER, '(/A)') 'Ellies PHASE FUNCTION '
     &      //'CALCULATION ? (Y/N)'
      READ(IU_STDIN, '(A)') CHAR_YN
      IF ( (CHAR_YN.EQ.'Y').OR.(CHAR_YN.EQ.'y') ) THEN
         L_ELLIE=.TRUE.
      ELSE
         L_ELLIE=.FALSE.
      ENDIF
         
      print*,L_ELLIE
! *************************************************************************

      IF ( (I_DISTRIBUTION.EQ.IP_EXTERNAL).AND.(N_SIZE.EQ.1)
     &     .AND.(I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE).OR.
     &     (L_ELLIE) ) THEN
         WRITE(IU_USER, '(/A)') 'IS THE PHASE FUNCTION TO BE '
     &      //'CALCULATED ? (Y/N)'
         READ(IU_STDIN, '(A)') CHAR_YN
         IF ( (CHAR_YN.EQ.'Y').OR.(CHAR_YN.EQ.'y') ) THEN
            L_PHASE=.TRUE.
            CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &         , 'SPECIFY NAME OF FILE CONTAINING SCATTERING ANGLES.')
            IF(IERR.NE.I_NORMAL) STOP
               L_DATA_REGION=.FALSE.
               N_ANGLE=0
200            READ(IU_FILE_IN, '(A)', END=210) LINE
                  IF (L_DATA_REGION) THEN
                     IF (LINE(1:4).NE.'*END') THEN
                        BACKSPACE(IU_FILE_IN)
                        N_ANGLE=N_ANGLE+1
                        IF (N_ANGLE.GT.NPD_SCATT_ANGLE) THEN
                           WRITE(IU_ERR, '(/A)') '*** ERROR: THERE ARE '
     &                        //'TOO MANY ANGLES IN THE INPUT FILE.'
                           IERR=I_ERR_FATAL
                           STOP
                        ENDIF
                        READ(IU_FILE_IN, *) ANGLE(N_ANGLE)
                     ELSE
                        L_DATA_REGION=.FALSE.
                     ENDIF
                  ELSE
                     IF (LINE(1:11).EQ.'*BEGIN_DATA') THEN
                        L_DATA_REGION=.TRUE.
                     ENDIF
                  ENDIF
!
               GOTO 200
!
210         CLOSE(IU_FILE_IN)
            DO J=1, N_ANGLE
               MU_ANGLE(J)=COS(ANGLE(J))
            ENDDO
         ELSE
            L_PHASE=.FALSE.
         ENDIF
      ELSE
         L_PHASE=.FALSE.
      ENDIF
!
!
!     READ IN THE NAME OF THE FILE FOR OUTPUT.
      WRITE(IU_USER, '(/A/)')
     &   'GIVE THE NAME OF THE FILE TO CONTAIN THE SCATTERING DATA.'
      READ(IU_STDIN, '(A)', IOSTAT=IOS) FILE_OUT
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(/A)')
     &      '*** ERROR: FILE NAME COULD NOT BE READ SUCCESSFULLY.'
         IERR=I_ERR_IO
         STOP
      ENDIF
!     CHECK FOR EXISTENCE OF THE FILE.
80    INQUIRE(FILE=FILE_OUT, EXIST=L_EXIST)
      IF (L_EXIST) THEN
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(/A)') 'THIS FILE ALREADY EXISTS.'
            WRITE(IU_USER, '(/A/)') 'DO YOU WISH TO OVERWRITE? (Y/N).'
            READ(IU_STDIN, '(A)') CHAR_YN
            IF ( (CHAR_YN.NE.'Y').AND.(CHAR_YN.NE.'y') ) THEN
               WRITE(IU_USER, '(/A/)')
     &            'PLEASE SPECIFY ANOTHER NAME OR !QUIT.'
               READ(IU_STDIN, '(A)') FILE_OUT
               IF ( (FILE_OUT(1:5).EQ.'!QUIT').OR.
     &              (FILE_OUT(1:5).EQ.'!quit').OR.
     &              (FILE_OUT(1:5).EQ.'!Quit') ) THEN
                  WRITE(IU_USER, '(/A/)') '*** PROGRAM TERMINATED.'
                  IERR=I_ERR_IO
                  STOP
               ELSE
                  GOTO 80
               ENDIF
            ENDIF
         ENDIF
      ENDIF
!
!     THE FILE MAY NOW SAFELY BE OPENED.
      OPEN(UNIT=IU_MIE, FILE=FILE_OUT, IOSTAT=IOS
     &   , STATUS='UNKNOWN'
     &   )
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(/, 3A)') '*** ERROR: FILE ', FILE_OUT
     &      , ' COULD NOT BE OPENED.'
         IERR=I_ERR_IO
         STOP
      ENDIF
!
!     OPEN A FILE FOR THE PHASE FUNCTION.
      IF (L_PHASE) THEN
!        FIND THE LENGTH OF THE OUTPUT FILENAME.
         LENGTH_OUT=LEN(FILE_OUT)
90       IF (FILE_OUT(LENGTH_OUT: LENGTH_OUT).EQ.' ') THEN
            LENGTH_OUT=LENGTH_OUT-1
            GOTO 90
         ENDIF
         FILE_PHASE(1: LENGTH_OUT+6)=FILE_OUT(1: LENGTH_OUT)//'.phase'
         WRITE(IU_STDOUT, '(/A)') 'PHASE FUNCTION WILL BE WRITTEN TO '
     &      //'FILE:'
         WRITE(IU_STDOUT, '(6X, A)') FILE_PHASE
         OPEN(UNIT=IU_PHASE, FILE=FILE_PHASE, IOSTAT=IOS
     &      , STATUS='UNKNOWN'
     &      )
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(/, 3A)') '*** ERROR: FILE ', FILE_PHASE
     &         , ' COULD NOT BE OPENED.'
            IERR=I_ERR_IO
            STOP
         ENDIF
      ENDIF
!
!     SET THE COMPONENT INDICES FOR THE FILE.
      WRITE(IU_USER, '(/A/)')
     &   'ENTER THE TYPE NUMBER FOR THE SCATTERER.'
      DO I=0, NPD_SCATTER_TYPE
         WRITE(IU_USER, '(3X, I3, A2, 3X, A20)')
     &      I, '.)', NAME_SCATTER_TYPE(I)
      ENDDO
      WRITE(IU_USER, '(/)')
110   READ(IU_STDIN, *, IOSTAT=IOS) I_SCATTER_TYPE
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(A)') '+++ ILLEGAL RESPONSE:'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RETYPE.'
            GOTO 110
         ENDIF
      ENDIF
      IF ( (I_SCATTER_TYPE.LT.0)
     &   .OR.(I_SCATTER_TYPE.GT.NPD_SCATTER_TYPE) ) THEN
         WRITE(IU_ERR, '(A)') '+++ RESPONSE IS OUT OF RANGE: '
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RETYPE.'
            GOTO 110
         ENDIF
      ENDIF
!
      IF (I_SCATTER_TYPE.EQ.IP_TYPE_UNASSIGNED) THEN
!        TYPE 0 IS UNASSIGNED TO ALLOW USE OF THE PROGRAM AS A PURE MIE
!        CODE FOR DIFFERENT SUBSTANCES WHEN THERE IS NO INTENTION OF
!        PROCEEDING TO ASSEMBLE A SPECTRAL FILE.
         I_COMPONENT=0
         NAME_COMPONENT='UNASSIGNED          '
      ELSE IF (I_SCATTER_TYPE.EQ.IP_TYPE_AEROSOL) THEN
         WRITE(IU_USER, '(/A)')
     &      'SPECIFY THE TYPE OF AEROSOL COMPONENT:'
         DO I=1, NPD_AEROSOL_COMPONENT
            WRITE(IU_USER, '(3X, I3, A2, 3X, A20)')
     &         I, '.)', NAME_AEROSOL_COMPONENT(I)
         ENDDO
         WRITE(IU_USER, '(/)')
111      READ(IU_STDIN, *, IOSTAT=IOS) I_COMPONENT
         IF ( (IOS.NE.0)
     &      .OR.(I_COMPONENT.GT.NPD_AEROSOL_COMPONENT) ) THEN
            WRITE(IU_ERR, '(A)') '+++ ILLEGAL RESPONSE:'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RETYPE.'
               GOTO 111
            ENDIF
         ENDIF
         IF ( (I_COMPONENT.LE.0)
     &      .OR.(I_COMPONENT.GT.NPD_AEROSOL_COMPONENT) ) THEN
            WRITE(IU_ERR, '(A)') '+++ RESPONSE IS OUT OF RANGE: '
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RETYPE.'
               GOTO 111
            ENDIF
         ENDIF
         NAME_COMPONENT=NAME_AEROSOL_COMPONENT(I_COMPONENT)
      ENDIF
!
!     CALCUALTE THE EFFECTIVE RADIUS AND VOLUME FRACTION.
!     WHEN THE EFFECTS OF HUMIDITY ARE INCLUDED THESE REFER TO THE
!     DRY DISTRIBUTION.
      IF ( (I_DISTRIBUTION.EQ.IP_EXTERNAL).AND.(N_SIZE.EQ.1) ) THEN
         IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
            RADIUS_EFF=RADIUS(N_SIZE)
            VOLUME_FRACTION=4.0D+00*PI*NUMBER(N_SIZE)
     &         *RADIUS(N_SIZE)**3/3.0D+00
            NUMBER_TOTAL=NUMBER(N_SIZE)
         ELSE IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_HEXCYL) THEN
            RADIUS_EFF=RADIUS(N_SIZE)
            VOLUME_FRACTION=SQRT(2.7D+01)*NUMBER(N_SIZE)
     &         *RADIUS(N_SIZE)**3*ASPECT(N_SIZE)
            NUMBER_TOTAL=NUMBER(N_SIZE)
         ENDIF
      ELSE
         CALL SIZE_INTEGRAL(IERR, RADIUS_0
     &      , 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
     &      , NUMBER_TOTAL, RADIUS_EFF, VOLUME_FRACTION)
         IF (IERR.NE.I_NORMAL) THEN
            IERR=I_ERR_FATAL
            STOP
         ENDIF
      ENDIF
      IF (L_HUMID) THEN
c         PRINT*,'SETTING DRY DISTRIBUTION PROPS'
         VOLUME_FRACTION_DRY=VOLUME_FRACTION
         RADIUS_EFF_DRY=RADIUS_EFF
         NUMBER_TOTAL_DRY=NUMBER_TOTAL
      ENDIF
!
      IF (L_PHASE) THEN
         WRITE(IU_PHASE, '(27X, A///)') 'SCATTERING PHASE FUNCTIONS'
      ENDIF
!
!     STORE THE PARAMETERS OF THE DRY DISTRIBUTION HERE.
!     AT PRESENT IT IS ASSUMED THAT A LOGNORMAL DISTRIBUTION IS USED.
!
      IF (L_HUMID) THEN
         IF (I_DISTRIBUTION.NE.IP_LOG_NORMAL) THEN
            WRITE(IU_ERR, '(/A)') '*** ERROR: AT PRESENT ONLY A LOG-'
     &         //'NORMAL DISTRIBUTION MAY BE USED FOR MOIST AEROSOLS.'
            STOP
         ENDIF
         LN_R0_LN_DRY = LN_R0_LN(1)
         LN_SIGMA_LN_DRY = LN_SIGMA_LN(1)
      ENDIF
!
!     BEGIN LOOP OVER HUMIDITY VALUES.
!
      DO 676 K = 1, N_HUMIDITIES
         
         IF (L_HUMID) THEN
c            PRINT*,K
!           ADJUST SIZE OF PARTICLES TO ALLOW FOR THE EFFECT OF 
!           HUMIDITY. NOTE THAT IT IS ASSUMED THAT BOTH THE DRY AND 
!           MOIST SIZE DISTRIBUTIONS ARE LOGNORMAL.
!
            ICHEM_TYPE = I_COMPONENT
            CALL GROW_PARTICLES(HUMIDITY(K),LN_R0_LN,LN_SIGMA_LN,
     &         ICHEM_TYPE,GROWTH_FACTOR,HUMIDITY(N_HUMIDITIES-1),IERR)
! !           PRINT*,LN_R0_LN,LN_SIGMA_LN,GROWTH_FACTOR
            IF (IERR.NE.I_NORMAL) THEN
               IERR=I_ERR_FATAL
               STOP
            ENDIF
!
            CALL SIZE_INTEGRAL(IERR, RADIUS_0
     &         , 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
     &         , NUMBER_TOTAL, RADIUS_EFF, VOLUME_FRACTION
     &         )
            IF (IERR.NE.I_NORMAL) THEN
               IERR=I_ERR_FATAL
               STOP
            ENDIF
!
         ENDIF
!
!        WRITE THE HEADERS TO THE FILE
!
         WRITE(IU_MIE, '(/, A13, I5)') '*FILE TYPE = ', I_OUTPUT_TYPE
         WRITE(IU_MIE, '(/, 16X, A, //)')
     &      'SINGLE SCATTERING PARAMETERS.'
         WRITE(IU_MIE, '(A23, I5, 3X, A1, A40, A1)')
     &      'SCATTERING ALGORITHM = ', I_SCATTER_ALGORITHM
     &      , '(', NAME_SCATTER_ALGORITHM(I_SCATTER_ALGORITHM), ')'
         WRITE(IU_MIE, '(A23, I5, 3X, A1, A20, A1)')
     &      'SCATTERING TYPE      = ', I_SCATTER_TYPE
     &      , '(', NAME_SCATTER_TYPE(I_SCATTER_TYPE), ')'
         IF (I_SCATTER_TYPE.EQ.IP_TYPE_AEROSOL) THEN
            WRITE(IU_MIE, '(A23, I5, 3X, A1, A20, A1)')
     &         'COMPONENT            = ', I_COMPONENT
     &         , '(', NAME_COMPONENT, ')'
         ENDIF
         IF (L_HUMID) THEN
            WRITE(IU_MIE,'(/,A11,F7.5)') 'HUMIDITY = ',HUMIDITY(K)
         ENDIF
         WRITE(IU_MIE, '(/, A18, 1X, 1PE12.5, 1X, A3)')
     &      'NUMBER DENSITY   =', NUMBER_TOTAL, 'm-3'
         WRITE(IU_MIE, '(A18, 1X, 1PE12.5, 1X, A1)')
     &      'EFFECTIVE RADIUS =', RADIUS_EFF, 'm'
         IF (L_HUMID) THEN
            WRITE(IU_MIE, '(A35, 1X, 1PE12.5, 1X, A1)')
     &         'EFFECTIVE RADIUS OF DRY COMPONENT ='
     &         , RADIUS_EFF_DRY, 'm'
         ENDIF
         WRITE(IU_MIE, '(A18, 1X, 1PE12.5)')
     &      'VOLUME FRACTION  =', VOLUME_FRACTION
         IF (L_HUMID) THEN
            WRITE(IU_MIE, '(A35, 1X, 1PE12.5)')
     &         'VOLUME FRACTION OF DRY COMPONENT  ='
     &         , VOLUME_FRACTION_DRY
         ENDIF
         WRITE(IU_MIE, '(/, 4X, A14, 6X, A16, 4X, A16, 4X, A9)')
     &      'Wavelength (m)', 'Absorption (m-1)'
     &      , 'Scattering (m-1)', 'Asymmetry'
!
!
!        LOOP OVER WAVELENGTHS.
!
         DO I=1, N_WAVELENGTH
!
!           FIND THE REFRACTIVE INDEX FOR THE CURRENT WAVELENGTH.
!
            N_RELATIVE=REFRACTIVE_INDEX(IERR, WAVELENGTH(I)
     &         , N_REFRACT, WAVELENGTH_REFRACT
     &         , RE_REFRACT, IM_REFRACT
     &         , D2_RE_REFRACT, D2_IM_REFRACT
     &         )
!
!           FOR MOIST AEROSOL, THE EFFECT OF DILUTION IS INCLUDED.
!           IT IS ASSUMED THAT ALL PARTICLES GROW BY THE SAME FACTOR.
!
            IF (L_HUMID) THEN
               N_RELATIVE_WATER = REFRACTIVE_INDEX(IERR, WAVELENGTH(I)
     &            , N_REFRACT_WATER, WAVELENGTH_REFRACT_WATER
     &            , RE_REFRACT_WATER, IM_REFRACT_WATER
     &            , D2_RE_REFRACT_WATER, D2_IM_REFRACT_WATER
     &            )
!
!              DILUTION PRODUCES A VOLUME-WEIGHTED AVERAGE OF THE
!              COMPLEX REFRACTIVE INDEX.
!
               DILUTION_FACTOR = 1.0D+00/(GROWTH_FACTOR)**3
               N_RELATIVE = N_RELATIVE*DILUTION_FACTOR +
     &                 N_RELATIVE_WATER*(1.0D+00 - DILUTION_FACTOR)
               print*,wavelength(i),n_relative
            ENDIF
!
!
            IF ( (I_DISTRIBUTION.EQ.IP_EXTERNAL).AND.
     &         (N_SIZE.EQ.1) ) THEN
               SIZE=2.0D+00*PI*RADIUS_0/WAVELENGTH(I)
               IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE) THEN
                  CALL MIE_SCATTER(IERR, SIZE, N_RELATIVE
     &               , Q_SCATTER, Q_EXT, ASYMMETRY
     &               , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     STOP
                  ENDIF
                  EXTINCTION=PI*NUMBER(N_SIZE)*Q_EXT*RADIUS(N_SIZE)**2
                  SCATTERING=PI*NUMBER(N_SIZE)
     &               *Q_SCATTER*RADIUS(N_SIZE)**2
               ELSE IF (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL)
     &            THEN
                  CALL TAKANO_LIOU_SCATTER(WAVELENGTH(I)
     &               , RADIUS(N_SIZE), ASPECT(N_SIZE)
     &               , N_RELATIVE
     &               , XS_SCATTER, XS_EXT, ASYMMETRY
     &               )
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     STOP
                  ENDIF
                  EXTINCTION=NUMBER(N_SIZE)*XS_EXT
                  SCATTERING=NUMBER(N_SIZE)*XS_SCATTER
               ENDIF
               
               
            ELSE IF ( (I_DISTRIBUTION.EQ.IP_EXTERNAL).AND.
     &          (N_SIZE.NE.1).AND.(L_ELLIE)) THEN

               DO ISIZE=1,N_SIZE
               SIZE=2.0D+00*PI*RADIUS(ISIZE)/WAVELENGTH(I)
               IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE) THEN
                  CALL MIE_SCATTER(IERR, SIZE, N_RELATIVE
     &               , Q_SCATTER, Q_EXT, ASYMMETRY
     &               , L_PHASE, N_ANGLE, MU_ANGLE, S1, S2)
                  IF (IERR.NE.I_NORMAL) THEN
                     IERR=I_ERR_FATAL
                     STOP
                  ENDIF
                  EXTINCTION=PI*NUMBER(ISIZE)*Q_EXT*RADIUS(ISIZE)**2
                  SCATTERING=PI*NUMBER(ISIZE)
     &               *Q_SCATTER*RADIUS(ISIZE)**2
               END IF

!   put integral in here********************
               

! ******************************************

               WRITE(IU_PHASE, '(/A13, 1PE12.5, 6X, A17, 1PE12.5)')
     &              'WAVELENGTH = ', WAVELENGTH(I)
     &              , 'SIZE PARAMETER = ', SIZE
               WRITE(IU_PHASE, '(6X, A19, 2(1X, 1PE12.5), /)')
     &              'REFRACTIVE INDEX = ', N_RELATIVE
               WRITE(IU_PHASE, '(9X, A5, 8X, A6, 16X, A2, 23X, A2)')
     &              'ANGLE', 'COSINE', 'S1', 'S2'
               DO J=1, N_ANGLE
                  WRITE(IU_PHASE, '(6X, 1PE10.3, 4X, 1PE10.3, 3X
     &                 , 2(1X, 1PE10.3), 3X, 2(1X, 1PE10.3))')
     &                 ANGLE(J), MU_ANGLE(J), S1(J), S2(J)
               ENDDO
               WRITE(IU_PHASE, '(/)')

               END DO
               L_ELLIE=.FALSE.
           

            
            ELSE

!              NOTICE THAT THE EFFECTIVE RADIUS IS NOW NORMALLY SUPPLIED
!              AS THE TYPICAL VALUE: THIS SHOULD BE MORE EFFICIENT. 
!              HOWEVER, IN THE CASE OF NON-SPHERICAL PARTICLES THE 
!              DISTRIBUTION MAY BE 0 AT THE EFFECTIVE RADIUS.
               IF (I_SHAPE_PARTICLE.EQ.IP_SHAPE_SPHERE) THEN
                  RADIUS_TYPICAL=RADIUS_EFF
               ELSE
                  NUMBER_0=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_0)
                  NUMBER_EFF=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_EFF)
                  IF (NUMBER_EFF*RADIUS_EFF**2.GT.1.0D-02*NUMBER_0
     &               *RADIUS_0**2) THEN
                     RADIUS_TYPICAL=RADIUS_EFF
                  ELSE
                     RADIUS_TYPICAL=RADIUS_0
                  ENDIF
               ENDIF
               IF ( (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_MIE).OR.
     &              (I_SCATTER_ALGORITHM.EQ.IP_TAKANO_LIOU_HEXCYL) )
     &            THEN
                  print*,RADIUS_TYPICAL
                  CALL SCATTER_INTEGRAL(IERR, RADIUS_TYPICAL
     &               , I_DISTRIBUTION, I_SHAPE_PARTICLE
     &               , I_SCATTER_ALGORITHM
     &               , 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
     &               , WAVELENGTH(I), N_RELATIVE
     &               , L_PHASE, N_ANGLE, MU_ANGLE
     &               , EXTINCTION, SCATTERING, ASYMMETRY
     &               , S1, S2)
!                  print*,S1,S2
               ELSE IF (I_SCATTER_ALGORITHM.EQ.IP_ALGORITHM_ADT) THEN
                  EXTINCTION=0.0D+00
                  SCATTERING=0.0D+00
                  DO J=1, N_MODE
                     SIZE=2.0D+00*PI*RM_MG(J)**BETA_MG(J)/WAVELENGTH(I)
                     CALL ADT_INTEGRAL(ALPHA_MG(J), SIZE, N_RELATIVE
     &                  , Q_EXT, Q_SCATTER)
                     EXTINCTION=EXTINCTION+WEIGHT_MODE(J)
     &                  +Q_EXT*7.5D-01*VOLUME_FRACTION/RADIUS_EFF
                     SCATTERING=SCATTERING+WEIGHT_MODE(J)
     &                  +Q_SCATTER*7.5D-01*VOLUME_FRACTION/RADIUS_EFF
                  ENDDO
!                 ADT IN ITS PRESENT FORM DOES NOT GIVE AN 
!                 ASYMMETRY FACTOR: ZERO IS PRINTED OUT AS A FLAG.
                  ASYMMETRY=0.0D+00
               ENDIF
               IF (IERR.NE.I_NORMAL) THEN
                  IERR=I_ERR_FATAL
                  STOP
               ENDIF
            ENDIF
!
            IF ( (L_PHASE).AND.(L_ELLIE)) THEN
               WRITE(IU_PHASE, '(/A13, 1PE12.5, 6X, A17, 1PE12.5)')
     &            'WAVELENGTH = ', WAVELENGTH(I)
     &            , 'SIZE PARAMETER = ', SIZE
               WRITE(IU_PHASE, '(6X, A19, 2(1X, 1PE12.5), /)')
     &            'REFRACTIVE INDEX = ', N_RELATIVE
               WRITE(IU_PHASE, '(9X, A5, 8X, A6, 16X, A2, 23X, A2)')
     &            'ANGLE', 'COSINE', 'S1', 'S2'
               DO J=1, N_ANGLE
                  WRITE(IU_PHASE, '(6X, 1PE10.3, 4X, 1PE10.3, 3X
     &               , 2(1X, 1PE10.3), 3X, 2(1X, 1PE10.3))')
     &               ANGLE(J), MU_ANGLE(J), S1(J), S2(J)
               ENDDO
               WRITE(IU_PHASE, '(/)')
            ENDIF
!
            ABSORPTION=EXTINCTION-SCATTERING
!
            WRITE (IU_MIE, '(4(4X, 1PE12.5, 4X))')
     &         WAVELENGTH(I), ABSORPTION, SCATTERING, ASYMMETRY
!
         ENDDO
!
         WRITE(IU_MIE, '(1X)')
!
!        REMEMBER TO RESTORE THE DRY VALUES OF THE DISTRIBUTION
!        PARAMETERS BEFORE GOING AROUND THE LOOP AGAIN (OR INDEED
!        LEAVING THE LOOP.)
!
         IF (L_HUMID) THEN
            LN_R0_LN(1) = LN_R0_LN_DRY
            LN_SIGMA_LN(1) = LN_SIGMA_LN_DRY
         ENDIF
!
!        END OF HUMIDITY LOOP.
!
676   CONTINUE
!
!
      CLOSE(IU_MIE)
!
!
      STOP
      END


