| 1 | subroutine conductionT(nz,z,dt,T,Tsurf,Tsurfp1,ti,rhoc,Fgeotherm,Fsurf) |
|---|
| 2 | !*********************************************************************** |
|---|
| 3 | ! conductionT: program to calculate the diffusion of temperature |
|---|
| 4 | ! into the ground with prescribed surface temperature |
|---|
| 5 | ! and variable thermal properties on irregular grid |
|---|
| 6 | ! Crank-Nicholson scheme, flux conservative |
|---|
| 7 | ! |
|---|
| 8 | ! Eqn: rhoc*T_t = (k*T_z)_z |
|---|
| 9 | ! BC (z=0): T=T(t) |
|---|
| 10 | ! BC (z=L): heat flux = Fgeotherm |
|---|
| 11 | ! |
|---|
| 12 | ! nz = number of grid points |
|---|
| 13 | ! dt = time step |
|---|
| 14 | ! T = vertical temperature profile [K] (in- and output) |
|---|
| 15 | ! Tsurf, Tsurfp1 = surface temperatures at times n and n+1 |
|---|
| 16 | ! ti = thermal inertia [J m^-2 K^-1 s^-1/2] VECTOR |
|---|
| 17 | ! rhoc = rho*c heat capacity per volume [J m^-3 K^-1] VECTOR |
|---|
| 18 | ! ti and rhoc are not allowed to vary in the layers immediately |
|---|
| 19 | ! adjacent to the surface or the bottom |
|---|
| 20 | ! Fgeotherm = geothermal heat flux at bottom boundary [W/m^2] |
|---|
| 21 | ! Fsurf = heat flux at surface [W/m^2] (output) |
|---|
| 22 | ! |
|---|
| 23 | ! Grid: surface is at z=0 |
|---|
| 24 | ! T(1) is at z(1); ...; T(i) is at z(i) |
|---|
| 25 | ! k(i) is midway between z(i-1) and z(i) |
|---|
| 26 | ! rhoc(i) is midway between z(i-1) and z(i) |
|---|
| 27 | !*********************************************************************** |
|---|
| 28 | implicit none |
|---|
| 29 | |
|---|
| 30 | integer, intent(IN) :: nz |
|---|
| 31 | real*8, intent(IN) :: z(nz), dt, Tsurf, Tsurfp1, ti(nz), rhoc(nz) |
|---|
| 32 | real*8, intent(IN) :: Fgeotherm |
|---|
| 33 | real*8, intent(INOUT) :: T(nz) |
|---|
| 34 | real*8, intent(OUT) :: Fsurf |
|---|
| 35 | integer i |
|---|
| 36 | real*8 alpha(nz), k(nz), gamma(nz), buf |
|---|
| 37 | real*8 a(nz), b(nz), c(nz), r(nz) |
|---|
| 38 | |
|---|
| 39 | ! set some constants |
|---|
| 40 | k(:) = ti(:)**2/rhoc(:) ! thermal conductivity |
|---|
| 41 | alpha(1) = k(2)*dt/rhoc(1)/(z(2)-z(1))/z(2) |
|---|
| 42 | gamma(1) = k(1)*dt/rhoc(1)/z(1)/z(2) |
|---|
| 43 | do i=2,nz-1 |
|---|
| 44 | buf = dt/(z(i+1)-z(i-1)) |
|---|
| 45 | alpha(i) = k(i+1)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i+1)-z(i)) |
|---|
| 46 | gamma(i) = k(i)*buf*2./(rhoc(i)+rhoc(i+1))/(z(i)-z(i-1)) |
|---|
| 47 | enddo |
|---|
| 48 | buf = dt/(z(nz)-z(nz-1))**2 |
|---|
| 49 | gamma(nz) = k(nz)*buf/(rhoc(nz)+rhoc(nz)) ! assumes rhoc(nz+1)=rhoc(nz) |
|---|
| 50 | |
|---|
| 51 | ! elements of tridiagonal matrix |
|---|
| 52 | a(:) = -gamma(:) ! a(1) is not used |
|---|
| 53 | b(:) = 1. + alpha(:) + gamma(:) |
|---|
| 54 | c(:) = -alpha(:) ! c(nz) is not used |
|---|
| 55 | b(nz) = 1. + gamma(nz) |
|---|
| 56 | |
|---|
| 57 | ! Set RHS |
|---|
| 58 | r(1)= alpha(1)*T(2) + (1.-alpha(1)-gamma(1))*T(1) + gamma(1)*(Tsurf+Tsurfp1) |
|---|
| 59 | do concurrent (i=2:nz-1) |
|---|
| 60 | r(i) = gamma(i)*T(i-1) + (1.-alpha(i)-gamma(i))*T(i) + alpha(i)*T(i+1) |
|---|
| 61 | enddo |
|---|
| 62 | r(nz) = gamma(nz)*T(nz-1) + (1.-gamma(nz))*T(nz) + & |
|---|
| 63 | & dt/rhoc(nz)*Fgeotherm/(z(nz)-z(nz-1)) ! assumes rhoc(nz+1)=rhoc(nz) |
|---|
| 64 | |
|---|
| 65 | ! Solve for T at n+1 |
|---|
| 66 | call tridag(a,b,c,r,T,nz) ! update by tridiagonal inversion |
|---|
| 67 | |
|---|
| 68 | Fsurf = -k(1)*(T(1)-Tsurfp1)/z(1) ! heat flux into surface |
|---|
| 69 | |
|---|
| 70 | end subroutine conductionT |
|---|
