subroutine conductionT(nz,z,dt,T,Tsurf,Tsurfp1,ti,rhoc, & Fgeotherm,Fsurf) C*********************************************************************** C conductionT: program to calculate the diffusion of temperature C into the ground with prescribed surface temperature C and variable thermal properties on irregular grid C Crank-Nicholson scheme, flux conservative C C Eqn: rhoc*T_t = (k*T_z)_z C BC (z=0): T=T(t) 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 heat capacity per volume [J m^-3] VECTOR 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] C Tsurf, Tsurfp1 = surface temperatures at times n and n+1 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) C rhoc(i) is midway between z(i-1) and z(i) C*********************************************************************** implicit none integer NMAX parameter (NMAX=2000) integer nz, i real*8 z(NMAX), dt, T(NMAX), Tsurf, Tsurfp1, ti(NMAX), rhoc(NMAX) real*8 alpha(NMAX), k(NMAX), gamma(NMAX), buf, Fgeotherm, Fsurf real*8 a(NMAX), b(NMAX), c(NMAX), r(NMAX) C set some constants do i=1,nz k(i) = ti(i)**2/rhoc(i) ! thermal conductivity enddo alpha(1) = k(2)*dt/rhoc(1)/(z(2)-z(1))/z(2) gamma(1) = k(1)*dt/rhoc(1)/z(1)/z(2) do i=2,nz-1 buf=dt/(z(i+1)-z(i-1)) alpha(i) = k(i+1)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) gamma(i) = k(i)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) enddo buf=dt/(z(nz)-z(nz-1))**2 gamma(nz) = k(nz)*buf/(rhoc(nz)+rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) 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) c(i) = -alpha(i) ! c(nz) is not used enddo b(nz) = 1. + gamma(nz) C Set RHS r(1)= alpha(1)*T(2) + (1.-alpha(1)-gamma(1))*T(1) & + gamma(1)*(Tsurf+Tsurfp1) 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 Fsurf = -k(1)*(T(1)-Tsurfp1)/z(1) ! heat flux into surface end