对于考虑重力源项的二维无黏欧拉方程如何利用二阶MUSCL进行空间离散，利用二阶runge kutta方法进行时间离散，请给出fortran代码示例

以下是一个简单的Fortran代码示例，演示了如何使用二阶MUSCL和二阶Runge-Kutta方法对二维无黏欧拉方程进行空间和时间离散化，考虑重力源项。

program main

implicit none

! Parameters
integer, parameter :: nx = 100, ny = 100 ! Number of grid points in x and y directions
real, parameter :: lx = 1.0, ly = 1.0 ! Domain size
real, parameter :: dx = lx/nx, dy = ly/ny ! Grid spacing
real, parameter :: cfl = 0.5 ! CFL number
real, parameter :: gamma = 1.4 ! Ratio of specific heats
real, parameter :: g = 9.81 ! Acceleration due to gravity

! Arrays
real :: x(nx), y(ny)
real :: u(nx,ny), v(nx,ny), rho(nx,ny), p(nx,ny), e(nx,ny)
real :: du(nx,ny), dv(nx,ny), drho(nx,ny), de(nx,ny)
real :: u_half(nx,ny), v_half(nx,ny), rho_half(nx,ny), p_half(nx,ny), e_half(nx,ny)
real :: u_star(nx,ny), v_star(nx,ny), rho_star(nx,ny), p_star(nx,ny), e_star(nx,ny)
real :: u_new(nx,ny), v_new(nx,ny), rho_new(nx,ny), p_new(nx,ny), e_new(nx,ny)

! Initialize arrays
do i = 1, nx
    x(i) = (i-0.5)*dx
end do
do j = 1, ny
    y(j) = (j-0.5)*dy
end do
do i = 1, nx
    do j = 1, ny
        rho(i,j) = 1.0 ! Density
        p(i,j) = 1.0 ! Pressure
        u(i,j) = 0.0 ! x-velocity
        v(i,j) = 0.0 ! y-velocity
        e(i,j) = p(i,j)/(gamma-1.0)/rho(i,j) + 0.5*(u(i,j)**2 + v(i,j)**2) ! Total energy
    end do
end do

! Time integration
do tstep = 1, 1000

    ! Compute time step
    dt = cfl*min(dx/u + dy/v + sqrt(gamma*p/rho)/sqrt(u**2+v**2))

    ! First stage of Runge-Kutta
    do i = 2, nx-1
        do j = 2, ny-1

            ! Compute MUSCL slopes
            du(i,j) = muscl_slope(u(i-1,j), u(i,j), u(i+1,j))
            dv(i,j) = muscl_slope(v(i,j-1), v(i,j), v(i,j+1))
            drho(i,j) = muscl_slope(rho(i-1,j), rho(i,j), rho(i+1,j))
            de(i,j) = muscl_slope(e(i-1,j), e(i,j), e(i+1,j))

            ! Compute half-time values
            u_half(i,j) = u(i,j) - 0.5*dt/dx*(u(i,j)-u(i-1,j)) - 0.5*dt/dy*(v(i,j)-v(i,j-1))
            v_half(i,j) = v(i,j) - 0.5*dt/dx*(u(i,j)-u(i-1,j)) - 0.5*dt/dy*(v(i,j)-v(i,j-1))
            rho_half(i,j) = rho(i,j) - 0.5*dt/dx*(rho(i,j)*u(i,j)-rho(i-1,j)*u(i-1,j)) - 0.5*dt/dy*(rho(i,j)*v(i,j)-rho(i,j-1)*v(i,j-1))
            p_half(i,j) = p(i,j) - 0.5*dt/dx*(u(i,j)*(p(i,j)+rho(i,j)*u(i,j))-u(i-1,j)*(p(i-1,j)+rho(i-1,j)*u(i-1,j))) - 0.5*dt/dy*(v(i,j)*(p(i,j)+rho(i,j)*v(i,j))-v(i,j-1)*(p(i,j-1)+rho(i,j-1)*v(i,j-1)))
            e_half(i,j) = e(i,j) - 0.5*dt/dx*(u(i,j)*(e(i,j)+p(i,j))-u(i-1,j)*(e(i-1,j)+p(i-1,j))) - 0.5*dt/dy*(v(i,j)*(e(i,j)+p(i,j))-v(i,j-1)*(e(i,j-1)+p(i,j-1)))

            ! Compute star values
            u_star(i,j) = u_half(i,j) - dt/dx*(p_half(i+1,j)-p_half(i,j)) / (0.5*(rho_half(i+1,j)+rho_half(i,j)))
            v_star(i,j) = v_half(i,j) - dt/dy*(p_half(i,j+1)-p_half(i,j)) / (0.5*(rho_half(i,j+1)+rho_half(i,j)))
            rho_star(i,j) = rho_half(i,j) - dt/dx*(rho_half(i+1,j)*u_star(i+1,j)-rho_half(i,j)*u_star(i,j)) - dt/dy*(rho_half(i,j+1)*v_star(i,j+1)-rho_half(i,j)*v_star(i,j))
            p_star(i,j) = p_half(i,j) - dt/dx*(u_star(i+1,j)*(p_half(i+1,j)+rho_half(i+1,j)*u_star(i+1,j))-u_star(i,j)*(p_half(i,j)+rho_half(i,j)*u_star(i,j))) - dt/dy*(v_star(i,j+1)*(p_half(i,j+1)+rho_half(i,j+1)*v_star(i,j+1))-v_star(i,j)*(p_half(i,j)+rho_half(i,j)*v_star(i,j)))
            e_star(i,j) = e_half(i,j) - dt/dx*(u_star(i+1,j)*(e_half(i+1,j)+p_half(i+1,j))-u_star(i,j)*(e_half(i,j)+p_half(i,j))) - dt/dy*(v_star(i,j+1)*(e_half(i,j+1)+p_half(i,j+1))-v_star(i,j)*(e_half(i,j)+p_half(i,j)))

        end do
    end do

    ! Second stage of Runge-Kutta
    do i = 2, nx-1
        do j = 2, ny-1

            ! Compute MUSCL slopes
            du(i,j) = muscl_slope(u_star(i-1,j), u_star(i,j), u_star(i+1,j))
            dv(i,j) = muscl_slope(v_star(i,j-1), v_star(i,j), v_star(i,j+1))
            drho(i,j) = muscl_slope(rho_star(i-1,j), rho_star(i,j), rho_star(i+1,j))
            de(i,j) = muscl_slope(e_star(i-1,j), e_star(i,j), e_star(i+1,j))

            ! Compute new values
            u_new(i,j) = 0.5*(u(i,j)+u_star(i,j)) - 0.5*dt/dx*(p_star(i,j+1)-p_star(i,j)) / (0.5*(rho_star(i,j+1)+rho_star(i,j))) - 0.5*dt/dy*(p_star(i+1,j)-p_star(i,j)) / (0.5*(rho_star(i+1,j)+rho_star(i,j)))
            v_new(i,j) = 0.5*(v(i,j)+v_star(i,j)) - 0.5*dt/dx*(p_star(i,j+1)-p_star(i,j)) / (0.5*(rho_star(i,j+1)+rho_star(i,j))) - 0.5*dt/dy*(p_star(i+1,j)-p_star(i,j)) / (0.5*(rho_star(i+1,j)+rho_star(i,j)))
            rho_new(i,j) = 0.5*(rho(i,j)+rho_star(i,j)) - 0.5*dt/dx*(rho_star(i,j+1)*u_star(i,j+1)-rho_star(i,j)*u_star(i,j)) - 0.5*dt/dy*(rho_star(i+1,j)*v_star(i+1,j)-rho_star(i,j)*v_star(i,j))
            p_new(i,j) = 0.5*(p(i,j)+p_star(i,j)) - 0.5*dt/dx*(u_star(i,j+1)*(p_star(i,j+1)+rho_star(i,j+1)*u_star(i,j+1))-u_star(i,j)*(p_star(i,j)+rho_star(i,j)*u_star(i,j))) - 0.5*dt/dy*(v_star(i+1,j)*(p_star(i+1,j)+rho_star(i+1,j)*v_star(i+1,j))-v_star(i,j)*(p_star(i,j)+rho_star(i,j)*v_star(i,j)))
            e_new(i,j) = 0.5*(e(i,j)+e_star(i,j)) - 0.5*dt/dx*(u_star(i,j+1)*(e_star(i,j+1)+p_star(i,j+1))-u_star(i,j)*(e_star(i,j)+p_star(i,j))) - 0.5*dt/dy*(v_star(i+1,j)*(e_star(i+1,j)+p_star(i+1,j))-v_star(i,j)*(e_star(i,j)+p_star(i,j)))

        end do
    end do

    ! Update arrays
    u = u_new
    v = v_new
    rho = rho_new
    p = p_new
    e = e_new

    ! Apply boundary conditions
    do j = 1, ny
        u(1,j) = 0.0
        v(1,j) = 0.0
        rho(1,j) = rho(2,j)
        p(1,j) = p(2,j)
        e(1,j) = e(2,j)
        u(nx,j) = 0.0
        v(nx,j) = 0.0
        rho(nx,j) = rho(nx-1,j)
        p(nx,j) = p(nx-1,j)
        e(nx,j) = e(nx-1,j)
    end do
    do i = 1, nx
        u(i,1) = 0.0
        v(i,1) = 0.0
        rho(i,1) = rho(i,2)
        p(i,1) = p(i,2)
        e(i,1) = e(i,2)
        u(i,ny) = 0.0
        v(i,ny) = 0.0
        rho(i,ny) = rho(i,ny-1)
        p(i,ny) = p(i,ny-1)
        e(i,ny) = e(i,ny-1)
    end do

end do

contains

! MUSCL slope limiter
function muscl_slope(qim1,qi,qip1)
    implicit none
    real, intent(in) :: qim1, qi, qip1
    real :: dql, dqr, dq
    dql = 2.0*(qi-qim1)
    dqr = 2.0*(qip1-qi)
    dq = min(abs(dql),abs(dqr))
    if (dql*dqr < 0.0) then
        muscl_slope = 0.0
    else
        muscl_slope = sign(dq,dql+dqr)
    end if
end function muscl_slope

end program main