!!!!!!!!!!!!!!! conduction.f90 for Fortran 90 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Purpose:      Three methods to simulate conduction through a multi-layered material
!               1. Half-layer scheme (standard method)
!               2. Modified half-layer scheme (standard method with additional surface node)
!               3. Interface scheme  method (proposed method)
!               4. Hgh resolution (Standard method witih 200 layers at 1 second timestep)
!               More information at https://www.theurbanist.com.au/2017/06/solving-the-heat-equation/
! Developer:    Mathew Lipson <m.lipson@unsw.edu.au>
! Revision:     27 Nov 2017
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program conduction
implicit none

integer, parameter :: dp = SELECTED_REAL_KIND(8)       ! double precision for energy closure check

type materialdata
  real(kind=dp), dimension(:,:), allocatable :: depth, volcp, lambda, nodetemp, nodetemp_prev
end type materialdata

real(kind=dp), dimension(:), allocatable :: temp_ext    ! external sol-air temperature (varying) [K]
real(kind=dp), dimension(:), allocatable :: temp_int    ! internal sol-air temperature (fixed )  [K]
real(kind=dp), dimension(:), allocatable :: storageflux ! total storage flux between timesteps   

real(kind=dp), parameter :: pi = 4 * atan (1.0_8) ! pi
real(kind=dp), parameter :: spinup = 10.          ! number of periods spin up to be discarded
integer, parameter :: ufull = 1                   ! number of tiles (for vectorisation)

real(kind=dp) :: ddt=1800.                        ! timestep length                         [s]
real(kind=dp) :: period=86400.                    ! sinusoidal period                       [s]
real(kind=dp) :: res_ext=0.04                     ! external surface resistance             [K m^2/W]  
real(kind=dp) :: res_int=0.13                     ! internal surface resistance             [K m^2/W]
integer :: nl = 4                                 ! number of layers in material            [-]
integer :: timesteps                              ! total timesteps in simulation           [-]
integer :: t                                      ! time loop integer
character(len=15) :: filename

integer, save :: conductmeth                      ! conduction method (0=half-layer, 1=interface)
type(materialdata), save :: material

print *, "Select conduction method (0=standard half-layer, 1=modified half-layer, 2=interface, 3=high-res): "
read *,conductmeth
write(filename, '(a,I0,a)') "output_",conductmeth,'.txt'
print *, "Writing to ", filename
open(1,file=filename,action="write",status="replace")
write(1,*) "timestep   temp_ext    storageflux nodetemp(0) nodetemp(1)"

! high-res standard half-layer run (exact)
if (conductmeth==3) then
  conductmeth=0
  nl = 200
  ddt = 1.
end if

allocate(material%depth(ufull,nl),material%volcp(ufull,nl),material%lambda(ufull,nl))
allocate(material%nodetemp(ufull,0:nl),material%nodetemp_prev(ufull,0:nl))
allocate(temp_ext(ufull),temp_int(ufull),storageflux(ufull))

! initialisation
timesteps = (int(spinup)+1)*(24*3600)/int(ddt)
material%nodetemp=290.                                            ! [K]
temp_ext = 290.                                                   ! [K]
temp_int = 290.                                                   ! [K]
! 4 layered wall
material%depth(1,:) =(/ 0.01, 0.04, 0.10, 0.05 /)                 ! [m]
material%volcp(1,:) =(/ 1.55E6, 1.55E6, 1.55E6, 0.29E6 /)         ! [J/m^3/K]
material%lambda(1,:)=(/ 0.9338, 0.9338, 0.9338, 0.0500 /)         ! [W/m/K]
! 200 layered wall (for high-res)
if (nl==200) then
  do t = 1,10
      material%depth(:,t)  = 0.001
      material%volcp(:,t)  = 1.55E6
      material%lambda(:,t) = 0.9338
  end do
  do t = 11,50
      material%depth(:,t)  = 0.001
      material%volcp(:,t)  = 1.55E6
      material%lambda(:,t) = 0.9338
  end do
  do t = 51,150
      material%depth(:,t)  = 0.001
      material%volcp(:,t)  = 1.55E6
      material%lambda(:,t) = 0.9338
  end do
  do t = 151,200
      material%depth(:,t)  = 0.001
      material%volcp(:,t)  = 0.29E6
      material%lambda(:,t) = 0.050
  end do
end if

! main loop
do t = 1, timesteps
  temp_ext = 290. + sin(2*pi*t*ddt/period)                        ! external temperature variation
  material%nodetemp_prev = material%nodetemp                      ! previous node temperatures
  call solvetridiag(temp_ext,temp_int,res_ext,res_int,ddt,     &       
                    material%nodetemp,                         &  ! layer temeperature
                    material%depth*material%volcp,             &  ! layer capacitance (cap)
                    material%depth/material%lambda)               ! layer resistance  (res)
  call energycheck(temp_ext,temp_int,res_ext,res_int,ddt,      &
                    material%nodetemp, material%nodetemp_prev, &  ! layer temeperature
                    material%depth*material%volcp,             &  ! layer capacitance (cap)
                    material%depth/material%lambda,            &  ! layer resistance  (res)
                    storageflux) 
  if ((t>=spinup*period/ddt) .and. MOD(t,1800/int(ddt))==0) then 
    write(1,'(I8,F12.4,F12.6,5F12.4)') t, temp_ext,storageflux, material%nodetemp(:,0:1)
  end if
end do
close(1)
deallocate(material%depth,material%volcp,material%lambda,material%nodetemp)
deallocate(material%nodetemp_prev,temp_ext,temp_int,storageflux)

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Tridiagonal matrix form:
!  [ ggB ggC         ] [ nodetemp ] = [ ggD ]
!  [ ggA ggB ggC     ] [ nodetemp ] = [ ggD ]
!  [     ggA ggB ggC ] [ nodetemp ] = [ ggD ]
!  [         ggA ggB ] [ nodetemp ] = [ ggD ]

!  ufull is number of systems to be calculated (for vectorisation)
!  nl is number of layers in system
!  conductmeth is conduction method (standard half-layer=0, modified half-layer=1, interface=2)

subroutine solvetridiag(temp_ext,temp_int,res_ext,res_int,ddt,nodetemp,cap,res)

implicit none

real(kind=dp), dimension(ufull),     intent(in)    :: temp_ext,temp_int  ! boundary temperatures
real(kind=dp), dimension(ufull,0:nl),intent(inout) :: nodetemp           ! temperature of each node
real(kind=dp), dimension(ufull,nl),  intent(in)    :: cap,res            ! layer capacitance & resistance
real(kind=dp), intent(in)                          :: ddt                ! timestep
real(kind=dp), intent(in)                          :: res_ext,res_int    ! aerodynamic surface resistance

! local variables
real(kind=dp), dimension(ufull,0:nl)               :: ggA,ggB,ggC,ggD    ! tridiagonal matrices
real(kind=dp), dimension(ufull)                    :: ggX                ! tridiagonal coefficient
integer sidx ! index start integer
integer k  ! loop integer

select case(conductmeth)
  case(0) !!!!!!!!! standard half-layer conduction !!!!!!!!!!!
    sidx = 1
    ggA(:,2:nl)   =-2./(res(:,1:nl-1) +res(:,2:nl))
    ggB(:,1)      = 2./(2.*res_ext+ res(:,1))  +2./(res(:,1)+res(:,2)) + cap(:,1)/ddt
    ggB(:,2:nl-1) = 2./(res(:,1:nl-2) +res(:,2:nl-1)) +2./(res(:,2:nl-1) +res(:,3:nl)) +cap(:,2:nl-1)/ddt
    ggB(:,nl)     = 2./(res(:,nl-1)+res(:,nl)) + 1./(0.5*res(:,nl)+res_int) + cap(:,nl)/ddt
    ggC(:,1:nl-1) =-2./(res(:,1:nl-1)+res(:,2:nl))
    ggD(:,1)      = temp_ext/(res_ext+0.5*res(:,1)) + nodetemp(:,1)*cap(:,1)/ddt
    ggD(:,2:nl-1) = nodetemp(:,2:nl-1)*cap(:,2:nl-1)/ddt
    ggD(:,nl)     = nodetemp(:,nl)*cap(:,nl)/ddt + temp_int/(0.5*res(:,nl) + res_int)
  case(1) !!!!!!!!! modified half-layer conduction !!!!!!!!!!!
    sidx = 0
    ggA(:,1)      =-2./res(:,1)
    ggA(:,2:nl)   =-2./(res(:,1:nl-1) +res(:,2:nl))
    ggB(:,0)      = 2./res(:,1) +1./res_ext
    ggB(:,1)      = 2./res(:,1) +2./(res(:,1)+res(:,2)) + cap(:,1)/ddt
    ggB(:,2:nl-1) = 2./(res(:,1:nl-2) +res(:,2:nl-1)) +2./(res(:,2:nl-1) +res(:,3:nl)) +cap(:,2:nl-1)/ddt
    ggB(:,nl)     = 2./(res(:,nl-1)+res(:,nl)) + 1./(0.5*res(:,nl)+res_int) + cap(:,nl)/ddt
    ggC(:,0)      =-2./res(:,1)
    ggC(:,1:nl-1) =-2./(res(:,1:nl-1)+res(:,2:nl))
    ggD(:,0)      = temp_ext/res_ext
    ggD(:,1:nl-1) = nodetemp(:,1:nl-1)*cap(:,1:nl-1)/ddt
    ggD(:,nl)     = nodetemp(:,nl)*cap(:,nl)/ddt + temp_int/(0.5*res(:,nl) + res_int)
  case(2) !!!!!!!!! interface conduction !!!!!!!!!!!
    sidx = 0
    ggA(:,1:nl)   = -1./res(:,1:nl)
    ggB(:,0)      =  1./res(:,1) +1./res_ext +0.5*cap(:,1)/ddt
    ggB(:,1:nl-1) =  1./res(:,1:nl-1) +1./res(:,2:nl) +0.5*(cap(:,1:nl-1) +cap(:,2:nl))/ddt
    ggB(:,nl)     =  1./res(:,nl) +1./res_int +0.5*cap(:,nl)/ddt
    ggC(:,0:nl-1) = -1./res(:,1:nl)
    ggD(:,0)      = nodetemp(:,0)*0.5*cap(:,1)/ddt +temp_ext/res_ext
    ggD(:,1:nl-1) = nodetemp(:,1:nl-1)*0.5*(cap(:,1:nl-1)+cap(:,2:nl))/ddt
    ggD(:,nl)     = nodetemp(:,nl)*0.5*cap(:,nl)/ddt +temp_int/res_int
  case DEFAULT
    write(6,*) "ERROR: Unknown conduction mode selected: ",conductmeth
    stop
end select
! tridiagonal solver (Thomas algorithm) for node temperatures
do k=sidx+1,nl
  ggX(:)   = ggA(:,k)/ggB(:,k-1)
  ggB(:,k) = ggB(:,k)-ggX(:)*ggC(:,k-1)
  ggD(:,k) = ggD(:,k)-ggX(:)*ggD(:,k-1)
end do
nodetemp(:,nl) = ggD(:,nl)/ggB(:,nl)
do k=nl-1,sidx,-1
  nodetemp(:,k) = (ggD(:,k) - ggC(:,k)*nodetemp(:,k+1))/ggB(:,k)
end do

end subroutine solvetridiag

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! conservation of energy: check for flux across boundaries matches change in heat storage

subroutine energycheck(temp_ext,temp_int,res_ext,res_int,ddt,nodetemp,nodetemp_prev,cap,res,storageflux)

implicit none

real(kind=dp), intent(in)                          :: res_ext,res_int        ! convective heat transfer
real(kind=dp), intent(in)                          :: ddt                    ! timestep
real(kind=dp), dimension(ufull),     intent(in)    :: temp_ext,temp_int      ! boundary temperatures
real(kind=dp), dimension(ufull,nl),  intent(in)    :: cap,res                ! layer capacitance & resistance
real(kind=dp), dimension(ufull,0:nl),intent(in)    :: nodetemp,nodetemp_prev ! temperature of each node
real(kind=dp), dimension(ufull),     intent(out)   :: storageflux            ! layer capacitance & resistance

! local variables
real(kind=dp), dimension(ufull) :: ggext,ggint,error

select case(conductmeth)
  case(0) !!!!!!!!! standard half-layer conduction !!!!!!!!!!!
    storageflux = sum((nodetemp(:,1:nl)-nodetemp_prev(:,1:nl))*cap(:,:), dim=2)/ddt
    ggext = (temp_ext-nodetemp(:,1))/(res_ext+0.5*res(:,1))
    ggint = (nodetemp(:,nl)-temp_int)/(0.5*res(:,nl)+res_int)
  case(1) !!!!!!!!! modified half-layer conduction !!!!!!!!!!!
    storageflux = sum((nodetemp(:,1:nl)-nodetemp_prev(:,1:nl))*cap(:,:), dim=2)/ddt
    ggext = (temp_ext-nodetemp(:,0))/res_ext
    ggint = (nodetemp(:,nl)-temp_int)/(0.5*res(:,nl)+res_int)
  case(2) !!!!!!!!! interface conduction !!!!!!!!!!!
    storageflux = sum((nodetemp(:,0:nl-1)-nodetemp_prev(:,0:nl-1))*0.5*cap(:,:), dim=2)/ddt   &
                 +sum((nodetemp(:,1:nl)  -nodetemp_prev(:,1:nl))  *0.5*cap(:,:), dim=2)/ddt
    ggext = (temp_ext-nodetemp(:,0))/res_ext
    ggint = (nodetemp(:,nl)-temp_int)/res_int
end select

error = (storageflux - (ggext - ggint))
! write(6,*) "closure error: ",error
if (any(abs(error)>1.E-6)) write(6,*) "closure error: ",error

end subroutine energycheck

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

end program conduction