!+ Program to calclate average transmissions from grey data.
!
! Method:
!	A spectral file is read in to specify the spectral bands
!	and a file of cross-sectional data is read.
!	Subroutines are then called to perform the averaging and
!	the final values are written.
!
! 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 GREY_AVERAGE
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/def_arch_io.icf'
      INCLUDE 'incl/def_pp_io.icf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/dimensions_arch.ucf'
      INCLUDE 'incl/arch_record_size.pcf'
      INCLUDE 'incl/ckd1_table_size.ecf'
      INCLUDE 'incl/ckd1_table_parm.ecf'
      INCLUDE 'incl/ckd1_parm.ecf'
      INCLUDE 'incl/record_size_trans.pcf'
      INCLUDE 'incl/gas_list.pcf'
      INCLUDE 'incl/continuum.pcf'
      INCLUDE 'incl/physical_constants_pp.ccf'
      INCLUDE 'incl/error.pcf'
!
!
!     DECLARATION OF VARIABLES.
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
      EXTERNAL
     &     LOCK_CODE
!
      CHARACTER
     &     FILE_SPECTRAL*80
!             NAME OF SPECTRAL FILE
     &   , CHAR_YN*1
!             CHARACTER RESPONSE VARIABLE
     &   , TRANS_TYPE*4
!             TYPE OF TRANSMISSION
     &   , DATE(NPD_ARCH_ENTRY)*9
!             DATES OF ARCHIVAL ENTRIES
      INTEGER
     &     IERR
!             ERROR FLAG
     &   , IOS
!             I/O ERROR FLAG
      INTEGER
     &     N_SOLAR_POINTS
!             SIZE OF SOLAR SPECTRUM
     &   , I_INDEX_GAS(NPD_GASES)
!             INDICIES OF GASES
     &   , I_SPECIES
!             SPECIES OF GAS
     &   , N_BAND_LIST
!             NUMBER OF BANDS IN LIST
     &   , BAND_LIST(NPD_BAND)
!             LIST OF BANDS FOR ARCHIVE
     &   , I_BAND
!             BAND BEING CONSIDERED
     &   , I_SET
!             SET OF DATA
     &   , I_PLACE
!             PLACE IN SET
      INTEGER
     &     I_VERSION
!             VERSION OF ARCHIVE
     &   , N_ENTRY
!             NUMBER OF ENTRIES IN THE ARCHIVE
     &   , I_TYPE_ENTRY(NPD_ARCH_ENTRY)
!             TYPES OF GASES
     &   , N_TRANS_RECORD(NPD_ARCH_ENTRY)
!             NUMBER OF RECORDS IN ENTRY
     &   , I1_RECORD(NPD_ARCH_ENTRY)
!             POINTER TO FIRST RECORD OF ENTRY
     &   , N_TRANS_ENTRY(NPD_ARCH_ENTRY)
!             NUMBER OF TRANSMISSIONS IN ENTRY
     &   , IDUMMY(NPD_ARCH_ENTRY)
!             DUMMY ARRAY
     &   , N_PER_ARCH_RECORD
!             NUMBER OF TRIPLES PER RECORD
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , K
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , I_BEGIN
!             INITIAL POINT JUST OUTSIDE BAND
     &   , I_WEIGHT
!             FLAG FOR METHOD OF WEIGHTING
     &   , N_INT_WEIGHT
!             NUMBER OF WEIGHTING INTERVALS
     &   , FIRST_BAND
!             FIRST BAND TO CONSIDER
     &   , LAST_BAND
!             LAST BAND TO CONSIDER
      INTEGER
     &     N_SET
!             NUMBER OF SETS
     &   , LENGTH_SET(NPD_SET)
!             LENGTHS OF EACH SET
     &   , BAND_SET(NPD_SET)
!             BAND FOR EACH SET
     &   , INDEX_SET(NPD_SET)
!             INDEX NUMBER OF GAS IN SET
      INTEGER
     &     N_PER_TRANS_RECORD
!             NUMBER OF ENTRIES PER TRANSMISSION RECORD
     &   , I_FIRST_RECORD
!             FIRST RECORD OF SET
     &   , I_REC
!             POINTER TO SCRATCH RECORD
      LOGICAL
     &     LDUMMY(NPD_ARCH_ENTRY)
!             DUMMY ARRAY
     &   , L_CONVOL_FRN
!             FLAG FOR CONVOLVING FOREIGN CONTINUUM
      REAL*8
     &     P_SET(NPD_SET)
!             PRESSURES IN SETS
     &   , T_SET(NPD_SET)
!             TEMPERATURES IN SETS
     &   , U_SET(NPD_AMOUNT)
!             AMOUNTS IN SETS
     &   , TRANS_SET(NPD_AMOUNT)
!             TRANSMISSIONS IN SETS
     &   , DENSITY_BROAD(NPD_AMOUNT)
!             BROADENING DENISTY
      REAL*8
     &     WAVE_NUMBER_LOW
!             LOW WAVENUMBER LIMIT OF BAND
     &   , WAVE_NUMBER_HIGH
!             HIGH WAVENUMBER LIMIT OF BAND
     &   , WEIGHT_COEFF(0: NPD_FREQ_ENTRY)
!             WEIGHTING COEFFICIENTS
     &   , WORK(0: NPD_FREQ_ENTRY)
!             WORK ARRAY
     &   , DIFFUSIVITY_FACTOR
!             DIFFUSIVITY FACTOR FOR FLUXES
      REAL*8
     &     SOLAR_WAVELENGTH(NPD_SOLAR_POINTS)
!             POINTS OF SOLAR SPECTRUM
     &   , SOLAR_IRRAD(NPD_SOLAR_POINTS)
!             SOLAR IRRADIANCES
!
      REAL*8
     &     EXCLUDE_SHORTER(NPD_EXCLUDE)
!             LOWER LIMITS OF EXCLUDED INTERVALS
     &   , EXCLUDE_LONGER(NPD_EXCLUDE)
!             UPPER LIMITS OF EXCLUDED INTERVALS
!
      REAL*8
     &     RDUMMY(NPD_ARCH_ENTRY)
!             DUMMY ARRAY
     &   , P_ENTRY(NPD_ARCH_ENTRY)
!             PRESSURE
     &   , T_ENTRY(NPD_ARCH_ENTRY)
!             TEMPERATURE
     &   , AMOUNT_ENTRY(NPD_ARCH_ENTRY)
!             AMOUNT
     &   , PARTIAL_P_ENTRY(NPD_ARCH_ENTRY)
!             PARTIAL PRESSURE
     &   , FREQ_ENTRY_LOW(NPD_ARCH_ENTRY)
!             LOW FREQUENCY LIMIT OF ENTRY
     &   , FREQ_ENTRY_HIGH(NPD_ARCH_ENTRY)
!             HIGH FREQUENCY LIMIT OF ENTRY
     &   , D_FREQ_ENTRY(NPD_ARCH_ENTRY)
!             INCREMENT IN FREQUENCY OF ENTRIES
      REAL*8
     &     FREQ_LOW(NPD_FREQ_ENTRY)
!             LOW FREQUENCIES OF INTERVALS
     &   , FREQ_HIGH(NPD_FREQ_ENTRY)
!             HIGH FREQUENCIES OF INTERVALS
     &   , TRANS(NPD_FREQ_ENTRY)
!             TRANSMITTANCES IN INTERVALS
     &   , WAVENUMBER(NPD_FREQ_ENTRY)
!             WAVENUMBERS OF MEAN POINTS
     &   , TRANS_AVE_ENTRY
!             AVERAGED TRANSMITTANCE
     &   , WAVENUMBER_TABLE_FRN(NPD_TABLE_FRN)
!             WAVENUMBERS OF ENTRIES IN TABLE
     &   , C_H2O_FRN_SI(NPD_TABLE_FRN)
!             SI VALUES OF COEFFICIENTS IN TABLE
     &   , D2_C_H2O_FRN_SI(NPD_TABLE_FRN)
!             2ND DERIVATIVES OF COEFFICIENTS
     &   , CONVERSION_FACTOR
!             CONVERSION FACTOR TO SI UNITS
     &   , WORK(NPD_TABLE_SELF_260+NPD_TABLE_SELF_296+NPD_TABLE_FRN)
!             WORKING VARIABLE
     &   , CONTINUUM_TRANS(NPD_WAVENUMBER_TRANS)
!             CONTINUUM TRANSMISSIONS
!
      INCLUDE 'dec_spec.finc'
!
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     IN_SPECTRUM, OPEN_ARCHIVE_IN
     &   , READ_ARCHIVE_HEADER_0, READ_ARCHIVE_HEADER_1
     &   , READ_ARCHIVE_ENTRY_0, READ_ARCHIVE_ENTRY_1
     &   , SPLINE_FIT, SELECT_WEIGHT_TRANS, OPEN_TRANS_OUT
     &   , FIND_GAS_BAND, CONTINUUM_FRN_TRANS, LOCATE_ENTRY
     &   , WEIGHTINGS, ORDER_SET, WRITE_TRANS, CONVERT_AMOUNT
!
!
!
!     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
!
!     READ IN THE SPECTRAL FILE.
      INCLUDE 'read_spec.finc'
!
      IF (IERR.NE.I_NORMAL) THEN
         IF (IERR.EQ.I_ERR_EXIST) THEN
            IERR=I_NORMAL
            GOTO 100
         ELSE
            STOP
         ENDIF
      ENDIF
!
!     SET LIMITS ON THE RANGE OF BANDS TO BE CONSIDERED.
      WRITE(IU_USER, '(/A)')
     &   'ENTER THE RANGE OF BANDS TO BE CONSIDERED.'
200   READ(IU_STDIN, *, IOSTAT=IOS) FIRST_BAND, LAST_BAND
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(A)')
     &      '+++ UNRECOGNIZED RESPONSE: PLEASE REENTER.'
         GOTO 200
      ENDIF
      IF ( (FIRST_BAND.LT.1).OR.(LAST_BAND.GT.N_BAND) ) THEN
         WRITE(IU_ERR, '(A)')
     &      '+++ BAND OUT OF RANGE: PLEASE REENTER.'
         GOTO 200
      ENDIF
!
!     FORM THE LIST OF SPECIES FOR EACH TYPE.
      DO I=1, NPD_GASES
         I_INDEX_GAS(I)=0
         DO J=1, N_ABSORB
            IF (TYPE_ABSORB(J).EQ.I) THEN
               I_INDEX_GAS(I)=J
            ENDIF
         ENDDO
      ENDDO
!
!
!     OPEN A FILE OF CROSS-SECTIONAL DATA.
      CALL OPEN_FILE_IN(IERR, IU_FILE_IN
     &   , 'GIVE THE NAME OF THE FILE OF CROSS-SECTIONAL DATA.')
      IF (IERR.NE.I_NORMAL) STOP
!
!     READ THE FILE OF CROSS-SECTIONS.
      CALL READ_X_SECTION(IU_FILE_IN
     &   , N_TRANS_RECORD(I), I1_RECORD(I)
     &   , N_PER_ARCH_RECORD, N_TRANS_ENTRY(I)
     &   , FREQ_LOW, FREQ_HIGH, TRANS
     &   )
!
!     SET THE AMOUNTS OF ABSORBERS TO BE USED.
      CALL SET_AMOUNT_GREY
!
!     SELECT THE WEIGHTING FUNCTION
      CALL SELECT_WEIGHT_TRANS(IERR, I_WEIGHT
     &   , N_SOLAR_POINTS, SOLAR_WAVELENGTH, SOLAR_IRRAD)
      IF (IERR.NE.I_NORMAL) STOP
!
!     PREPARE THE OUTPUT FILES.
      CALL OPEN_TRANS_OUT(IERR, .FALSE., N_PER_TRANS_RECORD
     &   , IU_TRANSMISSION_SUMMARY, IU_TRANSMISSION)
      IF (IERR.NE.I_NORMAL) STOP
!
!     CONVERT TO FLUX TRANSMITTANCES IF REQUIRED.
      WRITE(IU_USER, '(/A)')
     &   'CONVERT TO FLUX TRANSMITTANCES? (Y/N)'
1     READ(IU_STDIN, '(A)') CHAR_YN
      IF ( (CHAR_YN.EQ.'Y').OR.(CHAR_YN.EQ.'y') ) THEN
         TRANS_TYPE='FLUX'
         WRITE(IU_USER, '(/A)') 'ENTER DIFFUSIVITY FACTOR.'
2        READ(IU_STDIN, *, IOSTAT=IOS) DIFFUSIVITY_FACTOR
         IF (IOS.NE.0) THEN
            WRITE(IU_ERR, '(A)') '+++ INCORRECT VALUE.'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE-ENTER'
               GOTO 2
            ENDIF
         ENDIF
         IF (DIFFUSIVITY_FACTOR.LE.0.0D+00) THEN
            WRITE(IU_ERR, '(A)') '+++ INADMISSIBLE VALUE.'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE-ENTER'
               GOTO 2
            ENDIF
         ENDIF
      ELSE IF ( (CHAR_YN.EQ.'N').OR.(CHAR_YN.EQ.'n') ) THEN
         TRANS_TYPE='PATH'
      ELSE
         WRITE(IU_ERR, '(A)') '+++ UNRECOGNISED RESPONSE.'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_ERR, '(A)') 'PLEASE RE-ENTER.'
            GOTO 1
         ENDIF
      ENDIF
!
!     THE AVERAGED TRANSMISSIONS ARE WRITTEN TO A SCRATCH FILE
!     TO REDUCE STORAGE.
      OPEN(UNIT=IU_AVE_SCRATCH, RECL=16, ACCESS='DIRECT'
     &   , FORM='UNFORMATTED', STATUS='SCRATCH')
!
!
!     AVERAGE THE CONTENTS OF EACH ENTRY IF VALID BUILDING UP
!     UNORDERED SETS OF DATA.
      N_SET=0
      DO I=1, N_ENTRY
!
!        CHECK WHETHER THE ENTRY IS NEEDED AND DETERMINE WHICH
!        GASEOUS SPECIES AND BAND IS TO BE CONSIDERED.
         I_SPECIES=I_INDEX_GAS(I_TYPE_ENTRY(I))
         CALL FIND_GAS_BAND(I_SPECIES
     &      , I, FREQ_ENTRY_LOW(I), FREQ_ENTRY_HIGH(I)
     &      , FIRST_BAND, LAST_BAND
     &      , WAVE_LENGTH_SHORT, WAVE_LENGTH_LONG
     &      , N_BAND_ABSORB, INDEX_ABSORB
     &      , N_BAND_LIST, BAND_LIST, I_SPECIES
     &      )
!
!        READ THE ENTRY IF IT IS REQUIRED
         IF (N_BAND_LIST.GT.0) THEN
!           CHECK THE NUMBER OF FREQUENCIES IN THE ARCHIVE.
            IF (N_TRANS_ENTRY(I).GT.NPD_FREQ_ENTRY) THEN
               WRITE(IU_ERR, '(/A30, A24, /11X, A8, I5, A26, I5)')
     &            '*** ERROR: THERE ARE TOO MANY '
     &            , 'FREQUENCIES IN AN ENTRY.'
     &            , 'ENTRY = ', I, ': NUMBER OF FREQUENCIES = '
     &            , N_TRANS_ENTRY(I)
               IERR=I_ERR_FATAL
               STOP
            ENDIF
            IF (I_VERSION.EQ.0) THEN
               CALL READ_ARCHIVE_ENTRY_0(IU_ARCHIVE_IN
     &            , N_TRANS_RECORD(I), I1_RECORD(I)
     &            , N_PER_ARCH_RECORD, N_TRANS_ENTRY(I)
     &            , FREQ_LOW, FREQ_HIGH, TRANS
     &            )
            ELSE IF (I_VERSION.EQ.1) THEN
               CALL READ_ARCHIVE_ENTRY_1(IU_ARCHIVE_IN
     &            , N_TRANS_RECORD(I), I1_RECORD(I)
     &            , N_PER_ARCH_RECORD, N_TRANS_ENTRY(I)
     &            , FREQ_ENTRY_LOW(I), D_FREQ_ENTRY(I)
     &            , FREQ_LOW, FREQ_HIGH, TRANS
     &            )
            ENDIF
!           ALLOCATE THE TRANSMISSIONS 
!           TO THE MID-POINTS OF EACH INTERVAL.
            DO J=1, N_TRANS_ENTRY(I)
               WAVENUMBER(J)=0.5D+00*(FREQ_LOW(J)+FREQ_HIGH(J))
            ENDDO
         ENDIF
!
!        PROCESS EACH BAND IN THE LIST, FORMING A NEW SET OF DATA
!        IF NECESSARY.
         DO J=1, N_BAND_LIST
            CALL LOCATE_ENTRY(N_SET, LENGTH_SET
     &         , P_SET, T_SET, BAND_SET, INDEX_SET
     &         , P_ENTRY(I), T_ENTRY(I), I_SPECIES
     &         , BAND_LIST(J)
     &         , I_SET, I_PLACE
     &         )
            IF (I_PLACE.GT.NPD_AMOUNT) THEN
               WRITE(IU_ERR, '(/A)')
     &            '*** ERROR: A SET WILL BE TOO LONG: INCREASE '
     &            //'NPD_AMOUNT AND RECOMPILE.'
               IERR=I_ERR_FATAL
               STOP
            ENDIF
!
!           CALCULATE THE WEIGHTINGS. HERE WEIGHTING IS BY WAVENUMBER.
            I_BAND=BAND_LIST(J)
            WAVE_NUMBER_LOW=1.0D+00/WAVE_LENGTH_LONG(I_BAND)
            WAVE_NUMBER_HIGH=1.0D+00/WAVE_LENGTH_SHORT(I_BAND)
            DO K=1, N_BAND_EXCLUDE(I_BAND)
               EXCLUDE_SHORTER(K)
     &            =1.0D+00/WAVE_LENGTH_LONG(INDEX_EXCLUDE(K, I_BAND))
               EXCLUDE_LONGER(K)
     &            =1.0D+00/WAVE_LENGTH_SHORT(INDEX_EXCLUDE(K, I_BAND))
            ENDDO
            CALL WEIGHTINGS(IERR
     &         , I_WEIGHT, .TRUE.
     &         , N_SOLAR_POINTS, SOLAR_WAVELENGTH, SOLAR_IRRAD
     &         , T_ENTRY(I)
     &         , WAVE_NUMBER_LOW, WAVE_NUMBER_HIGH
     &         , N_BAND_EXCLUDE(I_BAND)
     &         , EXCLUDE_SHORTER, EXCLUDE_LONGER
     &         , N_TRANS_ENTRY(I), WAVENUMBER
     &         , I_BEGIN, N_INT_WEIGHT, WEIGHT_COEFF
     &         , WORK
     &         )
            IF (IERR.EQ.I_NORMAL) THEN
!
               TRANS_AVE_ENTRY=0.0D+00
               DO K=0, N_INT_WEIGHT
                  TRANS_AVE_ENTRY=TRANS_AVE_ENTRY
     &               +WEIGHT_COEFF(K)*TRANS(I_BEGIN+K)
               ENDDO
!              PUT THE AVERAGE IN A WELL-DEFINED PLACE IN 
!              THE SCRATCH FILE.
               I_REC=NPD_AMOUNT*(I_SET-1)+I_PLACE
               WRITE(IU_AVE_SCRATCH, REC=I_REC)
     &            AMOUNT_ENTRY(I), TRANS_AVE_ENTRY
!
            ENDIF
         ENDDO
      ENDDO
!
!
!     ORDER THE CASES WITHIN THE SETS BY ABSORBER AMOUNT.
      I_FIRST_RECORD=1
      DO I=1, N_SET
!
!        COLLECT ALL ENTRIES BELONGING TO A SET. THIS IS NOT
!        COMPUTATIONALLY EFFICIENT, BUT AVOIDS USING VERY LARGE AMOUNTS
!        OF STORAGE.
         DO J=1, LENGTH_SET(I)
            I_REC=NPD_AMOUNT*(I-1)+J
            READ(IU_AVE_SCRATCH, REC=I_REC) U_SET(J), TRANS_SET(J)
         ENDDO
!
         CALL ORDER_SET(.FALSE., LENGTH_SET(I)
     &      , DENSITY_BROAD, U_SET, TRANS_SET)
         IF (TRANS_TYPE.EQ.'FLUX') THEN
            CALL CONVERT_AMOUNT(IERR, LENGTH_SET(I)
     &         , DIFFUSIVITY_FACTOR, U_SET)
         ENDIF
         CALL WRITE_TRANS(TYPE_ABSORB(INDEX_SET(I))
     &      , TRANS_TYPE
     &      , I, LENGTH_SET(I), BAND_SET(I), INDEX_SET(I)
     &      , P_SET(I), T_SET(I)
     &      , U_SET, TRANS_SET
     &      , WAVE_LENGTH_SHORT(BAND_SET(I))
     &      , WAVE_LENGTH_LONG(BAND_SET(I))
     &      , N_PER_TRANS_RECORD, I_FIRST_RECORD
     &      )
         IF (IERR.NE.I_NORMAL) STOP
      ENDDO
!
      CLOSE(IU_AVE_SCRATCH)
!
!
      STOP
      END
