subroutine boundary
    use commondata
    implicit none
    integer :: jfs1
    real*8 :: rho, u, v, p
    
    if(FLG_prob==1) then    ! Rayleigh Taylor
        !  up stream
        rho = 1.d0;    p = 2.5d0
        u = 0.d0;    v = 0.d0
        v1(1,:,jn+1:jn+4) = rho
        v1(2,:,jn+1:jn+4) = rho*u
        v1(3,:,jn+1:jn+4) = rho*v
        v1(4,:,jn+1:jn+4) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
        
        !  down stream
        rho = 2.d0;    p = 1.d0
        u = 0.d0;    v = 0.d0        
        v1(1,:,-4:-1) = rho
        v1(2,:,-4:-1) = rho*u
        v1(3,:,-4:-1) = rho*v
        v1(4,:,-4:-1) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
        
        !  left wall
        v1(:,-4,:) = v1(:,4,:);    v1(2,-4,:) = -v1(2,4,:)
        v1(:,-3,:) = v1(:,3,:);    v1(2,-3,:) = -v1(2,3,:)
        v1(:,-2,:) = v1(:,2,:);    v1(2,-2,:) = -v1(2,2,:)
        v1(:,-1,:) = v1(:,1,:);    v1(2,-1,:) = -v1(2,1,:)
        
        !  right wall
        v1(:,im+4,:) = v1(:,im-4,:);    v1(2,im+4,:) = -v1(2,im-4,:)
        v1(:,im+3,:) = v1(:,im-3,:);    v1(2,im+3,:) = -v1(2,im-3,:)
        v1(:,im+2,:) = v1(:,im-2,:);    v1(2,im+2,:) = -v1(2,im-2,:)
        v1(:,im+1,:) = v1(:,im-1,:);    v1(2,im+1,:) = -v1(2,im-1,:)
        
        
    else if(FLG_prob==2) then    ! Double Mach Reflection
        !  left inlet boundary
        rho = 8.d0
        u = 8.25d0*sqrt(3.d0)*0.5d0
        v = -8.25d0*0.5d0
        p = 116.5d0
        v1(1,-4:-1,:) = rho
        v1(2,-4:-1,:) = rho*u
        v1(3,-4:-1,:) = rho*v
        v1(4,-4:-1,:) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
        
        !  right outlet boundary
        v1(:,im+1,:) = v1(:,im,:)
        v1(:,im+2,:) = v1(:,im,:)
        v1(:,im+3,:) = v1(:,im,:)
        v1(:,im+4,:) = v1(:,im,:)
        
        !  bottom reflective boundary
        do i = -4, im+4
            if( xr(i,0) <= (1.d0/6.d0) ) then
                rho = 8.d0
                u = 8.25d0*sqrt(3.d0)*0.5d0
                v = -8.25d0*0.5d0
                p = 116.5d0
                v1(1,i,-4:-1) = rho
                v1(2,i,-4:-1) = rho*u
                v1(3,i,-4:-1) = rho*v
                v1(4,i,-4:-1) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
            else
                v1(:,i,-4) = v1(:,i,4);    v1(3,i,-4) = -v1(3,i,4)
                v1(:,i,-3) = v1(:,i,3);    v1(3,i,-3) = -v1(3,i,3)
                v1(:,i,-2) = v1(:,i,2);    v1(3,i,-2) = -v1(3,i,2)
                v1(:,i,-1) = v1(:,i,1);    v1(3,i,-1) = -v1(3,i,1)
            end if
        end do
        
        !  top boundary
        do i = -4, im+4
            if( xr(i,jn) <= (1.d0/6.d0+(1.d0+20.d0*tn)/sqrt(3.d0)) ) then
                rho = 8.d0
                u = 8.25d0*sqrt(3.d0)*0.5d0
                v = -8.25d0*0.5d0
                p = 116.5d0
            else
                rho = 1.4d0
                u = 0.d0
                v = 0.d0
                p = 1.d0
            end if
            v1(1,i,jn+1:jn+4) = rho
            v1(2,i,jn+1:jn+4) = rho*u
            v1(3,i,jn+1:jn+4) = rho*v
            v1(4,i,jn+1:jn+4) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
        end do
        
        
    else if(FLG_prob==3) then    ! Forward Step
        !  left inlet boundary
        rho = 1.4d0
        u = 3.d0
        v = 0.d0
        p = 1.d0
        v1(1,-4:-1,:) = rho
        v1(2,-4:-1,:) = rho*u
        v1(3,-4:-1,:) = rho*v
        v1(4,-4:-1,:) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
        
        !  right outlet boundary
        v1(:,im+1,:) = v1(:,im,:)
        v1(:,im+2,:) = v1(:,im,:)
        v1(:,im+3,:) = v1(:,im,:)
        v1(:,im+4,:) = v1(:,im,:)
        
        !  top reflective boundary
        v1(:,:,jn+4) = v1(:,:,jn-3);    v1(3,:,jn+4) = -v1(3,:,jn-3)
        v1(:,:,jn+3) = v1(:,:,jn-2);    v1(3,:,jn+3) = -v1(3,:,jn-2)
        v1(:,:,jn+2) = v1(:,:,jn-1);    v1(3,:,jn+2) = -v1(3,:,jn-1)
        v1(:,:,jn+1) = v1(:,:,jn  );    v1(3,:,jn+1) = -v1(3,:,jn  )
        
        !  bottom reflective boundary
        do i = -4, im+4
            if( i<ifs ) then
                jfs1 = 0
            else
                jfs1 = jfs+1
            end if
            v1(:,i,jfs1-4) = v1(:,i,jfs1+3);    v1(3,i,jfs1-4) = -v1(3,i,jfs1+3)
            v1(:,i,jfs1-3) = v1(:,i,jfs1+2);    v1(3,i,jfs1-3) = -v1(3,i,jfs1+2)
            v1(:,i,jfs1-2) = v1(:,i,jfs1+1);    v1(3,i,jfs1-2) = -v1(3,i,jfs1+1)
            v1(:,i,jfs1-1) = v1(:,i,jfs1  );    v1(3,i,jfs1-1) = -v1(3,i,jfs1  )
        end do
        
        !  forward reflective boundary
        v1(:,ifs+3,-4:jfs) = v1(:,ifs-4,-4:jfs);    v1(2,ifs+3,-4:jfs) = -v1(2,ifs-4,-4:jfs)
        v1(:,ifs+2,-4:jfs) = v1(:,ifs-3,-4:jfs);    v1(2,ifs+2,-4:jfs) = -v1(2,ifs-3,-4:jfs)
        v1(:,ifs+1,-4:jfs) = v1(:,ifs-2,-4:jfs);    v1(2,ifs+1,-4:jfs) = -v1(2,ifs-2,-4:jfs)
        v1(:,ifs  ,-4:jfs) = v1(:,ifs-1,-4:jfs);    v1(2,ifs  ,-4:jfs) = -v1(2,ifs-1,-4:jfs)
        
        
    else if(FLG_prob==4) then    ! Isentropic Vortex
        !  periodic boundary condition
        do i = 1, 4
            v1(:,-i,:) = v1(:,im+1-i,:)
            v1(:,:,-i) = v1(:,:,jn+1-i)
            v1(:,im+i,:) = v1(:,i-1,:)
            v1(:,:,jn+i) = v1(:,:,i-1)
        end do
        do i = 1, 4
            do j = 1, 4
                v1(:,-i,-j) = v1(:,im+1-i,jn+1-j)
                v1(:,im+i,-j) = v1(:,i-1,jn+1-j)
                v1(:,-i,jn+j) = v1(:,im+1-i,j-1)
                v1(:,im+i,jn+j) = v1(:,i-1,j-1)
            end do
        end do
        !v1(:,-4,:) = v1(:,im-3,:);    v1(:,-3,:) = v1(:,im-2,:)
        !v1(:,-2,:) = v1(:,im-1,:);    v1(:,-1,:) = v1(:,im,:)
        !v1(:,:,-4) = v1(:,:,jn-3);    v1(:,:,-3) = v1(:,:,jn-2)
        !v1(:,:,-2) = v1(:,:,jn-1);    v1(:,:,-1) = v1(:,:,jn)
        !v1(:,im+4,:) = v1(:,-3,:);    v1(:,im+3,:) = v1(:,-2,:)
        !v1(:,im+2,:) = v1(:,-1,:);    v1(:,im+1,:) = v1(:,0,:)
        !v1(:,:,jn+4) = v1(:,:,-3);    v1(:,:,jn+3) = v1(:,:,-2)
        !v1(:,:,jn+2) = v1(:,:,-1);    v1(:,:,jn+1) = v1(:,:,0)
        
        
    else    !  2D Riemann -- free flow boundary
        !  left free flow
        v1(:,-4,0:jn) = v1(:,4,0:jn)!;    v1(2,-4,0:jn) = -v1(2,4,0:jn)
        v1(:,-3,0:jn) = v1(:,3,0:jn)!;    v1(2,-3,0:jn) = -v1(2,3,0:jn)
        v1(:,-2,0:jn) = v1(:,2,0:jn)!;    v1(2,-2,0:jn) = -v1(2,2,0:jn)
        v1(:,-1,0:jn) = v1(:,1,0:jn)!;    v1(2,-1,0:jn) = -v1(2,1,0:jn)
        !  right free flow
        v1(:,im+4,0:jn) = v1(:,im-4,0:jn)!;    v1(2,im+4,0:jn) = -v1(2,im-4,0:jn)
        v1(:,im+3,0:jn) = v1(:,im-3,0:jn)!;    v1(2,im+3,0:jn) = -v1(2,im-3,0:jn)
        v1(:,im+2,0:jn) = v1(:,im-2,0:jn)!;    v1(2,im+2,0:jn) = -v1(2,im-2,0:jn)
        v1(:,im+1,0:jn) = v1(:,im-1,0:jn)!;    v1(2,im+1,0:jn) = -v1(2,im-1,0:jn)
        !  bottom free flow
        v1(:,:,-4) = v1(:,:,4)!;    v1(3,:,-4) = -v1(3,:,4)
        v1(:,:,-3) = v1(:,:,3)!;    v1(3,:,-3) = -v1(3,:,3)
        v1(:,:,-2) = v1(:,:,2)!;    v1(3,:,-2) = -v1(3,:,2)
        v1(:,:,-1) = v1(:,:,1)!;    v1(3,:,-1) = -v1(3,:,1)
        !  left free flow
        v1(:,:,jn+4) = v1(:,:,jn-4)!;    v1(3,:,jn+4) = -v1(3,:,jn-4)
        v1(:,:,jn+3) = v1(:,:,jn-3)!;    v1(3,:,jn+3) = -v1(3,:,jn-3)
        v1(:,:,jn+2) = v1(:,:,jn-2)!;    v1(3,:,jn+2) = -v1(3,:,jn-2)
        v1(:,:,jn+1) = v1(:,:,jn-1)!;    v1(3,:,jn+1) = -v1(3,:,jn-1)
        
        
    end if
    
end subroutine boundary
    
    
subroutine corner_x
    use commondata
    implicit none
    do j = jfs-3, jfs
        var(:,ifs  ,j) = var(:,ifs-1,j);    var(2,ifs  ,j) = -var(2,ifs-1,j)
        var(:,ifs+1,j) = var(:,ifs-2,j);    var(2,ifs+1,j) = -var(2,ifs-2,j)
        var(:,ifs+2,j) = var(:,ifs-3,j);    var(2,ifs+2,j) = -var(2,ifs-3,j)
        var(:,ifs+3,j) = var(:,ifs-4,j);    var(2,ifs+3,j) = -var(2,ifs-4,j)
        v1(:,ifs  ,j) = v1(:,ifs-1,j);    v1(2,ifs  ,j) = -v1(2,ifs-1,j)
        v1(:,ifs+1,j) = v1(:,ifs-2,j);    v1(2,ifs+1,j) = -v1(2,ifs-2,j)
        v1(:,ifs+2,j) = v1(:,ifs-3,j);    v1(2,ifs+2,j) = -v1(2,ifs-3,j)
        v1(:,ifs+3,j) = v1(:,ifs-4,j);    v1(2,ifs+3,j) = -v1(2,ifs-4,j)
    end do
end subroutine corner_x
    
    
subroutine corner_y
    use commondata
    implicit none
    do i = ifs, ifs+3
        var(:,i,jfs  ) = var(:,i,jfs+1);    var(3,i,jfs  ) = -var(3,i,jfs+1)
        var(:,i,jfs-1) = var(:,i,jfs+2);    var(3,i,jfs-1) = -var(3,i,jfs+2)
        var(:,i,jfs-2) = var(:,i,jfs+3);    var(3,i,jfs-2) = -var(3,i,jfs+3)
        var(:,i,jfs-3) = var(:,i,jfs+4);    var(3,i,jfs-3) = -var(3,i,jfs+4)
        v1(:,i,jfs  ) = v1(:,i,jfs+1);    v1(3,i,jfs  ) = -v1(3,i,jfs+1)
        v1(:,i,jfs-1) = v1(:,i,jfs+2);    v1(3,i,jfs-1) = -v1(3,i,jfs+2)
        v1(:,i,jfs-2) = v1(:,i,jfs+3);    v1(3,i,jfs-2) = -v1(3,i,jfs+3)
        v1(:,i,jfs-3) = v1(:,i,jfs+4);    v1(3,i,jfs-3) = -v1(3,i,jfs+4)
    end do
end subroutine corner_y