! source file: /Users/nmengis/Documents/UVic_ESCM/2.10/updates/02_CE_permafrost_merge_CMIP6forcing/source/mtlm/soil_advect.F
      subroutine SOIL_ADVECT (NGND, NSOIL, HCAP_W, HCON_ROCK, DZE
     &,                       HC_B, TG, W_FL, A_FLUX, DF1_DT1
     &,                       DF1_DT2, DF2_DT1, DF2_DT2)

      integer NGND, NSOIL, N

      real DZE(NGND), HC_B(NGND), W_FL(0:NGND), HCON_ROCK, HCAP_W
      real A_FLUX(NGND), HCAP(NGND), TG(NGND), DF1_DT1(NGND)
      real DF1_DT2(NGND), DF2_DT1(NGND), DF2_DT2(NGND)

      A_FLUX(:) = 0.

!---------------------------------------------------------------------
! Surface soil layer (heat flux through top set by penmon)
!---------------------------------------------------------------------
      DF1_DT1(1) = 0.
      DF1_DT2(1) = 0.
      DF2_DT1(1) = 2.*HC_B(1)/(DZE(1) + DZE(2))
      DF2_DT2(1) = - DF2_DT1(1)
      if (W_FL(1) .lt. 0.)
     &  A_FLUX(1) = -HCAP_W*W_FL(1)*(TG(2) - TG(1))

!--------------------------------------------------------------------
! Middle soil layers
!--------------------------------------------------------------------
      do N=2,NSOIL-1
       DF1_DT1(N) = 2.*HC_B(N-1)/(DZE(N-1) + DZE(N))
       DF1_DT2(N) = - DF1_DT1(N)
       DF2_DT1(N) = 2.*HC_B(N)/(DZE(N) + DZE(N+1))
       DF2_DT2(N) = -DF2_DT1(N)
       if (W_FL(N-1) .gt. 0.)
     &   A_FLUX(N) = A_FLUX(N) + HCAP_W*W_FL(N-1)*(TG(N-1) - TG(N))
       if (W_FL(N) .lt. 0.)
     &   A_FLUX(N) = A_FLUX(N) - HCAP_W*W_FL(N)*(TG(N+1) - TG(N))
      enddo

!--------------------------------------------------------------------
! Bottom soil layer
!--------------------------------------------------------------------
      DF1_DT1(NSOIL) = 2*HC_B(NSOIL-1)/(DZE(NSOIL-1)+DZE(NSOIL))
      DF1_DT2(NSOIL) = - DF1_DT1(NSOIL)
      DF2_DT1(NSOIL) = 2*HC_B(NSOIL)/(DZE(NSOIL) + DZE(NSOIL+1))
      DF2_DT2(NSOIL) = -DF2_DT1(NSOIL)
      if (W_FL(NSOIL-1) .gt. 0.)
     &  A_FLUX(NSOIL) = HCAP_W*W_FL(NSOIL-1)*(TG(NSOIL-1) - TG(NSOIL))

!---------------------------------------------------------------------
! Bedrock - top layer
!---------------------------------------------------------------------
      N = NSOIL+1
!     This is already adjusted for bedrock / slab boundary
      DF1_DT1(N) = 2*HC_B(N-1)/(DZE(N-1) + DZE(N))
      DF1_DT2(N) = -DF1_DT1(N)
      DF2_DT1(N) = 2*HCON_ROCK/(DZE(N) + DZE(N+1))
      DF2_DT2(N) = -DF2_DT1(N)

!---------------------------------------------------------------------
! Bedrock - middle layers
!---------------------------------------------------------------------
      do N=NSOIL+2,NGND-1
        DF1_DT1(N) = 2*HCON_ROCK/(DZE(N-1) + DZE(N))
        DF1_DT2(N) = -DF1_DT1(N)
        DF2_DT1(N) = 2*HCON_ROCK/(DZE(N) + DZE(N+1))
        DF2_DT2(N) = -DF2_DT1(N)
      enddo

!---------------------------------------------------------------------
! Bedrock - bottom layer
!---------------------------------------------------------------------
      DF1_DT1(NGND) = 2*HCON_ROCK/(DZE(NGND-1) + DZE(NGND))
      DF1_DT2(NGND) = -DF1_DT1(N)
      DF2_DT1(NGND) = 0.  ! Geothermal heat flux is not T dependent
      DF2_DT2(NGND) = 0.

      return
      end