module mod_dt_hd

  use mod_variables
  use mod_physics_vars
  implicit none
  private

  public :: setdt

contains

! instantiate here for inlining in kernels
pure subroutine to_primitive(u)
  
  real(dp), intent(inout) :: u(nw_phys)

  ! Compute velocity from momentum
      u(iw_mom(1)) = u(iw_mom(1))/u(iw_rho)
      u(iw_mom(2)) = u(iw_mom(2))/u(iw_rho)
      u(iw_mom(3)) = u(iw_mom(3))/u(iw_rho)

  ! Compute pressure from energy
  u(iw_e) = gamma_1 * (u(iw_e) - 0.5_dp * u(iw_rho) * &
     sum(u(iw_mom(1:ndim))**2) )

end subroutine to_primitive
pure real(dp) function get_cmax(u, x, flux_dim) result(wC)
  
  real(dp), intent(in)  :: u(nw_phys)
  real(dp), intent(in)  :: x(1:ndim)
  integer, intent(in)   :: flux_dim

  wC = sqrt(hd_gamma * u(iw_e) / u(iw_rho)) + abs(u(iw_mom(flux_dim)))

end function get_cmax
  subroutine phys_get_dt(w, x, dx, dtnew)
  use mod_usr, only: gravity_field
  
    real(dp), intent(in)   :: w(nw_phys), x(1:ndim), dx(1:ndim)
    real(dp), intent(out)  :: dtnew
    ! .. local ..
    integer                :: idim
    real(dp)               :: field

    dtnew = huge(1.0d0)
    
    do idim = 1, ndim
       field = gravity_field(w, x, idim)
       field = max( abs(field), epsilon(1.0d0) )
       dtnew = min( dtnew, 1_dp / sqrt( field/dx(idim) ) )
    end do
    
  end subroutine phys_get_dt

  !>setdt  - set dt for all levels between levmin and levmax.
  !>         dtpar>0  --> use fixed dtpar for all level
  !>         dtpar<=0 --> determine CFL limited timestep
  subroutine setdt()
    use mod_global_parameters

    integer :: iigrid, igrid, idims, ix1,ix2,ix3, ifile
    double precision :: dtmin_mype, factor, dx1,dx2,dx3, dtmax

    double precision :: dxinv(1:ndim), cmaxtot, cmax, u(1:nw_phys), cmax_mype
    double precision :: xloc(1:ndim), qdtnew

    if (dtpar<=zero) then
       dtmin_mype=bigdouble

       
       do iigrid=1,igridstail_active; igrid=igrids_active(iigrid)

          dx1=rnode(rpdx1_,igrid);dx2=rnode(rpdx2_,igrid)
          dx3=rnode(rpdx3_,igrid);

          dxinv(1)=one/dx1;dxinv(2)=one/dx2;dxinv(3)=one/dx3;

          
          do ix3=ixMlo3,ixMhi3 
             do ix2=ixMlo2,ixMhi2 
                do ix1=ixMlo1,ixMhi1 
                   u = bg(1)%w(ix1, ix2, ix3, 1:nw_phys, igrid)
                   call to_primitive(u)
                   xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid)


                   cmaxtot = 0.0d0
                   
                   do idims = 1, ndim
                      cmax = get_cmax(u, xloc, idims)
                      cmaxtot = cmaxtot + cmax * dxinv(idims)
                   end do
                   dtmin_mype     = min( dtmin_mype, courantpar / cmaxtot )

                   u = bg(1)%w(ix1,ix2,ix3,1:nw_phys,igrid)
                   xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid)
                   call phys_get_dt(u, xloc, one/dxinv(1:ndim), qdtnew)
                   dtmin_mype = min( dtmin_mype, qdtnew )

                end do
             end do
          end do


       end do

    else
       dtmin_mype=dtpar
    end if

    
    if (need_global_cmax) then
       cmax_mype=-bigdouble

       
       do iigrid=1,igridstail_active; igrid=igrids_active(iigrid)

          
          do ix3=ixMlo3,ixMhi3 
             do ix2=ixMlo2,ixMhi2 
                do ix1=ixMlo1,ixMhi1

                   u(1:nw_phys) = bg(1)%w(ix1, ix2, ix3, 1:nw_phys, igrid)
                   call to_primitive(u)
                   xloc(1:ndim) = bgeo%x(ix1, ix2, ix3, 1:ndim, igrid)
                   
                   do idims = 1, ndim
                      cmax = get_cmax(u, xloc, idims)
                      cmax_mype = max( cmax_mype, cmax )
                   end do

                end do
             end do
          end do
       end do
    end if


    if (dtmin_mype<dtmin) then
       write(unitterm,*)"Error: Time step too small!", dtmin_mype
       write(unitterm,*)"on processor:", mype, "at time:", global_time,&
          " step:", it
       write(unitterm,*)"Lower limit of time step:", dtmin
       crash=.true.
    end if

    if (slowsteps>it-it_init+1) then
       factor=one-(one-dble(it-it_init+1)/dble(slowsteps))**2
       dtmin_mype=dtmin_mype*factor
    end if

    if(final_dt_reduction)then
       if(time_max-global_time<=dtmin) then
          final_dt_exit=.true.
       endif
       dtmin_mype=min(dtmin_mype,time_max-global_time)
    end if

    if (dtpar<=zero) then
       call MPI_ALLREDUCE(dtmin_mype,dt,1,MPI_DOUBLE_PRECISION,MPI_MIN, icomm,&
          ierrmpi)
    else
       dt=dtmin_mype
    end if

    if (need_global_cmax) then
      call MPI_ALLREDUCE(cmax_mype, cmax_global, 1, MPI_DOUBLE_PRECISION, MPI_MAX, icomm, &
             ierrmpi)
      
    end if

    if(any(dtsave(1:nfile)<bigdouble).or.any(tsave(isavet(1:nfile),&
       1:nfile)<bigdouble))then
       dtmax = minval(ceiling(global_time/dtsave(1:nfile))*dtsave(1:nfile))-&
          global_time
       do ifile=1,nfile
          dtmax = min(tsave(isavet(ifile),ifile)-global_time,dtmax)
       end do
       if(dtmax > smalldouble)then
          dt=min(dt,dtmax)
       else
          ! dtmax=0 means dtsave is divisible by global_time
          dt=min(dt,minval(dtsave(1:nfile)))
       end if
    end if

    if(mype==0) then
       if(any(dtsave(1:nfile)<dt)) then
          write(unitterm,*) 'Warning: timesteps: ',dt,&
             ' exceeding output intervals ', dtsave(1:nfile)
       endif
    endif

    
  end subroutine setdt
end module mod_dt_hd