!+ Subroutine to calculate coefficients in the two-stream equations.
!
! Method:
!	The basic two-stream coefficients in the differential equations
!	are calculated. These are then used to determine the 
!	transmission and reflection coefficients. Coefficients for
!	determining the solar or infra-red source terms are calculated.
!
! 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 TWO_COEFF(IERR
     &   , N_PROFILE, I_LAYER_FIRST, I_LAYER_LAST
     &   , I_2STREAM, L_IR_SOURCE_QUAD
     &   , ASYMMETRY, OMEGA, TAU
     &   , ISOLIR, SEC_0
     &   , TRANS, REFLECT, TRANS_0
     &   , SOURCE_COEFF
     &   )
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/dimensions_fixed.pcf'
      INCLUDE 'incl/spectral_region.pcf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/error.pcf'
!
!
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , I_LAYER_FIRST
!             FIRST LAYER TO CONSIDER
     &   , I_LAYER_LAST
!             LAST LAYER TO CONSIDER
     &   , ISOLIR
!             SPECTRAL REGION
     &   , I_2STREAM
!             TWO STREAM SCHEME
      LOGICAL	!, INTENT(IN)
     &     L_IR_SOURCE_QUAD
!               USE A QUADRATIC SOURCE FUNCTION
!
!     OPTICAL PROPERTIES OF LAYER:
      REAL*8	!, INTENT(IN)
     &     ASYMMETRY(NPD_PROFILE, NPD_LAYER)
!             ASYMMETRY FACTOR
     &   , OMEGA(NPD_PROFILE, NPD_LAYER)
!             ALBEDO OF SINGLE SCATTERING
     &   , TAU(NPD_PROFILE, NPD_LAYER)
!             OPTICAL DEPTH
!
!     SOLAR BEAM
      REAL*8	!, INTENT(IN)
     &     SEC_0(NPD_PROFILE)
!             SECANT OF ZENITH ANGLE
!
!
!     COEFFICIENTS IN THE TWO-STREAM EQUATIONS:
      REAL*8	!, INTENT(OUT)
     &     TRANS(NPD_PROFILE, NPD_LAYER)
!             DIFFUSE TRANSMISSION COEFFICIENT
     &   , REFLECT(NPD_PROFILE, NPD_LAYER)
!             DIFFUSE REFLECTION COEFFICIENT
     &   , TRANS_0(NPD_PROFILE, NPD_LAYER)
!             DIRECT TRANSMISSION COEFFICIENT
     &   , SOURCE_COEFF(NPD_PROFILE, NPD_LAYER, NPD_SOURCE_COEFF)
!             SOURCE COEFFICIENTS IN TWO-STREAM EQUATIONS
!
!
!     LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
     &   , K
!             LOOP VARIABLE
     &   , N_INDEX
!             NUMBER OF INDICES SATISFYING TEST
     &   , INDEX(NPD_PROFILE)
!             INDICES OF TESTED POINTS
!
!     COEFFICIENTS IN THE TWO-STREAM EQUATIONS:
      REAL*8
     &     LAMBDA(NPD_PROFILE, NPD_LAYER)
!             COEFFICIENTS IN TWO-STREAM EQUATIONS
     &   , SUM(NPD_PROFILE, NPD_LAYER)
!             SUM OF ALPHA_1 AND ALPHA_2
     &   , DIFF(NPD_PROFILE, NPD_LAYER)
!             DIFFERENCE OF ALPHA_1 AND ALPHA_2
     &   , GAMMA_UP(NPD_PROFILE, NPD_LAYER)
!             BASIC SOLAR COEFFICIENT FOR UPWARD RADIATION
     &   , GAMMA_DOWN(NPD_PROFILE, NPD_LAYER)
!             BASIC SOLAR COEFFICIENT FOR DOWNWARD RADIATION
!
      REAL*8
     &     TARGET
!             TARGET TO SEARCH FOR
!
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     WHENFLT, TWO_COEFF_BASIC, SOLAR_COEFFICIENT_BASIC
     &   , TRANS_SOURCE_COEFF
!
!     CRAY DIRECTIVES FOR THE WHOLE ROUTINE:
!     POINTS ARE NOT REPEATED IN THE INDEXING ARRAY, SO IT IS SAFE
!     TO VECTORIZE OVER INDIRECTLY ADDRESSED ARRAYS.
Cfpp$ NODEPCHK R
!
!
!
!     PERTURB THE SINGLE SCATTERING ALBEDO AWAY FROM 1 TO AVOID
!     LATER DIVISION BY 0.
      TARGET=1.0D+00-TOL_DIV
      DO I=I_LAYER_FIRST, I_LAYER_LAST
         CALL WHENFGT(N_PROFILE, OMEGA(1, I), 1, TARGET
     &      , INDEX, N_INDEX)
         DO K=1, N_INDEX
            OMEGA(INDEX(K), I)=TARGET
         ENDDO
      ENDDO
!
!     CALCULATE THE BASIC TWO-STREAM COEFFICIENTS.
      CALL TWO_COEFF_BASIC(IERR
     &   , N_PROFILE, I_LAYER_FIRST, I_LAYER_LAST
     &   , I_2STREAM
     &   , ASYMMETRY, OMEGA
     &   , SUM, DIFF
     &   )
      IF (IERR.NE.I_NORMAL) THEN
         RETURN
      ENDIF
!
!     LAMBDA IS NOW CALCULATED.
      DO I=I_LAYER_FIRST, I_LAYER_LAST
         DO L=1, N_PROFILE
            LAMBDA(L, I)=SQRT(SUM(L, I)*DIFF(L, I))
         ENDDO
      ENDDO
!
!
!     CALCULATE THE BASIC COEFFICIENTS FOR THE SOLAR SOURCE TERMS.
      IF (ISOLIR.EQ.IP_SOLAR) THEN
!        LAMBDA MAY BE PERTURBED BY THIS ROUTINE TO AVOID
!        ILL-CONDITIONING FOR THE SINGULAR ZENITH ANGLE.
         CALL SOLAR_COEFFICIENT_BASIC(IERR
     &      , N_PROFILE, I_LAYER_FIRST, I_LAYER_LAST
     &      , OMEGA, ASYMMETRY, SEC_0
     &      , I_2STREAM
     &      , SUM, DIFF, LAMBDA
     &      , GAMMA_UP, GAMMA_DOWN
     &      )
         IF (IERR.NE.I_NORMAL) RETURN
      ENDIF
!
!
!     DETERMINE THE TRANSMISSION AND REFLECTION COEFFICIENTS.
      CALL TRANS_SOURCE_COEFF(N_PROFILE, I_LAYER_FIRST, I_LAYER_LAST
     &   , ISOLIR, L_IR_SOURCE_QUAD
     &   , TAU, SUM, DIFF, LAMBDA, SEC_0
     &   , GAMMA_UP, GAMMA_DOWN
     &   , TRANS, REFLECT, TRANS_0, SOURCE_COEFF
     &   )
!
!
!
      RETURN
      END
