subroutine time_iteration(irk)
    use commondata
    implicit none
    integer :: irk
    
    vtemp = v1
    if(irk==1) then
        do j = 0, jn
            do i = 0, im
                v1(:,i,j) = v0(:,i,j) - dt/dx*( fx(:,i,j)-fx(:,i-1,j) ) &
                    - dt/dy*( fy(:,i,j)-fy(:,i,j-1) )
                if(FLG_prob==1) then
                    v1(3,i,j) = v1(3,i,j) + dt*vtemp(1,i,j)
                    v1(4,i,j) = v1(4,i,j) + dt*vtemp(3,i,j)
                end if
            end do
        end do
    else if(irk==2) then
        do j = 0, jn
            do i = 0, im
                v1(:,i,j) = 0.25d0*( 3.d0*v0(:,i,j)  + v1(:,i,j) &
                    - dt/dx*( fx(:,i,j)-fx(:,i-1,j) ) - dt/dy*( fy(:,i,j)-fy(:,i,j-1) ) )
                if(FLG_prob==1) then
                    v1(3,i,j) = v1(3,i,j) + 0.25d0*dt*vtemp(1,i,j)
                    v1(4,i,j) = v1(4,i,j) + 0.25d0*dt*vtemp(3,i,j)
                end if
            end do
        end do        
    else if(irk==3) then
        do j = 0, jn
            do i = 0, im
                v1(:,i,j) = 1.d0/3.d0*( v0(:,i,j) + 2.d0*v1(:,i,j) &
                    - 2.d0*dt/dx*( fx(:,i,j)-fx(:,i-1,j) ) - 2.d0*dt/dy*( fy(:,i,j)-fy(:,i,j-1) ) )
                if(FLG_prob==1) then
                    v1(3,i,j) = v1(3,i,j) + 2.d0/3.d0*dt*vtemp(1,i,j)
                    v1(4,i,j) = v1(4,i,j) + 2.d0/3.d0*dt*vtemp(3,i,j)
                end if
            end do
        end do
    end if
    
end subroutine time_iteration