subroutine timestep(s,par)
use params
use spaceparams

IMPLICIT NONE

type(spacepars)			    :: s
type(parameters)		    :: par

integer                     :: i
integer                     :: j
integer                     :: itheta

! Robert new time step criterion
if (par%t==0.0d0) then		! conservative estimate
	par%dt    = par%CFL*min(minval(s%xz(2:s%nx+1)-s%xz(1:s%nx)),minval(s%yz(2:s%ny+1)-s%yz(1:s%ny)))/(maxval(s%hh)*par%g)
else
	par%dt=huge(0.0d0)	! Seed dt
	do itheta=1,s%ntheta
		do j=2,s%ny
			do i=2,s%nx					
				! dx dy direction
				!	par%dt=min(par%dt,par%CFL*0.5d0/sqrt(&
				!	(par%g*s%hu(i,j)+max(s%uu(i,j)**2+s%vv(i,j)**2,s%ueu(i,j)**2+s%vev(i,j)**2))&
				!   *(1.d0/(s%xu(i)-s%xu(i-1))**2.d0+1.d0/(s%yv(j)-s%yv(j-1))**2.d0)))
				    par%dt=min(par%dt,par%CFL*0.5d0/sqrt(&
					(par%g*s%hh(i,j)+abs(s%u(i,j))+abs(s%v(i,j)))&
				   *(1.d0/(s%xu(i)-s%xu(i-1))**2.d0+1.d0/(s%yv(j)-s%yv(j-1))**2.d0)))

				! Theta direction
					par%dt=min(par%dt,par%CFL*s%dtheta/max(abs(s%ctheta(i,j,itheta)),tiny(0.0d0)))
			end do
		end do
	end do
	
end if

par%t=par%t+par%dt;
if(par%t>=par%tnext) then
    par%dt=par%dt-(par%t-par%tnext)
    par%t=par%tnext
    par%tnext=par%tnext+par%tint
end if  

end subroutine