!+ Subroutine to calculate two-stream coefficients in the regions.
!
! Method:
!	The coeffients for each region are determined and
!	averaged.
!
! 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_REGION(IERR
     &   , N_PROFILE, N_LAYER, N_CLOUD_TOP
     &   , I_2STREAM, L_IR_SOURCE_QUAD, N_SOURCE_COEFF
     &   , N_CLOUD_TYPE, FRAC_CLOUD
     &   , I_REGION_CLOUD, FRAC_REGION
     &   , ASYMMETRY_FREE, OMEGA_FREE, TAU_FREE
     &   , ASYMMETRY_CLOUD, OMEGA_CLOUD, TAU_CLOUD
     &   , 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/error.pcf'
      INCLUDE 'incl/cloud_region.pcf'
!
!
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOPMOST CLOUDY LAYER 
     &   , ISOLIR
!             SPECTRAL REGION
     &   , N_CLOUD_TYPE
!             NUMBER OF TYPES OF CLOUDS
     &   , I_2STREAM
!             TWO STREAM SCHEME
     &   , N_SOURCE_COEFF
!             NUMBER OF SOURCE COEFFICIENTS
!
      INTEGER	!, INTENT(IN)
     &     I_REGION_CLOUD(NPD_CLOUD_TYPE)
!             REGIONS IN WHICH TYPES OF CLOUDS FALL
!
      LOGICAL	!, INTENT(IN)
     &     L_IR_SOURCE_QUAD
!             USE A QUADRATIC SOURCE IN THE INFRA-RED
!
!     OPTICAL PROPERTIES OF LAYER:
      REAL*8	!, INTENT(IN)
     &     FRAC_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             FRACTIONS OF DIFFERENT TYPES OF CLOUDS
     &   , FRAC_REGION(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             FRACTIONS OF TOTAL CLOUD OCCUPIED BY EACH REGION
     &   , ASYMMETRY_FREE(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLEAR-SKY ASYMMETRY FACTOR
     &   , OMEGA_FREE(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLEAR-SKY ALBEDO OF SINGLE SCATTERING
     &   , TAU_FREE(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLEAR-SKY OPTICAL DEPTH
     &   , ASYMMETRY_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             ASYMMETRY FACTOR
     &   , OMEGA_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             ALBEDO OF SINGLE SCATTERING
     &   , TAU_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             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, NPD_REGION)
!             DIFFUSE TRANSMISSION COEFFICIENT
     &   , REFLECT(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             DIFFUSE REFLECTION COEFFICIENT
     &   , TRANS_0(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             DIRECT TRANSMISSION COEFFICIENT
     &   , SOURCE_COEFF(NPD_PROFILE, NPD_LAYER
     &      , NPD_SOURCE_COEFF, NPD_REGION)
!             SOURCE COEFFICIENTS IN TWO-STREAM EQUATIONS
!
!     LOCAL VARIABLES.
      INTEGER
     &     N_REGION
!             NUMBER OF REGIONS
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , K
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
     &   , I_REGION
!             LOOP VARIABLE OVER REGIONS
!
!     COEFFICIENTS IN THE TWO-STREAM EQUATIONS:
      REAL*8	!, INTENT(OUT)
     &     TRANS_TEMP(NPD_PROFILE, NPD_LAYER)
!             TEMPORARY DIFFUSE TRANSMISSION COEFFICIENT
     &   , REFLECT_TEMP(NPD_PROFILE, NPD_LAYER)
!             TEMPORARY DIFFUSE REFLECTION COEFFICIENT
     &   , TRANS_0_TEMP(NPD_PROFILE, NPD_LAYER)
!             TEMPORARY DIRECT TRANSMISSION COEFFICIENT
     &   , SOURCE_COEFF_TEMP(NPD_PROFILE, NPD_LAYER, NPD_SOURCE_COEFF)
!             TEMPORARY SOURCE COEFFICIENTS IN TWO-STREAM EQUATIONS
!
!     VARIABLES FOR GATHERING:
      INTEGER
     &     N_LIST
!             NUMBER OF POINTS IN LIST
     &   , L_LIST(NPD_PROFILE)
!             LIST OF COLLECTED POINTS
     &   , LL
!             LOOP VARIABLE
      REAL*8
     &     TARGET
!             TARGET FOR SEARCHING
      REAL*8
     &     TAU_GATHERED(NPD_PROFILE, NPD_LAYER)
!             GATHERED OPTICAL DEPTH
     &   , OMEGA_GATHERED(NPD_PROFILE, NPD_LAYER)
!             GATHERED ALEBDO OF SINGLE SCATTERING
     &   , ASYMMETRY_GATHERED(NPD_PROFILE, NPD_LAYER)
!             GATHERED ASYMMETRY
     &   , SEC_0_GATHERED(NPD_PROFILE)
!             GATHERED ASYMMETRY
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     TWO_COEFF, WHENFGT
!
!     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
!
!
!
!     FOR THE TRIPLE OVERLAP THE NUMBER OF REGIONS IS 3.
      N_REGION=3
!
!     DETERMINE THE OPTICAL PROPERTIES OF THE CLEAR-SKY REGIONS OF
!     THE LAYERS.
!
      CALL TWO_COEFF(IERR
     &   , N_PROFILE, 1, N_LAYER
     &   , I_2STREAM, L_IR_SOURCE_QUAD
     &   , ASYMMETRY_FREE, OMEGA_FREE, TAU_FREE
     &   , ISOLIR, SEC_0
     &   , TRANS(1, 1, IP_REGION_CLEAR)
     &   , REFLECT(1, 1, IP_REGION_CLEAR)
     &   , TRANS_0(1, 1, IP_REGION_CLEAR)
     &   , SOURCE_COEFF(1, 1, 1, IP_REGION_CLEAR)
     &   )
      IF (IERR.NE.I_NORMAL) RETURN
!
!
!     NOW DEAL WITH CLOUDS.
!
!     INITIALIZE THE FULL ARRAYS FOR CLOUDY REGIONS.
!
      DO I_REGION=1, N_REGION
         IF (I_REGION.NE.IP_REGION_CLEAR) THEN
            DO I=N_CLOUD_TOP, N_LAYER
               DO L=1, N_PROFILE
                  TRANS(L, I, I_REGION)=0.0D+00
                  REFLECT(L, I, I_REGION)=0.0D+00
               ENDDO
            ENDDO
            DO J=1, N_SOURCE_COEFF
               DO I=N_CLOUD_TOP, N_LAYER
                  DO L=1, N_PROFILE
                     SOURCE_COEFF(L, I, J, I_REGION)=0.0D+00
                  ENDDO
               ENDDO
            ENDDO
!
            IF (ISOLIR.EQ.IP_SOLAR) THEN
               DO I=N_CLOUD_TOP, N_LAYER
                  DO L=1, N_PROFILE
                     TRANS_0(L, I, I_REGION)=0.0D+00
                  ENDDO
               ENDDO
            ENDIF
!
         ENDIF
!
      ENDDO
!
!
!
!     CONSIDER EACH TYPE OF CLOUD IN TURN, CHECKING WHICH REGION IT
!     CONTRUBUTES TO AND FORM WEIGHTED SUMS OF CLOUD PROPERTIES.
!
      DO K=1, N_CLOUD_TYPE
!
!
!        SET THE REGION IN WHICH CLOUDS OF THIS TYPE ARE INCLUDED.
         I_REGION=I_REGION_CLOUD(K)
!
!        GATHER POINTS WHERE THERE IS CLOUD OF THE PRESENT TYPE.
         TARGET=0.0D+00
!
         DO I=N_CLOUD_TOP, N_LAYER
!
!           DETERMINE WHETHER THIS TYPE OF CLOUD EXISTS IN THIS ROW.
            CALL WHENFGT(N_PROFILE, FRAC_CLOUD(1, I, K), 1
     &         , TARGET, L_LIST, N_LIST)
!
!
            IF (N_LIST.GT.0) THEN
!
!              GATHER THE OPTICAL PROPERTIES. THOUGH WE CONSIDER ONLY
!              ONE LAYER AT A TIME THE LOWER ROUTINES WILL OPERATE ON
!              ARRAYS WITH VERTICAL STRUCTURE, SO THE GATHERED ARRAYS
!              ARE TWO-DIMENSIONAL.
!
               DO L=1, N_LIST
                  TAU_GATHERED(L, I)
     &              =TAU_CLOUD(L_LIST(L), I, K)
                  OMEGA_GATHERED(L, I)
     &              =OMEGA_CLOUD(L_LIST(L), I, K)
                  ASYMMETRY_GATHERED(L, I)
     &              =ASYMMETRY_CLOUD(L_LIST(L), I, K)
               ENDDO
               IF (ISOLIR.EQ.IP_SOLAR) THEN
                  DO L=1, N_LIST
                     SEC_0_GATHERED(L)=SEC_0(L_LIST(L))
                  ENDDO
               ENDIF
!
!
               CALL TWO_COEFF(IERR
     &            , N_LIST, I, I
     &            , I_2STREAM, L_IR_SOURCE_QUAD
     &            , ASYMMETRY_GATHERED, OMEGA_GATHERED
     &            , TAU_GATHERED
     &            , ISOLIR, SEC_0_GATHERED
     &            , TRANS_TEMP, REFLECT_TEMP, TRANS_0_TEMP
     &            , SOURCE_COEFF_TEMP
     &            )
               IF (IERR.NE.I_NORMAL) RETURN
!
!
               DO L=1, N_LIST
                  LL=L_LIST(L)
                  TRANS(LL, I, I_REGION)=TRANS(LL, I, I_REGION)
     &               +FRAC_CLOUD(LL, I, K)*TRANS_TEMP(L, I)
                  REFLECT(LL, I, I_REGION)=REFLECT(LL, I, I_REGION)
     &               +FRAC_CLOUD(LL, I, K)*REFLECT_TEMP(L, I)
               ENDDO
               DO J=1, N_SOURCE_COEFF
                  DO L=1, N_LIST
                     LL=L_LIST(L)
                     SOURCE_COEFF(LL, I, J, I_REGION)
     &                  =SOURCE_COEFF(LL, I, J, I_REGION)
     &                  +FRAC_CLOUD(LL, I, K)
     &                  *SOURCE_COEFF_TEMP(L, I, J)
                  ENDDO
               ENDDO
               IF (ISOLIR.EQ.IP_SOLAR) THEN
                  DO L=1, N_LIST
                     LL=L_LIST(L)
                     TRANS_0(LL, I, I_REGION)=TRANS_0(LL, I, I_REGION)
     &                  +FRAC_CLOUD(LL, I, K)*TRANS_0_TEMP(L, I)
                  ENDDO
               ENDIF
!
            ENDIF
!
         ENDDO
      ENDDO
!
!     FINALLY, SCALE THE WEIGHTED SUMS BY THE CLOUD FRACTIONS.
      DO I_REGION=1, N_REGION
         IF (I_REGION.NE.IP_REGION_CLEAR) THEN
            DO I=N_CLOUD_TOP, N_LAYER
               CALL WHENFGT(N_PROFILE, FRAC_REGION(1, I, I_REGION), 1
     &            , TARGET, L_LIST, N_LIST)
               DO L=1, N_LIST
                  LL=L_LIST(L)
                  TRANS(LL, I, I_REGION)=TRANS(LL, I, I_REGION)
     &               /FRAC_REGION(LL, I, I_REGION)
                  REFLECT(LL, I, I_REGION)=REFLECT(LL, I, I_REGION)
     &               /FRAC_REGION(LL, I, I_REGION)
               ENDDO
               DO J=1, N_SOURCE_COEFF
                  DO L=1, N_LIST
                     LL=L_LIST(L)
                     SOURCE_COEFF(LL, I, J, I_REGION)
     &                  =SOURCE_COEFF(LL, I, J, I_REGION)
     &                  /FRAC_REGION(LL, I, I_REGION)
                  ENDDO
               ENDDO
               IF (ISOLIR.EQ.IP_SOLAR) THEN
                  DO L=1, N_LIST
                     LL=L_LIST(L)
                     TRANS_0(LL, I, I_REGION)=TRANS_0(LL, I, I_REGION)
     &                  /FRAC_REGION(LL, I, I_REGION)
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDDO
!
!
!
      RETURN
      END
