subroutine time_step
#define DEBUG_LEVEL 1
!  Routine to take an explicit time step via Euler,
!  2nd or third order Adams Bashforth method. The differential
!  equations are as follows:
!
!    d/dt u   = rhs1
!    d/dt v   = rhs2
!    d/dt w   = rhs3
!    d/dt s1  = rhs4
!    d/dt s2  = rhs5
!
!  notation:
!    MM0   integer index identifying the time t    storage locations
!    MM1   integer index identifying the time t-dt storage locations
!    MM2   integer index identifying the time t-2*dt storage locations
!    rhs1(:,:,:,MM0),rhs1(:,:,:,MM1),rhs1(:,:,:,MM2)  rhs1 at t,t-dt,t-2*dt
!    etc for rhs2,rhs3,rhs4,rhs5
!    dt   time increment
!
!  outputs:
!    u_star(1:nx,1:ny,1:nz)  calculated values of u at time t+dt
!    v_star(1:nx,1:ny,1:nz)  calculated values of v at time t+dt
!    w_star(1:nx,1:ny,1:nz)  calculated values of w at time t+dt
!    s1(1:nx,1:ny,1:nz)      calculated values of s1 at time t+dt (overwrite)
!    s2(1:nx,1:ny,1:nz)      calculated values of s2 at time t+dt (overwrite)

 use grid_info, only: dt
 use dependent_variables
 use intermediate_variables
 use rhs_variables
 use counters_flags_etc
 use mpi_parameters
 implicit none

 include '../input/problem_size.h'
 real dt2,dt12,dt24
#if DEBUG_LEVEL >= 1
       if(myid==0) write(0,*) 'hello world from time_step'
#endif

 dt2=dt/2.
 dt12=dt/12.
 dt24=dt/24.

 if( step_flag == 'EULER') then
  u_star(:,:,:)= u(:,:,:)+dt*rhs1(:,:,:,MM0)
  v_star(:,:,:)= v(:,:,:)+dt*rhs2(:,:,:,MM0)
  w_star(:,:,:)= w(:,:,:)+dt*rhs3(:,:,:,MM0)
  s1(:,:,:)    =s1(:,:,:)+dt*rhs4(:,:,:,MM0)
  if(num_scalars>1) s2(:,:,:)=s2(:,:,:)+dt*rhs5(:,:,:,MM0)
  step_flag='AB2'
 elseif( step_flag == 'AB2' ) then
  u_star(:,:,:)= u(:,:,:)+dt2*(3.*rhs1(:,:,:,MM0)-rhs1(:,:,:,MM1))
  v_star(:,:,:)= v(:,:,:)+dt2*(3.*rhs2(:,:,:,MM0)-rhs2(:,:,:,MM1))
  w_star(:,:,:)= w(:,:,:)+dt2*(3.*rhs3(:,:,:,MM0)-rhs3(:,:,:,MM1))
  s1(:,:,:)    =s1(:,:,:)+dt2*(3.*rhs4(:,:,:,MM0)-rhs4(:,:,:,MM1))
  if(num_scalars>1) s2(:,:,:)=s2(:,:,:)+dt2*(3.*rhs5(:,:,:,MM0)-rhs5(:,:,:,MM1))
  if( AB_order >= 3 ) step_flag='AB3'
 elseif( step_flag == 'AB3' ) then
  u_star= u+dt12*(23.*rhs1(:,:,:,MM0)-16.*rhs1(:,:,:,MM1)+5.*rhs1(:,:,:,MM2))
  v_star= v+dt12*(23.*rhs2(:,:,:,MM0)-16.*rhs2(:,:,:,MM1)+5.*rhs2(:,:,:,MM2))
  w_star= w+dt12*(23.*rhs3(:,:,:,MM0)-16.*rhs3(:,:,:,MM1)+5.*rhs3(:,:,:,MM2))
  s1    =s1+dt12*(23.*rhs4(:,:,:,MM0)-16.*rhs4(:,:,:,MM1)+5.*rhs4(:,:,:,MM2))
  if(num_scalars>1) s2=s2+dt12*(23.*rhs5(:,:,:,MM0)-16.*rhs5(:,:,:,MM1)+5.*rhs5(:,:,:,MM2))
  if( AB_order == 4 ) step_flag='AB4'
 elseif( step_flag == 'AB4' ) then
  u_star= u+dt24*(55.*rhs1(:,:,:,MM0)-59.*rhs1(:,:,:,MM1)+37.*rhs1(:,:,:,MM2)-9.*rhs1(:,:,:,MM3))
  v_star= v+dt24*(55.*rhs2(:,:,:,MM0)-59.*rhs2(:,:,:,MM1)+37.*rhs2(:,:,:,MM2)-9.*rhs2(:,:,:,MM3))
  w_star= w+dt24*(55.*rhs3(:,:,:,MM0)-59.*rhs3(:,:,:,MM1)+37.*rhs3(:,:,:,MM2)-9.*rhs3(:,:,:,MM3))
  s1    =s1+dt24*(55.*rhs4(:,:,:,MM0)-59.*rhs4(:,:,:,MM1)+37.*rhs4(:,:,:,MM2)-9.*rhs4(:,:,:,MM3))
  if(num_scalars>1) s2=s2+dt24*(55.*rhs5(:,:,:,MM0)-59.*rhs5(:,:,:,MM1)+37.*rhs5(:,:,:,MM2)-9.*rhs5(:,:,:,MM3))
 endif

end subroutine time_step