!+ Subroutine to set the matrix equations for the net fluxes.
!
! Method:
!	Straightforward.
!
! 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 SET_MATRIX_NET(N_PROFILE, N_LAYER
     &   , TRANS, REFLECT
     &   , S_DOWN, S_UP
     &   , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR
     &   , FLUX_DIRECT_GROUND, FLUX_INC_DOWN
     &   , SOURCE_GROUND
     &   , A3, B
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/prec_machine.ucf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
      REAL*8	!, INTENT(IN)
     &     TRANS(NPD_PROFILE, NPD_LAYER)
!             TRANSMISSION COEFFICIENTS
     &   , REFLECT(NPD_PROFILE, NPD_LAYER)
!             REFLECT COEFFICIENTS
     &   , S_UP(NPD_PROFILE, NPD_LAYER)
!             UPWARD SOURCE FUNCTION
     &   , S_DOWN(NPD_PROFILE, NPD_LAYER)
!             DOWNWARD SOURCE FUNCTION
     &   , ALBEDO_SURFACE_DIFF(NPD_PROFILE)
!             DIFFUSE REFLECTION COEFFICIENTS
     &   , ALBEDO_SURFACE_DIR(NPD_PROFILE)
!             DIRECT REFLECTION COEFFICIENTS
     &   , FLUX_DIRECT_GROUND(NPD_PROFILE)
!             DIRECT FLUX AT SURFACE
     &   , FLUX_INC_DOWN(NPD_PROFILE)
!             INCIDENT DOWNWARD FLUX
     &   , SOURCE_GROUND(NPD_PROFILE)
!             SOURCE FROM GROUND
      REAL*8	!, INTENT(OUT)
     &     A3(NPD_PROFILE, 3, 2*NPD_LAYER+2)
!             TRIDIAGONAL MATRIX
     &   , B(NPD_PROFILE, 2*NPD_LAYER+2)
!             RHS OF EQUATIONS
!
!     LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
      REAL*8
     &     SIGMA(NPD_PROFILE, NPD_LAYER)
!             ARRAY SIGMA
     &   , DELTA(NPD_PROFILE, NPD_LAYER)
!             ARRAY DELTA
     &   , ZETA_TOP(NPD_PROFILE, NPD_LAYER)
!             2.0*ARRAY ZETA AT TOP OF LAYER
     &   , ZETA_BASE(NPD_PROFILE, NPD_LAYER)
!             2.0*ARRAY ZETA AT BASE OF LAYER
!
!
!     THE ARRAYS P AND Q MUST BE CALCULATED SEPARATELY SINCE THEY WILL
!     BE COUPLED BETWEEN LAYERS.
      DO I=1, N_LAYER
         DO L=1, N_PROFILE
            SIGMA(L, I)=2.0D+00/(1.0D+00-(TRANS(L, I)+REFLECT(L, I))
     &         *(1.0D+00-SQRT_TOL_MACHINE))-1.0D+00
            DELTA(L, I)=2.0D+00/(1.0D+00+(TRANS(L, I)-REFLECT(L, I))
     &         *(1.0D+00-SQRT_TOL_MACHINE))-1.0D+00
            ZETA_BASE(L, I)=SIGMA(L, I)*(S_DOWN(L, I)+S_UP(L, I))
     &         +DELTA(L, I)*(S_DOWN(L, I)-S_UP(L, I))
     &         +2.0D+00*S_DOWN(L, I)
            ZETA_TOP(L, I)=SIGMA(L, I)*(S_DOWN(L, I)+S_UP(L, I))
     &         -DELTA(L, I)*(S_DOWN(L, I)-S_UP(L, I))+2.0D+00*S_UP(L, I)
         ENDDO
      ENDDO
!
!     CODE THE EQUATIONS INTO THE MATRIX:
!     UPPER SURFACE:
      DO L=1, N_PROFILE
         A3(L, 2, 1)=-2.0D+00-SIGMA(L, 1)-DELTA(L, 1)
         A3(L, 3, 1)=SIGMA(L, 1)-DELTA(L, 1)
         B(L, 1)=ZETA_TOP(L, 1)-4.0D+00*FLUX_INC_DOWN(L)
      ENDDO
!     THE INTERIOR EQUATIONS:
      DO I=1, N_LAYER-1
         DO L=1, N_PROFILE
            A3(L, 1, I+1)=SIGMA(L, I)-DELTA(L, I)
            A3(L, 2, I+1)=-SIGMA(L, I)-DELTA(L, I)
     &          -SIGMA(L, I+1)-DELTA(L, I+1)
            A3(L, 3, I+1)=SIGMA(L, I+1)-DELTA(L, I+1)
            B(L, I+1)=ZETA_TOP(L, I+1)-ZETA_BASE(L, I)
         ENDDO
      ENDDO
!     LOWER BOUNDARY
      DO L=1, N_PROFILE
         A3(L, 1, N_LAYER+1)=(1.0D+00-ALBEDO_SURFACE_DIFF(L))
     &      *(SIGMA(L, N_LAYER)-DELTA(L, N_LAYER))
         A3(L, 2, N_LAYER+1)=(ALBEDO_SURFACE_DIFF(L)-1.0D+00)
     &      *(SIGMA(L, N_LAYER)+DELTA(L, N_LAYER))
     &      -2.0D+00*(1.0D+00+ALBEDO_SURFACE_DIFF(L))
         B(L, N_LAYER+1)
     &      =4.0D+00*((ALBEDO_SURFACE_DIR(L)-ALBEDO_SURFACE_DIFF(L))
     &      *FLUX_DIRECT_GROUND(L)+SOURCE_GROUND(L))
     &      -(1.0D+00-ALBEDO_SURFACE_DIFF(L))*ZETA_BASE(L, N_LAYER)
      ENDDO
!
!
      RETURN
      END
