!------------------------------------------------------------------------

! this is a fortran code for 2d WENO_TVD algorithm by Yu Wang, SPA, SYSU
! Contact wangy665@mail2.sysu.edu.cn if necessary

!------------------------------------------------------------------------


function step(value) ! step function, used in calculation of source term
    implicit none 
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: value,step
    if(value > 0) then
        step = 1
    else 
        step = 0
    end if
end function

function logic2int(logic)  !  boole  transformed into int 
    implicit none
    logical :: logic
    integer :: logic2int
    if(logic) then
        logic2int = 1
    else 
        logic2int = 0
    end if 
end function

Function minus_r(v0,v1,v2,n) ! WENO reconstruction
    implicit none 
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n,i,j
    real(kind = DP) :: minus_r(n)
    Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
    Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
    
    v_0 = 1/2.0_DP*v1 + 1/2.0_DP*v2
    v_1 = -1/2.0_DP*v0 + 3/2.0_DP*v1

    v(1,:) = v_0
    v(2,:) = v_1
    !weight
    d = (/2.0_DP/3,1.0_DP/3/)
    epsn = 1e-6_DP
    beta(1,:) = (v2 - v1)**2
    beta(2,:) = (v1 - v0)**2

    do i = 1,2
        do j = 1,n
            alpha(i,j) = d(i) / (epsn + beta(i,j))**2
        end do 
    end do 

    do i = 1,2
        do j = 1,n
            w(i,j) = alpha(i,j) / (alpha(1,j) + alpha(2,j))
        end do 
    end do
    minus_r = .0_DP
    do i = 1,n
        do j = 1,2
            minus_r(i) = minus_r(i) + v(j,i)*w(j,i)
        end do 
    end do
end function

Function plus_l(v0,v1,v2,n) ! WENO reconstruction
    implicit none 
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n,i,j
    real(kind = DP) :: plus_l(n)
    Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
    Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
    
    v_0 = 3/2.0_DP*v1 - 1/2.0_DP*v2
    v_1 = 1/2.0_DP*v0 + 1/2.0_DP*v1

    v(1,:) = v_0
    v(2,:) = v_1
    !weight
    d = (/1.0_DP/3,2.0_DP/3/)
    epsn = 1e-6_DP
    beta(1,:) = (v2 - v1)**2
    beta(2,:) = (v1 - v0)**2

    do i = 1,2
        do j = 1,n
            alpha(i,j) = d(i) / (epsn + beta(i,j))**2
        end do 
    end do 

    do i = 1,2
        do j = 1,n
            w(i,j) = alpha(i,j) / (alpha(1,j) + alpha(2,j))
        end do 
    end do
    plus_l = .0_DP
    do i = 1,n
        do j = 1,2
            plus_l(i) = plus_l(i) + v(j,i)*w(j,i)
        end do 
    end do
end function

Function absF(u,v,c,u_l,u_r,c_l,c_r)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: i
    real(kind = DP) :: absF(3,3)
    real(kind = DP), intent(IN) :: u,v,c,u_l,u_r,c_l,c_r
    real(kind = DP) :: diag(3,3) = 0, R(3,3),L(3,3),delta(3)

    delta(1) = max(.0_DP,u-c-(u_l-c_l),u_r-c_r-(u-c))
    delta(2) = max(.0_DP,u-u_l,u_r-u)
    delta(3) = max(.0_DP,u+c-(u_l+c_l),u_r+c_r-(u+c))

    diag(1,1) = abs(u-c)
    diag(2,2) = abs(u)
    diag(3,3) = abs(u+c)

    do i = 1,3
        if (diag(i,i) < delta(i)) then
            ! print *,'modifyF'
            diag(i,i) = 1.0_DP/2*(diag(i,i)**2 / delta(i) + delta(i))
        end if 
    end do

    R(1,:) = (/1.0_DP,.0_DP,1.0_DP/)
    R(2,:) = (/u-c,.0_DP,u+c/)
    R(3,:) = (/v,c,v/)

    L(1,:) = 1.0_DP /(2.0_DP*c)*(/u+c,-1.0_DP,.0_DP/)
    L(2,:) = (/-v/c,.0_DP,1.0_DP/c/)
    L(3,:) = 1.0_DP /(2.0_DP*c)*(/c-u,1.0_DP,.0_DP/)

    absF = matmul(matmul(R,diag),L)
end function

Function absG(u,v,c,v_l,v_r,c_l,c_r)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: i
    real(kind = DP) :: absG(3,3)
    real(kind = DP), intent(IN) :: u,v,c,v_l,v_r,c_l,c_r
    real(kind = DP) :: diag(3,3) = 0, R(3,3),L(3,3),delta(3)

    delta(1) = max(.0_DP,v-c-(v_l-c_l),v_r-c_r-(v-c))
    delta(2) = max(.0_DP,v-v_l,v_r-v)
    delta(3) = max(.0_DP,v+c-(v_l+c_l),v_r+c_r-(v+c))

    diag(1,1) = abs(v-c)
    diag(2,2) = abs(v)
    diag(3,3) = abs(v+c)

    do i = 1,3
        if (diag(i,i) < delta(i)) then
            ! print *,'modifyG'
            diag(i,i) = 1.0_DP/2*(diag(i,i)**2 / delta(i) + delta(i))
        end if 
    end do

    R(1,:) = (/1.0_DP,.0_DP,1.0_DP/)
    R(2,:) = (/u,-c,u/)
    R(3,:) = (/v-c,.0_DP,v+c/)

    L(1,:) = 1.0_DP /(2.0_DP*c)*(/v+c,.0_DP,-1.0_DP/)
    L(2,:) = (/u/c,-1.0_DP/c,.0_DP/)
    L(3,:) = 1.0_DP /(2.0_DP*c)*(/c-v,.0_DP,1.0_DP/)

    absG = matmul(matmul(R,diag),L)
end function

function func1(F,g)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: F(3),g,func1(3)
    func1(1) = F(2)
    func1(2) = F(2)**2 / F(1) + g*F(1)**2/2.0_DP
    func1(3) = F(3)*F(2)/F(1)
end function

function func2(G,g1)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: G(3),g1,func2(3)
    func2(1) = G(3)
    func2(2) = G(3)*G(2)/G(1)  
    func2(3) = G(3)**2 / G(1) + g1*G(1)**2/2.0_DP
end function

Function roe_bar(h_l,h_r,u_l,u_r,v_l,v_r,g)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real(kind = DP) :: roe_bar(4),h_l,h_r,u_l,u_r,v_l,v_r,g
    !h_bar
    roe_bar(1) = (h_l+h_r) / 2 
    !u_bar
    roe_bar(2) = (sqrt(h_l)*u_l + sqrt(h_r)*u_r) / (sqrt(h_l) + sqrt(h_r))
    ! v_bar
    roe_bar(3) = (sqrt(h_l)*v_l + sqrt(h_r)*v_r) / (sqrt(h_l) + sqrt(h_r))
    !c_bar
    roe_bar(4) = sqrt(roe_bar(1)*g)
end function

function flux1(Q02,Q01,Q0,Q1,Q2,g)
    implicit none
    interface
        Function minus_r(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: minus_r(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
            
        End Function minus_r

        Function plus_l(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: plus_l(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
        end function

        Function absF(u,v,c,u_l,u_r,c_l,c_r)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: i
            real(kind = DP) :: absF(3,3)
            real(kind = DP), intent(IN) :: u,v,c,u_l,u_r,c_l,c_r
            real(kind = DP) :: diag(3,3) = 0, R(3,3),L(3,3),delta(3)
        end function

        function func1(F,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: F(3),g,func1(3)
        end function

        function func2(G,g1)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: G(3),g1,func2(3)
        end function 

        Function roe_bar(h_l,h_r,u_l,u_r,v_l,v_r,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: roe_bar(4),h_l,h_r,u_l,u_r,v_l,v_r,g
        end function
    end interface
        integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
        real(kind = DP) :: flux1(3),Q02(3),Q01(3),Q0(3),Q1(3),Q2(3),g
        real(kind = DP) :: Q_minus_r(3),Q_plus_r(3),Q_minus_l(3),Q_plus_l(3)
        real(kind = DP) :: h_minus_r,u_minus_r,h_plus_r,u_plus_r,h_minus_l,u_minus_l,h_plus_l,u_plus_l
        real(kind = DP) :: v_minus_r,v_minus_l,v_plus_l,v_plus_r,c_minus_l,c_minus_r,c_plus_l,c_plus_r
        real(kind = DP) :: roe_bar_l(4),roe_bar_r(4)
        real(kind = DP) :: absa_l(3,3),absa_r(3,3)
        real(kind = DP) :: f_l(3),f_r(3)

        Q_minus_r = minus_r(Q01,Q0,Q1,3)
        Q_plus_r = plus_l(Q0,Q1,Q2,3)
        Q_minus_l = minus_r(Q02,Q01,Q0,3)
        Q_plus_l = plus_l(Q01,Q0,Q1,3)
        
        h_minus_r = Q_minus_r(1)
        u_minus_r = Q_minus_r(2) / Q_minus_r(1)
        v_minus_r = Q_minus_r(3) / Q_minus_r(1)
        c_minus_r = sqrt(h_minus_r*g)
        
        h_plus_r = Q_plus_r(1)
        u_plus_r = Q_plus_r(2) / Q_plus_r(1)
        v_plus_r = Q_plus_r(3) / Q_plus_r(1)
        c_plus_r = sqrt(h_plus_r*g)
        
        h_minus_l = Q_minus_l(1)
        u_minus_l = Q_minus_l(2) / Q_minus_l(1)
        v_minus_l = Q_minus_l(3) / Q_minus_l(1)
        c_minus_l = sqrt(h_minus_l*g)
        
        h_plus_l = Q_plus_l(1)
        u_plus_l = Q_plus_l(2) / Q_plus_l(1)
        v_plus_l = Q_plus_l(3) / Q_plus_l(1)
        c_plus_l = sqrt(h_plus_l*g)

        roe_bar_l = roe_bar(h_minus_l,h_plus_l,u_minus_l,u_plus_l,v_minus_l,v_plus_l,g)
        roe_bar_r = roe_bar(h_minus_r,h_plus_r,u_minus_r,u_plus_r,v_minus_r,v_plus_r,g)

        absa_l = absF(roe_bar_l(2),roe_bar_l(3),roe_bar_l(4),u_minus_l,u_plus_l,c_minus_l,c_plus_l)
        absa_r = absF(roe_bar_r(2),roe_bar_r(3),roe_bar_r(4),u_minus_r,u_plus_r,c_minus_r,c_plus_r)

        f_l = 1.0_DP/2*(func1(Q_minus_l,g)+func1(Q_plus_l,g)) - 1.0_DP/2*(matmul(absa_l,(Q_plus_l - Q_minus_l)))
        f_r = 1.0_DP/2*(func1(Q_minus_r,g)+func1(Q_plus_r,g)) - 1.0_DP/2*(matmul(absa_r,(Q_plus_r - Q_minus_r)))

        flux1 = -(f_r-f_l)
end function

function flux2(Q02,Q01,Q0,Q1,Q2,g)
    implicit none
    interface
        Function minus_r(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: minus_r(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
            
        End Function minus_r

        Function plus_l(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: plus_l(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
        end function

        Function absG(u,v,c,v_l,v_r,c_l,c_r)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: i
            real(kind = DP) :: absG(3,3)
            real(kind = DP), intent(IN) :: u,v,c,v_l,v_r,c_l,c_r
            real(kind = DP) :: diag(3,3) = 0, R(3,3),L(3,3),delta(3)
        end function

        function func1(F,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: F(3),g,func1(3)
        end function

        function func2(G,g1)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: G(3),g1,func2(3)
        end function 

        Function roe_bar(h_l,h_r,u_l,u_r,v_l,v_r,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: roe_bar(4),h_l,h_r,u_l,u_r,v_l,v_r,g
        end function
    end interface
        integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
        real(kind = DP) :: flux2(3),Q02(3),Q01(3),Q0(3),Q1(3),Q2(3),g
        real(kind = DP) :: Q_minus_r(3),Q_plus_r(3),Q_minus_l(3),Q_plus_l(3)
        real(kind = DP) :: h_minus_r,u_minus_r,h_plus_r,u_plus_r,h_minus_l,u_minus_l,h_plus_l,u_plus_l
        real(kind = DP) :: v_minus_r,v_minus_l,v_plus_l,v_plus_r,c_minus_l,c_minus_r,c_plus_l,c_plus_r
        real(kind = DP) :: roe_bar_l(4),roe_bar_r(4)
        real(kind = DP) :: absa_l(3,3),absa_r(3,3)
        real(kind = DP) :: f_l(3),f_r(3)

        Q_minus_r = minus_r(Q01,Q0,Q1,3)
        Q_plus_r = plus_l(Q0,Q1,Q2,3)
        Q_minus_l = minus_r(Q02,Q01,Q0,3)
        Q_plus_l = plus_l(Q01,Q0,Q1,3)
        
        h_minus_r = Q_minus_r(1)
        u_minus_r = Q_minus_r(2) / Q_minus_r(1)
        v_minus_r = Q_minus_r(3) / Q_minus_r(1)
        c_minus_r = sqrt(h_minus_r*g)
        
        h_plus_r = Q_plus_r(1)
        u_plus_r = Q_plus_r(2) / Q_plus_r(1)
        v_plus_r = Q_plus_r(3) / Q_plus_r(1)
        c_plus_r = sqrt(h_plus_r*g)
        
        h_minus_l = Q_minus_l(1)
        u_minus_l = Q_minus_l(2) / Q_minus_l(1)
        v_minus_l = Q_minus_l(3) / Q_minus_l(1)
        c_minus_l = sqrt(h_minus_l*g)
        
        h_plus_l = Q_plus_l(1)
        u_plus_l = Q_plus_l(2) / Q_plus_l(1)
        v_plus_l = Q_plus_l(3) / Q_plus_l(1)
        c_plus_l = sqrt(h_plus_l*g)

        roe_bar_l = roe_bar(h_minus_l,h_plus_l,u_minus_l,u_plus_l,v_minus_l,v_plus_l,g)
        roe_bar_r = roe_bar(h_minus_r,h_plus_r,u_minus_r,u_plus_r,v_minus_r,v_plus_r,g)

        absa_l = absG(roe_bar_l(2),roe_bar_l(3),roe_bar_l(4),v_minus_l,v_plus_l,c_minus_l,c_plus_l)
        absa_r = absG(roe_bar_r(2),roe_bar_r(3),roe_bar_r(4),v_minus_r,v_plus_r,c_minus_r,c_plus_r)

        f_l = 1.0_DP/2*(func2(Q_minus_l,g)+func2(Q_plus_l,g)) - 1.0_DP/2*(matmul(absa_l,(Q_plus_l - Q_minus_l)))
        f_r = 1.0_DP/2*(func2(Q_minus_r,g)+func2(Q_plus_r,g)) - 1.0_DP/2*(matmul(absa_r,(Q_plus_r - Q_minus_r)))

        flux2 = -(f_r-f_l)
end function

function Gussian(n_x,n_y)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n_x,n_y,i,j
    real(kind = DP) :: Gussian(n_y,n_x,3),x,y
    do i = 1,n_y
        do j = 1,n_x
            x = (j+.0_DP)/n_x-0.5_DP
            y = (i+.0_DP)/n_y-0.5_DP
            ! print *,x
            Gussian(i,j,1) = 0.3_DP*exp(-((x-0.2_DP)**2 + (y-0.2_DP)**2)/0.001_DP)+0.7_DP
            Gussian(i,j,2) = .0_DP
            Gussian(i,j,3) = .0_DP
        end do 
    end do 
end function

function east_flow(n_x,n_y)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    real, parameter :: PI = 3.1415926
    integer :: n_x,n_y,i,j
    real(kind = DP) :: east_flow(n_y,n_x,3),x,y,Long = 2*PI,lati = PI
    do i = 1,n_y
        do j = 1,n_x
            ! x = (j+.0_DP)/n_x-0.5_DP
            y = ((i+.0_DP)/n_y-0.5_DP)*lati
            ! print *,x
            east_flow(i,j,1) = 10000.0_DP - 400.0_DP * cos(y)
            east_flow(i,j,2) = .0_DP
            east_flow(i,j,3) = .0_DP
        end do 
    end do 
end function

function solitary(n_x,n_y)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n_x,n_y,i,j
    real(kind = DP) :: solitary(n_y,n_x,3),x,y,B = 0.395,h,u,v,cauchy,dcau
    do i = 1,n_y
        do j = 1,n_x
            x = ((j+.0_DP)/n_x-0.5_DP)*48.0_DP
            y = ((i+.0_DP)/n_y-0.5_DP)*24.0_DP
            cauchy = 0.772_DP*B**2*(1/cosh(B*x))**2
            dcau = cauchy*B*2.0_DP*(-tanh(B*x))
            h = cauchy*(3.0_DP+6.0_DP*y**2)/4.0_DP*exp(-0.5_DP*y**2)+1.0
            u = cauchy*(-9.0_DP+6.0_DP*y**2)/4.0_DP*exp(-0.5_DP*y**2)
            v = 2.0_DP*y*dcau*exp(-0.5_DP*y**2)
            ! print *,x
            solitary(i,j,1) = h
            solitary(i,j,2) = h*u
            solitary(i,j,3) = h*v
        end do 
    end do 
end function

function inerio(n_x,n_y)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n_x,n_y,i,j
    real(kind = DP) :: inerio(n_y,n_x,3),x,y,h,u,v,k = 0.5,w1
    w1 = -k/2-sqrt((k/2)**2+1.0_DP)
    do i = 1,n_y
        do j = 1,n_x
            x = ((j+.0_DP)/n_x-0.5_DP)*12.0_DP
            y = ((i+.0_DP)/n_y-0.5_DP)*12.0_DP
            h = 0.1_DP*exp(-1.0_DP/2*y**2)*2*y*cos(0.5_DP*x)+1.0_DP
            u = 0.1_DP*exp(-1.0_DP/2*y**2)*2*y*cos(0.5_DP*x)
            v = 0.1_DP*exp(-1.0_DP/2*y**2)*2*(w1+k)*sin(0.5_DP*x)
            ! print *,x
            inerio(i,j,1) = h
            inerio(i,j,2) = h*u
            inerio(i,j,3) = h*v
        end do 
    end do 
end function

function refraction(n_x,n_y,Lx,Ly)
    implicit none
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: n_x,n_y,i,j
    real(kind = DP) :: refraction(n_y,n_x,3),x,y,Lx,Ly
    do i = 1,n_y
        do j = 1,n_x
            x = ((j+.0_DP)/n_x) * Lx
            y = ((i+.0_DP)/n_y) * Ly
            refraction(i,j,1) = 40-(sqrt(3.0_DP)/2.0_DP*x-y*0.5_DP)*0.02_DP
            refraction(i,j,2) = .0_DP
            refraction(i,j,3) = .0_DP
        end do
    end do 
end function

subroutine evolve(T,n_x,n_y,x,y,tau_d,tau_r,heq,dt,Q_ini,Q_list)
    implicit none
    interface

        function step(value)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: value,step
        end function

        function logic2int(logic)
            implicit none
            logical :: logic
            integer :: logic2int
        end function

        Function minus_r(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: minus_r(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
            
        End Function minus_r

        Function plus_l(v0,v1,v2,n)
            implicit none 
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n,i,j
            real(kind = DP) :: plus_l(n)
            Real(kind = DP), intent(IN) :: v0(n),v1(n),v2(n)
            Real(kind = DP) :: v_0(n),v_1(n),v(2,n),d(2),epsn,beta(2,n),alpha(2,n),w(2,n)
        end function

        function func1(F,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: F(3),g,func1(3)
        end function

        function func2(G,g1)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: G(3),g1,func2(3)
        end function 

        Function roe_bar(h_l,h_r,u_l,u_r,v_l,v_r,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: roe_bar(4),h_l,h_r,u_l,u_r,v_l,v_r,g
        end function

        function flux1(Q02,Q01,Q0,Q1,Q2,g)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: flux1(3),Q02(3),Q01(3),Q0(3),Q1(3),Q2(3),g
            real(kind = DP) :: Q_minus_r(3),Q_plus_r(3),Q_minus_l(3),Q_plus_l(3)
            real(kind = DP) :: h_minus_r,u_minus_r,h_plus_r,u_plus_r,h_minus_l,u_minus_l,h_plus_l,u_plus_l
            real(kind = DP) :: v_minus_r,v_minus_l,v_plus_l,v_plus_r
            real(kind = DP) :: roe_bar_l(4),roe_bar_r(4)
            real(kind = DP) :: absa_l(3,3),absa_r(3,3)
            real(kind = DP) :: f_l(3),f_r(3)
        end function

        function flux2(Q02,Q01,Q0,Q1,Q2,g)
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real(kind = DP) :: flux2(3),Q02(3),Q01(3),Q0(3),Q1(3),Q2(3),g
            real(kind = DP) :: Q_minus_r(3),Q_plus_r(3),Q_minus_l(3),Q_plus_l(3)
            real(kind = DP) :: h_minus_r,u_minus_r,h_plus_r,u_plus_r,h_minus_l,u_minus_l,h_plus_l,u_plus_l
            real(kind = DP) :: v_minus_r,v_minus_l,v_plus_l,v_plus_r
            real(kind = DP) :: roe_bar_l(4),roe_bar_r(4)
            real(kind = DP) :: absa_l(3,3),absa_r(3,3)
            real(kind = DP) :: f_l(3),f_r(3)
        end function

        function Gussian(n_x,n_y)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n_x,n_y,i,j
            real(kind = DP) :: Gussian(n_y,n_x,3),x,y
        end function

        function east_flow(n_x,n_y)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            real, parameter :: PI = 3.1415926
            integer :: n_x,n_y,i,j
            real(kind = DP) :: east_flow(n_y,n_x,3),x,y,Long = 2*PI,lati = PI
        end function

        function solitary(n_x,n_y)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n_x,n_y,i,j
            real(kind = DP) :: solitary(n_y,n_x,3),x,y,B = 0.395,h,u,v,cauchy,dcau
        end function

        function inerio(n_x,n_y)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n_x,n_y,i,j
            real(kind = DP) :: inerio(n_y,n_x,3),x,y,h,u,v,k = 0.5,w1
        end function

        function refraction(n_x,n_y,Lx,Ly)
            implicit none
            integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
            integer :: n_x,n_y,i,j
            real(kind = DP) :: refraction(n_y,n_x,3),x,y,Lx,Ly
        end function
    end interface
    ! INPUT parameters: (T/steps,n_x,n_y)
    integer(kind = 4), intent(in) :: T,n_x,n_y
    real(kind = 8), intent(in) :: heq(n_y,n_x),Q_ini(n_y,n_x,3),tau_d,tau_r,x,y,dt

    ! OUTPUT varibles: 
    real(kind = 8), intent(out) :: Q_list(n_y,n_x,3,T)


    ! inner parameters:
    integer, parameter :: DP = selected_real_kind(r = 50,p = 14)
    integer :: i,j,m
    real(kind = DP), parameter :: beta = 1.0,PI = 3.1415926
    real(kind = DP) :: t_total(T),g = 1.0,cfl = 0.5,dx,dy,prelambda(n_y,n_x,2),lambda1,lambda2,qsource
    real(kind = DP) :: Q(n_y,n_x,3),rotation(n_y,n_x,3),source(n_y,n_x,3),drag(n_y,n_x,3)
    real(kind = DP) :: Q_1(n_y,n_x,3),Q_2(n_y,n_x,3),Q_new(n_y,n_x,3)
    dx = x / n_x
    dy = y / n_y
    ! dt = 0.15_DP
    lambda1 = dt/dx
    lambda2 = dt/dy
    Q = Q_ini
    Q_1(:,:,:) = Q(:,:,:)
    Q_2(:,:,:) = Q(:,:,:)
    Q_new(:,:,:) = Q(:,:,:)

    
    do m = 1,T
        Q_list(:,:,:,m) = Q
        print *,m
    !     do i = 1,n_y
    !         do j = 1,n_x
    !             prelambda(i,j,1) = cfl / (abs(Q(i,j,2)/Q(i,j,1))+sqrt(g*Q(i,j,1)))
    !             prelambda(i,j,2) = cfl / (abs(Q(i,j,3)/Q(i,j,1))+sqrt(g*Q(i,j,1)))
    !         end do 
    !     end do
    !     lambda = minval(prelambda)
        ! dt = lambda * 0.04
        ! t_total(m) = dt
            

        !coriolis & source term & drag term
        do i = 1,n_y
            do j = 1,n_x
                ! coriolis
                rotation(i,j,1) = .0_DP
                rotation(i,j,2) = beta*Q(i,j,3)*((i+.0_DP)/n_y-0.5_DP)*y
                rotation(i,j,3) = -beta*Q(i,j,2)*((i+.0_DP)/n_y-0.5_DP)*y
                ! source term 
                qsource = (heq(i,j)-Q(i,j,1))/tau_r
                source(i,j,1) = qsource
                source(i,j,2) = qsource*Q(i,j,2)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                source(i,j,3) = qsource*Q(i,j,3)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                ! drag term
                drag(i,j,1) = .0_DP
                drag(i,j,2) = -Q(i,j,2)/tau_d
                drag(i,j,3) = -Q(i,j,3)/tau_d
            end do
        end do

        ! 1st iteration
        do i = 1,n_y
            do j = 1,n_x
                Q_1(i,j,:) = Q(i,j,:) + lambda1 * &
                flux1(Q(i,j-2 + n_x * logic2int(j-2<=0),:),Q(i,j-1+ n_x * logic2int(j-1<=0),:),Q(i,j,:),&
                Q(i,j+1-n_x*logic2int(j+1>n_x),:),Q(i,j+2-n_x*logic2int(j+2>n_x),:),g)&
                + lambda2 * flux2(Q(i-2+n_y*logic2int(i-2<=0),j,:),Q(i-1+n_y*logic2int(i-1<=0),j,:),Q(i,j,:),&
                Q(i+1-n_y*logic2int(i+1>n_y),j,:),Q(i+2-n_y*logic2int(i+2>n_y),j,:),g)+ &
                dt * (rotation(i,j,:)+ source(i,j,:)+drag(i,j,:))
            end do
        end do
        


        !coriolis & source term & drag term
        do i = 1,n_y
            do j = 1,n_x
                ! coriolis
                rotation(i,j,1) = .0_DP
                rotation(i,j,2) = beta*Q(i,j,3)*((i+.0_DP)/n_y-0.5_DP)*y
                rotation(i,j,3) = -beta*Q(i,j,2)*((i+.0_DP)/n_y-0.5_DP)*y
                ! source term 
                qsource = (heq(i,j)-Q(i,j,1))/tau_r
                source(i,j,1) = qsource
                source(i,j,2) = qsource*Q(i,j,2)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                source(i,j,3) = qsource*Q(i,j,3)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                ! drag term
                drag(i,j,1) = .0_DP
                drag(i,j,2) = -Q(i,j,2)/tau_d
                drag(i,j,3) = -Q(i,j,3)/tau_d
            end do
        end do

        ! 2nd iteration
        do i = 1,n_y
            do j = 1,n_x
                Q_2(i,j,:) =0.75_DP*Q(i,j,:)+0.25_DP*Q_1(i,j,:)+0.25_DP*lambda1*&
                flux1(Q_1(i,j-2+n_x*logic2int(j-2<=0),:),Q_1(i,j-1+n_x*logic2int(j-1<=0),:),Q_1(i,j,:),&
                Q_1(i,j+1-n_x*logic2int(j+1>n_x),:),Q_1(i,j+2-n_x*logic2int(j+2>n_x),:),g)+0.25_DP*lambda2*&
                flux2(Q_1(i-2+n_y*logic2int(i-2<=0),j,:),Q_1(i-1+n_y*logic2int(i-1<=0),j,:),Q_1(i,j,:),&
                Q_1(i+1-n_y*logic2int(i+1>n_y),j,:),Q_1(i+2-n_y*logic2int(i+2>n_y),j,:),g)&
                +0.25_DP*dt * (rotation(i,j,:)+source(i,j,:)+drag(i,j,:))
            end do
        end do 

        !coriolis & source term & drag term
        do i = 1,n_y
            do j = 1,n_x
                ! coriolis
                rotation(i,j,1) = .0_DP
                rotation(i,j,2) = beta*Q(i,j,3)*((i+.0_DP)/n_y-0.5_DP)*y
                rotation(i,j,3) = -beta*Q(i,j,2)*((i+.0_DP)/n_y-0.5_DP)*y
                ! source term 
                qsource = (heq(i,j)-Q(i,j,1))/tau_r
                source(i,j,1) = qsource
                source(i,j,2) = qsource*Q(i,j,2)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                source(i,j,3) = qsource*Q(i,j,3)/Q(i,j,1)*(1-step(heq(i,j)-Q(i,j,1)))
                ! drag term
                drag(i,j,1) = .0_DP
                drag(i,j,2) = -Q(i,j,2)/tau_d
                drag(i,j,3) = -Q(i,j,3)/tau_d
            end do
        end do

        ! 3rd iteration
        do i = 1,n_y
            do j = 1,n_x
                Q_new(i,j,:) = 1.0_DP/3 * Q(i,j,:) + 2.0_DP/3*Q_2(i,j,:)+2.0_DP/3*lambda1*&
                flux1(Q_2(i,j-2+n_x*logic2int(j-2<=0),:),Q_2(i,j-1+n_x*logic2int(j-1<=0),:),Q_2(i,j,:),&
                Q_2(i,j+1-n_x*logic2int(j+1>n_x),:),Q_2(i,j+2-n_x*logic2int(j+2>n_x),:),g)+2.0_DP/3*lambda2*&
                flux2(Q_2(i-2+n_y*logic2int(i-2<=0),j,:),Q_2(i-1+n_y*logic2int(i-1<=0),j,:),Q_2(i,j,:),&
                Q_2(i+1-n_y*logic2int(i+1>n_y),j,:),Q_2(i+2-n_y*logic2int(i+2>n_y),j,:),g)&
                +2.0_DP/3*dt*(rotation(i,j,:)+source(i,j,:)+drag(i,j,:))
            end do 
        end do 
        
    ! !reflection
    ! Q_new(3,:,3) = -Q_new(3,:,3)
    ! Q_new(:,3,2) = -Q_new(:,3,2)
    ! Q_new(n_y-2,:,3) = -Q_new(n_y-2,:,3)
    ! Q_new(:,n_x-2,2) = -Q_new(:,n_x-2,2)
    ! Q_new(2,:,:) = Q_new(3,:,:)
    ! Q_new(1,:,:) = Q_new(2,:,:)
    ! Q_new(:,2,:) = Q_new(:,3,:)
    ! Q_new(:,1,:) = Q_new(:,2,:)
    ! Q_new(n_y-1,:,:) = Q_new(n_y-2,:,:)
    ! Q_new(n_y,:,:) = Q_new(n_y-1,:,:)
    ! Q_new(:,n_x-1,:) = Q_new(:,n_x-2,:)
    ! Q_new(:,n_x,:) = Q_new(:,n_x-1,:)

        


        Q(:,:,:) = Q_new(:,:,:)
    end do 

    ! open(15,file = 'F:\\exoplanet\\paper\\data\\20200321_refraction4.bin' &
    ! ,Access = 'stream' , Form = 'Unformatted')
    ! write(15) Q_list(:,:,1,:)
    ! Close( 15 )
    
    ! open(16,file = 'F:\\exoplanet\\paper\\data\\20200317_inerio_hu1.bin' &
    ! ,Access = 'stream' , Form = 'Unformatted')
    ! write(16) Q_list(:,:,2,:)
    ! Close( 16 )

    ! open(17,file = 'F:\\exoplanet\\paper\\data\\20200317_inerio_hv1.bin' &
    ! ,Access = 'stream' , Form = 'Unformatted')
    ! write(17) Q_list(:,:,3,:)
    ! Close( 17 )

    ! print*, Q(:,:,1)

end subroutine