diff --git a/ogstools/fe2vtu/__init__.py b/ogstools/fe2vtu/__init__.py
index 19435c26911f4f91e75071728d0f91973428332a..3cf5b46852bfed91e1a81ea3ea7194124a6d307e 100644
--- a/ogstools/fe2vtu/__init__.py
+++ b/ogstools/fe2vtu/__init__.py
@@ -12,6 +12,7 @@ from .fe2vtu import (
 from .tools import (
     get_specific_surface,
     helpFormat,
+    write_cell_boundary_conditions,
     write_point_boundary_conditions,
     write_xml,
 )
@@ -23,5 +24,6 @@ __all__ = [
     "helpFormat",
     "update_geo_mesh",
     "write_xml",
+    "write_cell_boundary_conditions",
     "write_point_boundary_conditions",
 ]
diff --git a/ogstools/fe2vtu/_cli.py b/ogstools/fe2vtu/_cli.py
index 45a6e59e59d2d49b85225260754e93a997534af4..c119bb78a2c559feeb1bc441ada2e09c8933d666 100644
--- a/ogstools/fe2vtu/_cli.py
+++ b/ogstools/fe2vtu/_cli.py
@@ -14,6 +14,7 @@ from ogstools.fe2vtu import (
     get_geo_mesh,
     helpFormat,
     update_geo_mesh,
+    write_cell_boundary_conditions,
     write_point_boundary_conditions,
 )
 
@@ -88,3 +89,4 @@ def cli():
         return
 
     write_point_boundary_conditions(args.output, mesh)
+    write_cell_boundary_conditions(args.output, mesh)
diff --git a/ogstools/fe2vtu/tools.py b/ogstools/fe2vtu/tools.py
index 269287ed7ea5472c00dee547224ddef80876b87e..9eb4683020576cf0da2055b295261aa17952c7e8 100644
--- a/ogstools/fe2vtu/tools.py
+++ b/ogstools/fe2vtu/tools.py
@@ -63,6 +63,8 @@ def write_xml(mesh_name: str, data: pv.DataSetAttributes, mesh_type: str):
         "2ND": "Neumann",
         "3RD": "Robin",
         "4TH": "NodalSourceTerm",
+        "P_IOFLOW": "Neumann?",
+        "P_SOUF": "Neumann?",
     }
     mesh_name = mesh_name.replace(".vtu", "")
     xml_meshes = ET.Element("meshes")
@@ -107,7 +109,7 @@ def write_point_boundary_conditions(mesh_name: str, mesh: pv.UnstructuredGrid):
     """
     Writes the point boundary condition of the mesh. It works by iterating all point data and looking for
     data arrays that include the string "_BC". Depending on what follows, it defines the boundary condition type.
-    This function also writes then the corresponding xml-files using "write_xml"
+    This function also writes then the corresponding xml-files using the function "write_xml"
 
     :param mesh_name: name of the mesh
     :type mesh_name: str
@@ -153,10 +155,20 @@ def write_point_boundary_conditions(mesh_name: str, mesh: pv.UnstructuredGrid):
     write_xml(mesh_name, mesh.point_data, "MeshNode")
 
 
-"""
-def write_cell_bc():
+def write_cell_boundary_conditions(mesh_name: str, mesh: pv.UnstructuredGrid):
+    """
+    Writes the cell boundary condition of the mesh. It works by iterating all cell data and looking for
+    data arrays that include the strings "P_SOUF" or "P_IOFLOW".
+    This function also writes then the corresponding xml-files using the function "write_xml".
+    +++WARNING+++: This function still in a experimental state since it is not clear how exactly this function will
+    be used in the future.
+
+    :param mesh_name: name of the mesh
+    :type mesh_name: str
+    :param mesh: mesh
+    :type mesh: pyvista.UnstructuredGrid
+    """
     BC_mesh = mesh.copy()
-    # print("no nan:", [x for x in mesh["P_BCFLOW_4TH"] if not np.isnan(x)])
     for cd in [
         cell_data
         for cell_data in BC_mesh.cell_data
@@ -166,16 +178,14 @@ def write_cell_bc():
     # Only cell data are needed
     BC_mesh.point_data.clear()
     # get the topsurface since there are the cells of interest
+    # TODO: Allow a generic definition of the normal vector for the filter condition.
     topsurf = get_specific_surface(
         BC_mesh.extract_surface(), lambda normals: normals[:, 2] > 0
     )
-    topsurf.save("topsurface_" + args.output)
-    # pv.save_meshio("topsurface_" + args.output, topsurf, file_format="vtu")
+    topsurf.save("topsurface_" + mesh_name)
     # create the xml-file
     write_xml(
-        "topsurface_" + args.output,
-        "Neumann",
+        "topsurface_" + mesh_name,
         topsurf.cell_data,
         "MeshElement",
     )
-"""