!+ Subroutine to solve the two-stream equations in a mixed column.
!
! Method:
!	The two-stream coefficients are calculated in clear regions
!	and in stratiform and convective clouds. From these
!	coefficients transmission and reflection coefficients are
!	determined. The coefficients for convective and stratiform
!	clouds are appropriately mixed to form single cloudy values
!	and an appropriate solver is called.
!
! 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 MIX_COLUMN(IERR
!			Atmospheric Properties
     &   , N_PROFILE, N_LAYER
!			Two-stream Scheme
     &   , I_2STREAM
!			Corrections to Two-stream Equations
     &   , L_2_STREAM_CORRECT, PLANCK_SOURCE, GROUND_EMISSION
!			Options for Solver
     &   , I_SOLVER, I_EXTRAPOLATE_MIX_NET, L_NET
!			Options for Equivalent Extinction
     &   , L_SCALE_SOLAR, ADJUST_SOLAR_KE
!			Spectral Region
     &   , ISOLIR
!			Infra-red Properties
     &   , DIFF_PLANCK
     &   , L_IR_SOURCE_QUAD, DIFF_PLANCK_2
!			Conditions at TOA
     &   , FLUX_INC_DOWN, FLUX_INC_DIRECT, SEC_0
!			Conditions at Surface
     &   , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
!			Clear-sky Single Scattering Properties
     &   , TAU_FREE, OMEGA_FREE, ASYMMETRY_FREE
!			Cloud Geometry
     &   , N_CLOUD_TOP
     &   , N_CLOUD_TYPE, FRAC_CLOUD
     &   , W_FREE, N_FREE_PROFILE, I_FREE_PROFILE
     &   , W_CLOUD, N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &   , CLOUD_OVERLAP
!			Cloudy Optical Properties
     &   , TAU_CLOUD, OMEGA_CLOUD, ASYMMETRY_CLOUD
!			Fluxes Calculated
     &   , FLUX_DIRECT, FLUX_TOTAL
!			Flags for Clear-sky Calculations
     &   , L_CLEAR, I_SOLVER_CLEAR
!			Clear-sky Fluxes Calculated
     &   , FLUX_DIRECT_CLEAR, FLUX_TOTAL_CLEAR
     &   )
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/dimensions_fixed.pcf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/solver.pcf'
      INCLUDE 'incl/spectral_region.pcf'
      INCLUDE 'incl/cloud_coeff_pointer.pcf'
!
!     DUMMY VARIABLES.
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOP CLOUDY LAYER
     &   , N_CLOUD_TYPE
!             NUMBER OF TYPES OF CLOUDS
     &   , N_FREE_PROFILE(NPD_LAYER)
!             NUMBER OF FREE PROFILES
     &   , I_FREE_PROFILE(NPD_PROFILE, NPD_LAYER)
!             INDICES OF FREE PROFILES
     &   , N_CLOUD_PROFILE(NPD_LAYER)
!             NUMBER OF CLOUDY PROFILES
     &   , I_CLOUD_PROFILE(NPD_PROFILE, NPD_LAYER)
!             INDICES OF CLOUDY PROFILES
     &   , ISOLIR
!             SPECTRAL REGION
     &   , I_2STREAM
!             TWO-STREAM SCHEME
     &   , I_SOLVER
!             SOLVER USED
     &   , I_EXTRAPOLATE_MIX_NET
!             METHOD FOR EXTRAPOLATING WITH A MIXED SOLVER
     &   , I_SOLVER_CLEAR
!             SOLVER FOR CLEAR-SKY FLUXES
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      LOGICAL	!, INTENT(IN)
     &     L_NET
!             CALCULATE NET FLUXES
     &   , L_CLEAR
!             CALCULATE CLEAR-SKY FLUXES
     &   , L_SCALE_SOLAR
!             FLAG TO SCALE SOLAR
     &   , L_IR_SOURCE_QUAD
!             USE QUADRATIC SOURCE TERM
     &   , L_2_STREAM_CORRECT
!             EDGE CORRECTION TO 2-STREAM
!
!     OPTICAL PROPERTIES:
      REAL*8	!, INTENT(IN)
     &     TAU_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE OPTICAL DEPTH
     &   , OMEGA_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE ALBEDO OF SINGLE SCATTERING
     &   , ASYMMETRY_FREE(NPD_PROFILE, NPD_LAYER)
!             CLEAR-SKY ASYMMETRY
     &   , TAU_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLOUDY OPTICAL DEPTH
     &   , OMEGA_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLOUDY ALBEDO OF SINGLE SCATTERING
     &   , ASYMMETRY_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             CLOUDY ASYMMETRY 
!
!     CLOUD GEOMETRY:
      REAL*8	!, INTENT(IN)
     &     W_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY FRACTIONS IN EACH LAYER
     &   , W_FREE(NPD_PROFILE, NPD_LAYER)
!             CLEAR SKY FRACTIONS IN EACH LAYER
     &   , FRAC_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             FRACTIONS OF DIFFERENT TYPES OF CLOUD
     &   , CLOUD_OVERLAP(NPD_PROFILE, 0: NPD_LAYER, NPD_OVERLAP_COEFF)
!             ENERGY TRANSFER COEFFICIENTS
      REAL*8	!, INTENT(IN)
     &     SEC_0(NPD_PROFILE)
!             SECANT OF SOLAR ZENITH ANGLE
     &   , ALBEDO_SURFACE_DIFF(NPD_PROFILE)
!             DIFFUSE ALBEDO
     &   , ALBEDO_SURFACE_DIR(NPD_PROFILE)
!             DIRECT ALBEDO
     &   , FLUX_INC_DOWN(NPD_PROFILE)
!             INCIDENT TOTAL FLUX
     &   , FLUX_INC_DIRECT(NPD_PROFILE)
!             INCIDENT DIRECT FLUX
     &   , DIFF_PLANCK(NPD_PROFILE, NPD_LAYER)
!             CHANGE IN PLANCK FUNCTION
     &   , SOURCE_GROUND(NPD_PROFILE)
!             FLUX FROM SURFACE
     &   , ADJUST_SOLAR_KE(NPD_PROFILE, NPD_LAYER)
!             ADJUSTMENT OF SOLAR BEAM WITH EQUIVALENT EXTINCTION
     &   , DIFF_PLANCK_2(NPD_PROFILE, NPD_LAYER)
!             2x2ND DIFFERENCE OF PLANCKIAN
     &   , PLANCK_SOURCE(NPD_PROFILE, 0: NPD_LAYER)
!             PLANCKIAN SOURCE FUNCTION
     &   , GROUND_EMISSION(NPD_PROFILE)
!             TOTAL FLUX EMITTED FROM GROUND
!
!     FLUXES CALCULATED
      REAL*8	!, INTENT(OUT)
     &     FLUX_DIRECT(NPD_PROFILE, 0: NPD_LAYER)
!             DIRECT FLUX
     &   , FLUX_TOTAL(NPD_PROFILE, 2*NPD_LAYER+2)
!             LONG FLUX VECTOR
     &   , FLUX_DIRECT_CLEAR(NPD_PROFILE, 0: NPD_LAYER)
!             CLEAR DIRECT FLUX
     &   , FLUX_TOTAL_CLEAR(NPD_PROFILE, 2*NPD_LAYER+2)
!             CLEAR TOTAL FLUX
!
!
!
!     LOCAL VARIABALES.
      INTEGER
     &     N_SOURCE_COEFF
!             NUMBER OF SOURCE COEFFICIENTS
     &   , N_EQUATION
!             NUMBER OF EQUATIONS
     &   , I
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
!
!
!     CLEAR-SKY COEFFICIENTS:
      REAL*8
     &     TRANS_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE TRANSMISSION OF LAYER
     &   , REFLECT_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE REFLECTANCE OF LAYER
     &   , TRANS_0_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE DIRECT TRANSMISSION OF LAYER
     &   , SOURCE_COEFF_FREE(NPD_PROFILE, NPD_LAYER, NPD_SOURCE_COEFF)
!             FREE SOURCE COEFFICIENTS
     &   , S_DOWN_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE DOWNWARD SOURCE
     &   , S_UP_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE UPWARD SOURCE
     &   , S_DOWN_CLEAR(NPD_PROFILE, NPD_LAYER)
!             CLEAR DOWNWARD SOURCE
     &   , S_UP_CLEAR(NPD_PROFILE, NPD_LAYER)
!             CLEAR UPWARD SOURCE
!
!     CLOUDY COEFFICIENTS:
      REAL*8
     &     TRANS_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY TRANSMISSION OF LAYER
     &   , REFLECT_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY REFLECTANCE OF LAYER
     &   , TRANS_0_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY DIRECT TRANSMISSION OF LAYER
     &   , SOURCE_COEFF_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_SOURCE_COEFF)
!             CLOUDY SOURCE COEFFICIENTS
     &   , S_DOWN_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY DOWNWARD SOURCE
     &   , S_UP_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY UPWARD SOURCE
!
!     SOURCE FUNCTIONS AT THE CROUND
      REAL*8
     &     SOURCE_GROUND_FREE(NPD_PROFILE)
!             SOURCE FROM GROUND UNDER CLEAR SKIES
     &   , SOURCE_GROUND_CLOUD(NPD_PROFILE)
!             SOURCE FROM GROUND UNDER CLOUDY SKIES
     &   , FLUX_DIRECT_GROUND_CLOUD(NPD_PROFILE)
!             DIRECT FLUX AT GROUND UNDER CLOUDY SKIES
!
!     VARIABLES FOR DEALING WITH ILL-CONDITIONING:
      LOGICAL
     &     L_SET_MIX_EXTRAPOLATE
!             FUNCTION FOR CHOOSING METHOD OF EXTRAPOLATION
      REAL*8
     &     TOL_CLOUD
!             TOLERANCE FOR RESTORING CONDITIONING
     &   , FLUX_TOTAL_2(NPD_PROFILE, 2*NPD_LAYER+2)
!             SECOND TOTAL FLUX ARRAY
!
!
!     NUMERICAL ARRAYS:
      REAL*8
     &     A5(NPD_PROFILE, 5, 2*NPD_LAYER+2)
!             PENTADIAGONAL MATRIX
     &   , B(NPD_PROFILE, 2*NPD_LAYER+2)
!             RHS OF MATRIX EQUATION
     &   , WORK(NPD_PROFILE)
!             WORKING ARRAY FOR SOLVER
!
!     FUNCTIONS CALLED:
      INTEGER
     &     SET_N_SOURCE_COEFF
!             FUNCTION TO SET NUMBER OF SOURCE COEFFICIENTS
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     TWO_COEFF, TWO_COEFF_CLOUD, IR_SOURCE, MIXED_SOLAR_SOURCE
     &   , L_SET_MIX_EXTRAPOLATE, MIX_NET_MATRIX, BAND_SOLVER
     &   , MIX_SUM_MATRIX, MIX_COLUMN_FULL, MIX_APP_SCAT
     &   , CLEAR_SUPPLEMENT
!
!
!
!     CALCULATE THE TRANSMISSION AND REFLECTION COEFFICIENTS AND
!     SOURCE TERMS FOR THE CLEAR AND CLOUDY PARTS OF THE COLUMN
!
!     SET THE NUMBER OF SOURCE COEFFICIENTS FOR THE APPROXIMATION
      N_SOURCE_COEFF=SET_N_SOURCE_COEFF(ISOLIR, L_IR_SOURCE_QUAD)
!
      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_FREE, REFLECT_FREE, TRANS_0_FREE
     &   , SOURCE_COEFF_FREE
     &   )
      IF (IERR.NE.I_NORMAL) RETURN
!
!
!     INFRA-RED SOURCE TERMS DEPEND ONLY ON THE LAYER AND MAY BE
!     CALCULATED NOW. SOLAR TERMS DEPEND ON CONDITIONS IN CLOUD
!     IN OVERLYING LAYERS AND MUST BE CALCULATED LATER.
!
      IF (ISOLIR.EQ.IP_INFRA_RED) THEN
!
         CALL IR_SOURCE(N_PROFILE, 1, N_LAYER
     &      , SOURCE_COEFF_FREE, DIFF_PLANCK
     &      , L_IR_SOURCE_QUAD, DIFF_PLANCK_2
     &      , L_2_STREAM_CORRECT, PLANCK_SOURCE
     &      , GROUND_EMISSION, N_LAYER
     &      , TAU_FREE, TRANS_FREE
     &      , S_DOWN_FREE, S_UP_FREE
     &      )
!
!        IF A CLEAR-SKY CALCULATION IS REQUIRED THESE SOURCE TERMS MUST
!        BE STORED.
         IF (L_CLEAR) THEN
            DO I=1, N_LAYER
               DO L=1, N_PROFILE
                  S_DOWN_CLEAR(L, I)=S_DOWN_FREE(L, I)
                  S_UP_CLEAR(L, I)=S_UP_FREE(L, I)
               ENDDO
            ENDDO
         ENDIF
!
!        SCALE THE SOURCES BY THE CLEAR-SKY FRACTIONS IN THE CLOUDY 
!        LAYERS. IN HIGHER LAYERS THE CLEAR-SKY FRACTION IS 1.
         DO I=N_CLOUD_TOP, N_LAYER
            DO L=1, N_PROFILE
               S_DOWN_FREE(L, I)=W_FREE(L, I)*S_DOWN_FREE(L, I)
               S_UP_FREE(L, I)=W_FREE(L, I)*S_UP_FREE(L, I)
            ENDDO
         ENDDO
!
      ENDIF
!
!
!
!     REPEAT THE CALCULATION FOR CLOUDY REGIONS.
!
!
      CALL TWO_COEFF_CLOUD(IERR
     &   , N_PROFILE, N_CLOUD_TOP, N_LAYER
     &   , I_2STREAM, L_IR_SOURCE_QUAD, N_SOURCE_COEFF
     &   , N_CLOUD_TYPE, FRAC_CLOUD
     &   , ASYMMETRY_CLOUD, OMEGA_CLOUD, TAU_CLOUD
     &   , ISOLIR, SEC_0
     &   , TRANS_CLOUD, REFLECT_CLOUD, TRANS_0_CLOUD
     &   , SOURCE_COEFF_CLOUD
     &   )
      IF (IERR.NE.I_NORMAL) RETURN
!
!
      IF (ISOLIR.EQ.IP_INFRA_RED) THEN
!
!        EDGE CORRECTIONS FOR THE TWO-STREAM EQUATIONS DO NOT
!        REALLY FIT WITH THIS METHOD OF TREATING CLOUDS. OPTICAL
!        DEPTHS AND TRANSMISSIONS MUST BE PASSED TO THE SUBROUTINE
!        TO FILL THE ARGUMENT LIST, BUT IT IS NOT INTENDED THAT
!        THESE ARRAYS WILL BE USED.
!
         CALL IR_SOURCE(N_PROFILE, N_CLOUD_TOP, N_LAYER
     &      , SOURCE_COEFF_CLOUD, DIFF_PLANCK
     &      , L_IR_SOURCE_QUAD, DIFF_PLANCK_2
     &      , L_2_STREAM_CORRECT, PLANCK_SOURCE
     &      , GROUND_EMISSION, N_LAYER
     &      , TAU_CLOUD, TRANS_CLOUD
     &      , S_DOWN_CLOUD, S_UP_CLOUD
     &      )
!
!
         DO I=N_CLOUD_TOP, N_LAYER
            DO L=1, N_PROFILE
               S_DOWN_CLOUD(L, I)=W_CLOUD(L, I)*S_DOWN_CLOUD(L, I)
               S_UP_CLOUD(L, I)=W_CLOUD(L, I)*S_UP_CLOUD(L, I)
            ENDDO
         ENDDO
!
      ENDIF
!
!
!     CALCULATE THE APPROPRIATE SOURCE TERMS FOR THE SOLAR: CLOUDY
!     AND CLEAR PROPERTIES ARE BOTH NEEDED HERE.
!
      IF (ISOLIR.EQ.IP_SOLAR) THEN
!
         CALL MIXED_SOLAR_SOURCE(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , FLUX_INC_DIRECT
     &      , L_SCALE_SOLAR, ADJUST_SOLAR_KE
     &      , TRANS_0_FREE, SOURCE_COEFF_FREE
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCC)
     &      , TRANS_0_CLOUD, SOURCE_COEFF_CLOUD
     &      , FLUX_DIRECT
     &      , FLUX_DIRECT_GROUND_CLOUD
     &      , S_UP_FREE, S_DOWN_FREE
     &      , S_UP_CLOUD, S_DOWN_CLOUD
     &   )
      ENDIF
!
!
!
!     FORMULATE THE MATRIX EQUATION FOR THE FLUXES.
!
      IF (I_SOLVER.EQ.IP_SOLVER_MIX_NET) THEN
!
!        TO AVOID NUMERICAL INSTABILITY WHEN SOLVING FOR NET FLUXES
!        ALONE A SMALL AMOUNT OF ABSORPTION MUST BE INTRODUCED. THIS
!        IS THE ROLE OF TOL_CLOUD: IT IS DEFINED FROM THE TOLERANCE
!        AVAILABLE ON THE MACHINE. FOR SOME APPLICATIONS THIS
!        ABSORPTION MAY BE TOO LARGE AND AN EXTRAPOLATION IS USED TO
!        INFER THE LIMITING VALUE.
!
         TOL_CLOUD=4.0D+01*SQRT_TOL_MACHINE
!
         CALL MIX_NET_MATRIX(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , TRANS_FREE, REFLECT_FREE
     &      , TRANS_CLOUD, REFLECT_CLOUD
     &      , S_UP_FREE, S_DOWN_FREE, S_UP_CLOUD, S_DOWN_CLOUD
     &      , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &      , FLUX_DIRECT(1, N_LAYER)
     &      , FLUX_INC_DOWN
     &      , N_FREE_PROFILE, I_FREE_PROFILE
     &      , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GP)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BP)
     &      , TOL_CLOUD
     &      , A5, B
     &      )
!
         N_EQUATION=N_LAYER+1
!
!        SOLVE THE SET OF EQUATIONS PRODUCED.
         CALL BAND_SOLVER(N_PROFILE, N_EQUATION
     &      , 2, 2
     &      , A5, B
     &      , FLUX_TOTAL
     &      , NPD_PROFILE, 2*NPD_LAYER+2
     &      , WORK
     &      )
!
         IF (L_SET_MIX_EXTRAPOLATE(I_EXTRAPOLATE_MIX_NET)) THEN
!
!           RECALCULATE THE FLUXES WITH THE PERTURBATION DOUBLED.
            TOL_CLOUD=2.0D+00*TOL_CLOUD
!
            CALL MIX_NET_MATRIX(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &         , TRANS_FREE, REFLECT_FREE
     &         , TRANS_CLOUD, REFLECT_CLOUD
     &         , S_UP_FREE, S_DOWN_FREE, S_UP_CLOUD, S_DOWN_CLOUD
     &         , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &         , FLUX_DIRECT(1, N_LAYER)
     &         , FLUX_INC_DOWN
     &         , N_FREE_PROFILE, I_FREE_PROFILE
     &         , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &         , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GM)
     &         , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GP)
     &         , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BM)
     &         , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BP)
     &         , TOL_CLOUD
     &         , A5, B
     &         )
!
            N_EQUATION=N_LAYER+1
!
!           SOLVE THE SET OF EQUATIONS PRODUCED.
            CALL BAND_SOLVER(N_PROFILE, N_EQUATION
     &         , 2, 2
     &         , A5, B
     &         , FLUX_TOTAL_2
     &         , NPD_PROFILE, 2*NPD_LAYER+2
     &         , WORK
     &         )
!
!           EXTRAPOLATE TO THE LIMITING VALUE.
            DO I=1, N_LAYER+1
               DO L=1, N_PROFILE
                  FLUX_TOTAL(L, I)=2.0D+00*FLUX_TOTAL(L, I)
     &               -FLUX_TOTAL_2(L, I)
               ENDDO
            ENDDO
!
         ENDIF
!
      ELSE IF ( (I_SOLVER.EQ.IP_SOLVER_MIX_APP_SCAT).OR.
     &          (I_SOLVER.EQ.IP_SOLVER_MIX_NET_APP_SCAT) ) THEN
!
         CALL MIX_APP_SCAT(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , TRANS_FREE, REFLECT_FREE, S_DOWN_FREE, S_UP_FREE
     &      , TRANS_CLOUD, REFLECT_CLOUD
     &      , S_DOWN_CLOUD, S_UP_CLOUD
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BCF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BCC)
     &      , L_NET
     &      , FLUX_INC_DOWN
     &      , SOURCE_GROUND, ALBEDO_SURFACE_DIFF
     &      , FLUX_TOTAL
     &      )
!
      ELSE IF ( (I_SOLVER.EQ.IP_SOLVER_MIX_DIRECT).OR.
     &          (I_SOLVER.EQ.IP_SOLVER_MIX_DIRECT_NET) ) THEN
!
!        SET THE PARTITIONED SOURCE FUNCTIONS AT THE GROUND.
         IF (ISOLIR.EQ.IP_SOLAR) THEN
            DO L=1, N_PROFILE
               SOURCE_GROUND_FREE(L)=(ALBEDO_SURFACE_DIR(L)
     &            -ALBEDO_SURFACE_DIFF(L))
     &            *(FLUX_DIRECT(L, N_LAYER)
     &            -FLUX_DIRECT_GROUND_CLOUD(L))
               SOURCE_GROUND_CLOUD(L)=(ALBEDO_SURFACE_DIR(L)
     &            -ALBEDO_SURFACE_DIFF(L))
     &            *FLUX_DIRECT_GROUND_CLOUD(L)
            ENDDO
         ELSE 
            DO L=1, N_PROFILE
               SOURCE_GROUND_FREE(L)
     &            =CLOUD_OVERLAP(L, N_LAYER, IP_CLOVLP_BFF)
     &            *SOURCE_GROUND(L)
               SOURCE_GROUND_CLOUD(L)
     &            =CLOUD_OVERLAP(L, N_LAYER, IP_CLOVLP_BCF)
     &            *SOURCE_GROUND(L)
            ENDDO
         ENDIF
!
         CALL SOLVER_MIX_DIRECT(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , TRANS_FREE, REFLECT_FREE, S_DOWN_FREE, S_UP_FREE
     &      , TRANS_CLOUD, REFLECT_CLOUD
     &      , S_DOWN_CLOUD, S_UP_CLOUD
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GCC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFC)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BCF)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BCC)
     &      , L_NET
     &      , FLUX_INC_DOWN
     &      , SOURCE_GROUND_FREE, SOURCE_GROUND_CLOUD
     &      , ALBEDO_SURFACE_DIFF
     &      , FLUX_TOTAL
     &      )
!
      ELSE IF (I_SOLVER.EQ.IP_SOLVER_MIX_FULL) THEN
!
         CALL MIX_SUM_MATRIX(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , TRANS_FREE, REFLECT_FREE
     &      , TRANS_CLOUD, REFLECT_CLOUD
     &      , S_UP_FREE, S_DOWN_FREE, S_UP_CLOUD, S_DOWN_CLOUD
     &      , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &      , FLUX_DIRECT(1, N_LAYER)
     &      , FLUX_INC_DOWN
     &      , N_FREE_PROFILE, I_FREE_PROFILE
     &      , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GFP)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BFP)
     &      , A5, B
     &      )
!
         N_EQUATION=2*N_LAYER+2
!
         CALL BAND_SOLVER(N_PROFILE, N_EQUATION
     &      , 2, 2
     &      , A5, B
     &      , FLUX_TOTAL
     &      , NPD_PROFILE, 2*NPD_LAYER+2
     &      , WORK
     &      )
!
      ELSE IF (I_SOLVER.EQ.IP_SOLVER_MIX_11) THEN
!
         CALL MIX_COLUMN_FULL(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &      , TRANS_FREE, REFLECT_FREE, S_DOWN_FREE, S_UP_FREE
     &      , TRANS_CLOUD, REFLECT_CLOUD
     &      , S_DOWN_CLOUD, S_UP_CLOUD
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_GP)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BM)
     &      , CLOUD_OVERLAP(1, 0, IP_CLOVLP_BP)
     &      , FLUX_INC_DOWN
     &      , SOURCE_GROUND, ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR
     &      , FLUX_DIRECT(1, N_LAYER)
     &      , FLUX_TOTAL
     &      )
!
      ELSE
!
         WRITE(IU_ERR, '(/A)') 
     &      '***ERROR: THE SOLVER SPECIFIED IS NOT VALID HERE.'
         IERR=I_ERR_FATAL
         RETURN
!
      ENDIF
!
!
!
      IF (L_CLEAR) THEN
!
         CALL CLEAR_SUPPLEMENT(IERR, N_PROFILE, N_LAYER, I_SOLVER_CLEAR
     &      , TRANS_FREE, REFLECT_FREE, TRANS_0_FREE, SOURCE_COEFF_FREE
     &      , ISOLIR, FLUX_INC_DIRECT, FLUX_INC_DOWN
     &      , S_DOWN_CLEAR, S_UP_CLEAR
     &      , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR
     &      , SOURCE_GROUND
     &      , L_SCALE_SOLAR, ADJUST_SOLAR_KE
     &      , FLUX_DIRECT_CLEAR, FLUX_TOTAL_CLEAR
     &      )
      ENDIF
!
!
!
      RETURN
      END
