program main
    use commondata
    implicit none
    
    integer :: irk, npointer
    real*8 :: time_start, time_finish, time_cost
    real*8 :: dtx, dty, wavolcx, wavolcy
    real*8 :: fluxp, fluxm, fluyp, fluym
    
    !open(404, file="error.dat")
    FLG_continue = 0
    if(FLG_continue==0) then
        tn = 0.d0
        call initial
        call boundary
    else if(FLG_continue==1) then
        open(111, file="Euler.dat")
        read(111,*) tn
        do k = 1, kn
            do i = -4, im+4
                do j = -4, jn+4
                    read(111,*) v1(k,i,j)
                end do
            end do
        end do
        close(111)
    end if
    var(1,:,:) = v1(1,:,:)
    var(2,:,:) = v1(2,:,:)/v1(1,:,:)
    var(3,:,:) = v1(3,:,:)/v1(1,:,:)
    var(4,:,:) = (gamma-1.d0)*( v1(4,:,:) - 0.5d0*var(1,:,:)*( var(2,:,:)**2 + var(3,:,:)**2 ) )
    var(5,:,:) = sqrt( gamma*var(4,:,:)/var(1,:,:) )
    
    nrk = 3
    npointer = 0
    call CPU_TIME(time_start)
    do while( tn < t_end - eps )
        v0 = v1
        do irk = 1, nrk
            !  Flux splitting
            if(FLG_split==1) then       !  Local L-F
                call flux_splitting_LLF(wavolcx, wavolcy)
            else if(FLG_split==2) then  !  Global L-F
                call flux_splitting_GLF(wavolcx, wavolcy)
            else if(FLG_split==3) then  !  Steger Warming
                call flux_splitting_SW(1, fxp, fxm, wavolcx)
                call flux_splitting_SW(2, fyp, fym, wavolcy)
            else if(FLG_split==4) then  !  VanLeer
                call flux_splitting_VL(1, fxp, fxm, wavolcx)
                call flux_splitting_VL(2, fyp, fym, wavolcy)                
            end if
            if(irk == 1) then
                dtx = dx/wavolcx
                dty = dy/wavolcy
                dt = CFL*dtx*dty/(dtx+dty)
                if( (tn+dt)>t_end ) dt = t_end - tn
                npointer = npointer + 1
                tn  = tn + dt
                write(*,*) npointer, tn, dt
            end if
            do j = -1, jn
                do i = -1, im
                    do k = 1, kn
                        if(FLG_order==5) then   ! 5th order
                            if(FLG_recon==1) then       !  WENOJS_5
                                call weno5_p( fxp(k,i-2:i+2,j), fluxp )
                                call weno5_m( fxm(k,i-1:i+3,j), fluxm )
                                call weno5_p( fyp(k,i,j-2:j+2), fluyp )
                                call weno5_m( fym(k,i,j-1:j+3), fluym )
                            else if(FLG_recon==2) then  !  WENOz_5
                                call wenoz5_p( fxp(k,i-2:i+2,j), fluxp )
                                call wenoz5_m( fxm(k,i-1:i+3,j), fluxm )
                                call wenoz5_p( fyp(k,i,j-2:j+2), fluyp )
                                call wenoz5_m( fym(k,i,j-1:j+3), fluym )
                            else if(FLG_recon==3) then  !  WENOzn_5
                                call wenozn5_p( fxp(k,i-2:i+2,j), fluxp )
                                call wenozn5_m( fxm(k,i-1:i+3,j), fluxm )
                                call wenozn5_p( fyp(k,i,j-2:j+2), fluyp )
                                call wenozn5_m( fym(k,i,j-1:j+3), fluym )
                            end if
                        else if(FLG_order==7) then   ! 7th order
                            !if(FLG_recon==1) then       !  WENOJS_7
                            !    call weno7_p( fxp(k,i-3:i+3,j), fluxp )
                            !    call weno7_m( fxm(k,i-2:i+4,j), fluxm )
                            !    call weno7_p( fyp(k,i,j-3:j+3), fluyp )
                            !    call weno7_m( fym(k,i,j-2:j+4), fluym )
                            !else if(FLG_recon==2) then  !  WENOz_7
                            !    call wenoz7_p( fxp(k,i-3:i+3,j), fluxp )
                            !    call wenoz7_m( fxm(k,i-2:i+4,j), fluxm )
                            !    call wenoz7_p( fyp(k,i,j-3:j+3), fluyp )
                            !    call wenoz7_m( fym(k,i,j-2:j+4), fluym )
                            !else if(FLG_recon==3) then  !  WENOzn_7
                            !    call wenozn7_p( fxp(k,i-3:i+3,j), fluxp )
                            !    call wenozn7_m( fxm(k,i-2:i+4,j), fluxm )
                            !    call wenozn7_p( fyp(k,i,j-3:j+3), fluyp )
                            !    call wenozn7_m( fym(k,i,j-2:j+4), fluym )
                            !end if
                            if(FLG_recon==1) then       !  WENOJS_7
                                call weno7balsara_p( fxp(k,i-3:i+3,j), fluxp )
                                call weno7balsara_m( fxm(k,i-2:i+4,j), fluxm )
                                call weno7balsara_p( fyp(k,i,j-3:j+3), fluyp )
                                call weno7balsara_m( fym(k,i,j-2:j+4), fluym )
                            else if(FLG_recon==2) then  !  WENOz_7
                                call wenoz7balsara_p( fxp(k,i-3:i+3,j), fluxp )
                                call wenoz7balsara_m( fxm(k,i-2:i+4,j), fluxm )
                                call wenoz7balsara_p( fyp(k,i,j-3:j+3), fluyp )
                                call wenoz7balsara_m( fym(k,i,j-2:j+4), fluym )
                            else if(FLG_recon==3) then  !  WENOzn_7
                                call wenozn7balsara_p( fxp(k,i-3:i+3,j), fluxp )
                                call wenozn7balsara_m( fxm(k,i-2:i+4,j), fluxm )
                                call wenozn7balsara_p( fyp(k,i,j-3:j+3), fluyp )
                                call wenozn7balsara_m( fym(k,i,j-2:j+4), fluym )
                            end if
                        end if
                        fx(k,i,j) = fluxp + fluxm
                        fy(k,i,j) = fluyp + fluym
                    end do
                end do
            end do
            call time_iteration(irk)
            call boundary
            var(1,:,:) = v1(1,:,:)
            var(2,:,:) = v1(2,:,:)/v1(1,:,:)
            var(3,:,:) = v1(3,:,:)/v1(1,:,:)
            var(4,:,:) = (gamma-1.d0)*( v1(4,:,:) - 0.5d0*var(1,:,:)*( var(2,:,:)**2 + var(3,:,:)**2 ) )
            where( var(4,:,:) <  epsilon )  var(4,:,:) = epsilon
            !var(4,:,:) = max( var(4,:,:), epsilon )
            var(5,:,:) = sqrt( gamma*var(4,:,:)/var(1,:,:) )
            
            !do i = -4, im+4
            !    do j = -4, jn+4
            !        if( var(1,i,j)<0 .or. var(4,i,j)<0 ) then
            !            write(404,*) i, j, var(1,i,j), var(4,i,j)
            !            pause
            !        end if
            !    end do
            !end do
                    
        end do
    end do
    call CPU_TIME(time_finish)
    time_cost = time_finish - time_start
    write(*,*) "Computation finished."
    write(*,*) "The time cost by computation is", time_cost, "second."
    call output
    
end program main