! **********************************************************************
! *                                                                    *
! *     1-D Half-space cooling model for plate thermal evolution       *
! *                                                                    * 
! **********************************************************************

    program plcm

      implicit none

      real(8):: rho,  Cp,   ak,   D
      real(8):: age,  T0,   TM
      real(8):: zk,   z,    zmax, dz,   dzk,  t
      real(8):: hpl,  si_ex
      integer:: mz,   j,    ndeg
      character(40):: apara
      character(30):: fname

      real(8),allocatable:: Tmp(:)

!  *****  read parameters  *****

!  Open file

      open ( 10, file = 'pl_para2_w.dat' )

!  read file 

      read (10,*) apara
      read (10,*) rho
      write (6,*) apara,' = ',rho

      read (10,*) apara
      read (10,*) Cp
      write (6,*) apara,' = ',Cp

      read (10,*) apara
      read (10,*) D
      write (6,*) apara,' = ',D

      read (10,*) apara
      read (10,*) age
      write (6,*) apara,' = ',age

      read (10,*) apara
      read (10,*) T0
      write (6,*) apara,' = ',T0

      read (10,*) apara
      read (10,*) TM
      write (6,*) apara,' = ',TM

      read (10,*) apara
      read (10,*) mz
      write (6,*) apara,' = ',mz

      read (10,*) apara
      read (10,*) zmax
      write (6,*) apara,' = ',zmax

      read (10,*) apara
      read (10,*) fname
      write (6,*) apara,' = ',fname

!  convert unit km -> m, Ma -> sec

      t    = age  * 1.0d6 * 365.2422d0 * 24.0d0 * 3600.0d0
      hpl  = zmax * 1.0d3

!  depth interval

      dz   = hpl / dble( mz )
      dzk  = dz  / 1.0d3

!  thermal conductivity

      ak = rho * Cp * D

!  allocate dimension variables

      allocate ( Tmp( 0:mz ) )

!  maximum degree of Fourier series

        ndeg = mz / 2

!
!  ******  calculate temperature  ******
!

!  at the surface

      Tmp(0)  = T0
      Tmp(mz) = TM

!  at z = dz to zmax

      do j = 1,mz-1
        z = dz * dfloat( j )

!  Fouries series of saw-tooth function and decay with time
 
        call     fourexp( si_ex, z, t,  D, hpl, ndeg ) 

!  composite the steady-state and varying temperatures

        Tmp(j) = T0 + ( TM-T0 ) * ( z/hpl + si_ex ) 

!        write (6,*) z(j),Tmp(j),si_ex

      end do
!
!   *****  Output results to a file
!

!  Open file

      open ( 20, file = fname )

!  write file

      do j = 0,mz
        zk = dzk * dfloat( j )
        write (20,*) zk, Tmp( j )
      end do

!  close files 

      close (10)
      close (20)

      stop

    end program plcm