!+ Subroutine to calculate the solar beam and P.I. in Schmetz's scheme.
!
! Method:
!	The solar beam is propagated down through the amosphere and
!	the terms of the particular integral are determined using
!	standard formulae.
!
! 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 S_DIRECT_SOURCE(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &   , TAU_FREE, TAU_CLOUD
     &   , GAMMA_1_FREE, GAMMA_2_FREE
     &   , TRANS_FREE, REFLECT_FREE
     &   , GAMMA_1_CLOUD, GAMMA_2_CLOUD
     &   , TRANS_CLOUD, REFLECT_CLOUD
     &   , MU_0
     &   , W_CLOUD
     &   , L_SCALE_SOLAR, SCALE_SOLAR_FREE, SCALE_SOLAR_CLOUD
     &   , FLUX_DIRECT, FLUX_INC_DIRECT
     &   , S_UP_FREE, S_DOWN_FREE, S_UP_CLOUD, S_DOWN_CLOUD
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES
      INCLUDE 'incl/dimensions_field.ucf'
!
!     DECLARATION OF DUMMY ARGUMENTS.
      INTEGER 	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOP CLOUDY LAYER
      LOGICAL	!, INTENT(IN)
     &     L_SCALE_SOLAR
!             APPLY SCALING TO SOLAR FLUX
      REAL*8 	!, INTENT(IN)
     &     TAU_FREE(NPD_PROFILE, NPD_LAYER)
!             OPTICAL DEPTH
     &   , TAU_CLOUD(NPD_PROFILE, NPD_LAYER)
!             OPTICAL DEPTH
     &   , TRANS_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE TRANSMISSION COEFFICIENT
     &   , REFLECT_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE REFLECTION COEFFICIENT
     &   , TRANS_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY TRANSMISSION COEFFICIENT
     &   , REFLECT_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY REFLECTION COEFFICIENT
     &   , GAMMA_1_FREE(NPD_PROFILE, NPD_LAYER)
!             COEFFICIENTS IN SOLAR SOURCE
     &   , GAMMA_2_FREE(NPD_PROFILE, NPD_LAYER)
!
     &   , GAMMA_1_CLOUD(NPD_PROFILE, NPD_LAYER)
!
     &   , GAMMA_2_CLOUD(NPD_PROFILE, NPD_LAYER)
!
     &   , W_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUD AMOUNT
     &   , MU_0(NPD_PROFILE)
!             COSINES OF ZENITH ANGLES
     &   , FLUX_INC_DIRECT(NPD_PROFILE)
!             INCIDENT DIRECT FLUXES
     &   , SCALE_SOLAR_FREE(NPD_PROFILE, NPD_LAYER)
!             SCALING OF FREE SOLAR BEAM
     &   , SCALE_SOLAR_CLOUD(NPD_PROFILE, NPD_LAYER)
!             SCALING OF CLOUDY SOLAR BEAM
      REAL*8 	!, INTENT(OUT)
     &     FLUX_DIRECT(NPD_PROFILE, 0: NPD_LAYER)
!             DIRECT FLUXES
     &   , S_UP_FREE(NPD_PROFILE, NPD_LAYER)
!             UPWARD SOURCE FUNCTION
     &   , S_DOWN_FREE(NPD_PROFILE, NPD_LAYER)
!             DOWNWARD SOURCE FUNCTION
     &   , S_UP_CLOUD(NPD_PROFILE, NPD_LAYER)
!             UPWARD SOURCE FUNCTION
     &   , S_DOWN_CLOUD(NPD_PROFILE, NPD_LAYER)
!             DOWNWARD SOURCE FUNCTION
!
!
!     LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
      REAL*8
     &     FLUX_DIRECT_TOP(NPD_PROFILE)
!             DIRECT FLUX AT TOP OF LAYER
     &   , FLUX_DIRECT_CLOUD(NPD_PROFILE)
!             DIRECT FLUX AT BASE OF CLOUD
!
!
!     CALCULATE THE INTENSITY OF THE DIRECT BEAM JUST ABOVE EACH
!     INTERFACE.
!
!     SURFACE CONDITION:
      DO L=1, N_PROFILE
         FLUX_DIRECT(L, 0)=FLUX_INC_DIRECT(L)
      ENDDO
!
!     CLEAR-SKY REGION ABOVE THE CLOUDS.
      DO I=1, N_CLOUD_TOP-1
         IF (L_SCALE_SOLAR) THEN
            DO L=1, N_PROFILE
               FLUX_DIRECT(L, I)
     &            =FLUX_DIRECT(L, I-1)*EXP(-TAU_FREE(L, I)
     &            *SCALE_SOLAR_FREE(L, I)/MU_0(L))
            ENDDO
         ELSE
            DO L=1, N_PROFILE
               FLUX_DIRECT(L, I)
     &            =FLUX_DIRECT(L, I-1)*EXP(-TAU_FREE(L, I)/MU_0(L))
            ENDDO
         ENDIF
         DO L=1, N_PROFILE
            S_UP_FREE(L, I)=(GAMMA_1_FREE(L, I)-GAMMA_2_FREE(L, I)
     &         *REFLECT_FREE(L, I))*FLUX_DIRECT(L, I-1)
     &         -GAMMA_1_FREE(L, I)*TRANS_FREE(L, I)*FLUX_DIRECT(L, I)
            S_DOWN_FREE(L, I)=(GAMMA_2_FREE(L, I)-GAMMA_1_FREE(L, I)
     &         *REFLECT_FREE(L, I))*FLUX_DIRECT(L, I)
     &         -GAMMA_2_FREE(L, I)*TRANS_FREE(L, I)*FLUX_DIRECT(L, I-1)
         ENDDO
      ENDDO
!
!     CLEAR AND CLOUDY REGION. FLUX_DIRECT_CLOUD IS USED TO HOLD THE
!     FLUX TRANSMITTED THROUGH THE CLOUD IGNORING ITS AERIAL COVERAGE.
      DO I=N_CLOUD_TOP, N_LAYER
         DO L=1, N_PROFILE
            IF (L_SCALE_SOLAR) THEN
               FLUX_DIRECT_CLOUD(L)=FLUX_DIRECT(L, I-1)
     &            *EXP(-TAU_CLOUD(L, I)*SCALE_SOLAR_CLOUD(L, I)
     &            /MU_0(L))
               FLUX_DIRECT_TOP(L)=(1.0D+00-W_CLOUD(L, I))
     &            *FLUX_DIRECT(L, I-1)+W_CLOUD(L, I)
     &            *FLUX_DIRECT_CLOUD(L)
               FLUX_DIRECT(L, I)=FLUX_DIRECT_TOP(L)
     &            *EXP(-TAU_FREE(L, I)*SCALE_SOLAR_FREE(L, I)/MU_0(L))
            ELSE
               FLUX_DIRECT_CLOUD(L)=FLUX_DIRECT(L, I-1)
     &            *EXP(-TAU_CLOUD(L, I)/MU_0(L))
               FLUX_DIRECT_TOP(L)=(1.0D+00-W_CLOUD(L, I))
     &            *FLUX_DIRECT(L, I-1)+W_CLOUD(L, I)
     &            *FLUX_DIRECT_CLOUD(L)
               FLUX_DIRECT(L, I)=FLUX_DIRECT_TOP(L)
     &            *EXP(-TAU_FREE(L, I)/MU_0(L))
            ENDIF
            S_UP_FREE(L, I)=(GAMMA_1_FREE(L, I)-GAMMA_2_FREE(L, I)
     &         *REFLECT_FREE(L, I))*FLUX_DIRECT_TOP(L)
     &         -GAMMA_1_FREE(L, I)*TRANS_FREE(L, I)*FLUX_DIRECT(L, I)
            S_DOWN_FREE(L, I)=(GAMMA_2_FREE(L, I)-GAMMA_1_FREE(L, I)
     &         *REFLECT_FREE(L, I))*FLUX_DIRECT(L, I)
     &         -GAMMA_2_FREE(L, I)*TRANS_FREE(L, I)*FLUX_DIRECT_TOP(L)
            S_UP_CLOUD(L, I)=(GAMMA_1_CLOUD(L, I)-GAMMA_2_CLOUD(L, I)
     &         *REFLECT_CLOUD(L, I))*FLUX_DIRECT(L, I-1)
     &         -GAMMA_1_CLOUD(L, I)*TRANS_CLOUD(L, I)
     &         *FLUX_DIRECT_CLOUD(L)
            S_DOWN_CLOUD(L, I)=(GAMMA_2_CLOUD(L, I)-GAMMA_1_CLOUD(L, I)
     &         *REFLECT_CLOUD(L, I))*FLUX_DIRECT_CLOUD(L)
     &         -GAMMA_2_CLOUD(L, I)*TRANS_CLOUD(L, I)
     &         *FLUX_DIRECT(L, I-1)
         ENDDO
      ENDDO
!
!
      RETURN
      END
