s_compute_cylindrical_geometry_source_flux()
!$acc routine seq
! This subroutine computes the cylindrical geometry source fluxes
#:if (NORM_DIR == 2)
    if (cyl_coord) then
        if (model_eqns == 3) then
            !Substituting the advective flux into the inviscid geometrical source flux
            !$acc loop seq
            do i = 1, E_idx
                flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
            end do
            !$acc loop seq
            do i = intxb, intxe
                flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
            end do
            ! Recalculating the radial momentum geometric source flux
            flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) = &
                flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) - p_Star
        else
            ! Substituting the advective flux into the inviscid geometrical source flux
            !$acc loop seq
            do i = 1, E_idx
                flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
            end do
            ! Recalculating the radial momentum geometric source flux
            flux_gsrc_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = &
                xi_M*(rho_L*(vel_L(dir_idx(1))* &
                                vel_L(dir_idx(1)) + &
                                s_M*(xi_L*(dir_flg(dir_idx(1))*s_S + &
                                        (1._wp - dir_flg(dir_idx(1)))* &
                                        vel_L(dir_idx(1))) - vel_L(dir_idx(1))))) &
                + xi_P*(rho_R*(vel_R(dir_idx(1))* &
                                vel_R(dir_idx(1)) + &
                                s_P*(xi_R*(dir_flg(dir_idx(1))*s_S + &
                                            (1._wp - dir_flg(dir_idx(1)))* &
                                            vel_R(dir_idx(1))) - vel_R(dir_idx(1)))))
        end if
        ! Geometrical source of the void fraction(s) is zero
        !$acc loop seq
        do i = advxb, advxe
            flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
        end do
    end if
#:endif
#:if (NORM_DIR == 3)
    if (grid_geometry == 3) then
        !$acc loop seq
        do i = 1, sys_size
            flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
        end do
        if (model_eqns == 3) then
            flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) = &
                flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) - p_Star
        else
            flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb + 1) = &
                -xi_M*(rho_L*(vel_L(dir_idx(1))* &
                                vel_L(dir_idx(1)) + &
                                s_M*(xi_L*(dir_flg(dir_idx(1))*s_S + &
                                            (1._wp - dir_flg(dir_idx(1)))* &
                                            vel_L(dir_idx(1))) - vel_L(dir_idx(1))))) &
                - xi_P*(rho_R*(vel_R(dir_idx(1))* &
                                vel_R(dir_idx(1)) + &
                                s_P*(xi_R*(dir_flg(dir_idx(1))*s_S + &
                                            (1._wp - dir_flg(dir_idx(1)))* &
                                            vel_R(dir_idx(1))) - vel_R(dir_idx(1)))))
        end if
        flux_gsrc_rs${XYZ}$_vf(j, k, l, momxe) = flux_rs${XYZ}$_vf(j, k, l, momxb + 1)
    end if
#:endif