!+ Subroutine to set the elements of the matrix for the net fluxes.
!
! Method:
!	The method is somewhat involved. It is described in a 
!	separate document.
!
! 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_NET_MATRIX(N_PROFILE, N_LAYER, N_CLOUD_TOP
     &   , T_FREE, R_FREE, T_CLOUD, R_CLOUD
     &   , S_UP_FREE, S_DOWN_FREE, S_UP_CLOUD, S_DOWN_CLOUD
     &   , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &   , FLUX_DIRECT_GROUND
     &   , FLUX_INC_DOWN
     &   , N_FREE_PROFILE, I_FREE_PROFILE
     &   , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &   , G_M, G_P
     &   , B_M, B_P
     &   , TOL_CLOUD
     &   , P, V
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
!
!     DUMMY VARIABLES
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOPMOST CLOUDY LAYER
     &   , 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
!
      REAL*8	!, INTENT(IN)
     &     T_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY TRANSMITTANCE
     &   , R_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY REFLECTANCE
     &   , T_FREE(NPD_PROFILE, NPD_LAYER)
!             CLOUD-FREE TRANSMITTANCE
     &   , R_FREE(NPD_PROFILE, NPD_LAYER)
!             CLOUD-FREE REFLEXTANCE
     &   , S_UP_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY UPWARD SOURCE FUNCTION
     &   , S_UP_FREE(NPD_PROFILE, NPD_LAYER)
!             CLOUD-FREE UPWARD SOURCE FUNCTION
     &   , S_DOWN_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY DOWNWARD SOURCE FUNCTION
     &   , S_DOWN_FREE(NPD_PROFILE, NPD_LAYER)
!             CLOUD-FREE DOWNWARD SOURCE FUNCTION
      REAL*8	!, INTENT(IN)
     &     ALBEDO_SURFACE_DIR(NPD_PROFILE)
!             DIRECT SURFACE RELFECTANCE
     &   , ALBEDO_SURFACE_DIFF(NPD_PROFILE)
!             DIFFUSE SURFACE REFLECTANCE
     &   , SOURCE_GROUND(NPD_PROFILE)
!             GROUND SOURCE FUNCTION
     &   , B_M(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , B_P(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , G_M(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , G_P(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , FLUX_DIRECT_GROUND(NPD_PROFILE)
!             DIRECT FLUX AT GROUND
     &   , FLUX_INC_DOWN(NPD_PROFILE)
!             DOWNWARD INCIDENT FLUX
     &   , TOL_CLOUD
!             TOLERANCE FOR RESTORING CONDITIONING WITH CLOUDS
      REAL*8	!, INTENT(OUT)
     &     P(NPD_PROFILE, 5, 2*NPD_LAYER+2)
!             BANDED MATRIX OF EQUATIONS
     &   , V(NPD_PROFILE, 2*NPD_LAYER+2)
!             RIGHT-HAND SIDES OF EQUATIONS
!
!     LOCAL VARIABLES.
      INTEGER
     &     L
!             LOOP VARIABLE
     &   , I
!             LOOP VARIABLE
     &   , LL
!             LOOP VARIABLE
     &   , N_INDEX
!             NUMBER OF INDICES SATISFYING TEST
     &   , INDEX(NPD_PROFILE)
!             INDICES SATISFYING TEST
      REAL*8
     &     DET
!             DETERMINANT IN SETS OF EQUATIONS
     &   , RATIO_DOWN(NPD_PROFILE, 0: NPD_LAYER)
!             RATIO OF TERMS FOR ELIMINATION
     &   , RATIO_UP(NPD_PROFILE, 0: NPD_LAYER)
!             RATIO OF TERMS FOR ELIMINATION
      REAL*8
     &     ZETA_TOP_CLOUD
!             J-SOURCE IN CLOUD
     &   , ZETA_BASE_CLOUD
!             J-SOURCE IN CLOUD
      REAL*8
     &     SIGMA_FREE(NPD_PROFILE, NPD_LAYER)
     &   , DELTA_FREE(NPD_PROFILE, NPD_LAYER)
     &   , SIGMA_CLOUD
     &   , DELTA_CLOUD
     &   , D_SIGMA
     &   , D_DELTA
     &   , PERTURBATION
     &   , SIGMA_MEAN
     &   , DELTA_MEAN
     &   , SIGMA_H
     &   , DELTA_K
     &   , U_PP(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , U_MP(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , U_PM(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , U_MM(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , P_0_P(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , P_0_M(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , P_1_M(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , P_1_P(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , Q_0_P(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , Q_0_M(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , Q_1_M(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , Q_1_P(NPD_PROFILE, 0: NPD_LAYER+1)
      REAL*8
     &     M_M_M1(NPD_PROFILE, 0: NPD_LAYER)
     &   , M_M_0(NPD_PROFILE, 0: NPD_LAYER)
     &   , M_M_P1(NPD_PROFILE, 0: NPD_LAYER)
     &   , H_M_M1(NPD_PROFILE, 0: NPD_LAYER)
     &   , H_M_0(NPD_PROFILE, 0: NPD_LAYER)
     &   , H_M_P1(NPD_PROFILE, 0: NPD_LAYER)
     &   , M_P_M1(NPD_PROFILE, NPD_LAYER)
     &   , M_P_0(NPD_PROFILE, NPD_LAYER)
     &   , M_P_P1(NPD_PROFILE, NPD_LAYER)
     &   , H_P_M1(NPD_PROFILE, NPD_LAYER)
     &   , H_P_0(NPD_PROFILE, NPD_LAYER)
     &   , H_P_P1(NPD_PROFILE, NPD_LAYER)
     &   , W_M(NPD_PROFILE, 0: NPD_LAYER)
     &   , W_P(NPD_PROFILE, NPD_LAYER)
      REAL*8
     &     F_P_M2(NPD_PROFILE, 0: NPD_LAYER)
     &   , F_P_M1(NPD_PROFILE, 0: NPD_LAYER)
     &   , F_P_0(NPD_PROFILE, 0: NPD_LAYER)
     &   , F_P_P1(NPD_PROFILE, 0: NPD_LAYER)
     &   , Q_P(NPD_PROFILE, 0: NPD_LAYER)
     &   , F_M_M2(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , F_M_M1(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , F_M_0(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , F_M_P1(NPD_PROFILE, 0: NPD_LAYER+1)
     &   , Q_M(NPD_PROFILE, 0: NPD_LAYER+1)
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     WHENFGT
!
!
!     SET UP THE EXPRESSIONS FOR J IN TERMS OF N IN CLOUDLESS LAYERS.
!
      DO I=1, N_LAYER
         DO L=1, N_PROFILE
            SIGMA_FREE(L, I)=2.0D+00/(1.0D+00-(T_FREE(L, I)
     &         +R_FREE(L, I))*(1.0D+00-SQRT_TOL_MACHINE))-1.0D+00
            DELTA_FREE(L, I)=2.0D+00/(1.0D+00+(T_FREE(L, I)
     &         -R_FREE(L, I))*(1.0D+00-SQRT_TOL_MACHINE))-1.0D+00
!
!           PERTURB SIGMA AND DELTA TO AVOID LATER ILL-CONDITIONING.
            PERTURBATION=TOL_MACHINE/(SIGMA_FREE(L, I)-DELTA_FREE(L, I)
     &         +SQRT_TOL_MACHINE)
            SIGMA_FREE(L, I)=SIGMA_FREE(L, I)+PERTURBATION
            DELTA_FREE(L, I)=DELTA_FREE(L, I)-PERTURBATION
!
            F_M_M2(L, I)=0.0D+00
            F_M_M1(L, I)=0.5D+00*(SIGMA_FREE(L, I)+DELTA_FREE(L, I))
            F_M_0(L, I)=-0.5D+00*(SIGMA_FREE(L, I)-DELTA_FREE(L, I))
            F_M_P1(L, I)=0.0D+00
            Q_M(L, I)=0.5D+00*SIGMA_FREE(L, I)*(S_UP_FREE(L, I)
     &         +S_DOWN_FREE(L, I))-0.5D+00*DELTA_FREE(L, I)
     &         *(S_DOWN_FREE(L, I)-S_UP_FREE(L, I))+S_UP_FREE(L, I)
            F_P_M2(L, I)=0.0D+00
            F_P_M1(L, I)=0.5D+00*(SIGMA_FREE(L, I)-DELTA_FREE(L, I))
            F_P_0(L, I)=-0.5D+00*(SIGMA_FREE(L, I)+DELTA_FREE(L, I))
            F_P_P1(L, I)=0.0D+00
            Q_P(L, I)=0.5D+00*SIGMA_FREE(L, I)*(S_UP_FREE(L, I)
     &         +S_DOWN_FREE(L, I))+0.5D+00*DELTA_FREE(L, I)
     &         *(S_DOWN_FREE(L, I)-S_UP_FREE(L, I))+S_DOWN_FREE(L, I)
         ENDDO
      ENDDO
!
!     SET UP THE EQUATIONS FOR THE DIFFERENCES IN THE CLEAR AND CLOUDY
!     NET AND SUMMED FLUXES AND SOLVE THEM. THESE SHOULD BE THOUGHT OF
!     AS EQUATIONS IN THE I'TH LAYER.
!
      DO I=N_CLOUD_TOP, N_LAYER
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
            SIGMA_CLOUD=2.0D+00/(1.0D+00-(T_CLOUD(L, I)+R_CLOUD(L, I))
     &         *(1.0D+00-TOL_CLOUD))-1.0D+00
            DELTA_CLOUD=2.0D+00/(1.0D+00+(T_CLOUD(L, I)-R_CLOUD(L, I))
     &         *(1.0D+00-TOL_CLOUD))-1.0D+00
!
            D_SIGMA=SIGMA_FREE(L, I)-SIGMA_CLOUD
            D_DELTA=DELTA_FREE(L, I)-DELTA_CLOUD
!
!           PERTURB AWAY FROM ILL-CONDITIONING:
            PERTURBATION=TOL_CLOUD*SIGN(TOL_CLOUD, D_SIGMA)
     &         /(D_SIGMA+SIGN(TOL_CLOUD, D_SIGMA))
            SIGMA_CLOUD=SIGMA_CLOUD-PERTURBATION
            D_SIGMA=D_SIGMA+PERTURBATION
            PERTURBATION=TOL_CLOUD*SIGN(TOL_CLOUD, D_DELTA)
     &         /(D_DELTA+SIGN(TOL_CLOUD, D_DELTA))
            DELTA_CLOUD=DELTA_CLOUD-PERTURBATION
            D_DELTA=D_DELTA+PERTURBATION

            SIGMA_MEAN=0.5D+00*(SIGMA_FREE(L, I)+SIGMA_CLOUD)
            DELTA_MEAN=0.5D+00*(DELTA_FREE(L, I)+DELTA_CLOUD)
            ZETA_TOP_CLOUD=0.5D+00*SIGMA_CLOUD*(S_UP_CLOUD(L, I)
     &         +S_DOWN_CLOUD(L, I))-0.5D+00*DELTA_CLOUD
     &         *(S_DOWN_CLOUD(L, I)-S_UP_CLOUD(L, I))+S_UP_CLOUD(L, I)
            ZETA_BASE_CLOUD=0.5D+00*SIGMA_CLOUD*(S_UP_CLOUD(L, I)
     &         +S_DOWN_CLOUD(L, I))+0.5D+00*DELTA_CLOUD
     &         *(S_DOWN_CLOUD(L, I)-S_UP_CLOUD(L, I))+S_DOWN_CLOUD(L, I)
!
            SIGMA_H=(Q_M(L, I)+ZETA_TOP_CLOUD+Q_P(L, I)
     &         +ZETA_BASE_CLOUD)/D_SIGMA
            DELTA_K=(Q_P(L, I)+ZETA_BASE_CLOUD-Q_M(L, I)
     &         -ZETA_TOP_CLOUD)/D_DELTA
            U_PP(L, I)=Q_P(L, I)-ZETA_BASE_CLOUD
     &         +((1.0D+00-DELTA_MEAN)*DELTA_K
     &         +(1.0D+00-SIGMA_MEAN)*SIGMA_H)
            U_PM(L, I)=Q_P(L, I)-ZETA_BASE_CLOUD
     &         -((1.0D+00+DELTA_MEAN)*DELTA_K
     &         +(1.0D+00+SIGMA_MEAN)*SIGMA_H)
            U_MP(L, I)=Q_M(L, I)-ZETA_TOP_CLOUD
     &         +((1.0D+00+DELTA_MEAN)*DELTA_K
     &         -(1.0D+00+SIGMA_MEAN)*SIGMA_H)
            U_MM(L, I)=Q_M(L, I)-ZETA_TOP_CLOUD
     &         -((1.0D+00-DELTA_MEAN)*DELTA_K
     &         -(1.0D+00-SIGMA_MEAN)*SIGMA_H)
!
            P_0_P(L, I)=SIGMA_FREE(L, I)*(SIGMA_CLOUD-1.0D+00)/D_SIGMA
     &         +DELTA_CLOUD*(DELTA_FREE(L, I)-1.0D+00)/D_DELTA
            P_0_M(L, I)=-SIGMA_FREE(L, I)*(SIGMA_CLOUD+1.0D+00)/D_SIGMA
     &         -DELTA_CLOUD*(DELTA_FREE(L, I)+1.0D+00)/D_DELTA
            P_1_P(L, I)=SIGMA_CLOUD*(1.0D+00-SIGMA_FREE(L, I))/D_SIGMA
     &         -DELTA_CLOUD*(1.0D+00-DELTA_FREE(L, I))/D_DELTA
            P_1_M(L, I)=SIGMA_CLOUD*(SIGMA_FREE(L, I)+1.0D+00)/D_SIGMA
     &         -DELTA_CLOUD*(DELTA_FREE(L, I)+1.0D+00)/D_DELTA
            Q_0_P(L, I)=(SIGMA_FREE(L, I)-1.0D+00)/D_SIGMA
     &         +(DELTA_CLOUD-1.0D+00)/D_DELTA
            Q_0_M(L, I)=-(SIGMA_FREE(L, I)+1.0D+00)/D_SIGMA
     &         -(DELTA_CLOUD+1.0D+00)/D_DELTA
            Q_1_P(L, I)=(SIGMA_CLOUD-1.0D+00)/D_SIGMA
     &         -(DELTA_CLOUD-1.0D+00)/D_DELTA
            Q_1_M(L, I)=-(SIGMA_CLOUD+1.0D+00)/D_SIGMA
     &         +(DELTA_CLOUD+1.0D+00)/D_DELTA
!
         ENDDO
!
         DO LL=1, N_FREE_PROFILE(I)
            L=I_FREE_PROFILE(LL, I)
            P_0_P(L, I)=1.0D+00
            P_1_P(L, I)=0.0D+00
            Q_0_P(L, I)=1.0D+00
            Q_1_P(L, I)=0.0D+00
            U_PP(L, I)=0.0D+00
            U_MM(L, I)=0.0D+00
         ENDDO
!
      ENDDO
!
!     COMPLETE THE ASSIGNMENT OF COEFFICIENTS IN THE LAYERS ADJACENT TO
!     THE CLOUDY PART OF THE COLUMN.
      DO L=1, N_PROFILE
         P_0_P(L, N_CLOUD_TOP-1)=1.0D+00
         P_1_P(L, N_CLOUD_TOP-1)=0.0D+00
         Q_0_P(L, N_CLOUD_TOP-1)=1.0D+00
         Q_1_P(L, N_CLOUD_TOP-1)=0.0D+00
         U_PP(L, N_CLOUD_TOP-1)=0.0D+00
      ENDDO
      DO L=1, N_PROFILE
         P_0_P(L, N_LAYER+1)=1.0D+00
         P_1_P(L, N_LAYER+1)=0.0D+00
         Q_0_P(L, N_LAYER+1)=1.0D+00
         Q_1_P(L, N_LAYER+1)=0.0D+00
         U_MM(L, N_LAYER+1)=0.0D+00
      ENDDO
!
!     FORM THE BASIC EQUATIONS IN N FOR THE TOTAL DIFFERENCE AND SUMMED
!     FLUXES. THESE ARE EQUATIONS AT BOUNDARIES.
!
!     WE INITIALIZE THE RATIOS TO 0 AND ADJUST THEM 
!     WHEN THERE ARE CLOUDS.
      DO I=N_CLOUD_TOP-1, N_LAYER
         DO L=1, N_PROFILE
            RATIO_UP(L, I)=0.0D+00
            RATIO_DOWN(L, I)=0.0D+00
         ENDDO
      ENDDO
!
      DO I=N_CLOUD_TOP, N_LAYER
!
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
!
!           CONDITION ON UPWARD FLUXES:
            M_P_M1(L, I)=2.0D+00*P_1_M(L, I)
            M_P_0(L, I)=2.0D+00*P_0_M(L, I)-B_P(L, I)
     &         -B_M(L, I)*P_0_P(L, I+1)
            M_P_P1(L, I)=-B_M(L, I)*P_1_P(L, I+1)
            H_P_M1(L, I)=2.0D+00*Q_1_M(L, I)
            H_P_0(L, I)=2.0D+00*Q_0_M(L, I)+B_P(L, I)
     &         +B_M(L, I)*Q_0_P(L, I+1)
            H_P_P1(L, I)=B_M(L, I)*Q_1_P(L, I+1)
            W_P(L, I)=-2.0D+00*U_PM(L, I)+B_M(L, I)*U_MM(L, I+1)
            RATIO_UP(L, I)=ABS(H_P_P1(L, I))
!
         ENDDO
      ENDDO
!
      DO I=N_CLOUD_TOP-1, N_LAYER-1
         DO LL=1, N_CLOUD_PROFILE(I+1)
            L=I_CLOUD_PROFILE(LL, I+1)
!
!           CONDITION ON DOWNWARD FLUXES:
            M_M_M1(L, I)=G_M(L, I)*P_1_P(L, I)
            M_M_0(L, I)=-2.0D+00*P_0_M(L, I+1)+G_P(L, I)
     &         +G_M(L, I)*P_0_P(L, I)
            M_M_P1(L, I)=-2.0D+00*P_1_M(L, I+1)
            H_M_M1(L, I)=G_M(L, I)*Q_1_P(L, I)
            H_M_0(L, I)=2.0D+00*Q_0_M(L, I+1)+G_P(L, I)
     &         +G_M(L, I)*Q_0_P(L, I)
            H_M_P1(L, I)=2.0D+00*Q_1_M(L, I+1)
            W_M(L, I)=-2.0D+00*U_MP(L, I+1)+G_M(L, I)*U_PP(L, I)
            RATIO_DOWN(L, I)=ABS(H_M_M1(L, I))
!
         ENDDO
      ENDDO
!
!     ELIMINATE THE VALUES OF J AT OTHER LAYERS:
!     THESE OPERATIONS ARE NEEDED ONLY WHEN THERE IS CLOUD BOTH ABOVE
!     AND BELOW THE LEVEL. THE ARRAYS FOR TESTING SHOULD BE O(1), SO
!     WE DO NOT NEED TO SCALE THEIR ENTRIES.
      DO I=N_CLOUD_TOP, N_LAYER-1
!
!        DOWNWARD FLUXES:
         CALL WHENFGT(N_PROFILE, RATIO_DOWN(1, I), 1, TOL_DIV
     &      , INDEX, N_INDEX)
!        INDIRECT ADDRESSING CREATES A POTENTIAL DATA DEPENDENCY
!        SINCE THE COMPILER CANNOT KNOW THAT NO ELEMENT IS REPEATED
!        IN INDEX: WE THEREFORE NEED A DIRECTIVE TO FORCE VECTORIZATION.
         DO LL=1, N_INDEX
            L=INDEX(LL)
            RATIO_DOWN(L, I)=H_M_M1(L, I)/H_P_M1(L, I)
!           THIS ASSIGNMENT IS NEEDED LATER.
            H_M_M1(L, I)=0.0D+00
            H_M_0(L, I)=H_M_0(L, I)-RATIO_DOWN(L, I)*H_P_0(L, I)
            H_M_P1(L, I)=H_M_P1(L, I)-RATIO_DOWN(L, I)*H_P_P1(L, I)
            M_M_M1(L, I)=M_M_M1(L, I)-RATIO_DOWN(L, I)*M_P_M1(L, I)
            M_M_0(L, I)=M_M_0(L, I)-RATIO_DOWN(L, I)*M_P_0(L, I)
            M_M_P1(L, I)=M_M_P1(L, I)-RATIO_DOWN(L, I)*M_P_P1(L, I)
            W_M(L, I)=W_M(L, I)-RATIO_DOWN(L, I)*W_P(L, I)
         ENDDO
!
!        UPWARD FLUXES:
         CALL WHENFGT(N_PROFILE, RATIO_UP(1, I), 1, TOL_DIV
     &      , INDEX, N_INDEX)
!        INDIRECT ADDRESSING CREATES A POTENTIAL DATA DEPENDENCY
!        SINCE THE COMPILER CANNOT KNOW THAT NO ELEMENT IS REPEATED
!        IN INDEX: WE THEREFORE NEED A DIRECTIVE TO FORCE VECTORIZATION.
         DO LL=1, N_INDEX
            L=INDEX(LL)
            RATIO_UP(L, I)=H_P_P1(L, I)/H_M_P1(L, I)
            H_P_M1(L, I)=H_P_M1(L, I)-RATIO_UP(L, I)*H_M_M1(L, I)
            H_P_0(L, I)=H_P_0(L, I)-RATIO_UP(L, I)*H_M_0(L, I)
            M_P_M1(L, I)=M_P_M1(L, I)-RATIO_UP(L, I)*M_M_M1(L, I)
            M_P_0(L, I)=M_P_0(L, I)-RATIO_UP(L, I)*M_M_0(L, I)
            M_P_P1(L, I)=M_P_P1(L, I)-RATIO_UP(L, I)*M_M_P1(L, I)
            W_P(L, I)=W_P(L, I)-RATIO_UP(L, I)*W_M(L, I)
         ENDDO
      ENDDO
!
!     SOLVE THESE PAIRS OF EQUATIONS TO FIND EXPRESSIONS FOR J AT
!     THE BOUNDARIES OF THE LAYERS: THESE EXPRESSIONS ARE REQUIRED ONLY
!     WHEN THE LAYER CONTAINS CLOUD.
      DO I=N_CLOUD_TOP, N_LAYER
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
            DET=H_P_M1(L, I)*H_M_P1(L, I-1)-H_P_0(L, I)*H_M_0(L, I-1)
            F_P_M2(L, I)=-H_P_M1(L, I)*M_M_M1(L, I-1)/DET
            F_P_M1(L, I)=(H_M_0(L, I-1)*M_P_M1(L, I)
     &         -H_P_M1(L, I)*M_M_0(L, I-1))/DET
            F_P_0(L, I)=(H_M_0(L, I-1)*M_P_0(L, I)
     &         -H_P_M1(L, I)*M_M_P1(L, I-1))/DET
            F_P_P1(L, I)=H_M_0(L, I-1)*M_P_P1(L, I)/DET
            Q_P(L, I)=(H_M_0(L, I-1)*W_P(L, I)
     &         -H_P_M1(L, I)*W_M(L, I-1))/DET
            F_M_M2(L, I)=H_P_0(L, I)*M_M_M1(L, I-1)/DET
            F_M_M1(L, I)=(H_P_0(L, I)*M_M_0(L, I-1)
     &         -H_M_P1(L, I-1)*M_P_M1(L, I))/DET
            F_M_0(L, I)=(H_P_0(L, I)*M_M_P1(L, I-1)
     &         -H_M_P1(L, I-1)*M_P_0(L, I))/DET
            F_M_P1(L, I)=-H_M_P1(L, I-1)*M_P_P1(L, I)/DET
            Q_M(L, I)=(H_P_0(L, I)*W_M(L, I-1)
     &         -H_M_P1(L, I-1)*W_P(L, I))/DET
         ENDDO
      ENDDO
!
!
!     SET THE BOUNDARY CONDITIONS:
      DO L=1, N_PROFILE
         F_P_M2(L, 0)=0.0D+00
         F_P_M1(L, 0)=0.0D+00
         F_P_0(L, 0)=-1.0D+00
         F_P_P1(L, 0)=0.0D+00
         Q_P(L, 0)=2.0*FLUX_INC_DOWN(L)
         F_M_M2(L, N_LAYER+1)=0.0D+00
         F_M_M1(L, N_LAYER+1)=(1.0D+00+ALBEDO_SURFACE_DIFF(L))
     &      /(1.0D+00-ALBEDO_SURFACE_DIFF(L))
         F_M_0(L, N_LAYER+1)=0.0D+00
         F_M_P1(L, N_LAYER+1)=0.0D+00
         Q_M(L, N_LAYER+1)=2.0D+00*(SOURCE_GROUND(L)
     &      +(ALBEDO_SURFACE_DIR(L)-ALBEDO_SURFACE_DIFF(L))
     &      *FLUX_DIRECT_GROUND(L))
     &      /(1.0D+00-ALBEDO_SURFACE_DIFF(L))
      ENDDO
!
!     FORMULATE THE FINAL SET OF EQUATIONS AT THE BOUNDARIES 
!     OF THE LAYERS.
      DO I=0, N_LAYER
         DO L=1, N_PROFILE
            P(L, 5, I+1)=F_P_M2(L, I)
            P(L, 4, I+1)=F_P_M1(L, I)-F_M_M2(L, I+1)
            P(L, 3, I+1)=F_P_0(L, I)-F_M_M1(L, I+1)
            P(L, 2, I+1)=F_P_P1(L, I)-F_M_0(L, I+1)
            P(L, 1, I+1)=-F_M_P1(L, I+1)
            V(L, I+1)=Q_M(L, I+1)-Q_P(L, I)
         ENDDO
      ENDDO
!
!
!
      RETURN
      END
