!-----------------------------------------------------------------------
!   Heat conduction: Euler -- Centered Difference
!
!   coded by Shigeo Yoden on Sept. 12, 2000
!-----------------------------------------------------------------------
program cond_eu

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: jmax=1000, mmax=10
  real,    parameter :: xmax=0.5, tmax=1.0
  real,    parameter :: rkappa=0.14
  real, dimension(0:jmax) :: x, u, du
  real, dimension(0:jmax,0:mmax) :: v
  character :: ct*30, cl*6, csgi*1

!!!    write(*,*) 'dx=?;   dt=? '
!!!    read (*,*)  dx, dt
    dx = 5.e-3;  dt = 5.e-5

    imax = xmax/dx
    nmax = tmax/dt;    ndel = nmax/mmax
    if(imax > jmax) stop 'error: dx must have a larger value'

!-- 初期条件 ----
    u = 100
    u(0) = 0;  u(imax) = 0      ! 実験(1)：両端が 0 度
!!!    u(0) = 0                    ! 実験(2)：左端が 0 度

    do i=0, imax
      x(i) = i*dx
      v(i,0) = u(i)
    end do
    mm = 0

!-- 時間積分 ----
    do n=1, nmax
      call drv(u, du)
      u = u + du

      if(mod(n,ndel) == 0) then
        mm = mm + 1
        v(0:imax,mm) = u
      end if
    end do

!-- グラフ化 ----
    call DclOpenGraphics()
    call graph
    call DclCloseGraphics

!-- 内部手続き ----
  contains
    subroutine drv(u, du)         ! 時間微分を計算する内部サブルーチン
      real, dimension(0:imax) :: u, du

      du(0)    = 0                ! 境界条件(x=0)
      du(imax) = 0                ! 境界条件(x=xmax)

      do i=1, imax-1
        du(i) = rkappa*(u(i+1) - 2*u(i) + u(i-1))/dx**2 *dt
      end do
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('x', 'temperature', 'cm', '|'//csgi(4)//'"C')
      call DclDrawScaledGraph(x(0:imax), v(0:imax,0))

      call DclSetParm('ENABLE_LINE_LABELING', .true.)
      call DclSetParm('LINE_CYCLE_LENGTH', 45.)
      call DclSetLineTextSize(0.018)
      cl = 't=?.?s'
      do m=1,mm
          write(cl(3:5),'(f3.1)') ndel*m*dt
        call DclSetLineText(cl)
          if(m == 6) call DclSetParm('ENABLE_LINE_LABELING', .false.)
        call DclDrawLine(x(0:imax), v(0:imax,m))
      end do

        ct = 'dt=??.??E???s'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('t', ct, position=-0.8)
      call DclDrawTitle('t', ' Euler-Centered Difference', position=-1.)
      call DclDrawTitle('t', 'Heat Conduction Equation:', position=-1.)
     end subroutine graph
end program cond_eu
!-----------------------------------------------------------------------
!   Heat conduction: Runge Kutta(4th) -- Centered Difference
!
!   coded by Satoshi Sakai for lectures on Geophysical Fluid Dynamics
!   modified by Shigeo Yoden on Sept. 12, 2000
!-----------------------------------------------------------------------
program cond_rk

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: jmax=1000, mmax=10
  real,    parameter :: xmax=0.5, tmax=1.0
  real,    parameter :: rkappa=0.14
  real, dimension(0:jmax) :: x, u, du1, du2, du3, du4
  real, dimension(0:jmax,0:mmax) :: v
  character :: ct*30, cl*6, csgi*1

!!!    write(*,*) 'dx=?;   dt=? '
!!!    read (*,*)  dx, dt
    dx = 5.e-3;  dt = 1.e-4

    imax = xmax/dx
    nmax = tmax/dt;    ndel = nmax/mmax
    if(imax > jmax) stop 'error: dx must have a larger value'

!-- 初期条件 ----
    u = 100
    u(0) = 0;  u(imax) = 0      ! 実験(1)：両端が 0 度
!!    u(0) = 0                    ! 実験(2)：左端が 0 度

    do i=0, imax
      x(i) = i*dx
      v(i,0) = u(i)
    end do
    mm = 0

!-- 時間積分 ----
    do n=1, nmax
      call drv(u      , du1)
      call drv(u+du1/2, du2)
      call drv(u+du2/2, du3)
      call drv(u+du3  , du4)
      u = u + (du1 + 2*(du2+du3) + du4) / 6

      if(mod(n,ndel) == 0) then
        mm = mm + 1
        v(0:imax,mm) = u
      end if
    end do

!-- グラフ化 ----
    call DclOpenGraphics()
    call graph
    call DclCloseGraphics

!-- 内部手続き ----
  contains
    subroutine drv(u, du)         ! 時間微分を計算する内部サブルーチン
      real, dimension(0:imax) :: u, du

      du(0)    = 0                ! 境界条件(x=0)
      du(imax) = 0                ! 境界条件(x=xmax)

      do i=1, imax-1
        du(i) = rkappa*(u(i+1) - 2*u(i) + u(i-1))/dx**2 *dt
      end do
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('x', 'temperature', 'cm', '|'//csgi(4)//'"C')
      call DclDrawScaledGraph(x(0:imax), v(0:imax,0))

      call DclSetParm('ENABLE_LINE_LABELING', .true.)
      call DclSetParm('LINE_CYCLE_LENGTH', 45.)
      call DclSetLineTextSize(0.018)
      cl = 't=?.?s'
      do m=1,mm
          write(cl(3:5),'(f3.1)') ndel*m*dt
        call DclSetLineText(cl)
          if(m == 6) call DclSetParm('ENABLE_LINE_LABELING', .false.)
        call DclDrawLine(x(0:imax), v(0:imax,m))
      end do

        ct = 'dt=??.??E???s'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('t', ct, position=-0.8)
      call DclDrawTitle('t', ' 4th Runge Kutta-Center Difference', &
                                                          position=-1.)
      call DclDrawTitle('t', 'Heat Conduction Equation:', position=-1.)
    end subroutine graph
end program cond_rk
!-----------------------------------------------------------------------
!   Heat conduction: an implicit method (Euler -- Centered Difference)
!
!   coded by Shigeo Yoden on Sept. 12, 2000
!-----------------------------------------------------------------------
program cond_imp

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: jmax=1000, mmax=10
  real,    parameter :: xmax=0.5, tmax=1.0
  real,    parameter :: rkappa=0.14
  real, dimension(0:jmax) :: x, u
  real, dimension(0:jmax,0:mmax) :: v
  real, dimension(1:jmax-1) :: a, b, c
  character :: ct*30, cl*6, csgi*1

!!!    write(*,*) 'dx=?;   dt=? '
!!!    read (*,*)  dx, dt
    dx = 5.e-3;  dt = 1.e-3

    imax = xmax/dx
    nmax = tmax/dt;    ndel = nmax/mmax
    if(imax > jmax) stop 'error: dx must have a larger value'

!-- 初期条件 ----
    u = 100
    u(0) = 0;  u(imax) = 0      ! 実験(1)：両端が 0 度

    do i=0, imax
      x(i) = i*dx
      v(i,0) = u(i)
    end do
    mm = 0

!-- 時間積分 ----
    r = dt*rkappa/dx**2
    a = - r
    b = 1 + 2*r
    c = - r
    do n=1, nmax
      call tridiag(u(1))

      if(mod(n,ndel) == 0) then
        mm = mm + 1
        v(0:imax,mm) = u
      end if
    end do

!-- グラフ化 ----
    call DclOpenGraphics()
    call graph
    call DclCloseGraphics

!-- 内部手続き ----
  contains
    subroutine tridiag(f)         ! 3重対角連立方程式を解く内部サブルーチン
!
!-- Durran(1999) pp.440;  Thomas tridiagonal algorism
! Solves a standard tridiagonal system
!
! Definition of the variables
!   jmax = dimension of all the following arrays  <==> imax-1
!   a    = sub(lower) diagonal
!   b    = center diagonal
!   c    = super(upper) diagonal
!   f    = right hand side
!   q    = work array provided by calling program
!
!   a(1) and c(jmax) need not be initialized
! The output is in f; a, b, c are unchanged
! 両端の境界条件は 0 のみ ?

  real, dimension(1:imax-1) :: f, q

      c(imax-1) = 0.

! Forward elimination sweep
      q(1) = - c(1)/b(1)
      f(1) =   f(1)/b(1)
      do j=2, imax-1
        p = 1./(b(j) + a(j)*q(j-1))
        q(j) = - c(j)*p
        f(j) = (f(j) - a(j)*f(j-1))*p
      end do

! Backword pass
      do j=imax-2, 1, -1
        f(j) = f(j) + q(j)*f(j+1)
      end do
    end subroutine tridiag

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('x', 'temperature', 'cm', '|'//csgi(4)//'"C')
      call DclDrawScaledGraph(x(0:imax), v(0:imax,0))

      call DclSetParm('ENABLE_LINE_LABELING', .true.)
      call DclSetParm('LINE_CYCLE_LENGTH', 45.)
      call DclSetLineTextSize(0.018)
      cl = 't=?.?s'
      do m=1,mm
          write(cl(3:5),'(f3.1)') ndel*m*dt
        call DclSetLineText(cl)
          if(m == 6) call DclSetParm('ENABLE_LINE_LABELING', .false.)
        call DclDrawLine(x(0:imax), v(0:imax,m))
      end do

        ct = 'dt=??.??E???s'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('t', ct, position=-0.8)
      call DclDrawTitle('t', ' implicit method(Euler-C.D.)', position=-1.)
      call DclDrawTitle('t', 'Heat Conduction Equation:', position=-1.)
    end subroutine graph
end program cond_imp
!-----------------------------------------------------------------------
!   Advection equation: Forward -- Upstream
!
!   coded by Shigeo Yoden on Sept. 12, 2000
!-----------------------------------------------------------------------
program adv_ftus

  use dcl        ! DCL Fortran90 (dcl-5.1/dcl-f90)

  integer, parameter :: jmax=1000, mmax=5
  real,    parameter :: xmax=1.0, tmax=1.0
  real,    parameter :: c=1.0, pi=3.141593
  real, dimension(0:jmax) :: x, u, du
  real, dimension(0:jmax,0:mmax) :: v
  character :: ct*30, cl*6

!!!    write(*,*) 'dx=?;   dt=? '
!!!    read (*,*)  dx, dt
    dx = 1.e-2;  dt = 5.e-3

    imax = xmax/dx
    nmax = tmax/dt;    ndel = nmax/mmax
    if(imax > jmax) stop 'error: dx must have a larger value'

!-- 初期条件 ----
!!!    write(*,*) 'initial condition # ?  1:cos  2:step  3:delta'
!!!    read (*,*) init
    init = 2

    if(init == 1) then
      do i=0, imax
        u(i) = cos(2*pi*i/imax)
      end do
    else if(init == 2) then
      u = - 1
      u(imax/4:imax/2) = 1
    else
      u = -1
      u(imax/4) = 1
    end if

    do i=0, imax
      x(i) = i*dx
      v(i,0) = u(i)
    end do
    mm = 0

!-- 時間積分 ----
    do n=1, nmax
      call drv(u, du)
      u = u + du

      if(mod(n,ndel) == 0) then
        mm = mm + 1
        v(0:imax,mm) = u
      end if
    end do

!-- グラフ化 ----
    call DclOpenGraphics()
    call graph
    call DclCloseGraphics

!-- 内部手続き ----
  contains
    subroutine drv(u, du)         ! 時間微分を計算する内部サブルーチン
      real, dimension(0:imax) :: u, du

      do i=1, imax
        du(i) = - c*(u(i) - u(i-1))/dx *dt
      end do
      du(0)    = du(imax)         ! 周期境界条件: x=0 <=> x=xmax
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('x', 'u(x,t)')
      call DclDrawScaledGraph(x(0:imax), v(0:imax,0))

      call DclSetParm('ENABLE_LINE_LABELING', .true.)
      call DclSetParm('LINE_CYCLE_LENGTH', 50.)
      call DclSetLineTextSize(0.02)
      cl = 't=?.?s'
      do m=1,mm
        call DclSetLineType(mod(m,4)+1)
          write(cl(3:5),'(f3.1)') ndel*m*dt
        call DclSetLineText(cl)
        call DclDrawLine(x(0:imax), v(0:imax,m))
      end do

        ct = 'dt=??.??E???s'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('t', ct, position=-0.8)
      call DclDrawTitle('t', ' forward and upstream', position=-1.)
      call DclDrawTitle('t', 'Advection Equation:', position=-1.)
    end subroutine graph
end program adv_ftus