!+ Subroutine to calculate total upward and downward fluxes.
!
! Method:
!	The equations for the clear and cloudy upward and downward 
!	fluxes can be re-arranged, as discussed in a separate document,
!	to give equations for the total fluxes. The equations for
!	these fluxes are set.
!
! 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_SUM_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
     &   , GM_FM, GM_FP, BT_FM, BT_FP
     &   , 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
!             TOP LCOUDY 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
     &   , BT_FM(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , BT_FP(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , GM_FM(NPD_PROFILE, 0: NPD_LAYER)
!             ENERGY TRANSFER FUNCTION
     &   , GM_FP(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
      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
      REAL*8
     &     DET(NPD_PROFILE, 0: NPD_LAYER)
!             DETERMINANTS FOR DIFF. EQNS.
     &   , DEL_M_MM(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_M_0M(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_M_PP(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_M_2P(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_P_MM(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_P_0M(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_P_PP(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , DEL_P_2P(NPD_PROFILE, 0: NPD_LAYER)
!             COEFFICIENTS FOR DIFFERENCE EQNS.
     &   , KSI_M(NPD_PROFILE, 0: NPD_LAYER)
!             DIFFERENCE SOURCE TERMS
     &   , KSI_P(NPD_PROFILE, 0: NPD_LAYER)
!             DIFFERENCE SOURCE TERMS
      REAL*8
     &     ETA_M_MM(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , ETA_M_0M(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , ETA_M_PP(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , ETA_P_0M(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , ETA_P_PP(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , ETA_P_2P(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , DELTA_M_M(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , DELTA_M_P(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , DELTA_P_M(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , DELTA_P_P(NPD_PROFILE)
!             COEFFICIENT IN DIFFERENCE EQN.
     &   , SIGMA_M(NPD_PROFILE)
!             SOURCE TERM IN DIFFERENCE EQN.
     &   , SIGMA_P(NPD_PROFILE)
!             SOURCE TERM IN DIFFERENCE EQN.
      REAL*8
     &     H_M_MM
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , H_M_PP
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , H_P_MM
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , H_P_PP
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , D_M_MM
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , D_M_PP
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , D_P_MM
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , D_P_PP
!             COEFFICIENT IN SUMMED EQUATIONS
     &   , S_M
!             SOURCE TERM IN SUMMED EQNS.
     &   , S_P
!             SOURCE TERM IN SUMMED EQNS.
      REAL*8
     &     U_M_M2M
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_M_MM
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_M_0M
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_M_0P
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_M_PP
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_M_2P
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_M2M
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_MM
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_0M
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_0P
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_PP
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , U_P_2P
!             COEFFICIENT IN HEPTADIAGONAL EQN.
     &   , Q_M
!             SOURCE TERM IN HEPTADIAGONAL EQN.
     &   , Q_P
!             SOURCE TERM IN HEPTADIAGONAL EQN.
      REAL*8
     &     DIFF_T(NPD_PROFILE, NPD_LAYER)
!             DIFFERENCES IN TRANSMISSION
     &   , DIFF_R(NPD_PROFILE, NPD_LAYER)
!             DIFFERENCES IN REFLECTION
     &   , FACTOR(NPD_PROFILE, NPD_LAYER)
!             FACTOR IN ELIMINATING EQUATIONS
!
!
!     PRECALCULATE ARRAYS FOR ELIMINATION AND PERTURB SPECIAL CASES
!     AWAY FROM SINGULARITIES.
      DO I=N_CLOUD_TOP, N_LAYER
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
            DIFF_T(L, I)=T_FREE(L, I)-T_CLOUD(L, I)
            IF (DIFF_T(L, I).LT.TOL_DIV*T_FREE(L, I)) THEN
               T_CLOUD(L, I)=(1.0D+00-TOL_DIV)*T_CLOUD(L, I)
               DIFF_T(L, I)=T_FREE(L, I)-T_CLOUD(L, I)
            ENDIF
            DIFF_R(L, I)=R_FREE(L, I)-R_CLOUD(L, I)
            FACTOR(L, I)=DIFF_R(L, I)/DIFF_T(L, I)
            IF (ABS(1.0D+00-ABS(FACTOR(L, I))).LT.TOL_DIV) THEN
               R_CLOUD(L, I)=R_CLOUD(L, I)-2.0D+00*TOL_DIV*DIFF_T(L, I)
               FACTOR(L, I)=FACTOR(L, I)-2.0D+00*TOL_DIV
            ENDIF
         ENDDO
      ENDDO
!
!     SET UP THE EQUATIONS FOR THE DIFFERENCES IN THE CLEAR AND CLOUDY
!     FLUXES AND SOLVE THEM. THESE SHOULD BE THOUGHT OF AS EQUATIONS
!     AT THE I'TH LEVEL.

!     COEFFICIENTS AT THE TOP OF THE CLOUDY REGION 
!     (DOWNWARD TERMS ONLY):
      DO L=1, N_PROFILE
         DET(L, N_CLOUD_TOP-1)=1.0D+00
         DEL_M_MM(L, N_CLOUD_TOP-1)=0.0D+00
         DEL_M_0M(L, N_CLOUD_TOP-1)=1.0D+00
         DEL_M_PP(L, N_CLOUD_TOP-1)=0.0D+00
         DEL_M_2P(L, N_CLOUD_TOP-1)=0.0D+00
         KSI_M(L, N_CLOUD_TOP-1)=0.0D+00
      ENDDO
!
!     INTERMEDIATE LEVELS:
      DO I=N_CLOUD_TOP, N_LAYER-1
!
!        DOWNWARD FLUX:
         DO LL=1, N_FREE_PROFILE(I)
            L=I_FREE_PROFILE(LL, I)
!           DO NOT CANCEL TRANSMITTANCES SO THAT DET IS ALWAYS O(T^2)
!           FOR CONFORMITY IN TEST.
            ETA_M_MM(L)=0
            ETA_M_0M(L)=T_FREE(L, I)
            ETA_M_PP(L)=0
            SIGMA_M(L)=0
            DELTA_M_M(L)=T_FREE(L, I)
            DELTA_M_P(L)=0
         ENDDO
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
            ETA_M_MM(L)=2.0D+00*T_CLOUD(L, I)*T_FREE(L, I)
            ETA_M_0M(L)=-(T_FREE(L, I)+T_CLOUD(L, I))
            ETA_M_PP(L)=R_FREE(L, I)*T_CLOUD(L, I)*BT_FP(L, I)
     &         +R_CLOUD(L, I)*T_FREE(L, I)*(2.0D+00-BT_FP(L, I))
            SIGMA_M(L)=2.0D+00*(T_CLOUD(L, I)*S_DOWN_FREE(L, I)
     &         +T_FREE(L, I)*S_DOWN_CLOUD(L, I))
            DELTA_M_M(L)=-DIFF_T(L, I)
            DELTA_M_P(L)=-BT_FM(L, I)*(R_FREE(L, I)*T_CLOUD(L, I)
     &         -R_CLOUD(L, I)*T_FREE(L, I))
         ENDDO
!
!        UPWARD FLUX:
         DO LL=1, N_FREE_PROFILE(I+1)
            L=I_FREE_PROFILE(LL, I+1)
            ETA_P_0M(L)=0
            ETA_P_PP(L)=T_FREE(L, I+1)
            ETA_P_2P(L)=0
            SIGMA_P(L)=0
            DELTA_P_M(L)=0
            DELTA_P_P(L)=T_FREE(L, I+1)
         ENDDO
         DO LL=1, N_CLOUD_PROFILE(I+1)
            L=I_CLOUD_PROFILE(LL, I+1)
            ETA_P_0M(L)=T_CLOUD(L, I+1)*R_FREE(L, I+1)*GM_FP(L, I)
     &         +T_FREE(L, I+1)*R_CLOUD(L, I+1)*(2.0D+00-GM_FP(L, I))
            ETA_P_PP(L)=-(T_FREE(L, I+1)+T_CLOUD(L, I+1))
            ETA_P_2P(L)=2.0D+00*T_FREE(L, I+1)*T_CLOUD(L, I+1)
            SIGMA_P(L)=2.0D+00*(T_CLOUD(L, I+1)*S_UP_FREE(L, I+1)
     &         +T_FREE(L, I+1)*S_UP_CLOUD(L, I+1))
            DELTA_P_M(L)=-GM_FM(L, I)
     &         *(T_CLOUD(L, I+1)*R_FREE(L, I+1)
     &         -T_FREE(L, I+1)*R_CLOUD(L, I+1))
            DELTA_P_P(L)=-DIFF_T(L, I+1)
         ENDDO
!
!        SOLVE THE EQUATIONS:
!        THESE COEFFICIENTS SHOULD ALL BE DIVIDED BY DET; WE CAN
!        SAVE OPERATIONS BY DIVIDING THE D-VARIABLES BELOW BY DET
!        INSTEAD OF DIVIDING THESE.
         DO L=1, N_PROFILE
            DET(L, I)=DELTA_P_P(L)*DELTA_M_M(L)
     &         -DELTA_P_M(L)*DELTA_M_P(L)
            DEL_M_MM(L, I)=DELTA_P_P(L)*ETA_M_MM(L)
            DEL_M_0M(L, I)=(DELTA_P_P(L)*ETA_M_0M(L)
     &         -DELTA_M_P(L)*ETA_P_0M(L))
            DEL_M_PP(L, I)=(DELTA_P_P(L)*ETA_M_PP(L)
     &         -DELTA_M_P(L)*ETA_P_PP(L))
            DEL_M_2P(L, I)=-DELTA_M_P(L)*ETA_P_2P(L)
            KSI_M(L, I)=(DELTA_P_P(L)*SIGMA_M(L)
     &         -DELTA_M_P(L)*SIGMA_P(L))
            DEL_P_MM(L, I)=-DELTA_P_M(L)*ETA_M_MM(L)
            DEL_P_0M(L, I)=(ETA_P_0M(L)*DELTA_M_M(L)
     &         -ETA_M_0M(L)*DELTA_P_M(L))
            DEL_P_PP(L, I)=(ETA_P_PP(L)*DELTA_M_M(L)
     &         -ETA_M_PP(L)*DELTA_P_M(L))
            DEL_P_2P(L, I)=DELTA_M_M(L)*ETA_P_2P(L)
            KSI_P(L, I)=(DELTA_M_M(L)*SIGMA_P(L)
     &         -DELTA_P_M(L)*SIGMA_M(L))
!
         ENDDO
!        DEAL WITH ILL-CONDITIONED CASES.
         DO L=1, N_PROFILE
            IF (DET(L, I).LT.(TOL_DIV*T_FREE(L, I)*T_FREE(L, I+1))) THEN
!              FREAK COMBINATION OF OPTICAL PROPERTIES AND GEOMETRY:
!              HERE PERTURB THE UPPER CLOUDY TRANSMISSION.
               IF ( ((T_CLOUD(L, I)-T_FREE(L, I))
     &            * (T_CLOUD(L, I+1)-T_FREE(L, I+1)))
     &            .GT.(TOL_DIV*T_FREE(L, I)*T_FREE(L, I+1))) THEN
                  T_CLOUD(L, I)=(1.0D+00-1.0D+02*TOL_DIV)*T_CLOUD(L, I)
                  ETA_M_MM(L)=2.0D+00*T_CLOUD(L, I)*T_FREE(L, I)
                  ETA_M_0M(L)=-(T_FREE(L, I)+T_CLOUD(L, I))
                  ETA_M_PP(L)=R_FREE(L, I)*T_CLOUD(L, I)*BT_FP(L, I)
     &               +R_CLOUD(L, I)*T_FREE(L, I)*(2.0D+00-BT_FP(L, I))
                  SIGMA_M(L)=2.0D+00*(T_CLOUD(L, I)*S_DOWN_FREE(L, I)
     &               +T_FREE(L, I)*S_DOWN_CLOUD(L, I))
                  DELTA_M_M(L)=T_CLOUD(L, I)-T_FREE(L, I)
                  DELTA_M_P(L)=-BT_FM(L, I)*(R_FREE(L, I)*T_CLOUD(L, I)
     &               -R_CLOUD(L, I)*T_FREE(L, I))
               ENDIF
!
!              RESOLVE THE EQUATIONS:
               DET(L, I)=DELTA_P_P(L)*DELTA_M_M(L)
     &            -DELTA_P_M(L)*DELTA_M_P(L)
               DEL_M_MM(L, I)=DELTA_P_P(L)*ETA_M_MM(L)
               DEL_M_0M(L, I)=(DELTA_P_P(L)*ETA_M_0M(L)
     &            -DELTA_M_P(L)*ETA_P_0M(L))
               DEL_M_PP(L, I)=(DELTA_P_P(L)*ETA_M_PP(L)
     &            -DELTA_M_P(L)*ETA_P_PP(L))
               DEL_M_2P(L, I)=-DELTA_M_P(L)*ETA_P_2P(L)
               KSI_M(L, I)=(DELTA_P_P(L)*SIGMA_M(L)
     &            -DELTA_M_P(L)*SIGMA_P(L))
               DEL_P_MM(L, I)=-DELTA_P_M(L)*ETA_M_MM(L)
               DEL_P_0M(L, I)=(ETA_P_0M(L)*DELTA_M_M(L)
     &            -ETA_M_0M(L)*DELTA_P_M(L))
               DEL_P_PP(L, I)=(ETA_P_PP(L)*DELTA_M_M(L)
     &            -ETA_M_PP(L)*DELTA_P_M(L))
               DEL_P_2P(L, I)=DELTA_M_M(L)*ETA_P_2P(L)
               KSI_P(L, I)=(DELTA_M_M(L)*SIGMA_P(L)
     &            -DELTA_P_M(L)*SIGMA_M(L))
            ENDIF
         ENDDO
      ENDDO
!
!     THE BOTTOM LEVEL (UPWARD TERMS ONLY):
      DO L=1, N_PROFILE
         DET(L, N_LAYER)=1.0D+00
         DEL_P_MM(L, N_LAYER)=0.0D+00
         DEL_P_0M(L, N_LAYER)=0.0D+00
         DEL_P_PP(L, N_LAYER)=1.0D+00
         DEL_P_2P(L, N_LAYER)=0.0D+00
         KSI_P(L, N_LAYER)=0.0D+00
      ENDDO
!
!     ASSEMBLE THE SUMMED EQUATIONS WORKING OVER LAYERS.
!
!     REGION ABOVE CLOUDS:
      DO I=1, N_CLOUD_TOP-1
         DO L=1, N_PROFILE
!           UPWARD EQUATION:
            P(L, 5, 2*I)=0.0D+00
            P(L, 4, 2*I)=1.0D+00
            P(L, 3, 2*I)=-R_FREE(L, I)
            P(L, 2, 2*I)=-T_FREE(L, I)
            P(L, 1, 2*I)=0.0D+00
            V(L, 2*I)=S_UP_FREE(L, I)
!           DOWNWARD EQUATION:
            P(L, 5, 2*I+1)=0.0D+00
            P(L, 4, 2*I+1)=-T_FREE(L, I)
            P(L, 3, 2*I+1)=-R_FREE(L, I)
            P(L, 2, 2*I+1)=1.0D+00
            P(L, 1, 2*I+1)=0.0D+00
            V(L, 2*I+1)=S_DOWN_FREE(L, I)
         ENDDO
      ENDDO
!
!     REGION WHERE CLOUDS ARE POSSIBLE:
      DO I=N_CLOUD_TOP, N_LAYER
         DO LL=1, N_CLOUD_PROFILE(I)
            L=I_CLOUD_PROFILE(LL, I)
!           SUMMED DOWNWARD EQUATIONS:
            D_M_MM=GM_FM(L, I-1)*DIFF_T(L, I)/DET(L, I-1)
            H_M_MM=T_FREE(L, I)*GM_FP(L, I-1)
     &          +T_CLOUD(L, I)*(2.0D+00-GM_FP(L, I-1))
            H_M_PP=R_FREE(L, I)*BT_FP(L, I)+R_CLOUD(L, I)
     &          *(2.0D+00-BT_FP(L, I))
            D_M_PP=BT_FM(L, I)*DIFF_R(L, I)/DET(L, I)
            S_M=2.0D+00*(S_DOWN_FREE(L, I)+S_DOWN_CLOUD(L, I))
!
!           SUMMED UPWARD EQUATIONS:
            D_P_MM=GM_FM(L, I-1)*DIFF_R(L, I)/DET(L, I-1)
            H_P_MM=R_FREE(L, I)*GM_FP(L, I-1)
     &         +R_CLOUD(L, I)*(2.0D+00-GM_FP(L, I-1))
            H_P_PP=T_FREE(L, I)*BT_FP(L, I)+T_CLOUD(L, I)
     &         *(2.0D+00-BT_FP(L, I))
            D_P_PP=BT_FM(L, I)*DIFF_T(L, I)/DET(L, I)
            S_P=2.0D+00*(S_UP_FREE(L, I)+S_UP_CLOUD(L, I))
!
!           ESTABLISH THE COEFFICIENTS OF
!           THE PRIMITIVE HEPTADIAGONAL EQUATIONS:
!           EQUATION FOR FLUX_UP(I-1):
            U_P_M2M=D_P_MM*DEL_M_MM(L, I-1)
            U_P_MM=D_P_MM*DEL_M_0M(L, I-1)+H_P_MM+D_P_PP*DEL_P_MM(L, I)
            U_P_0M=D_P_PP*DEL_P_0M(L, I)
            U_P_0P=D_P_MM*DEL_M_PP(L, I-1)-2.0D+00
            U_P_PP=D_P_MM*DEL_M_2P(L, I-1)+H_P_PP+D_P_PP*DEL_P_PP(L, I)
            U_P_2P=D_P_PP*DEL_P_2P(L, I)
            Q_P=-S_P-D_P_MM*KSI_M(L, I-1)-D_P_PP*KSI_P(L, I)
!           EQUATION FOR FLUX_DOWN(I-1):
            U_M_M2M=D_M_MM*DEL_M_MM(L, I-1)
            U_M_MM=D_M_MM*DEL_M_0M(L, I-1)+H_M_MM+D_M_PP*DEL_P_MM(L, I)
            U_M_0M=-2.0D+00+D_M_PP*DEL_P_0M(L, I)
            U_M_0P=D_M_MM*DEL_M_PP(L, I-1)
            U_M_PP=D_M_MM*DEL_M_2P(L, I-1)+H_M_PP+D_M_PP*DEL_P_PP(L, I)
            U_M_2P=D_M_PP*DEL_P_2P(L, I)
            Q_M=-S_M-D_M_MM*KSI_M(L, I-1)-D_M_PP*KSI_P(L, I)
!
!
!           ASSIGN THE COEFFICIENTS OF THE PENTADIAGONAL MATRIX:
!
!           DOWNWARD EQUATION:
            P(L, 5, 2*I)=U_M_M2M-FACTOR(L, I)*U_P_M2M
            P(L, 4, 2*I)=U_M_0P-FACTOR(L, I)*U_P_0P
            P(L, 3, 2*I)=U_M_MM-FACTOR(L, I)*U_P_MM
            P(L, 2, 2*I)=U_M_PP-FACTOR(L, I)*U_P_PP
            P(L, 1, 2*I)=U_M_0M-FACTOR(L, I)*U_P_0M
            V(L, 2*I)=Q_M-FACTOR(L, I)*Q_P
!           UPWARD EQUATION:
            P(L, 5, 2*I+1)=U_P_0P-FACTOR(L, I)*U_M_0P
            P(L, 4, 2*I+1)=U_P_MM-FACTOR(L, I)*U_M_MM
            P(L, 3, 2*I+1)=U_P_PP-FACTOR(L, I)*U_M_PP
            P(L, 2, 2*I+1)=U_P_0M-FACTOR(L, I)*U_M_0M
            P(L, 1, 2*I+1)=U_P_2P-FACTOR(L, I)*U_M_2P
            V(L, 2*I+1)=Q_P-FACTOR(L, I)*Q_M
!
         ENDDO
!
!        CLOUDLESS PROFILES:
         DO LL=1, N_FREE_PROFILE(I)
            L=I_FREE_PROFILE(LL, I)
!
!           ASSIGN THE COEFFICIENTS OF THE PENTADIAGONAL MATRIX:
!           UPWARD EQUATION:
            P(L, 5, 2*I)=0.0D+00
            P(L, 4, 2*I)=1.0D+00
            P(L, 3, 2*I)=-R_FREE(L, I)
            P(L, 2, 2*I)=-T_FREE(L, I)
            P(L, 1, 2*I)=0.0D+00
            V(L, 2*I)=S_UP_FREE(L, I)
!
!           DOWNWARD EQUATION:
            P(L, 5, 2*I+1)=0.0D+00
            P(L, 4, 2*I+1)=-T_FREE(L, I)
            P(L, 3, 2*I+1)=-R_FREE(L, I)
            P(L, 2, 2*I+1)=1.0D+00
            P(L, 1, 2*I+1)=0.0D+00
            V(L, 2*I+1)=S_DOWN_FREE(L, I)
         ENDDO
!
      ENDDO
!
!     SET THE EQUATIONS AT THE BOUNDARIES.
      DO L=1, N_PROFILE
!        TOP BOUNDARY CONDITION:
         P(L, 5, 1)=0.0D+00
         P(L, 4, 1)=0.0D+00
         P(L, 3, 1)=0.0D+00
         P(L, 2, 1)=1.0D+00
         P(L, 1, 1)=0.0D+00
         V(L, 1)=FLUX_INC_DOWN(L)
!        BOTTOM BOUNDARY CONDITION:
!        A PERTURBATION IS REQUIRED IN CASES WHERE THE ALBEDO IS 0.
         P(L, 5, 2*N_LAYER+2)=0.0D+00
         P(L, 4, 2*N_LAYER+2)=-1.0D+00
         P(L, 3, 2*N_LAYER+2)=ALBEDO_SURFACE_DIFF(L)
     &       +TOL_MACHINE/(ALBEDO_SURFACE_DIFF(L)+SQRT_TOL_MACHINE)
         P(L, 2, 2*N_LAYER+2)=0.0D+00
         P(L, 1, 2*N_LAYER+2)=0.0D+00
         V(L, 2*N_LAYER+2)
     &      =-(ALBEDO_SURFACE_DIR(L)-ALBEDO_SURFACE_DIFF(L))
     &      *FLUX_DIRECT_GROUND(L)
     &      -SOURCE_GROUND(L)
      ENDDO
!
!
      RETURN
      END
