! $Id$
!
module Timestep
!
  use Cdata
!
  implicit none
!
  private
!
  include 'timestep.h'
!
  ! Parameters for adaptive time stepping
  real, parameter :: safety           = 0.9
  real, parameter :: dt_decrease      = -0.25
  real, parameter :: dt_increase      = -0.20
  real            :: errcon, dt_next
  real, dimension(mvar) :: farraymin
  logical :: fixed_dt=.false.
!
  contains
!
!***********************************************************************
    subroutine initialize_timestep
!
      use Messages, only: warning
      use General, only: rtoa
!
      ldt = .false.
!
      if (dt0>0.) then
        dt=dt0
      elseif (dt0<0.) then
        fixed_dt=.true.
        dt=-dt0
      else
        if (dt==0) then
          call warning('initialize_timestep','dt=0 not appropriate for Runge-Kutta-Fehlberg'// &
                     'set to dt_epsi='//trim(rtoa(dt_epsi)))
          dt=dt_epsi
        endif
      endif
!
      if (eps_rkf0/=0.) eps_rkf=eps_rkf0
      dt_next=dt

    endsubroutine initialize_timestep
!***********************************************************************
    subroutine time_step(f,df,p)
!
!  Runge-Kutta-Fehlberg accurate to 5th order
!
!  22-jun-06/tony: coded
!
      use Messages
!!      use Particles_main
!!      use Shear, only: advance_shear
!
      real, dimension (mx,my,mz,mfarray) :: f
      real, dimension (mx,my,mz,mvar) :: df
      type (pencil_case) :: p
!
      real :: errmax, tnew
      real :: dt_temp, dt_next, dt_did
      integer :: j,i
!
      ! General error condition
      errcon = (5.0/safety)**(1.0/dt_increase)
!
      if (itorder/=5) then
        call warning('time_step','itorder set to 5 for Runge-Kutta-Fehlberg')
        itorder=5
      endif
!
!  dt_beta_ts may be needed in other modules (like Dustdensity) for fixed dt
!
!      if (.not. ldt) dt_beta_ts=dt*beta_ts
!
      if (lshear.or.lparticles) call fatal_error("time_step","Shear and particles are not"// &
                                                 " yet supported by the adaptive rkf scheme")
!
      lfirst=.true.
      ! Time step to try after previous time
      dt = dt_next
      do i=1,10
        ! Do a Runge-Kutta step
        call rkck(f, df, p, errmax)
        ! Step succeeded so exit
        if (errmax <= safety) exit
        ! Step didn't succeed so decrease the time step
!        print*,"Decreasing"
        dt_temp = safety*dt*(errmax**dt_decrease)
        ! Don't decrease the time step by more than a factor of ten
        dt = sign(max(abs(dt_temp), 0.1*abs(dt)), dt)
        ! New time
        tnew = t+dt
        if (tnew == t) then
          ! Guard against infinitesimal time steps
          print*, 'WARNING: Timestep underflow in rkqs()'
        endif
      enddo
!
!      print*,"errmax, errcon", errmax,errcon
      if (errmax > errcon) then
        ! Increase the time step
        dt_next = safety*dt*(errmax**dt_increase)
      else
        ! But not by more than a factor of 5
        dt_next = 5.0*dt
      endif
!
      ! Time step that was actually performed
      dt_did = dt
!
      if (ip<=6) print*,'TIMESTEP: iproc,dt=',iproc_world,dt  !(all have same dt?)
      ! Increase time
      t = t+dt
!
!  Time evolution of grid variables
!  (do this loop in pencils, for cache efficiency)
!
      do j=1,mvar; do n=n1,n2; do m=m1,m2
        f(l1:l2,m,n,j)=f(l1:l2,m,n,j)+df(l1:l2,m,n,j)
      enddo; enddo; enddo
!
    endsubroutine time_step
!***********************************************************************
    subroutine rkck(f, df, p, errmax)
    ! Explicit fifth order Runge--Kutta--Fehlberg time stepping
      use Cdata
      use Mpicomm, only: mpiallreduce_max,MPI_COMM_WORLD
      use Messages
      use Equ
    ! RK parameters by Cash and Karp
      real, parameter :: b21      = 0.2
      real, parameter :: b31      = 0.075
      real, parameter :: b32      = 0.225
      real, parameter :: b41      = 0.3
      real, parameter :: b42      = -0.9
      real, parameter :: b43      = 1.2
      real, parameter :: b51      = -11.0 / 54.0
      real, parameter :: b52      = 2.5
      real, parameter :: b53      = -70.0 / 27.0
      real, parameter :: b54      = 35.0 / 27.0
      real, parameter :: b61      = 1631.0 / 55296.0
      real, parameter :: b62      = 175.0 / 512.0
      real, parameter :: b63      = 575.0 / 13824.0
      real, parameter :: b64      = 44275.0 / 110592.0
      real, parameter :: b65      = 253.0 / 4096.0
      real, parameter :: c1       = 37.0 / 378.0
      real, parameter :: c2       = 0.0
      real, parameter :: c3       = 250.0 / 621.0
      real, parameter :: c4       = 125.0 / 594.0
      real, parameter :: c5       = 0.0
      real, parameter :: c6       = 512.0 / 1771.0
      real, parameter :: dc1      = c1 - 2825.0 / 27648.0
      real, parameter :: dc2      = c2 - 0.0
      real, parameter :: dc3      = c3 - 18575.0 / 48384.0
      real, parameter :: dc4      = c4 - 13525.0 / 55296.0
      real, parameter :: dc5      = c5 - 277.0 / 14336.0
      real, parameter :: dc6      = c6 - 0.25
!
      real, dimension (mx,my,mz,mfarray), intent(inout) :: f
      real, dimension (mx,my,mz,mvar), intent(out) :: df
      type (pencil_case), intent(inout) :: p
      real, dimension(mx,my,mz,mfarray,5) :: k
      real, save, dimension(mx,my,mz,mfarray) :: tmp1
      real, dimension(mx,my,mz,mfarray) :: tmp2
      real, dimension(nx) :: scal, err
      real, intent(inout) :: errmax
      real :: errmaxs
      integer :: j
      logical, save :: first_call=.true.
!
      if (ny /= 1 .or. nz /= 1) then
        call fatal_error("rkck", "timestep_rkf_1d only works for the 1D case")
      endif
!
      do j=1,mvar
        farraymin(j) = max(dt_ratio*maxval(abs(f(l1:l2,m1:m2,n1:n2,j))),dt_epsi)
      enddo
      if (lroot.and.it==1) print*,"farraymin",farraymin
!
      errmax=0.
!
      if (first_call) then
        ! Initialize tmp1 to arbitrary value /= 0, so chemistry.f90 (which
        ! operates in the ghost zones where it most probably shouldn't)
        ! doesn't divide by 0.
        tmp1 = real(0.577215664901532860606512090082402431042159335)
        first_call = .false.
      endif
!
      tmp2(:,m1:m2,n1:n2,:) = 0.
      call pde(f,tmp2,p)
      k(:,m1:m2,n1:n2,:,1) = tmp2(:,m1:m2,n1:n2,:)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          k(l1:l2,m,n,j,1) = dt*k(l1:l2,m,n,j,1)
!
      enddo; enddo; enddo
!
      lfirst=.false.
!
      tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
          + b21*k(:,m1:m2,n1:n2,:,1)
      tmp2(:,m1:m2,n1:n2,:) = 0.
      call pde(tmp1,tmp2,p)
      k(:,m1:m2,n1:n2,:,2) = tmp2(:,m1:m2,n1:n2,:)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          k(l1:l2,m,n,j,2) = dt*k(l1:l2,m,n,j,2)
      enddo; enddo; enddo
!
      tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
          + b31*k(:,m1:m2,n1:n2,:,1) &
          + b32*k(:,m1:m2,n1:n2,:,2)
      tmp2(:,m1:m2,n1:n2,:) = 0.
      call pde(tmp1,tmp2,p)
      k(:,m1:m2,n1:n2,:,3) = tmp2(:,m1:m2,n1:n2,:)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          k(l1:l2,m,n,j,3) = dt*k(l1:l2,m,n,j,3)
      enddo; enddo; enddo
!
      tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
          + b41*k(:,m1:m2,n1:n2,:,1) &
          + b42*k(:,m1:m2,n1:n2,:,2) &
          + b43*k(:,m1:m2,n1:n2,:,3)
      tmp2(:,m1:m2,n1:n2,:) = 0.
      call pde(tmp1,tmp2,p)
      k(:,m1:m2,n1:n2,:,4) = tmp2(:,m1:m2,n1:n2,:)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          k(l1:l2,m,n,j,4) = dt*k(l1:l2,m,n,j,4)
      enddo; enddo; enddo
!
      tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
          + b51*k(:,m1:m2,n1:n2,:,1) &
          + b52*k(:,m1:m2,n1:n2,:,2) &
          + b53*k(:,m1:m2,n1:n2,:,3) &
          + b54*k(:,m1:m2,n1:n2,:,4)
      tmp2(:,m1:m2,n1:n2,:) = 0.
      call pde(tmp1,tmp2,p)
      k(:,m1:m2,n1:n2,:,5) = tmp2(:,m1:m2,n1:n2,:)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          k(l1:l2,m,n,j,5) = dt*k(l1:l2,m,n,j,5)
      enddo; enddo; enddo
!
      errmaxs=0.
      tmp1(:,m1:m2,n1:n2,:) = f(:,m1:m2,n1:n2,:) &
          + b61*k(:,m1:m2,n1:n2,:,1) &
          + b62*k(:,m1:m2,n1:n2,:,2) &
          + b63*k(:,m1:m2,n1:n2,:,3) &
          + b64*k(:,m1:m2,n1:n2,:,4) &
          + b65*k(:,m1:m2,n1:n2,:,5)
      df(:,m1:m2,n1:n2,:)=0.
      call pde(tmp1,df,p)
      do j=1,mvar; do n=n1,n2; do m=m1,m2
          df(l1:l2,m,n,j) = dt*df(l1:l2,m,n,j)
!
          err = dc1*k(l1:l2,m,n,j,1) + dc2*k(l1:l2,m,n,j,2) + &
                dc3*k(l1:l2,m,n,j,3) + dc4*k(l1:l2,m,n,j,4) + &
                dc5*k(l1:l2,m,n,j,5) + dc6*df(l1:l2,m,n,j)
!
          df(l1:l2,m,n,j) = c1*k(l1:l2,m,n,j,1) + c2*k(l1:l2,m,n,j,2) + &
                            c3*k(l1:l2,m,n,j,3) + c4*k(l1:l2,m,n,j,4) + &
                            c5*k(l1:l2,m,n,j,5) + c6*df(l1:l2,m,n,j)
!
          ! Get the maximum error over the whole field
          !
          select case (timestep_scaling(j))
          case ('per_var_err')
            !
            ! Per variable error
            !
            scal=  ( &
                 sqrt(f(l1:l2,m,n,1)**2+f(l1:l2,m,n,2)**2)  + &
                 sqrt(k(l1:l2,m,n,1,1)**2 + k(l1:l2,m,n,2,1)**2) + &
                 1e-30)
            errmaxs = max(maxval(abs(err/scal)),errmaxs)
            !scal=  ( &
            !     abs(f(l1:l2,m,n,j))  + abs(k(l1:l2,m,n,j,1)) + 1e-30)
            !errmaxs = max(maxval(abs(err/scal)),errmaxs)
          case ('cons_frac_err')
            !
            ! Constant fractional error
            !
            errmaxs = max(maxval(abs(err/f(l1:l2,m,n,j))),errmaxs)
          case ('cons_err')
            !
            ! Constant error
            !
            scal = max(abs(f(l1:l2,n,m,j)), 1e-8)
            errmaxs = max(maxval(abs(err/scal)),errmaxs)
            !
          case ('rel_err')
            !
            ! Relative error to f constrained with farraymin floor
            !
            scal=max(abs(f(l1:l2,m,n,j)+df(l1:l2,m,n,j)),farraymin(j))
            errmaxs = max(maxval(abs(err/scal)),errmaxs)
            !
          case ('none')
            !
            ! No error check
            !
            errmaxs = 0
            !
          endselect
          !
        enddo; enddo; enddo
        !
        ! Divide your maximum error by the required accuracy
        !
        errmaxs=errmaxs/eps_rkf
        !
      call mpiallreduce_max(errmaxs,errmax,MPI_COMM_WORLD)
!
    endsubroutine rkck
!***********************************************************************
    subroutine pushpars2c(p_par)

    use Messages, only: fatal_error

    integer, parameter :: n_pars=0
    integer(KIND=ikind8), dimension(:) :: p_par

    call fatal_error('timestep_rkf_1d','alpha_ts, beta_ts not defined')

    endsubroutine pushpars2c
!***********************************************************************
endmodule Timestep
