!+ Program to average scattering properties across bands.
!
! Method:
!	A file containing blocks of monochromatic single 
!	scattering properties is read in. These monochromatic
!	values are averaged across the bands given in a 
!	spectral file. The averaged values may be written to
!	a file or fitted using a recognized parametrization.
!
! 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_AVERAGE
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/scatter_pp.pcf'
      INCLUDE 'incl/method_weight.pcf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/def_sct_io.icf'
      INCLUDE 'incl/file_type.pcf'
      INCLUDE 'incl/error.pcf'
!
!     DECLARATION OF VARIABLES.
      CHARACTER
     &     FILE_SPECTRAL*80
!             NAME OF SPECTRAL FILE
     &   , NAME_TYPE*20
!             NAME OF TYPE OF SCATTERER
     &   , NAME_COMPONENT*20
!             NAME OF COMPONEMT
      INTEGER
     &     IERR
!             ERROR FLAG
     &   , N_SOLAR_POINTS
!             NUMBER OF POINTS IN SOLAR SPECT.
     &   , N_WAVELENGTH
!             NUMBER OF WAVELENGTHS IN FILE OF SCATTERING DATA
     &   , N_INT_WEIGHT
!             NUMBER OF WEIGHTING POINTS
     &   , I_WEIGHT
!             TYPE OF WEIGHTING FUNCTION
     &   , I_METHOD_WEIGHT
!             METHOD OF WEIGHTING
     &   , N_BLOCK
!             NUMBER OF DATA BLOCK
     &   , N_BLOCK_P1
!             NUMBER OF DATA BLOCK PLUS 1
     &   , I_BEGIN
!             INITIAL POINT IN BAND
     &   , I_BAND
!             INDEX OF BAND
     &   , I_SCATTER_TYPE(NPD_MIE_BLOCK)
!             TYPE OF SCATTERER IN EACH BLOCK
     &   , I_COMPONENT_BLOCK(NPD_MIE_BLOCK)
!             COMPONENTS IN EACH BLOCK
     &   , I
!             LOOP VARIABLE
     &   , K
!             LOOP VARIABLE
     &   , IPB
!             OFFSET LOOP VARIABLE
      INTEGER
     &     I_INPUT_TYPE
!             TYPE OF INPUT FILE
     &   , I_OUTPUT_TYPE
!             TYPE OF OUTPUT FILE
      LOGICAL
     &     L_END
!             END OF FILE FLAG
     &   , L_WRITE
!             WRITING FLAG
     &   , L_WAVE_NUMBER
!             FLAG FOR WEIGHTING BY WAVENUMBER
      REAL*8
     &     SOLAR_WAVELENGTH(NPD_SOLAR_POINTS)
!             WAVELENGTHS OF SOLAR RADIATION
     &   , SOLAR_IRRAD(NPD_SOLAR_POINTS)
!             SOLAR IRRADIANCES
     &   , NUMBER_TOTAL(NPD_MIE_BLOCK)
!             TOTAL NUMBER DENSITY
     &   , RADIUS_EFF(NPD_MIE_BLOCK)
!             EFFECTIVE RADIUS
     &   , VOLUME_FRACTION(NPD_MIE_BLOCK)
!             VOLUME FRACTION OF SPHERES
     &   , RADIUS_EFF_DRY(NPD_MIE_BLOCK)
!             EFFECTIVE RADIUS OF DRY PARTICLES
     &   , VOLUME_FRACTION_DRY(NPD_MIE_BLOCK)
!             VOLUME FRACTION OF DRY PARTICLES
     &   , WAVELENGTH(NPD_WAVELENGTH_SCAT)
!             WAVELENGTHS OF SCATTERING DATA
     &   , T_WEIGHT
!             TEMPERATURE FOR WEIGHTING
     &   , WEIGHT_COEFF(0: NPD_WAVELENGTH_SCAT)
!             COEFFICIENTS FOR WEIGHTING
     &   , WORK(0: NPD_WAVELENGTH_SCAT)
!             WORK ARRAY
!
      REAL*8
     &     HUMIDITY
!             HUMIDITY
      REAL*8
     &     ABSORPTION(NPD_WAVELENGTH_SCAT)
!             ABSORPTION
     &   , SCATTERING(NPD_WAVELENGTH_SCAT)
!             ABSORPTION
     &   , ASYMMETRY(NPD_WAVELENGTH_SCAT)
!             ABSORPTION
     &   , ABSORPTION_BAND(NPD_BAND, NPD_MIE_BLOCK)
!             BAND_AVERAGED ABSORPTION
     &   , SCATTERING_BAND(NPD_BAND, NPD_MIE_BLOCK)
!             BAND_AVERAGED ABSORPTION
     &   , ASYMMETRY_BAND(NPD_BAND, NPD_MIE_BLOCK)
!             BAND_AVERAGED ABSORPTION
!     VARIABLES FOR THICK WEIGHTING:
      REAL*8
     &     PHI
!             SQUARE ROOT OF (1-OMEGA)
     &   , PSI
!             SQUARE ROOT OF (1-G*OMEGA)
     &   , REFLECT_INF
!             MONOCHROMATIC REFLECTION COEFFICIENT FOR AN INFINITE CLOUD
     &   , EXTINCTION_TOT
!             TOTAL MONOCHROMATIC EXTINCTION COEFFICIENT
     &   , MEAN_REFLECT_INF
!             MEAN INFINITE REFLECTION COEFFICIENT
     &   , MEAN_EXTINCTION
!             MEAN TOTAL EXTINCTION
     &   , MEAN_SCATTERING
!             MEAN SCATTERING
     &   , MEAN_SCATTERING_ASYMMETRY
!             MEAN OF PRODUCT OF SCATTERING AND ASYMMETRY
     &   , OMEGA_EFF
!             EFFECTIVE ALBEDO OF SINGLE SCATTERING
!
      REAL*8
     &     EXCLUDE_SHORTER(NPD_EXCLUDE)
!             LOWER LIMITS OF EXCLUDED INTERVALS
     &   , EXCLUDE_LONGER(NPD_EXCLUDE)
!             UPPER LIMITS OF EXCLUDED INTERVALS
!
!     DECLARE VARIABLES FOR THE SPECTRAL FILE.
      INCLUDE 'dec_spec.finc'
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     IN_SPECTRUM, OPEN_FILE_IN, SELECT_WEIGHT_SCATTER
     &   , OPEN_AVERAGE, READ_SCATTER_BLOCK, WEIGHTINGS
     &   , WRITE_AVERAGE, CLOUD_FIT, ICE_CLOUD_FIT
!
!
      DATA IERR/I_NORMAL/, L_END/.FALSE./, L_WAVE_NUMBER/.FALSE./
!
!
!
!     OBTAIN THE BASIC SPECTRAL DATA.
!     REQUEST THE NAME OF THE FILE CONTAINING SPECTRAL DATA.
100   WRITE(IU_USER, '(/A/)')
     &   'GIVE THE NAME OF THE FILE CONTAINING THE SPECTRAL DATA.'
      READ(IU_STDIN, '(A)') FILE_SPECTRAL
      write(*,*)'opened spectral file'
!
      INCLUDE 'read_spec.finc'
      IF (IERR.NE.I_NORMAL) THEN
         IF (IERR.EQ.I_ERR_EXIST) THEN
!           RESET ERROR FLAG.
            IERR=I_NORMAL
            GOTO 100
         ELSE
            STOP
         ENDIF
      ENDIF
     
!
!     OPEN THE SCATTERING FILES.
      CALL OPEN_FILE_IN(IERR, IU_MIE
     &   , 'GIVE THE NAME OF THE FILE OF SCATTERING '
     &   //'DATA TO BE AVERAGED.')
      IF (IERR.NE.I_NORMAL) STOP
      write(*,*)'opened scattering file'
!
!     SELECT THE WEIGHTING FUNCTION.
      CALL SELECT_WEIGHT_SCATTER(IERR, I_WEIGHT, T_WEIGHT
     &   , N_SOLAR_POINTS, SOLAR_WAVELENGTH, SOLAR_IRRAD
     &   , I_METHOD_WEIGHT)
      IF (IERR.NE.I_NORMAL) STOP
      write(*,*)'selected weights'
!
!     DETERMINE WHETHER THE AVERAGES ACROSS EACH BAND SHOULD BE WRITTEN
!     TO A FILE.
      CALL OPEN_AVERAGE(IERR, L_WRITE)
      IF (IERR.NE.I_NORMAL) STOP
      write(*,*)'opened average file'
!
!     EACH BLOCK OF DATA IN THE SCATTERING FILE IS READ AND AVERAGED
!     SEPARATELY.
      N_BLOCK=0
!
1        N_BLOCK_P1=N_BLOCK+1
         CALL READ_SCATTER_BLOCK(IERR, L_END
     &      , I_SCATTER_TYPE(N_BLOCK_P1), I_COMPONENT_BLOCK(N_BLOCK_P1)
     &      , I_INPUT_TYPE
     &      , NAME_TYPE, NAME_COMPONENT
     &      , HUMIDITY
     &      , NUMBER_TOTAL(N_BLOCK_P1)
     &      , RADIUS_EFF(N_BLOCK_P1), VOLUME_FRACTION(N_BLOCK_P1)
     &      , RADIUS_EFF_DRY(N_BLOCK_P1)
     &      , VOLUME_FRACTION_DRY(N_BLOCK_P1)
     &      , N_WAVELENGTH, WAVELENGTH
     &      , ABSORPTION, SCATTERING, ASYMMETRY)
         IF (IERR.NE.I_NORMAL) STOP
         write(*,*)N_BLOCK,'read data'
         write(*,*)L_END
!
         IF (.NOT.L_END) THEN
            N_BLOCK=N_BLOCK_P1
            write(*,*)N_BLOCK,'started averaging'
!           AVERAGE ACROSS EACH BAND.
            DO I_BAND=1, N_BAND
!              CALCULATE THE WEIGHTINGS.
               DO K=1, N_BAND_EXCLUDE(I_BAND)
                  EXCLUDE_SHORTER(K)
     &               =WAVE_LENGTH_SHORT(INDEX_EXCLUDE(K, I_BAND))
                  EXCLUDE_LONGER(K)
     &               =WAVE_LENGTH_LONG(INDEX_EXCLUDE(K, I_BAND))
               ENDDO
               CALL WEIGHTINGS(IERR
     &            , I_WEIGHT, L_WAVE_NUMBER
     &            , N_SOLAR_POINTS, SOLAR_WAVELENGTH, SOLAR_IRRAD
     &            , T_WEIGHT
     &            , WAVE_LENGTH_SHORT(I_BAND), WAVE_LENGTH_LONG(I_BAND)
     &            , N_BAND_EXCLUDE(I_BAND)
     &            , EXCLUDE_SHORTER, EXCLUDE_LONGER
     &            , N_WAVELENGTH, WAVELENGTH
     &            , I_BEGIN, N_INT_WEIGHT, WEIGHT_COEFF
     &            , WORK
     &            )
               IF (IERR.NE.I_NORMAL) STOP
!

               print*,weight_coeff

!              CALCULATE THE MEAN SCATTERING.
               MEAN_SCATTERING=0.0D+00
               MEAN_SCATTERING_ASYMMETRY=0.0D+00
               DO K=0, N_INT_WEIGHT
                  MEAN_SCATTERING=MEAN_SCATTERING+WEIGHT_COEFF(K)
     &               *SCATTERING(I_BEGIN+K)
                  MEAN_SCATTERING_ASYMMETRY=MEAN_SCATTERING_ASYMMETRY
     &               +WEIGHT_COEFF(K)*SCATTERING(I_BEGIN+K)
     &               *ASYMMETRY(I_BEGIN+K)
               ENDDO
               ASYMMETRY_BAND(I_BAND, N_BLOCK)
     &            =MEAN_SCATTERING_ASYMMETRY/MEAN_SCATTERING
!              THE ABSORPTION MAY BE WEIGHTED IN DIFFERENT WAYS FOR
!              DIFFERENT PURPOSES.
               IF (I_METHOD_WEIGHT.EQ.IP_WEIGHT_THIN) THEN
                  SCATTERING_BAND(I_BAND, N_BLOCK)
     &               =MEAN_SCATTERING
                  ABSORPTION_BAND(I_BAND, N_BLOCK)=0.0D+00
                  DO K=0, N_INT_WEIGHT
                     ABSORPTION_BAND(I_BAND, N_BLOCK)
     &                  =ABSORPTION_BAND(I_BAND, N_BLOCK)
     &                  +WEIGHT_COEFF(K)*ABSORPTION(I_BEGIN+K)
                  ENDDO
               ELSE IF (I_METHOD_WEIGHT.EQ.IP_WEIGHT_THICK) THEN
!                 THE REFLECTION COEFFICIENT OF AN INFINITE CLOUD IS
!                 CALCULATED AND WEIGHTED, AND THE TOTAL EXTINCTION IS
!                 WEIGHTED AS OPPOSED TO THE INDIVIDUAL COEFFICIENTS.
                  MEAN_EXTINCTION=0.0D+00
                  MEAN_REFLECT_INF=0.0D+00
                  DO K=0, N_INT_WEIGHT
                     IPB=K+I_BEGIN
                     EXTINCTION_TOT=ABSORPTION(IPB)+SCATTERING(IPB)
                     PHI=SQRT(ABSORPTION(IPB)/EXTINCTION_TOT)
                     PSI=SQRT((ABSORPTION(IPB)
     &                  +SCATTERING(IPB)
     &                  *(1.0D+00-ASYMMETRY(IPB)))
     &                  /EXTINCTION_TOT)
                     REFLECT_INF=(PSI-PHI)/(PSI+PHI)
                     MEAN_EXTINCTION=MEAN_EXTINCTION
     &                  +WEIGHT_COEFF(K)*EXTINCTION_TOT
                     MEAN_REFLECT_INF=MEAN_REFLECT_INF
     &                  +WEIGHT_COEFF(K)*REFLECT_INF
                  ENDDO
!                 DETERMINE THE MEAN VALUE OF OMEGA.
                  OMEGA_EFF=4.0D+00*MEAN_REFLECT_INF
     &               /((1.0D+00+MEAN_REFLECT_INF)**2
     &               -ASYMMETRY_BAND(I_BAND, N_BLOCK)
     &               *(1.0D+00-MEAN_REFLECT_INF)**2)
                  ABSORPTION_BAND(I_BAND, N_BLOCK)=(1.0D+00-OMEGA_EFF)
     &               *MEAN_EXTINCTION
                  SCATTERING_BAND(I_BAND, N_BLOCK)=OMEGA_EFF
     &               *MEAN_EXTINCTION
               ENDIF
            ENDDO
!
            IF (L_WRITE) THEN
               write(*,*)'writing'
!              SET THE OUTPUT TYPE FOR THE FILE OF RESULTS.
               IF (I_INPUT_TYPE.EQ.IT_FILE_MIE_DRY) THEN
                  I_OUTPUT_TYPE=IT_FILE_AVE_MIE_DRY
               ELSE IF (I_INPUT_TYPE.EQ.IT_FILE_MIE_HUMID) THEN
                  I_OUTPUT_TYPE=IT_FILE_AVE_MIE_HUMID
               ENDIF
               CALL WRITE_AVERAGE(N_BAND, N_BLOCK
     &            , I_SCATTER_TYPE(N_BLOCK), I_COMPONENT_BLOCK(N_BLOCK)
     &            , I_OUTPUT_TYPE, HUMIDITY
     &            , NAME_TYPE, NAME_COMPONENT
     &            , NUMBER_TOTAL(N_BLOCK)
     &            , RADIUS_EFF(N_BLOCK)
     &            , VOLUME_FRACTION(N_BLOCK)
     &            , RADIUS_EFF_DRY(N_BLOCK)
     &            , VOLUME_FRACTION_DRY(N_BLOCK)
     &            , ABSORPTION_BAND(1, N_BLOCK)
     &            , SCATTERING_BAND(1, N_BLOCK)
     &            , ASYMMETRY_BAND(1, N_BLOCK))
            ENDIF
!
!           READ IN THE NEXT BLOCK OF DATA
            GOTO 1
         ENDIF
!
!     WHERE A FIT IS REQUIRED, A CONJUGATE GRADIENT SCHEME WILL PERFORM
!     IT. THE DATA ARE CHECKED TO ENSURE THEAT THE EY ARE ALL FOR 
!     THE SAME TYPE OF SCATTERER.
      DO I=2, N_BLOCK
         IF (I_SCATTER_TYPE(I).NE.I_SCATTER_TYPE(1)) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR: BLOCKS DO NOT ALL REPRESENT THE SAME '
     &         //'TYPE OF SCATTERER.'
            WRITE(IU_ERR, '(/A)')
     &         'FITTING WILL NOT BE PERFORMED.'
            IERR=I_ERR_FATAL
            STOP
          ENDIF
      ENDDO
!
      IF (I_SCATTER_TYPE(N_BLOCK).EQ.IP_TYPE_DROPLET) THEN
         CALL CLOUD_FIT(IERR, N_BAND, N_BLOCK
     &      , VOLUME_FRACTION, RADIUS_EFF
     &      , SCATTERING_BAND, ABSORPTION_BAND, ASYMMETRY_BAND
     &      )
      ELSE IF (I_SCATTER_TYPE(N_BLOCK).EQ.IP_TYPE_ICE) THEN
         CALL ICE_CLOUD_FIT(IERR, N_BAND, N_BLOCK
     &      , VOLUME_FRACTION, RADIUS_EFF
     &      , SCATTERING_BAND, ABSORPTION_BAND, ASYMMETRY_BAND
     &      )
      ENDIF
!
!
      STOP
      END
