!----------------------------------------------------------------------
!   Euler method: dx/dt = 1 - x**2
!
!   coded by Satoshi Sakai for lectures on Geophysical Fluid Dynamics
!   modified by Shigeo Yoden on Sept. 10, 2000
!----------------------------------------------------------------------
program euler01

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

  integer, parameter :: mmax=200
  real, parameter :: tmax=2.
  real, dimension(0:mmax) :: t, t0, x, x0
  character :: ct*30

!!!    write(*,*) 'dt=?'
!!!    read (*,*)  dt
    dt = 0.2
 
    nmax = tmax/dt
    if(nmax > mmax) stop 'error: dt must have a larger value'

!-- 初期条件 ----
    t(0) = 0
    x(0) = 0

!-- 時間積分 ----
    do n=0, nmax-1
      dx = (1 - x(n)**2) *dt
      x(n+1) = x(n) + dx
      t(n+1) = t(n) + dt
    end do

!-- 解析解 ----
    do m=0, mmax
      t0(m) = tmax/mmax *m
      x0(m) = tanh(t0(m))
    end do

!-- グラフ化 ----
    call DclOpenGraphics()
    call DclNewFrame
      
    call DclSetTitle('time', 'x')
    call DclDrawScaledGraph(t0, x0)
    call DclDrawLine  (t(0:nmax), x(0:nmax), type=3)
    call DclDrawMarker(t(0:nmax), x(0:nmax), type=3)

      ct = 'dt=xx.xxExxx'
      write(ct(4:12),'(e9.2)') dt
    call DclDrawTitle('T', ct, position=-0.8)
    call DclDrawTitle('T', 'Euler method', position=-1.)

    call DclCloseGraphics
end program euler01
!----------------------------------------------------------------------
!   Improved Euler(Heun) method: dx/dt = 1 - x**2
!
!   coded by Shigeo Yoden on Sept. 12, 2000
!----------------------------------------------------------------------
program heun01

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

  integer, parameter :: mmax=200
  real, parameter :: tmax=2.
  real, dimension(0:mmax) :: t, t0, x, x0
  character :: ct*30
  
!!!    write(*,*) 'dt=?'
!!!    read (*,*)  dt
    dt = 0.2

    nmax = tmax/dt
    if(nmax > mmax) stop 'dt must have a larger value'

!-- 初期条件 ----
    t(0) = 0
    x(0) = 0

!-- 時間積分 ----
    do n=0, nmax-1
      call drv(x(n)    , dx1)
      call drv(x(n)+dx1, dx2)
      x(n+1) = x(n) + (dx1+dx2)/2
      t(n+1) = t(n) + dt
    end do

!-- 解析解 ----
    do m=0, mmax
      t0(m) = tmax/mmax *m
      x0(m) = tanh(t0(m))
    end do

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

!-- 内部手続き ----
  contains
    subroutine drv(x, dx)         ! 時間微分を計算する内部サブルーチン
      real :: x, dx
      dx = (1 - x**2) *dt
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('time', 'x')
      call DclDrawScaledGraph(t0, x0)
      call DclDrawLine  (t(0:nmax), x(0:nmax), type=3)
      call DclDrawMarker(t(0:nmax), x(0:nmax), type=3)

        ct = 'dt=xx.xxExxx'
        write(ct(4:12),'(e9.2)') dt
    call DclDrawTitle('T', ct, position=-0.8)
    call DclDrawTitle('T', 'Modified Euler method', position=-1.)
    end subroutine graph
end program heun01
!----------------------------------------------------------------------
!   4-th Order Runge-Kutta method: dx/dt = 1 - x**2
!
!   coded by Satoshi Sakai for lectures on Geophysical Fluid Dynamics
!   modified by Shigeo Yoden on Sept. 10, 2000
!----------------------------------------------------------------------
program rk01

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

  integer, parameter :: mmax=200
  real, parameter :: tmax=2.
  real, dimension(0:mmax) :: t, t0, x, x0
  character :: ct*30
  
!!!    write(*,*) 'dt=?'
!!!    read (*,*)  dt
    dt = 0.2

    nmax = tmax/dt
    if(nmax > mmax) stop 'error: dt must have a larger value'

!-- 初期条件 ----
    t(0) = 0
    x(0) = 0

!-- 時間積分 ----
    do n=0, nmax-1
      call drv(x(n)      , dx1)
      call drv(x(n)+dx1/2, dx2)
      call drv(x(n)+dx2/2, dx3)
      call drv(x(n)+dx3  , dx4)
      x(n+1) = x(n) + (dx1 + 2*(dx2 + dx3) + dx4) / 6
      t(n+1) = t(n) + dt
    end do

!-- 解析解 ----
    do m=0, mmax
      t0(m) = tmax/mmax *m
      x0(m) = tanh(t0(m))
    end do

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

!-- 内部手続き ----
  contains
    subroutine drv(x, dx)         ! 時間微分を計算する内部サブルーチン
      real :: x, dx
      dx = (1 - x**2) *dt
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('time', 'x')
      call DclDrawScaledGraph(t0, x0)
      call DclDrawLine  (t(0:nmax), x(0:nmax), type=3)
      call DclDrawMarker(t(0:nmax), x(0:nmax), type=3)

        ct = 'dt=xx.xxExxx'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('T', ct, position=-0.8)
      call DclDrawTitle('T', '4th-Order Runge-Kutta method', position=-1.)
    end subroutine graph
end program rk01
!----------------------------------------------------------------------
!   4-th Order Runge-Kutta method: Satellite Orbit
!
!   coded by Shigeo Yoden on Sept. 12, 2000
!----------------------------------------------------------------------
program rk02

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

  integer, parameter :: mmax=2000
  real, parameter :: tmax=16.
  real, dimension(0:mmax) :: t, x, y, u, v 
  real, dimension(4) :: z, dz1, dz2, dz3, dz4
  character :: ct*30

!!!    write(*,*) 'dt=?'
!!!    read (*,*)  dt
    dt = 0.05

    nmax = tmax/dt
    if(nmax > mmax) stop 'error: dt must have a larger value'

!-- 初期条件 ----
    t(0)=0
    x(0)=3  ;  z(1)=x(0)
    y(0)=0  ;  z(2)=y(0)
    u(0)=0.3;  z(3)=u(0)
    v(0)=0.2;  z(4)=v(0)

!-- 時間積分 ----
    do n=0, nmax-1
      call drv(z      , dz1)
      call drv(z+dz1/2, dz2)
      call drv(z+dz2/2, dz3)
      call drv(z+dz3  , dz4)
      z = z + (dz1 + 2*(dz2 + dz3) + dz4) / 6
      t(n+1) = t(n) + dt
      x(n+1) = z(1);  y(n+1) = z(2)
      u(n+1) = z(3);  v(n+1) = z(4)
    end do

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

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

      fxy = z(1)**2 + z(2)**2
      fxy = sqrt( fxy**3 )
      dz(1) = z(3)       *dt
      dz(2) = z(4)       *dt
      dz(3) = - z(1)/fxy *dt
      dz(4) = - z(2)/fxy *dt
    end subroutine drv

    subroutine graph              ! グラフを描く内部サブルーチン
      call DclNewFrame
      
      call DclSetTitle('x', 'y')
      call DclScalingPoint   (y(0:nmax), x(0:nmax)) ! 縦軸と横軸の範囲を同一に
      call DclDrawScaledGraph(x(0:nmax), y(0:nmax))
      call DclDrawMarker(x(0:0), y(0:0), type=230)

        ct = 'dt=xx.xxExxx'
        write(ct(4:12),'(e9.2)') dt
      call DclDrawTitle('T', ct, position=-0.8)
      call DclDrawTitle('T', 'Satellite Orbit(RK-4)', position=-1.)
    end subroutine graph
end program rk02