!
!  [Auto-generated file, so think twice before editing]
!
!  A second-order timestepping module similar to RKC (Runge-Kutta-Chebyshev).
!  The schemes used here are all second-order (p=2) accurate Runge-Kutta
!  schemes of stage number (number of substeps) s > 2 that trade order for
!  extended stability interval.
!    For this file, s=40, so we have a 2nd-order, 40-step Runge-Kutta
!  scheme with a critical Courant number of ~1044.8 as compared to 2.513 for
!  any p=s=3 Runge-Kutta scheme (like the Williamson scheme in timestep.f90).
!  Here the Courant number is
!    Cou = c nu dt / dx^2 ,
!  where
!    c = 272/45 = 6.04444
!  for 6th-order central finite differences in space.
!
!  This scheme uses 5N array slots (as opposed to 2N for the Williamson
!  scheme in timestep.f90), irrespective of s.
!  [TODO: it currently uses more, but this should be fixed...]

module Timestep

    use Cparam
    use Cdata
    use Equ

    implicit none

    private

    include 'timestep.h'

contains

!***********************************************************************
    subroutine initialize_timestep
    endsubroutine initialize_timestep
!***********************************************************************
    subroutine time_step(f,df,p)
    !
    !  Long-time-step Runge--Kutta--Chebyshev stepping, accurate to second
    !  order.
    !
    !  18-aug-08/perl: generated
    !

        use Mpicomm, only: mpiallreduce_max,MPI_COMM_WORLD

        real, dimension (mx,my,mz,mfarray) :: f
        ! real, dimension (mx,my,mz,mvar) :: fn_target, fn1_target
        real, dimension (mx,my,mz,mvar) :: df, dfn
        type (pencil_case) :: p
        real :: dt1, dt1_local
        real, save :: dt1_last=0.0
        double precision :: t0
        integer :: iv

        ! Use pointers for cheaply flipping fn and fn1 after each substep
        ! target :: f, df
        ! target :: fn_target, fn1_target
        ! real, pointer :: f0(:,:,:,:), df0(:,:,:,:)
        ! real, pointer :: fn(:,:,:,:), fn1(:,:,:,:)
        real, dimension(mx,my,mz,mvar) :: f0, df0, fn, fn1

        ! f0  => f(:,:,:,1:mvar)
        f0  = f(:,:,:,1:mvar)
        ! fn  => fn_target;
        ! fn1 => fn1_target;

        ! Step n = 0:
        lfirst = .true.
        t0 = t
        df0 = 0.
        call pde(f,df0,p)
        !
        ! In the first time substep we need to calculate timestep dt.
        ! Only do this on root processor, then broadcast dt to all others.
        !
        if (ldt) then
            dt1_local=maxval(dt1_max(1:nx))

            ! Timestep growth limiter
            if (real(ddt) > 0.) dt1_local=max(dt1_local,dt1_last)
            call mpiallreduce_max(dt1_local,dt1,MPI_COMM_WORLD)
            dt=1.0/dt1
            ! Timestep growth limiter
            if (ddt/=0.) dt1_last=dt1_local/ddt
        endif
        !
        ! IMPLEMENT ME:
        ! What do we need to do with dt_beta_ts?
        ! if (ldt) dt_beta_ts=dt*beta_ts
        !
        if (ip<=6) print*, 'TIMESTEP: iproc,dt=',iproc_world,dt  ! same dt everywhere?

        lfirst = .false.

        ! Step n = 1:
        fn1 = f0
        fn = f0 + 0.000478510350458706*dt*df0
        t = t0 + 0.000478510350458706*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 2:
        fn1 = -1*fn1 &
              + 2.00019230769231*fn &
              + -0.000192307692307692*f0 &
              + 0.00382881900883236*dt*dfn &
              + -0.00287170628669371*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.00191422544427731*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 3:
        fn1 = -1.18491934033202*fn1 &
              + 2.37006654976796*fn &
              + -0.185147209435933*f0 &
              + 0.00453684169419635*dt*dfn &
              + -0.00340241318471333*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0051042740348718*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 4:
        fn1 = -1.24932722343328*fn1 &
              + 2.10891544854153*fn &
              + 0.140411774891754*f0 &
              + 0.00403693961142769*dt*dfn &
              + -0.00284027278217515*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.00956965523565877*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 5:
        fn1 = -1.07922511744088*fn1 &
              + 2.04737021857798*fn &
              + 0.0318548988628987*f0 &
              + 0.00391912816625756*dt*dfn &
              + -0.00269341187986087*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0153096826732083*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 6:
        fn1 = -1.03610718224988*fn1 &
              + 2.0246586256007*fn &
              + 0.0114485566491798*f0 &
              + 0.00387565305710038*dt*dfn &
              + -0.00263387110185681*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0223234745134432*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 7:
        fn1 = -1.01931052070667*fn1 &
              + 2.01417965180629*fn &
              + 0.00513086890038122*f0 &
              + 0.00385559393883322*dt*dfn &
              + -0.00260380749872705*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0306099539452179*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 8:
        fn1 = -1.01122971861282*fn1 &
              + 2.00860771144968*fn &
              + 0.00262200716314171*f0 &
              + 0.00384492798882859*dt*dfn &
              + -0.00258630690932264*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0401678497703005*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 9:
        fn1 = -1.00678529919261*fn1 &
              + 2.00532720315342*fn &
              + 0.00145809603918911*f0 &
              + 0.00383864835637774*dt*dfn &
              + -0.00257498145280966*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0509956970991234*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 10:
        fn1 = -1.00409312955981*fn1 &
              + 2.00323664307064*fn &
              + 0.000856486489173866*f0 &
              + 0.00383464655307447*dt*dfn &
              + -0.00256699548586237*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0630918381515577*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 11:
        fn1 = -1.00233448501582*fn1 &
              + 2.00181492163199*fn &
              + 0.000519563383824295*f0 &
              + 0.00383192505772194*dt*dfn &
              + -0.00256093915097348*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0764544231618548*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 12:
        fn1 = -1.00111179478493*fn1 &
              + 2.0007930098222*fn &
              + 0.000318784962728245*f0 &
              + 0.00382996888813383*dt*dfn &
              + -0.00255604960538671*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.0910814113867886*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 13:
        fn1 = -1.00021525216958*fn1 &
              + 2.00002220238262*fn &
              + 0.000193049786956464*f0 &
              + 0.00382849338892036*dt*dfn &
              + -0.00255188602860517*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.106970572215923*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 14:
        fn1 = -0.999526697737571*fn1 &
              + 1.99941565140714*fn &
              + 0.000111046330430768*f0 &
              + 0.00382733231360984*dt*dfn &
              + -0.00254817991056254*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.124119486382828*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 15:
        fn1 = -0.998976003137801*fn1 &
              + 1.99892027941372*fn &
              + 5.57237240841453e-05*f0 &
              + 0.00382638405993564*dt*dfn &
              + -0.00254476049940758*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.142525547275952*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 16:
        fn1 = -0.99851970518724*fn1 &
              + 1.99850238696316*fn &
              + 1.73182240811886e-05*f0 &
              + 0.00382558411957381*dt*dfn &
              + -0.00254151525505059*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.16218596234777*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 17:
        fn1 = -0.998129868175228*fn1 &
              + 1.99813987107032*fn &
              + -1.00028950899976e-05*f0 &
              + 0.00382489018242778*dt*dfn &
              + -0.00253836774254408*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.183097754620724*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 18:
        fn1 = -0.997787970010627*fn1 &
              + 1.9978178224392*fn &
              + -2.985242857474e-05*f0 &
              + 0.00382427370874381*dt*dfn &
              + -0.00253526471506265*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.205257764288362*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 19:
        fn1 = -0.997481396066635*fn1 &
              + 1.9975259354053*fn &
              + -4.45393386597833e-05*f0 &
              + 0.00382371497115665*dt*dfn &
              + -0.00253216827471412*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.228662650410007*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 20:
        fn1 = -0.997201351522887*fn1 &
              + 1.99725693179725*fn &
              + -5.5580274360411e-05*f0 &
              + 0.00382320003760552*dt*dfn &
              + -0.00252905095601984*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.253308892697181*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 21:
        fn1 = -0.996941577073741*fn1 &
              + 1.99700557289402*fn &
              + -6.39958202792111e-05*f0 &
              + 0.00382271887999732*dt*dfn &
              + -0.00252589255233084*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.279192793389938*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 22:
        fn1 = -0.996697535885486*fn1 &
              + 1.99676802355418*fn &
              + -7.04876686956548e-05*f0 &
              + 0.00382226415700673*dt*dfn &
              + -0.00252267801450401*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.306310479221137*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 23:
        fn1 = -0.996465885656639*fn1 &
              + 1.99654143329167*fn &
              + -7.55476350310585e-05*f0 &
              + 0.0038218304121608*dt*dfn &
              + -0.00251939602766347*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.33465790346667*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 24:
        fn1 = -0.996244127955706*fn1 &
              + 1.99632365441664*fn &
              + -7.95264609337592e-05*f0 &
              + 0.00382141353429699*dt*dfn &
              + -0.00251603802744033*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.364230848079495*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 25:
        fn1 = -0.996030370522874*fn1 &
              + 1.99611304876287*fn &
              + -8.26782399957655e-05*f0 &
              + 0.00382101038759584*dt*dfn &
              + -0.00251259750739641*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.395024925905329*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 26:
        fn1 = -0.995823163157098*fn1 &
              + 1.99590835285766*fn &
              + -8.5189700563595e-05*f0 &
              + 0.00382061855348575*dt*dfn &
              + -0.00250906952325496*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.427035582977741*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 27:
        fn1 = -0.995621382501872*fn1 &
              + 1.99570858237387*fn &
              + -8.71998719977567e-05*f0 &
              + 0.00382023614774266*dt*dfn &
              + -0.00250545033257436*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.460258100890312*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 28:
        fn1 = -0.995424149919027*fn1 &
              + 1.99551296343908*fn &
              + -8.881352005202e-05*f0 &
              + 0.0038198616890003*dt*dfn &
              + -0.00250173712918076*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.494687599243495*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 29:
        fn1 = -0.995230772125691*fn1 &
              + 1.99532088259601*fn &
              + -9.01104703180183e-05*f0 &
              + 0.00381949400296314*dt*dfn &
              + -0.00249792784490238*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.530319038163722*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 30:
        fn1 = -0.995040697730879*fn1 &
              + 1.99513184990291*fn &
              + -9.11521720349449e-05*f0 &
              + 0.00381913215177221*dt*dfn &
              + -0.00249402099977041*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.567147220892262*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 31:
        fn1 = -0.99485348503399*fn1 &
              + 1.9949454714149*fn &
              + -9.19863809070311e-05*f0 &
              + 0.00381877538132819*dt*dfn &
              + -0.00249001558757054*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.605166796441267*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 32:
        fn1 = -0.99466877790394*fn1 &
              + 1.99476142844527*fn &
              + -9.26505413303426e-05*f0 &
              + 0.00381842308159289*dt*dfn &
              + -0.00248591098748475*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.644372262314417*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 33:
        fn1 = -0.994486287526112*fn1 &
              + 1.99457946178344*fn &
              + -9.3174257327099e-05*f0 &
              + 0.0038180747563788*dt*dfn &
              + -0.00248170689520068*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.684757967289498*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 34:
        fn1 = -0.994305778457965*fn1 &
              + 1.994399359575*fn &
              + -9.35811170338285e-05*f0 &
              + 0.00381773000014883*dt*dfn &
              + -0.00247740326869543*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.726318114260251*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 35:
        fn1 = -0.994127057881371*fn1 &
              + 1.99422094793469*fn &
              + -9.38900533153558e-05*f0 &
              + 0.00381738848004739*dt*dfn &
              + -0.0024730002851864*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.769046763134735*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 36:
        fn1 = -0.993949967249846*fn1 &
              + 1.99404408361776*fn &
              + -9.41163679188661e-05*f0 &
              + 0.00381704992187175*dt*dfn &
              + -0.00246849830665584*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.81293783378748*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 37:
        fn1 = -0.993774375746278*fn1 &
              + 1.9938686482554*fn &
              + -9.42725091174407e-05*f0 &
              + 0.00381671409903728*dt*dfn &
              + -0.00246389785201274*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.85798510906262*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 38:
        fn1 = -0.993600175121108*fn1 &
              + 1.99369454378814*fn &
              + -9.43686670321194e-05*f0 &
              + 0.0038163808238361*dt*dfn &
              + -0.0024591995744337*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.904182237825207*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 39:
        fn1 = -0.993427275591491*fn1 &
              + 1.9935216888244*fn &
              + -9.44132329048022e-05*f0 &
              + 0.00381604994046633*dt*dfn &
              + -0.00245440424277436*dt*df0
        call swap(fn, fn1)

        t = t0 + 0.951522738057872*dt
        f(:,:,:,1:mvar) = fn
        dfn = 0.
        call pde(f,dfn,p)

        ! Step n = 40:
        fn1 = -0.993255602562126*fn1 &
              + 1.99335001571812*fn &
              + -9.44131559979028e-05*f0 &
              + 0.0038157213194382*dt*dfn &
              + -0.0024495127262026*dt*df0
        call swap(fn, fn1)

        ! Done: last fn is the updated f:
        ! Increase time
        t = t0 + dt
        ! need explicit loop here, as f(:,:,:,1:mvar) = fn
        ! causes a `stack smashing' exception
        do iv=1,mvar
          f(:,:,:,iv) = fn(:,:,:,iv)
        enddo

    endsubroutine time_step
!***********************************************************************

    subroutine swap(a, b)
    !
    ! Swap two pointers
    !
    !  18-aug-08/perl: generated
    !

!        real, pointer :: a(:,:,:,:), b(:,:,:,:), tmp(:,:,:,:)

!        tmp => a
!        a   => b
!        b   => tmp

        real :: a(:,:,:,:), b(:,:,:,:), tmp(size(a,1), size(a,2), size(a,3), size(a,4))

        tmp = a
        a   = b
        b   = tmp



    endsubroutine swap
!***********************************************************************
    subroutine pushpars2c(p_par)
        
    use Messages, only: fatal_error

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

    call fatal_error('time_step_RKC-40','alpha_ts, beta_ts not defined')

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