subroutine conductionQ(nz,z,dt,Qn,Qnp1,T,ti,rhoc,emiss,Tsurf,
     &     Fgeotherm,Fsurf)
C***********************************************************************
C   conductionQ:  program to calculate the diffusion of temperature 
C                 into the ground and thermal emission at the surface 
C                 with variable thermal properties on irregular grid
C   Crank-Nicholson scheme, flux conservative
C                           uses Samar's radiation formula
C   Eqn: rhoc*T_t = (k*T_z)_z 
C   BC (z=0): Q(t) + kT_z = em*sig*T^4
C   BC (z=L): heat flux = Fgeotherm
C
C   ti = thermal inertia [J m^-2 K^-1 s^-1/2]  VECTOR
C   rhoc = rho*c  VECTOR where rho=density [kg m^-3] and 
C                              c=specific heat [J K^-1 kg^-1]
C   ti and rhoc are not allowed to vary in the layers immediately adjacent
C               to the surface or the bottom
C   T = vertical temperature profile [K] (output)
C   Qn,Qnp1 = net solar insolation at time steps n and n+1 [Watts m^-2]
C   emiss = emissivity
C   Tsurf = surface Temperature [K]  (output)
C
C   Grid: surface is at z=0
C         T(1) is at z(1); ...; T(i) is at z(i)
C         k(i) is midway between z(i-1) and z(i); same for rhoc(i)
C
C   orginally written by Samar Khatiwala, 2001
C   extended to variable thermal properties
C         and irregular grid by Norbert Schorghofer
C***********************************************************************

      implicit none
      integer NMAX
      real*8 sigSB
      parameter (NMAX=2000)
      parameter (sigSB=5.67051d-8)

      integer nz, i
      real*8 z(NMAX), dt, Qn, Qnp1, T(NMAX), ti(NMAX),rhoc(NMAX)
      real*8 emiss, Tsurf, Fgeotherm, Fsurf
      real*8 k(NMAX), k1, alpha(NMAX), gamma(NMAX), Tr
      real*8 a(NMAX), b(NMAX), c(NMAX), r(NMAX)
      real*8 arad, brad, ann, annp1, bn, buf, dz, beta


C     set some constants
      do i=1,nz
         k(i) = ti(i)**2/rhoc(i) ! thermal conductivity
      enddo
      dz=2.*z(1)
      beta = dt/rhoc(1)/(2.*dz**2)   ! assumes rhoc(0)=rhoc(1)
      alpha(1) = beta*k(2)
      gamma(1) = beta*k(1)
      do i=2,nz-1
         buf = dt/(z(i+1)-z(i-1))
         alpha(i) = 2.*k(i+1)*buf/(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i))
         gamma(i) = 2.*k(i)*buf/(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1))
      enddo
      buf=dt/(z(nz)-z(nz-1))**2
      gamma(nz) = k(nz)*buf/(2*rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz)

      k1=k(1)/dz

C     elements of tridiagonal matrix
      do i=1,nz
         a(i) = -gamma(i)   !  a(1) is not used
         b(i) = 1. + alpha(i) + gamma(i) !  b(1) has to be reset at every timestep
         c(i) = -alpha(i)   !  c(nz) is not used
      enddo
      b(nz) = 1. + gamma(nz)
      
C     Emission
      Tr = Tsurf                !   'reference' temperature
      arad = -3.*emiss*sigSB*Tr**4
      brad = 2.*emiss*sigSB*Tr**3
      ann = (Qn-arad)/(k1+brad)
      annp1 = (Qnp1-arad)/(k1+brad)
      bn = (k1-brad)/(k1+brad)
      b(1) = 1. + alpha(1) + gamma(1) - gamma(1)*bn

C     Set RHS         
      r(1) = gamma(1)*(annp1+ann) + 
     &     (1.-alpha(1)-gamma(1)+gamma(1)*bn)*T(1) + alpha(1)*T(2)
      do i=2,nz-1
         r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i)
     &        + alpha(i)*T(i+1)
      enddo
      r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz)
     &     + dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz)

C     Solve for T at n+1
      call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion
      
      Tsurf = 0.5*(annp1 + bn*T(1) + T(1)) ! (T0+T1)/2
      Fsurf = - k(1)*(T(1)-Tsurf)/z(1) ! heat flux into surface
      end