! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/common/gaussian.F
      SUBROUTINE GAUSSIAN (NLEVS, A, B, C, D, XMIN, XMAX, X)

      IMPLICIT NONE

!-----------------------------------------------------------------------
! Description:
!     Solves a tridiagnonal matrix equation of the form:
!             A(n) X(n-1) + B(n) X(n) + C(n) X(n+1) = D(n)
!     by Gausian elimination, assuming boundary conditions:
!             A(1) = 0.0    at the top
!             C(N) = 0.0    at the bottom.
!-----------------------------------------------------------------------

      INTEGER NLEVS    ! IN Number of levels for Gaussian elimination
      INTEGER N        ! WORK Loop counters.

      REAL A(NLEVS)    ! IN Matrix elements corresponding to X(n-1).
      REAL B(NLEVS)    ! IN Matrix elements corresponding to X(n).
      REAL C(NLEVS)    ! IN Matrix elements corresponding to X(n+1).
      REAL D(NLEVS)    ! IN Matrix elements corresponding to the RHS
      REAL XMIN        ! IN Minimum permitted value of X.
      REAL XMAX        ! IN Maximum permitted value of X.
      REAL X(NLEVS)    ! OUT Solution.
!     WORK Transformed matrix elements
      REAL ADASH(NLEVS), BDASH(NLEVS), CDASH(NLEVS), DDASH(NLEVS)

!-----------------------------------------------------------------------
! By default set the implicit increment to the explicit increment
! (for when denominators vanish).
!-----------------------------------------------------------------------
      DO N=1,NLEVS
        X(N) = D(N)
      ENDDO

!-----------------------------------------------------------------------
! Upward Sweep: eliminate "C" elements by replacing nth equation with:
!                  B'(n+1)*Eq(n)-C(n)*Eq'(n+1)
! where "'" denotes a previously tranformed equation. The resulting
! equations take the form:
!                A'(n) X(n-1) + B'(n) X(n) = D'(n)
! (NB. The bottom boundary condition implies that the NLEV equation does
!  not need transforming.)
!-----------------------------------------------------------------------
      ADASH(NLEVS) = A(NLEVS)
      BDASH(NLEVS) = B(NLEVS)
      DDASH(NLEVS) = D(NLEVS)

      DO N=NLEVS-1,1,-1
        ADASH(N) = BDASH(N+1)*A(N)
        BDASH(N) = BDASH(N+1)*B(N)-C(N)*ADASH(N+1)
        DDASH(N) = BDASH(N+1)*D(N)-C(N)*DDASH(N+1)
      ENDDO

!-----------------------------------------------------------------------
! Top boundary condition: A(1) = 0.0, allows X(1) to be diagnosed
!-----------------------------------------------------------------------
      IF (BDASH(1) .NE. 0.) X(1) = DDASH(1)/BDASH(1)
      X(1) = MAX(X(1),XMIN)
      X(1) = MIN(X(1),XMAX)

!-----------------------------------------------------------------------
! Downward Sweep: calculate X(n) from X(n-1):
!                 X(n) = (D'(n) - A'(n) X(n-1)) / B'(n)
!-----------------------------------------------------------------------
      DO N=2,NLEVS
        IF (BDASH(N) .NE. 0.) THEN
          X(N) = (DDASH(N)-ADASH(N)*X(N-1))/BDASH(N)
        ENDIF
        X(N) = MAX(X(N),XMIN)
        X(N) = MIN(X(N),XMAX)
      ENDDO

      RETURN
      END