diff --git a/cherab/solps/formats/balance.py b/cherab/solps/formats/balance.py
index 51cdd79..3a3d2da 100644
--- a/cherab/solps/formats/balance.py
+++ b/cherab/solps/formats/balance.py
@@ -41,16 +41,7 @@ def load_solps_from_balance(balance_filename):
     fhandle = netcdf.netcdf_file(balance_filename, 'r')
 
     # Load SOLPS mesh geometry
-    cr_x = fhandle.variables['crx'].data.copy()
-    cr_z = fhandle.variables['cry'].data.copy()
-    vol = fhandle.variables['vol'].data.copy()
-
-    # Re-arrange the array dimensions in the way CHERAB expects...
-    cr_x = np.moveaxis(cr_x, 0, -1)
-    cr_z = np.moveaxis(cr_z, 0, -1)
-
-    # Create the SOLPS mesh
-    mesh = SOLPSMesh(cr_x, cr_z, vol)
+    mesh = load_mesh_from_netcdf(fhandle)
 
     # Load each plasma species in simulation
 
@@ -58,7 +49,7 @@ def load_solps_from_balance(balance_filename):
     n_species = len(fhandle.variables['am'].data)
     for i in range(n_species):
 
-        # Extract the nuclear charge    
+        # Extract the nuclear charge
         if fhandle.variables['species'].data[i, 1] == b'D':
             zn = 1
         if fhandle.variables['species'].data[i, 1] == b'C':
@@ -75,7 +66,7 @@ def load_solps_from_balance(balance_filename):
         isotope = lookup_isotope(zn, number=am)
         species = prefer_element(isotope)  # Prefer Element over Isotope if the mass number is the same
 
-        # If we only need to populate species_list, there is probably a faster way to do this...        
+        # If we only need to populate species_list, there is probably a faster way to do this...
         species_list.append((species, charge))
 
     sim = SOLPSSimulation(mesh, species_list)
@@ -114,24 +105,24 @@ def load_solps_from_balance(balance_filename):
     # Calculate the total radiated power
     if eirene_run:
         # Total radiated power from B2, not including neutrals
-        b2_ploss = np.sum(fhandle.variables['b2stel_she_bal'].data, axis=0) / vol
+        b2_ploss = np.sum(fhandle.variables['b2stel_she_bal'].data, axis=0) / mesh.vol
 
         # Electron energy loss due to interactions with neutrals
         if 'eirene_mc_eael_she_bal' in fhandle.variables.keys():
-            eirene_ecoolrate = np.sum(fhandle.variables['eirene_mc_eael_she_bal'].data, axis=0) / vol
+            eirene_ecoolrate = np.sum(fhandle.variables['eirene_mc_eael_she_bal'].data, axis=0) / mesh.vol
 
         # Ionisation rate from EIRENE, needed to calculate the energy loss to overcome the ionisation potential of atoms
         if 'eirene_mc_papl_sna_bal' in fhandle.variables.keys():
-            eirene_potential_loss = rydberg_energy * np.sum(fhandle.variables['eirene_mc_papl_sna_bal'].data, axis=(0))[1, :, :] * el_charge / vol
+            eirene_potential_loss = rydberg_energy * np.sum(fhandle.variables['eirene_mc_papl_sna_bal'].data, axis=(0))[1, :, :] * el_charge / mesh.vol
 
         # This will be negative (energy sink); multiply by -1
         sim.total_radiation = -1.0 * (b2_ploss + (eirene_ecoolrate - eirene_potential_loss))
 
     else:
         # Total radiated power from B2, not including neutrals
-        b2_ploss = np.sum(fhandle.variables['b2stel_she_bal'].data, axis=0) / vol
+        b2_ploss = np.sum(fhandle.variables['b2stel_she_bal'].data, axis=0) / mesh.vol
 
-        potential_loss = np.sum(fhandle.variables['b2stel_sna_ion_bal'].data, axis=0) / vol
+        potential_loss = np.sum(fhandle.variables['b2stel_sna_ion_bal'].data, axis=0) / mesh.vol
 
         # Save total radiated power to the simulation object
         sim.total_radiation = rydberg_energy * el_charge * potential_loss - b2_ploss
@@ -139,3 +130,41 @@ def load_solps_from_balance(balance_filename):
     fhandle.close()
 
     return sim
+
+
+def load_mesh_from_netcdf(fhandle):
+
+    # Load SOLPS mesh geometry
+    r = fhandle.variables['crx'].data.copy()
+    z = fhandle.variables['cry'].data.copy()
+    vol = fhandle.variables['vol'].data.copy()
+
+    # Re-arrange the array dimensions in the way CHERAB expects...
+    r = np.moveaxis(r, 0, -1)
+    z = np.moveaxis(z, 0, -1)
+
+    # Loading neighbouring cell indices
+    neighbix = np.zeros(r.shape, dtype=np.int)
+    neighbiy = np.zeros(r.shape, dtype=np.int)
+
+    neighbix[:, :, 0] = fhandle.variables['leftix'].data.copy().astype(np.int)  # poloidal prev.
+    neighbix[:, :, 1] = fhandle.variables['bottomix'].data.copy().astype(np.int)  # radial prev.
+    neighbix[:, :, 2] = fhandle.variables['rightix'].data.copy().astype(np.int)  # poloidal next
+    neighbix[:, :, 3] = fhandle.variables['topix'].data.copy().astype(np.int)  # radial next
+
+    neighbiy[:, :, 0] = fhandle.variables['leftiy'].data.copy().astype(np.int)
+    neighbiy[:, :, 1] = fhandle.variables['bottomiy'].data.copy().astype(np.int)
+    neighbiy[:, :, 2] = fhandle.variables['rightiy'].data.copy().astype(np.int)
+    neighbiy[:, :, 3] = fhandle.variables['topiy'].data.copy().astype(np.int)
+
+    # In SOLPS cell indexing starts with -1 (guarding cell), but in SOLPSMesh -1 means no neighbour.
+    if neighbix.min() < -1 or neighbiy.min() < -1:
+        neighbix += 1
+        neighbiy += 1
+    neighbix[neighbix == r.shape[0]] = -1
+    neighbiy[neighbiy == r.shape[1]] = -1
+
+    # Create the SOLPS mesh
+    mesh = SOLPSMesh(r, z, vol, neighbix, neighbiy)
+
+    return mesh
diff --git a/cherab/solps/formats/mdsplus.py b/cherab/solps/formats/mdsplus.py
index d802ed5..ca0d0b8 100644
--- a/cherab/solps/formats/mdsplus.py
+++ b/cherab/solps/formats/mdsplus.py
@@ -22,7 +22,7 @@
 
 from cherab.core.atomic.elements import lookup_isotope
 from cherab.solps.mesh_geometry import SOLPSMesh
-from cherab.solps.solps_plasma import SOLPSSimulation, prefer_element
+from cherab.solps.solps_plasma import SOLPSSimulation, prefer_element, b2_flux_to_velocity, eirene_flux_to_velocity
 
 from matplotlib import pyplot as plt
 
@@ -44,19 +44,22 @@ def load_solps_from_mdsplus(mds_server, ref_number):
     conn.openTree('solps', ref_number)
 
     # Load SOLPS mesh geometry and lookup arrays
-    mesh = load_mesh_from_mdsplus(conn)
+    mesh = load_mesh_from_mdsplus(conn, mdsExceptions)
 
-    # Load each plasma species
+    # Load each plasma species in simulation
     ns = conn.get('\SOLPS::TOP.IDENT.NS').data()  # Number of species
     zn = conn.get('\SOLPS::TOP.SNAPSHOT.GRID.ZN').data().astype(np.int)  # Nuclear charge
     am = np.round(conn.get('\SOLPS::TOP.SNAPSHOT.GRID.AM').data()).astype(np.int)  # Atomic mass number
     charge = conn.get('\SOLPS::TOP.SNAPSHOT.GRID.ZA').data().astype(np.int)   # Ionisation/charge
 
     species_list = []
+    neutral_indx = []
     for i in range(ns):
         isotope = lookup_isotope(zn[i], number=am[i])
         species = prefer_element(isotope)  # Prefer Element over Isotope if the mass number is the same
         species_list.append((species, charge[i]))
+        if charge[i] == 0:
+            neutral_indx.append(i)
 
     sim = SOLPSSimulation(mesh, species_list)
     ni = mesh.nx
@@ -75,55 +78,59 @@ def load_solps_from_mdsplus(mds_server, ref_number):
     sim.ion_temperature = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.TI').data(), 0, 1)
 
     # Load species density
-    sim.species_density = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.NA').data(), 0, 2)
+    species_density = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.NA').data(), 0, 2)
 
-    # Load the neutral atom density from Eirene if available
-    try:
-        dab2 = conn.get('\SOLPS::TOP.SNAPSHOT.DAB2').data()
-        if isinstance(dab2, np.ndarray):
-            # Replace the species densities
-            neutral_densities = np.swapaxes(dab2, 0, 2)
-
-            neutral_i = 0  # counter for neutral atoms
-            for k, sp in enumerate(sim.species_list):
-                charge = sp[1]
-                if charge == 0:
-                    sim.species_density[:, :, k] = neutral_densities[:, :, neutral_i]
-                    neutral_i += 1
-    except mdsExceptions.TreeNNF:
-        pass
-
-    # Load the neutral atom temperature from Eirene if available
-    try:
-        tab2 = conn.get('\SOLPS::TOP.SNAPSHOT.TAB2').data()
-        if isinstance(tab2, np.ndarray):
-            sim.neutral_temperature = np.swapaxes(tab2, 0, 2)
-    except mdsExceptions.TreeNNF:
-        pass
+    # Load parallel velocity
+    parallel_velocity = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.UA').data(), 0, 2)
 
-    # TODO: Eirene data (TOP.SNAPSHOT.PFLA, TOP.SNAPSHOT.RFLA) should be used for neutral atoms.
-    velocities = np.zeros((ni, nj, len(sim.species_list), 3))
-    velocities[:, :, :, 0] = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.UA').data(), 0, 2)
+    # Load poloidal and radial particle fluxes for velocity calculation
+    poloidal_flux = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.FNAX').data(), 0, 2)
+    radial_flux = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.FNAY').data(), 0, 2)
 
-    ################################################
-    # Calculate the species' velocity distribution #
+    # B2 fluxes are defined between cells, so correcting array shapes if needed
+    if poloidal_flux.shape[0] == ni - 1:
+        poloidal_flux = np.vstack((np.zeros((1, nj, ns)), poloidal_flux))
 
-    # calculate field component ratios for velocity conversion
-    bplane2 = sim.b_field[:, :, 0]**2 + sim.b_field[:, :, 2]**2
-    parallel_to_toroidal_ratio = sim.b_field[:, :, 0] * sim.b_field[:, :, 2] / bplane2
+    if radial_flux.shape[1] == nj - 1:
+        radial_flux = np.hstack((np.zeros((ni, 1, ns)), radial_flux))
 
-    # Calculate toroidal velocity components
-    velocities[:, :, :, 2] = velocities[:, :, :, 0] * parallel_to_toroidal_ratio[:, :, None]
+    # Obtaining velocity from B2 flux
+    velocities_cartesian = b2_flux_to_velocity(mesh, species_density, poloidal_flux, radial_flux, parallel_velocity, sim.b_field_cartesian)
 
-    # Radial velocity is obtained from radial particle flux
-    radial_particle_flux = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.FNAY').data(), 0, 2)  # radial particle flux
-    radial_area = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.SY').data(), 0, 1)  # radial contact area
-    for k, sp in enumerate(sim.species_list):
-        i, j = np.where(sim.species_density[:, :-1, k] > 0)  # radial_area array corresponds to [:, 1:] in mesh, so maybe [:, 1:, k]
-        velocities[i, j, k, 1] = radial_particle_flux[i, j, k] / radial_area[i, j] / sim.species_density[i, j, k]
+    # Obtaining additional data from EIRENE and replacing data for neutrals
 
-    sim.velocities = velocities
-    # sim.velocities_cartesian is created authomatically
+    b2_standalone = False
+    try:
+        # Replace the species densities
+        neutral_density = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.DAB2').data(), 0, 2)
+        species_density[:, :, neutral_indx] = neutral_density
+    except (mdsExceptions.TreeNNF, np.AxisError):
+        print("Warning! This is B2 stand-alone simulation.")
+        b2_standalone = True
+
+    if not b2_standalone:
+        # Obtaining neutral atom velocity from EIRENE flux
+        # Note that if the output for fluxes was turned off, PFLA and RFLA' are all zeros
+        try:
+            neutral_poloidal_flux = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.PFLA').data(), 0, 2)
+            neutral_radial_flux = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.RFLA').data(), 0, 2)
+
+            if np.any(neutral_poloidal_flux) or np.any(neutral_radial_flux):
+                neutral_velocities_cartesian = eirene_flux_to_velocity(mesh, neutral_density, neutral_poloidal_flux, neutral_radial_flux,
+                                                                       parallel_velocity[:, :, neutral_indx], sim.b_field_cartesian)
+
+                velocities_cartesian[:, :, neutral_indx, :] = neutral_velocities_cartesian
+        except (mdsExceptions.TreeNNF, np.AxisError):
+            pass
+
+        # Obtaining neutral temperatures
+        try:
+            sim.neutral_temperature = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.TAB2').data(), 0, 2)
+        except (mdsExceptions.TreeNNF, np.AxisError):
+            pass
+
+    sim.species_density = species_density
+    sim.velocities_cartesian = velocities_cartesian  # this also updates sim.velocities
 
     ###############################
     # Load extra data from server #
@@ -131,22 +138,25 @@ def load_solps_from_mdsplus(mds_server, ref_number):
 
     ####################
     # Integrated power #
-    linerad = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.RQRAD').data(), 0, 2)
-    linerad = np.sum(linerad, axis=2)
-    brmrad = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.RQBRM').data(), 0, 2)
-    brmrad = np.sum(brmrad, axis=2)
+    try:
+        linerad = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.RQRAD').data(), 0, 2)
+        linerad = np.sum(linerad, axis=2)
+    except (mdsExceptions.TreeNNF, np.AxisError):
+        linerad = 0
 
-    # need to cope with fact that neurad may not be present!!!
     try:
-        neurad = conn.get('\SOLPS::TOP.SNAPSHOT.ENEUTRAD').data()
-        if isinstance(neurad, np.ndarray):
-            if len(neurad.shape) == 3:
-                neurad = np.swapaxes(np.abs(np.sum(neurad, axis=0)), 0, 1)
-            else:
-                neurad = np.swapaxes(np.abs(neurad), 0, 1)
+        brmrad = np.swapaxes(conn.get('\SOLPS::TOP.SNAPSHOT.RQBRM').data(), 0, 2)
+        brmrad = np.sum(brmrad, axis=2)
+    except (mdsExceptions.TreeNNF, np.AxisError):
+        brmrad = 0
+
+    try:
+        eneutrad = conn.get('\SOLPS::TOP.SNAPSHOT.ENEUTRAD').data()
+        if np.ndim(eneutrad) == 3:  # this will not return error if eneutrad is not np.ndarray
+            neurad = np.swapaxes(np.abs(np.sum(eneutrad, axis=0)), 0, 1)
         else:
-            neurad = 0
-    except mdsExceptions.TreeNNF:
+            neurad = np.swapaxes(np.abs(eneutrad), 0, 1)
+    except (mdsExceptions.TreeNNF, np.AxisError):
         neurad = 0
 
     sim.total_radiation = (linerad + brmrad + neurad) / mesh.vol
@@ -154,30 +164,51 @@ def load_solps_from_mdsplus(mds_server, ref_number):
     return sim
 
 
-def load_mesh_from_mdsplus(mds_connection):
+def load_mesh_from_mdsplus(mds_connection, mdsExceptions):
     """
     Load the SOLPS mesh geometry for a given MDSplus connection.
 
     :param mds_connection: MDSplus connection object. Already set to the SOLPS tree with pulse #ID.
+    :param mdsExceptions: MDSplus mdsExceptions module for error handling.
     """
 
     # Load the R, Z coordinates of the cell vertices, original coordinates are (4, 38, 98)
     # re-arrange axes (4, 38, 98) => (98, 38, 4)
-    x = np.swapaxes(mds_connection.get('\TOP.SNAPSHOT.GRID:R').data(), 0, 2)
+    r = np.swapaxes(mds_connection.get('\TOP.SNAPSHOT.GRID:R').data(), 0, 2)
     z = np.swapaxes(mds_connection.get('\TOP.SNAPSHOT.GRID:Z').data(), 0, 2)
 
     vol = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.VOL').data(), 0, 1)
 
+    # Loading neighbouring cell indices
+    neighbix = np.zeros(r.shape, dtype=np.int)
+    neighbiy = np.zeros(r.shape, dtype=np.int)
+
+    neighbix[:, :, 0] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:LEFTIX').data().astype(np.int), 0, 1)
+    neighbix[:, :, 1] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:BOTTOMIX').data().astype(np.int), 0, 1)
+    neighbix[:, :, 2] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:RIGHTIX').data().astype(np.int), 0, 1)
+    neighbix[:, :, 3] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:TOPIX').data().astype(np.int), 0, 1)
+
+    neighbiy[:, :, 0] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:LEFTIY').data().astype(np.int), 0, 1)
+    neighbiy[:, :, 1] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:BOTTOMIY').data().astype(np.int), 0, 1)
+    neighbiy[:, :, 2] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:RIGHTIY').data().astype(np.int), 0, 1)
+    neighbiy[:, :, 3] = np.swapaxes(mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:TOPIY').data().astype(np.int), 0, 1)
+
+    neighbix[neighbix == r.shape[0]] = -1
+    neighbiy[neighbiy == r.shape[1]] = -1
+
     # build mesh object
-    mesh = SOLPSMesh(x, z, vol)
+    mesh = SOLPSMesh(r, z, vol, neighbix, neighbiy)
 
     #############################
     # Add additional parameters #
     #############################
 
     # add the vessel geometry
-    vessel = mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:VESSEL').data()
-    if isinstance(vessel, np.ndarray):
-        mesh.vessel = vessel
+    try:
+        vessel = mds_connection.get('\SOLPS::TOP.SNAPSHOT.GRID:VESSEL').data()
+        if isinstance(vessel, np.ndarray):
+            mesh.vessel = vessel
+    except mdsExceptions.TreeNNF:
+        pass
 
     return mesh
diff --git a/cherab/solps/formats/raw_simulation_files.py b/cherab/solps/formats/raw_simulation_files.py
index 58c4bcf..2d6033e 100644
--- a/cherab/solps/formats/raw_simulation_files.py
+++ b/cherab/solps/formats/raw_simulation_files.py
@@ -28,7 +28,7 @@
 from cherab.solps.eirene import load_fort44_file
 from cherab.solps.b2.parse_b2_block_file import load_b2f_file
 from cherab.solps.mesh_geometry import SOLPSMesh
-from cherab.solps.solps_plasma import SOLPSSimulation, prefer_element
+from cherab.solps.solps_plasma import SOLPSSimulation, prefer_element, b2_flux_to_velocity, eirene_flux_to_velocity
 
 
 # Code based on script by Felix Reimold (2016)
@@ -59,7 +59,7 @@ def load_solps_from_raw_output(simulation_path, debug=False):
         RuntimeError("No B2 b2fstate file found in SOLPS output directory.")
 
     if not(os.path.isfile(eirene_fort44_file)):
-        print("Warning! No EIRENE fort.44 file found in SOLPS output directory. Assuming B2.5 stand-alone simulation.")
+        print("Warning! No EIRENE fort.44 file found in SOLPS output directory. Assuming B2 stand-alone simulation.")
         b2_standalone = True
     else:
         # Load data for neutral species from EIRENE output file
@@ -69,7 +69,8 @@ def load_solps_from_raw_output(simulation_path, debug=False):
     # Load SOLPS mesh geometry
     _, _, geom_data_dict = load_b2f_file(mesh_file_path, debug=debug)  # geom_data_dict is needed also for magnetic field
 
-    mesh = SOLPSMesh(geom_data_dict['crx'], geom_data_dict['cry'], geom_data_dict['vol'])
+    mesh = create_mesh_from_geom_data(geom_data_dict)
+
     ni = mesh.nx
     nj = mesh.ny
 
@@ -77,6 +78,7 @@ def load_solps_from_raw_output(simulation_path, debug=False):
 
     # Load each plasma species in simulation
     species_list = []
+    neutral_indx = []
     for i in range(len(sim_info_dict['zn'])):
 
         zn = int(sim_info_dict['zn'][i])  # Nuclear charge
@@ -85,10 +87,12 @@ def load_solps_from_raw_output(simulation_path, debug=False):
         isotope = lookup_isotope(zn, number=am)
         species = prefer_element(isotope)  # Prefer Element over Isotope if the mass number is the same
         species_list.append((species, charge))
+        if charge == 0:  # updating neutral index
+            neutral_indx.append(i)
 
     sim = SOLPSSimulation(mesh, species_list)
 
-    # Load magnetic field    
+    # Load magnetic field
     sim.b_field = geom_data_dict['bb'][:, :, :3]
     # sim.b_field_cartesian is created authomatically
 
@@ -100,53 +104,48 @@ def load_solps_from_raw_output(simulation_path, debug=False):
     sim.ion_temperature = mesh_data_dict['ti'] / elementary_charge
 
     # Load species density
-    sim.species_density = mesh_data_dict['na']
+    species_density = mesh_data_dict['na']
 
-    if not b2_standalone:
-        # Replacing B2 neutral densities with EIRENE ones
-        da_raw_data = eirene.da
-        neutral_i = 0  # counter for neutral atoms
-        for k, sp in enumerate(sim.species_list):
-            charge = sp[1]
-            if charge == 0:
-                sim.species_density[1:-1, 1:-1, k] = da_raw_data[:, :, neutral_i]
-                neutral_i += 1
+    # Load parallel velocity
+    parallel_velocity = mesh_data_dict['ua']
 
-    # TODO: Eirene data (TOP.SNAPSHOT.PFLA, TOP.SNAPSHOT.RFLA) should be used for neutral atoms.
-    velocities = np.zeros((ni, nj, len(sim.species_list), 3))
-    velocities[:, :, :, 0] = mesh_data_dict['ua']
+    # Load poloidal and radial particle fluxes for velocity calculation
+    poloidal_flux = mesh_data_dict['fna'][:, :, ::2]
+    radial_flux = mesh_data_dict['fna'][:, :, 1::2]
 
-    ################################################
-    # Calculate the species' velocity distribution #
+    # Obtaining velocity from B2 flux
+    velocities_cartesian = b2_flux_to_velocity(mesh, species_density, poloidal_flux, radial_flux, parallel_velocity, sim.b_field_cartesian)
 
-    # calculate field component ratios for velocity conversion
-    bplane2 = sim.b_field[:, :, 0]**2 + sim.b_field[:, :, 2]**2
-    parallel_to_toroidal_ratio = sim.b_field[:, :, 0] * sim.b_field[:, :, 2] / bplane2
+    if not b2_standalone:
+        # Obtaining additional data from EIRENE and replacing data for neutrals
+        # Note EIRENE data grid is slightly smaller than SOLPS grid, for example (98, 38) => (96, 36)
+        # Need to pad EIRENE data to fit inside larger B2 array
 
-    # Calculate toroidal velocity component
-    velocities[:, :, :, 2] = velocities[:, :, :, 0] * parallel_to_toroidal_ratio[:, :, None]
+        neutral_density = np.zeros((ni, nj, len(neutral_indx)))
+        neutral_density[1:-1, 1:-1, :] = eirene.da
+        species_density[:, :, neutral_indx] = neutral_density
 
-    # Radial velocity is obtained from radial particle flux
-    radial_particle_flux = mesh_data_dict['fna'][:, :, 1::2]
+        # Obtaining neutral atom velocity from EIRENE flux
+        # Note that if the output for fluxes was turned off, eirene.ppa and eirene.rpa are all zeros
+        if np.any(eirene.ppa) or np.any(eirene.rpa):
+            neutral_poloidal_flux = np.zeros((ni, nj, len(neutral_indx)))
+            neutral_poloidal_flux[1:-1, 1:-1, :] = eirene.ppa
 
-    vec_r = mesh.r[:, :, 1] - mesh.r[:, :, 0]
-    vec_z = mesh.z[:, :, 1] - mesh.z[:, :, 0]
-    radial_area = np.pi * (mesh.r[:, :, 1] + mesh.r[:, :, 0]) * np.sqrt(vec_r**2 + vec_z**2)
+            neutral_radial_flux = np.zeros((ni, nj, len(neutral_indx)))
+            neutral_radial_flux[1:-1, 1:-1, :] = eirene.rpa
 
-    for k, sp in enumerate(sim.species_list):
-        i, j = np.where(sim.species_density[:, :, k] > 0)
-        velocities[i, j, k, 1] = radial_particle_flux[i, j, k] / radial_area[i, j] / sim.species_density[i, j, k]
+            neutral_parallel_velocity = np.zeros((ni, nj, len(neutral_indx)))  # must be zero outside EIRENE grid
+            neutral_parallel_velocity[1:-1, 1:-1, :] = parallel_velocity[1:-1, 1:-1, neutral_indx]
 
-    sim.velocities = velocities
-    # sim.velocities_cartesian is created authomatically
+            neutral_velocities_cartesian = eirene_flux_to_velocity(mesh, neutral_density, neutral_poloidal_flux, neutral_radial_flux,
+                                                                   neutral_parallel_velocity, sim.b_field_cartesian)
 
-    if not b2_standalone:
-        # Note EIRENE data grid is slightly smaller than SOLPS grid, for example (98, 38) => (96, 36)
-        # Need to pad EIRENE data to fit inside larger B2 array
+            velocities_cartesian[:, :, neutral_indx, :] = neutral_velocities_cartesian
 
         # Obtaining neutral temperatures
         ta = np.zeros((ni, nj, eirene.ta.shape[2]))
         ta[1:-1, 1:-1, :] = eirene.ta
+        # extrapolating
         for i in (0, -1):
             ta[i, 1:-1, :] = eirene.ta[i, :, :]
             ta[1:-1, i, :] = eirene.ta[:, i, :]
@@ -161,4 +160,37 @@ def load_solps_from_raw_output(simulation_path, debug=False):
 
         sim.eirene_simulation = eirene
 
+    sim.species_density = species_density
+    sim.velocities_cartesian = velocities_cartesian  # this also updates sim.velocities
+
     return sim
+
+def create_mesh_from_geom_data(geom_data):
+
+    r = geom_data['crx']
+    z = geom_data['cry']
+    vol = geom_data['vol']
+
+    # Loading neighbouring cell indices
+    neighbix = np.zeros(r.shape, dtype=np.int)
+    neighbiy = np.zeros(r.shape, dtype=np.int)
+
+    neighbix[:, :, 0] = geom_data['leftix'].astype(np.int)  # poloidal prev.
+    neighbix[:, :, 1] = geom_data['bottomix'].astype(np.int)  # radial prev.
+    neighbix[:, :, 2] = geom_data['rightix'].astype(np.int)  # poloidal next
+    neighbix[:, :, 3] = geom_data['topix'].astype(np.int)  # radial next
+
+    neighbiy[:, :, 0] = geom_data['leftiy'].astype(np.int)
+    neighbiy[:, :, 1] = geom_data['bottomiy'].astype(np.int)
+    neighbiy[:, :, 2] = geom_data['rightiy'].astype(np.int)
+    neighbiy[:, :, 3] = geom_data['topiy'].astype(np.int)
+
+    # In SOLPS cell indexing starts with -1 (guarding cell), but in SOLPSMesh -1 means no neighbour.
+    neighbix += 1
+    neighbiy += 1
+    neighbix[neighbix == r.shape[0]] = -1
+    neighbiy[neighbiy == r.shape[1]] = -1
+
+    mesh = SOLPSMesh(r, z, vol, neighbix, neighbiy)
+
+    return mesh
diff --git a/cherab/solps/mesh_geometry.py b/cherab/solps/mesh_geometry.py
index 52e1d2a..c9422de 100755
--- a/cherab/solps/mesh_geometry.py
+++ b/cherab/solps/mesh_geometry.py
@@ -43,9 +43,17 @@ class SOLPSMesh:
     :param ndarray r: Array of cell vertex r coordinates, must be 3 dimensional. Example shape is (98 x 32 x 4).
     :param ndarray z: Array of cell vertex z coordinates, must be 3 dimensional. Example shape is (98 x 32 x 4).
     :param ndarray vol: Array of cell volumes. Example shape is (98 x 32).
+    :param ndarray neighbix: Array of poloidal indeces of neighbouring cells in order: left, bottom, right, top,
+                             must be 3 dimensional. Example shape is (98 x 32 x 4).
+                             In SOLPS notation: left/right - poloidal prev./next, bottom/top - radial prev./next.
+                             Cell indexing starts with 0 and -1 means no neighbour.
+    :param ndarray neighbix: Array of radial indeces of neighbouring cells in order: left, bottom, right, top,
+                             must be 3 dimensional. Example shape is (98 x 32 x 4).
     """
 
-    def __init__(self, r, z, vol):
+    # TODO Make neighbix and neighbix optional in the future, as they can be reconstructed with _tri_index_loopup
+
+    def __init__(self, r, z, vol, neighbix, neighbiy):
 
         if r.shape != z.shape:
             raise ValueError('Shape of r array: %s mismatch the shape of z array: %s.' % (r.shape, z.shape))
@@ -53,6 +61,12 @@ def __init__(self, r, z, vol):
         if vol.shape != r.shape[:-1]:
             raise ValueError('Shape of vol array: %s mismatch the grid dimentions: %s.' % (vol.shape, r.shape[:-1]))
 
+        if neighbix.shape != r.shape:
+            raise ValueError('Shape of neighbix array must be %s, but it is  %s.' % (r.shape, neighbix.shape))
+
+        if neighbiy.shape != r.shape:
+            raise ValueError('Shape of neighbix array must be %s, but it is  %s.' % (r.shape, neighbiy.shape))
+
         self._cr = r.sum(2) / 4.
         self._cz = z.sum(2) / 4.
 
@@ -64,6 +78,9 @@ def __init__(self, r, z, vol):
 
         self._vol = vol
 
+        self._neighbix = neighbix.astype(np.int)
+        self._neighbiy = neighbiy.astype(np.int)
+
         self.vessel = None
 
         # Calculating poloidal basis vector
@@ -71,18 +88,25 @@ def __init__(self, r, z, vol):
         vec_r = r[:, :, 1] - r[:, :, 0]
         vec_z = z[:, :, 1] - z[:, :, 0]
         vec_magn = np.sqrt(vec_r**2 + vec_z**2)
-        self._poloidal_basis_vector[:, :, 0] = vec_r / vec_magn
-        self._poloidal_basis_vector[:, :, 1] = vec_z / vec_magn
+        self._poloidal_basis_vector[:, :, 0] = np.divide(vec_r, vec_magn, out=np.zeros((self._nx, self._ny)), where=(vec_magn > 0))
+        self._poloidal_basis_vector[:, :, 1] = np.divide(vec_z, vec_magn, out=np.zeros((self._nx, self._ny)), where=(vec_magn > 0))
+
+        # Calculating radial contact areas
+        self._radial_area = np.pi * (r[:, :, 1] + r[:, :, 0]) * vec_magn[:, :]
 
         # Calculating radial basis vector
         self._radial_basis_vector = np.zeros((self._nx, self._ny, 2))
         vec_r = r[:, :, 2] - r[:, :, 0]
         vec_z = z[:, :, 2] - z[:, :, 0]
         vec_magn = np.sqrt(vec_r**2 + vec_z**2)
-        self._radial_basis_vector[:, :, 0] = vec_r / vec_magn
-        self._radial_basis_vector[:, :, 1] = vec_z / vec_magn
+        self._radial_basis_vector[:, :, 0] = np.divide(vec_r, vec_magn, out=np.zeros((self._nx, self._ny)), where=(vec_magn > 0))
+        self._radial_basis_vector[:, :, 1] = np.divide(vec_z, vec_magn, out=np.zeros((self._nx, self._ny)), where=(vec_magn > 0))
+
+        # Calculating poloidal contact areas
+        self._poloidal_area = np.pi * (r[:, :, 2] + r[:, :, 0]) * vec_magn[:, :]
 
         # For convertion from Cartesian to poloidal
+        # TODO Make it work with trianle cells
         self._inv_det = 1. / (self._poloidal_basis_vector[:, :, 0] * self._radial_basis_vector[:, :, 1] -
                               self._poloidal_basis_vector[:, :, 1] * self._radial_basis_vector[:, :, 0])
 
@@ -161,6 +185,26 @@ def vol(self):
         """Volume/area of each grid cell."""
         return self._vol
 
+    @property
+    def neighbix(self):
+        """Poloidal indeces of neighbouring cells in order: left, bottom, right, top."""
+        return self._neighbix
+
+    @property
+    def neighbiy(self):
+        """Radial indeces of neighbouring cells in order: left, bottom, right, top."""
+        return self._neighbiy
+
+    @property
+    def radial_area(self):
+        """Radial contact area."""
+        return self._radial_area
+
+    @property
+    def poloidal_area(self):
+        """Poloidal contact area."""
+        return self._poloidal_area
+
     @property
     def vertex_coordinates(self):
         """RZ-coordinates of unique vertices."""
diff --git a/cherab/solps/solps_plasma.py b/cherab/solps/solps_plasma.py
index 82d0a36..70fe625 100755
--- a/cherab/solps/solps_plasma.py
+++ b/cherab/solps/solps_plasma.py
@@ -56,7 +56,7 @@ def __init__(self, mesh, species_list):
 
         if not len(species_list):
             raise ValueError('Argument "species_list" must contain at least one species.')
-        self._species_list = tuple(species_list)  # adding species is not allowed
+        self._species_list = tuple(species_list)  # adding additional species is not allowed
 
         self._electron_temperature = None
         self._electron_density = None
@@ -167,7 +167,7 @@ def velocities(self):
     def velocities(self, value):
         _check_array("velocities", value, (self.mesh.nx, self.mesh.ny, len(self.species_list), 3))
 
-        # Converting to cartesian system
+        # Converting to Cartesian coordinates
         self._velocities_cartesian = np.zeros(value.shape)
         self._velocities_cartesian[:, :, :, 2] = value[:, :, :, 2]
         for k in range(value.shape[2]):
@@ -183,7 +183,7 @@ def velocities_cartesian(self):
     def velocities_cartesian(self, value):
         _check_array("velocities_cartesian", value, (self.mesh.nx, self.mesh.ny, len(self.species_list), 3))
 
-        # Converting to poloidal system
+        # Converting to poloidal coordinates
         self._velocities = np.zeros(value.shape)
         self._velocities[:, :, :, 2] = value[:, :, :, 2]
         for k in range(value.shape[2]):
@@ -588,3 +588,138 @@ def prefer_element(isotope):
         return isotope.element
 
     return isotope
+
+
+def b2_flux_to_velocity(mesh, density, poloidal_flux, radial_flux, parallel_velocity, b_field_cartesian):
+    """
+    Calculates velocities of neutral particles using B2 particle fluxes defined at cell faces.
+
+    :param SOLPSMesh mesh: SOLPS simulation mesh.
+    :param ndarray density: Density of atoms in m-3. Must be 3 dimensiona array of
+                            shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray poloidal_flux: Poloidal flux of atoms in s-1. Must be a 3 dimensional array of
+                                  shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray radial_flux: Radial flux of atoms in s-1. Must be a 3 dimensional array of
+                                shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray parallel_velocity: Parallel velocity of atoms in m/s. Must be a 3 dimensional
+                                      array of shape (mesh.nx, mesh.ny, num_atoms).
+                                      Parallel velocity is a velocity projection on magnetic
+                                      field direction.
+    :param ndarray b_field_cartesian: Magnetic field in Cartesian (R, Z, phi) coordinates.
+                                      Must be a 3 dimensional array of shape (mesh.nx, mesh.ny, 3).
+
+    :return: Velocities of atoms in (R, Z, phi) coordinates as a 4-dimensional ndarray of
+             shape (mesh.nx, mesh.ny, num_atoms, 3)
+    """
+
+    nx = mesh.nx
+    ny = mesh.ny
+    ns = density.shape[2]
+
+    _check_array('density', poloidal_flux, (nx, ny, ns))
+    _check_array('poloidal_flux', poloidal_flux, (nx, ny, ns))
+    _check_array('radial_flux', radial_flux, (nx, ny, ns))
+    _check_array('parallel_velocity', parallel_velocity, (nx, ny, ns))
+    _check_array('b_field_cartesian', b_field_cartesian, (nx, ny, 3))
+
+    poloidal_area = mesh.poloidal_area[:, :, None]
+    radial_area = mesh.radial_area[:, :, None]
+    leftix = mesh.neighbix[:, :, 0]  # poloidal prev.
+    leftiy = mesh.neighbiy[:, :, 0]
+    bottomix = mesh.neighbix[:, :, 1]  # radial prev.
+    bottomiy = mesh.neighbiy[:, :, 1]
+    rightix = mesh.neighbix[:, :, 2]   # poloidal next.
+    rightiy = mesh.neighbiy[:, :, 2]
+    topix = mesh.neighbix[:, :, 3]  # radial next.
+    topiy = mesh.neighbiy[:, :, 3]
+
+    # Converting s-1 --> m-2 s-1
+    poloidal_flux = np.divide(poloidal_flux, poloidal_area, out=np.zeros_like(poloidal_flux), where=poloidal_area > 0)
+    radial_flux = np.divide(radial_flux, radial_area, out=np.zeros_like(radial_flux), where=radial_area > 0)
+
+    # Obtaining left velocity
+    dens_neighb = density[leftix, leftiy, :]  # density in the left neighbouring cell
+    has_neighbour = ((leftix > -1) * (leftiy > -1))[:, :, None]  # check if has left neighbour
+    neg_flux = (poloidal_flux < 0) * (density > 0)  # will use density in this cell if flux is negative
+    pos_flux = (poloidal_flux > 0) * (dens_neighb > 0) * has_neighbour  # will use density in neighbouring cell if flux is positive
+    velocity_left = np.divide(poloidal_flux, density, out=np.zeros((nx, ny, ns)), where=neg_flux)
+    velocity_left = np.divide(poloidal_flux, dens_neighb, out=velocity_left, where=pos_flux)
+    velocity_left = velocity_left[:, :, :, None] * mesh.poloidal_basis_vector[:, :, None, :]  # to vector in Cartesian
+
+    # Obtaining bottom velocity
+    dens_neighb = density[bottomix, bottomiy, :]
+    has_neighbour = ((bottomix > -1) * (bottomiy > -1))[:, :, None]
+    neg_flux = (radial_flux < 0) * (density > 0)
+    pos_flux = (poloidal_flux > 0) * (dens_neighb > 0) * has_neighbour
+    velocity_bottom = np.divide(radial_flux, density, out=np.zeros((nx, ny, ns)), where=neg_flux)
+    velocity_bottom = np.divide(radial_flux, dens_neighb, out=velocity_bottom, where=pos_flux)
+    velocity_bottom = velocity_bottom[:, :, :, None] * mesh.radial_basis_vector[:, :, None, :]  # to Cartesian
+
+    # Obtaining right and top velocities
+    velocity_right = velocity_left[rightix, rightiy, :, :]
+    velocity_right[(rightix < 0) + (rightiy < 0)] = 0
+
+    velocity_top = velocity_bottom[topix, topiy, :, :]
+    velocity_top[(topix < 0) + (topiy < 0)] = 0
+
+    vcart = np.zeros((nx, ny, ns, 3))  # velocities in Cartesian coordinates
+
+    # Projection of velocity on RZ-plane
+    vcart[:, :, :, :2] = 0.25 * (velocity_bottom + velocity_left + velocity_top + velocity_right)
+
+    # Obtaining toroidal velocity
+    b = b_field_cartesian[:, :, None, :]
+    bmagn = np.sqrt((b * b).sum(3))
+    vcart[:, :, :, 2] = (parallel_velocity * bmagn - vcart[:, :, :, 0] * b[:, :, :, 0] - vcart[:, :, :, 1] * b[:, :, :, 1]) / b[:, :, :, 2]
+
+    return vcart
+
+
+def eirene_flux_to_velocity(mesh, density, poloidal_flux, radial_flux, parallel_velocity, b_field_cartesian):
+    """
+    Calculates velocities of neutral particles using Eirene particle fluxes defined at cell centre.
+
+    :param SOLPSMesh mesh: SOLPS simulation mesh.
+    :param ndarray density: Density of atoms in m-3. Must be 3 dimensiona array of
+                            shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray poloidal_flux: Poloidal flux of atoms in m-2 s-1. Must be a 3 dimensional array of
+                                  shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray radial_flux: Radial flux of atoms in m-2 s-1. Must be a 3 dimensional array of
+                                shape (mesh.nx, mesh.ny, num_atoms).
+    :param ndarray parallel_velocity: Parallel velocity of atoms in m/s. Must be a 3 dimensional
+                                      array of shape (mesh.nx, mesh.ny, num_atoms).
+                                      Parallel velocity is a velocity projection on magnetic
+                                      field direction.
+    :param ndarray b_field_cartesian: Magnetic field in Cartesian (R, Z, phi) coordinates.
+                                      Must be a 3 dimensional array of shape (mesh.nx, mesh.ny, 3).
+
+    :return: Velocities of atoms in (R, Z, phi) coordinates as a 4-dimensional ndarray of
+             shape (mesh.nx, mesh.ny, num_atoms, 3)
+    """
+
+    nx = mesh.nx
+    ny = mesh.ny
+    ns = density.shape[2]
+
+    _check_array('density', poloidal_flux, (nx, ny, ns))
+    _check_array('poloidal_flux', poloidal_flux, (nx, ny, ns))
+    _check_array('radial_flux', radial_flux, (nx, ny, ns))
+    _check_array('parallel_velocity', parallel_velocity, (nx, ny, ns))
+    _check_array('b_field_cartesian', b_field_cartesian, (nx, ny, 3))
+
+    # Obtaining velocity
+    poloidal_velocity = np.divide(poloidal_flux, density, out=np.zeros((nx, ny, ns)), where=(density > 0))
+    radial_velocity = np.divide(radial_flux, density, out=np.zeros((nx, ny, ns)), where=(density > 0))
+
+    vcart = np.zeros((nx, ny, ns, 3))  # velocities in Cartesian coordinates
+
+    # Projection of velocity on RZ-plane
+    vcart[:, :, :, :2] = (poloidal_velocity[:, :, :, None] * mesh.poloidal_basis_vector[:, :, None, :] +
+                          radial_velocity[:, :, :, None] * mesh.radial_basis_vector[:, :, None, :])  # to vector in Cartesian
+
+    # Obtaining toroidal velocity
+    b = b_field_cartesian[:, :, None, :]
+    bmagn = np.sqrt((b * b).sum(3))
+    vcart[:, :, :, 2] = (parallel_velocity * bmagn - (vcart[:, :, :, 0] * b[:, :, :, 0] + vcart[:, :, :, 1] * b[:, :, :, 1])) / b[:, :, :, 2]
+
+    return vcart