!+ Subroutine to set the matrix for net fluxes in Schemtz's scheme.
!
! 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 S_SET_MATRIX_NET(N_PROFILE, N_LAYER
     &   , TRANS_FREE, REFLECT_FREE, TRANS_CLOUD, REFLECT_CLOUD
     &   , W_CLOUD, N_CLOUD_TOP
     &   , S_DOWN_FREE, S_UP_FREE, S_DOWN_CLOUD, S_UP_CLOUD
     &   , SOURCE_GROUND, ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR
     &   , FLUX_INC_DOWN, FLUX_DIRECT_GROUND
     &   , A3, B
     &   )
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/spectral_region.pcf'
!
!     DECLARATION OF DUMMY ARGUMENTS.
      INTEGER 	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOP CLOUDY LAYER
      REAL*8 	!, INTENT(IN)
     &     TRANS_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE TRANSMISSION
     &   , REFLECT_FREE(NPD_PROFILE, NPD_LAYER)
!             FREE REFLECTION
     &   , TRANS_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY TRANSMISSION
     &   , REFLECT_CLOUD (NPD_PROFILE, NPD_LAYER)
!             CLOUDY REFLECTION
     &   , W_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUD AMOUNT
     &   , S_UP_FREE(NPD_PROFILE, NPD_LAYER)
!             UPWARD SOURCE FUNCTION
     &   , S_DOWN_FREE(NPD_PROFILE, NPD_LAYER)
!             DOWNWARD SOURCE FUNCTION
     &   , S_UP_CLOUD(NPD_PROFILE, NPD_LAYER)
!             UPWARD SOURCE FUNCTION
     &   , S_DOWN_CLOUD(NPD_PROFILE, NPD_LAYER)
!             DOWNWARD SOURCE FUNCTION
      REAL*8	!, INTENT(IN)
     &     SOURCE_GROUND(NPD_PROFILE)
!             GROUND SOURCE FUNCTION
     &   , ALBEDO_SURFACE_DIFF(NPD_PROFILE)
!             DIFFUSE ALBEDO OF SURFACE
     &   , ALBEDO_SURFACE_DIR(NPD_PROFILE)
!             DIRECT ALBEDO OF SURFACE
     &   , FLUX_INC_DOWN(NPD_PROFILE)
!             INCIDENT DIFUSE FLUX
     &   , FLUX_DIRECT_GROUND(NPD_PROFILE)
!             DIRECT FLUX AT GROUND
      REAL*8	!, INTENT(OUT)
     &     A3(NPD_PROFILE, 3, NPD_LAYER+1)
!             TRIDIAGONAL MATRIX
     &   , B(NPD_PROFILE, NPD_LAYER+1)
!             RIGHT-HAND SIDES
!
!     DECLARATION OF LOCAL VARIABLES.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , L
!             LOOP VARIABLE
      REAL*8
     &     P
!             WORKING ARRAY
     &   , Q
!             WORKING ARRAY
     &   , D
!             WORKING ARRAY
     &   , TE(NPD_PROFILE, NPD_LAYER)
!             EFFECTIVE TRANSMISSION
     &   , RP(NPD_PROFILE, NPD_LAYER)
!             EFFECTIVE UPWARD REFLECTION
     &   , RM(NPD_PROFILE, NPD_LAYER)
!             EFFECTIVE DOWNWARD REFLECTION
     &   , SP(NPD_PROFILE, NPD_LAYER)
!             EFFECTIVE UPWARD SOURCE FUNCTION
     &   , SM(NPD_PROFILE, NPD_LAYER)
!             EFFECTIVE DOWNWARD SOURCE FUNCTION
     &   , QP(NPD_PROFILE, NPD_LAYER)
!             ARRAY Q-PRIMED +
     &   , QM(NPD_PROFILE, NPD_LAYER)
!             ARRAY Q-PRIMED -
     &   , PE(NPD_PROFILE, NPD_LAYER)
!             ARRAY P-PRIMED
!
!
!     PRECALCULATION OF FIELDS:
!     TOP CLEAR REGION:
      DO I=1, N_CLOUD_TOP-1
         DO L=1, N_PROFILE
            QP(L, I)=TRANS_FREE(L, I)**2-REFLECT_FREE(L, I)**2+1.0D+00
            PE(L, I)=(REFLECT_FREE(L, I)-1.0D+00)**2-TRANS_FREE(L, I)**2
         ENDDO
      ENDDO
!
!     CLEAR AND CLOUDY REGION:
      DO I=N_CLOUD_TOP, N_LAYER
         DO L=1, N_PROFILE
            P=W_CLOUD(L, I)*REFLECT_CLOUD(L, I)
            Q=1.0D+00-W_CLOUD(L, I)*(1.0D+00-TRANS_CLOUD(L, I))
            D=1.0D+00-P*REFLECT_FREE(L, I)
            TE(L, I)=Q*TRANS_FREE(L, I)/(1.0D+00-P*REFLECT_FREE(L, I))
            RP(L, I)=(REFLECT_FREE(L, I)*(Q**2-P**2)+P)
     &         /(1.0D+00-P*REFLECT_FREE(L, I))
            RM(L, I)=(P*(TRANS_FREE(L, I)**2-REFLECT_FREE(L, I)**2)
     &         +REFLECT_FREE(L, I))
     &         /(1.0D+00-P*REFLECT_FREE(L, I))
            SP(L, I)=Q*((REFLECT_FREE(L, I)*S_DOWN_CLOUD(L, I)
     &         +S_UP_FREE(L, I))/(1.0D+00-P*REFLECT_FREE(L, I)))
     &         +S_UP_CLOUD(L, I)
            SM(L, I)=TRANS_FREE(L, I)*((P*S_UP_FREE(L, I)
     &         +S_DOWN_CLOUD(L, I))/(1.0D+00-P*REFLECT_FREE(L, I)))
     &         +S_DOWN_FREE(L, I)
            QM(L, I)=TE(L, I)**2-(RM(L, I)-1.0D+00)*(RP(L, I)+1.0D+00)
            QP(L, I)=TE(L, I)**2-(RM(L, I)+1.0D+00)*(RP(L, I)-1.0D+00)
            PE(L, I)=(RM(L, I)-1.0D+00)*(RP(L, I)-1.0D+00)-TE(L, I)**2
         ENDDO
      ENDDO
!
!     CODE THE EQUATIONS INTO THE MATRIX:
!     UPPER SURFACE:
      DO L=1, N_PROFILE
         A3(L, 2, 1)=0.5D+00*(PE(L, 1)+QM(L, 1))
         A3(L, 3, 1)=-TE(L, 1)
         B(L, 1)=PE(L, 1)*FLUX_INC_DOWN(L)-TE(L, 1)*SM(L, 1)
     &      +(RM(L, 1)-1.0D+00)*SP(L, 1)
      ENDDO
!     THE INTERIOR EQUATIONS:
      DO I=1, N_CLOUD_TOP-1
         DO L=1, N_PROFILE
            A3(L, 1, I+1)=TRANS_FREE(L, I)*PE(L, I)
            A3(L, 2, I+1)=-0.5D+00*(QP(L, I+1)*PE(L, I)
     &         +PE(L, I+1)*QP(L, I))
            A3(L, 3, I+1)=TRANS_FREE(L, I)*PE(L, I)
            B(L, I+1)=PE(L, I)*(TRANS_FREE(L, I+1)*S_DOWN_FREE(L, I)
     &         -(REFLECT_FREE(L, I+1)-1.0D+00)*S_UP_FREE(L, I+1))
     &         +PE(L, I+1)*((REFLECT_FREE(L, I)-1.0D+00)
     &         *S_DOWN_FREE(L, I)-TRANS_FREE(L, I)*S_UP_FREE(L, I))
         ENDDO
      ENDDO
      DO I=N_CLOUD_TOP, N_LAYER-1
         DO L=1, N_PROFILE
            A3(L, 1, I+1)=TE(L, I)*PE(L, I+1)
            A3(L, 2, I+1)=-0.5D+00
     &         *(QM(L, I+1)*PE(L, I)+PE(L, I+1)*QP(L, I))
            A3(L, 3, I+1)=TE(L, I+1)*PE(L, I)
            B(L, I+1)=PE(L, I)*(TE(L, I+1)*SM(L, I+1)
     &         -(RM(L, I+1)-1.0D+00)*SP(L, I+1))
     &         +PE(L, I+1)*((RP(L, I)-1.0D+00)*SM(L, I)
     &         -TE(L, I)*SP(L, I))
         ENDDO
      ENDDO
!     LOWER BOUNDARY
      DO L=1, N_PROFILE
         A3(L, 1, N_LAYER+1)=(1.0D+00-ALBEDO_SURFACE_DIFF(L))
     &      *TE(L, N_LAYER)
         A3(L, 2, N_LAYER+1)=-0.5D+00*((1.0D+00
     &      -ALBEDO_SURFACE_DIFF(L))*QM(L, N_LAYER)
     &      +(1.0D+00+ALBEDO_SURFACE_DIFF(L))*PE(L, N_LAYER))
         B(L, N_LAYER+1)=PE(L, N_LAYER)*(ALBEDO_SURFACE_DIR(L)
     &      *FLUX_DIRECT_GROUND(L)+SOURCE_GROUND(L))
     &      +(1.0D+00-ALBEDO_SURFACE_DIFF(L))
     &      *((RP(L, N_LAYER)-1.0D+00)*SM(L, N_LAYER)
     &      -TE(L, N_LAYER)*SP(L, N_LAYER))
      ENDDO
!
!
      RETURN
      END
