!+ Subroutine to solve for fluxes within a multiple mixed column.
!
!
! Method:
!   The fluxes in a multiple mixed column are determined by a method
!   which amounts to Gaussian elimination.
!
! 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 SOLVER_MIX_DIRECT_MULT
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/dimensions_fixed.ucf'
      INCLUDE 'incl/cloud_region.pcf'
!
!
!
!     DUMMY ARGUMENTS.
      INTEGER   !, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOPMOST CLOUDY LAYER
     &   , N_REGION
!             NUMBER OF REGIONS WITHIN COLUMN
!
!     PROPERTIES OF LAYERS.
      REAL*8	!, INTENT(IN)
     &     T(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             TRANSMISSION COEFFICIENTS
     &   , R(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             REFLECTION COEFFICIENTS
     &   , S_UP(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             UPWARD SOURCE TERM 
     &   , S_DOWN(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             DOWNWARD SOURCE TERM 
!
!     ENERGY TRANSFER COEFFICIENTS:
      REAL*8	!, INTENT(IN)
     &     U(NPD_PROFILE, NPD_LAYER, NPD_REGION, NPD_REGION)
!             ENERGY TRANSFER COEFFICIENTS FOR UPWARD RADIATION
     &   , V(NPD_PROFILE, NPD_LAYER, NPD_REGION, NPD_REGION)
!             ENERGY TRANSFER COEFFICIENTS FOR DOWNWARD RADIATION
!
!     SURFACE FIELDS:
      REAL*8	!, INTENT(IN)
     &     ALBEDO_SURFACE_DIFF(NPD_PROFILE)
!             ALBEDO OF DIFFUSE RADIATION
     &   , SOURCE_GROUND(NPD_PROFILE, NPD_REGION)
!             SOURCE OF UPWARD FLUX AT GROUND
!
!     FIELDS AT THE TOP OF THE ATMOSPHERE:
      REAL*8	!, INTENT(IN)
     &     FLUX_INC_DOWN(NPD_PROFILE)
!             INCIDENT DOWNWARD FLUX
!
!     CALCULATED FLUX
      REAL*8	!, INTENT(OUT)
     &     FLUX_TOTAL(NPD_PROFILE, 2*NPD_LAYER+2)
!             TOTAL FLUX
!
!
!     LOCAL VARIABLES:
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
     &   , IJ
!             LOOP VARIABLE OVER REGIONS
     &   , IK
!             LOOP VARIABLE OVER REGIONS
     &   , IL
!             LOOP VARIABLE OVER REGIONS
!
      REAL*8
     &     SOURCE_LOWER(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             SOURCE TERMS FROM LOWER LAYERS
     &   , ALBEDO_EFF(NPD_PROFILE, NPD_LAYER, NPD_REGION, NPD_REGION)
!             EFFECTIVE ALBEDOS OF LOWER LAYERS
     &   , BETA(NPD_PROFILE, NPD_LAYER, NPD_REGION, NPD_REGION)
!             COEFFICIENTS OF UPWARD FLUXES
     &   , GAMMA(NPD_PROFILE, NPD_LAYER, NPD_REGION, NPD_REGION)
!             COEFFICIENTS OF DOWNWARD FLUXES
     &   , H(NPD_PROFILE, NPD_LAYER, NPD_REGION)
!             SOURCE TERMS FOR UPWARD FLUXES
     &   , ZETA(NPD_REGION, NPD_REGION)
!             TEMPORARY COEFFICIENTS FOR UPWARD FLUXES
     &   , LAMBDA(NPD_PROFILE)
!             VARIABLES USED IN ELIMINATION
!
!
!
!     INITIALIZE THE EFFECTIVE ALBEDOS AT THE SURFACE.
      DO IJ=1, N_REGION
         DO IK=1, N_REGION
            DO L=1, N_PROFILE
               ALBEDO_EFF(L, N_LAYER+1, IJ, IK)=V(L, N_LAYER, IJ, IK)
     &            *ALBEDO_SURFACE_DIFF(L)
            ENDDO
         ENDDO
         DO L=1, N_PROFILE
            SOURCE_LOWER(L, N_LAYER, IJ)=SOURCE_GROUND(L, IJ)
         ENDDO
      ENDDO
!
!
!     ELIMINATE UPWARDS THROUGH THE CLOUDY REGIONS.
      DO I=N_LAYER, N_CLOUD_TOP, -1
!
!        FORM THE SUB-MATRIX FOR THE UPWARD FLUXES AT THE LEVEL.
         DO IJ=1, N_REGION
            DO IL=1, N_REGION
               DO L=1, N_PROFILE
                  BETA(L, I, IJ, IL)=0.0D+00
                  GAMMA(L, I, IJ, IL)=0.0D+00
               ENDDO
            ENDDO
            DO L=1, N_PROFILE
               H(L, I, IJ)=SOURCE_LOWER(L, I)
               BETA(L, I, IJ, IJ)=1.0D+00
            ENDDO
            DO IL=1, N_REGION
               DO IK=1, N_REGION
                  DO L=1, N_PROFILE
                     BETA(L, I, IJ, IL)=BETA(L, I, IJ, IL)
     &                  -ALBEDO_EFF(L, I+1, IJ, IK)*R(L, I, IK)
     &                  *U(L, I, IK, IL)
                     GAMMA(L, I, IJ, IL)=GAMMA(L, I, IJ, IL)
     &                  +ALBEDO_EFF(L, I+1, IJ, IK)*T(L, I, IK)
     &                  *V(L, I-1, IK, IL)
                  ENDDO
               ENDDO
            ENDDO
            DO IK=1, N_REGION
               DO L=1, N_PROFILE
                  H(L, I, IJ)=H(L, I, IJ)+ALBEDO_EFF(L, I+1, IJ, IK)
     &                *S_UP(L, I, IK)
               ENDDO
            ENDDO
         ENDDO
!
!
!        PERFORM GAUSSIAN ELIMINATION ON THE SUBMATRIX, STORING
!        THE VALUES REQUIRED FOR PERMANENT USE.
         DO IJ=N_REGION, 2, -1
            DO IL=IJ-1, 1, -1
               DO L=1, N_PROFILE
                  LAMBDA(L)=BETA(L, I, IL, IJ)/BETA(L, I, IJ, IJ)
                  H(L, I, IL)=H(L, I, IL)-LAMBDA(L)*H(L, I, IJ)
               ENDDO
               DO IK=1, IJ-1
                  DO L=1, N_PROFILE
                     BETA(L, I, IL, IK)=BETA(L, I, IL, IK)
     &                  -LAMBDA(L)*BETA(L, I, IJ, IK)
                  ENDDO
               ENDDO
               DO IK=1, IJ-1
                  DO L=1, N_PROFILE
                     GAMMA(L, I, IL, IK)=GAMMA(L, I, IL, IK)
     &                  -LAMBDA(L)*GAMMA(L, I, IJ, IK)
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
!
!
!        PROCEED TO FORMULATE THE EFFECTIVE ALBEDOS AT THE NEXT
!        LEVEL UP.
!
!        INITIALIZE USING THE EQUATIONS FOR THE UPWARD FLUXES AT THAT
!        LEVEL.
         DO IJ=1, N_REGION
            DO IK=1, N_REGION
               DO L=1, N_PROFILE
                  ALBEDO_EFF(L, I, IJ, IK)
     &               =R(L, I, IJ)*V(L, I-1, IJ, IK)
                  ZETA(L, IJ, IK)=T(L, I, IJ)*U(L, I, IJ, IK)
                  SOURCE_LOWER(L, I, IJ)=S_UP(L, I, IJ)
               ENDDO
            ENDDO
         ENDDO
!
!        ELIMINATE ALL THE TERMS INVOLVING ZETA USING THE BETAS
!        AVAILABLE FROM THE PREVIOUS ELIMINATION.
         DO IJ=N_REGION, 1, -1
            DO IL=1, N_REGION
               DO L=1, N_PROFILE
                  LAMBDA(L)=ZETA(L, IL, IJ)/BETA(L, I, IJ, IJ)
                  SOURCE_LOWER(L, I, IJ)=SOURCE_LOWER(L, I, IJ)
     &               -LAMBDA(L)*H(L, I, IJ)
                  DO IK=1, IJ
                     ZETA(L, I, IL, IK)=ZETA(L, I, IL, IK)
     &                  -LAMBDA(L)*BETA(L, I, IJ, IK)
                  ENDDO
                  DO IK=1, N_REGION
                     ALBEDO_EFF(L, I, IL, IK)=ALBEDO_EFF(L, I, IL, IK)
     &                  -LAMBDA(L)*GAMMA(L, I, IJ, IK)
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
         
      ENDDO
!
!
!     CONTINUE THROUGH THE REGION FREE OF CLOUDS.
      DO I=N_CLOUD-1, 1, -1
         DO L=1, N_PROFILE
            BETA(L, I, IP_REGION_CLEAR, IP_REGION_CLEAR)
     &         =1.0D+00-ALPHA(L, I+1, IP_REGION_CLEAR
     &         , IP_REGION_CLEAR)*R(L, I, IP_REGION_CLEAR)
            ALBEDO_EFF(L, I, IP_REGION_CLEAR, IP_REGION_CLEAR)
     &         =R(L, I, IP_REGION_CLEAR)+T(L, I, IP_REGION_CLEAR)**2
     &         *ALPHA(L, I+1, IP_REGION_CLEAR, IP_REGION_CLEAR)
     &         /BETA(L, I, IP_REGION_CLEAR, IP_REGION_CLEAR)
            SOURCE_LOWER(L, I, IP_REGION_CLEAR)
     &         =S_UP(L, I, IP_REGION_CLEAR)+ T(L, I, IP_REGION_CLEAR)
     &         /BETA(L, I, IP_REGION_CLEAR, IP_REGION_CLEAR)
         ENDDO
      ENDDO
!
!
!
!     BACKWARD SUBSTITUTION:
!     
!     INITIALIZE:
      DO L=1, N_PROFILE
         FLUX_TOTAL(L, 2)=FLUX_INC_DOWN(L)
         FLUX_TOTAL(L, 1)=ALPHA(L, 1, IP_REGION_CLEAR, IP_REGION_CLEAR)
     &      *FLUX_TOTAL(L, 2)+SOURCE_LOWER(L, 1, IP_REGION_CLEAR)
      ENDDO
!
      DO I=1, N_CLOUD_TOP-1
         DO L=1, N_PROFILE
            FLUX_TOTAL(L, 2*I+1)=(SOURCE_LOWER(L, IP_REGION_CLEAR)
     &         +ALPHA_EFF(L, I+1, IP_REGION_CLEAR, IP_REGION_CLEAR)
     &         *T(L, I, IP_REGION_CLEAR)*FLUX_TOTAL(L, 2*I))
     &         /BETA(L, I, IP_REGION_CLEAR, IP_REGION_CLEAR)
            FLUX_TOTAL(L, 2*I+2)=S_DOWN(L, I, IP_REGION_CLEAR)
     &         +T(L, I, IP_REGION_CLEAR)*FLUX_TOTAL(L, 2*I)
     &         +R(L, I, IP_REGION_CLEAR)*FLUX_TOTAL(L, 2*I+1)
         ENDDO
      ENDDO
!
!     INITIALIZE FOR BACK-SUBSTITUTION IN CLOUDY REGIONS:
      DO IJ=1, N_REGION
         DO L=1, N_PROFILE
            FLUX_UP(L, N_CLOUD_TOP-1, IJ)=0.0D+00
         ENDDO
      ENDDO
      DO L=1, N_PROFILE
         FLUX_DOWN(L, N_CLOUD_TOP-1, IP_REGION_CLEAR)
     &      =FLUX_TOTAL(L, N_CLOUD_TOP-1)
      ENDDO
!
      DO I=N_CLOUD_TOP, N_LAYER
         DO IJ=1, N_REGION
!
!           UPWARD FLUXES:
            DO L=1, N_PROFILE
               FLUX_UP(L, I, IJ)=H(L, I, IJ)
            ENDDO
            DO IL=1, N_REGION
               FLUX_UP(L, I, IJ)=FLUX_UP(L, I, IJ)
     &            +GAMMA(L, I, IJ, IL)*FLUX_DOWN(L, I-1, IL)
            ENDDO
            DO IL=1, IJ-1
               FLUX_UP(L, I, IJ)=FLUX_UP(L, I, IJ)
     &            -BETA(L, I, IJ, IL)*FLUX_UP(L, I-1, IL)
            ENDDO
            DO L=1, N_PROFILE
               FLUX_UP(L, I, IJ)=FLUX_UP(L, I, IJ)/BETA(L, I, IJ, IJ)
            ENDDO
!
!           DOWNWARD FLUXES:
            DO L=1, N_PROFILE
               FLUX_DOWN(L, I, IJ)=SOURCE_DOWN(L, I, IJ)
            ENDDO
            DO IK=1, N_REGION
               FLUX_DOWN(L, I, IJ)=FLUX_DOWN(L, I, IJ)
     &            +T(L, I, IJ)*V(L, I-1, IJ, IK)*FLUX_DOWN(L, I-1, IK)
     &            +R(L, I, IJ)*U(L, I, IJ, IK)*FLUX_UP(L, I, IK)
            ENDDO
         ENDDO
      ENDDO
!
!
!
      RETURN
      END
