Repository: pyccel/pyccel (language: python)

Issue:
Strides not implemented for Fortran preventing passing slices to functions
I have the pyccelized function

```python
def fun(x : 'int[:]'):
    for i in range(len(x)):
    	print(x[i], ' ')
    print()
```

in `file.py` which just prints the entries of the integer input array `x`.

Running the code 

```python
import numpy as np
from file import fun

A = np.ones((5, 5))
nz = np.nonzero(A)

print('Array passed to pyccelized function : ', nz[0])
print('-----------------------------------------------')

# pass array nz[0]
print(nz[0].flags)
fun(nz[0])
print('-----------------------------------------------')

# pass a copy of the array nz[0]
print(nz[0].copy().flags)
fun(nz[0].copy())
```

reveals that passing the output `nz[0]` of `np.nonzero(A)` to the pyccelized function `fun` gives a wrong order because `nz[0]` is neither F nor C contiguous. Using a copy `nz[0].copy()` instead (which is both F and C contiguous) fixes the problem.

Output:
```
Array passed to pyccelized function :  [0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4]
-----------------------------------------------
  C_CONTIGUOUS : False
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

0 0 0 1 0 2 0 3 0 4 1 0 1 1 1 2 1 3 1 4 2 0 2 1 2 
-----------------------------------------------
  C_CONTIGUOUS : True
  F_CONTIGUOUS : True
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False

0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4
```

**Language**
Fortran


Discussion / hints:
Hi,
Thanks for your interest in pyccel, and thanks for opening this issue. I thought we already had an issue for this but it seems we forgot to open it! The problem arises because we don't currently pass the strides to Fortran. In C this is done so your example works when translating to C

Fix patch:
diff --git a/.pylintrc b/.pylintrc
index 002b237a5a..a4900f9ead 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -16,6 +16,8 @@ ignore=CVS
 # regex matches against base names, not paths.
 ignore-patterns=
 
+ignore-paths=obsolete/*
+
 # Python code to execute, usually for sys.path manipulation such as
 # pygtk.require().
 #init-hook=
@@ -405,7 +407,7 @@ name-group=
 
 # Regular expression which should only match function or class names that do
 # not require a docstring.
-no-docstring-rgx=^(_visit_)|(_print_)|(__init__)|(__new__)
+no-docstring-rgx=^(_visit_)|(_print_)|(__.*__)
 
 # List of decorators that produce properties, such as abc.abstractproperty. Add
 # to this list to register other decorators that produce valid properties.
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 95519d79ef..a173dcb006 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,9 +8,14 @@ All notable changes to this project will be documented in this file.
 ### Fixed
 
 -   #1302 : Raise error message in case of empty class
+-   #929 : Allow optional variables when compiling with intel or nvidia.
+-   #1117 : Allow non-contiguous arrays to be passed to Fortran code.
 
 ### Changed
+
 -   #1355 : Remove unused `BasicParser` arguments.
+-   \[INTERNALS\] Re-write bind-c wrapping stage (#1388)
+
 ## \[1.7.4\] - 2023-05-02
 
 ### Added
diff --git a/ci_tools/list_docs_tovalidate.py b/ci_tools/list_docs_tovalidate.py
index 53de77845b..4e66f9e953 100644
--- a/ci_tools/list_docs_tovalidate.py
+++ b/ci_tools/list_docs_tovalidate.py
@@ -31,7 +31,10 @@ def should_ignore(name):
     if 'SyntaxParser._visit_' in name:
         return True
     #ignore _print_ methods in the codegen.printing module
-    if '_print_' in name and 'pyccel.codegen.printing' in name:
+    if 'Printer._print_' in name:
+        return True
+    #ignore _wrap_ methods in the codegen.wrapper module
+    if 'Wrapper._wrap_' in name:
         return True
     return False
 
diff --git a/developer_docs/overview.md b/developer_docs/overview.md
index 84e99df63f..41cbc8e539 100644
--- a/developer_docs/overview.md
+++ b/developer_docs/overview.md
@@ -9,7 +9,7 @@ Pyccel's development is split into 4 main stages:
 -   [Syntactic Stage](#Syntactic-stage) (for more details see [syntactic stage](syntactic_stage.md))
 -   [Semantic Stage](#Semantic-stage) (for more details see [semantic stage](semantic_stage.md))
 -   [Code Generation Stage](#Code-generation-stage) (for more details see [codegen stage](codegen_stage.md))
--   [Wrapping Stage](#Wrapping-stage)
+-   [Wrapping Stage](#Wrapping-stage) (for more details see [wrapping stage](wrapping_stage.md))
 -   [Compilation Stage](#Compilation-stage)
 
 ### Syntactic Stage
diff --git a/developer_docs/wrapper_stage.md b/developer_docs/wrapper_stage.md
new file mode 100644
index 0000000000..5560f09dd5
--- /dev/null
+++ b/developer_docs/wrapper_stage.md
@@ -0,0 +1,60 @@
+# Wrapper stage
+
+Python is written in C. In order to create a module which can be imported from Python, we therefore have to print a wrapper using the [Python-C API](https://docs.python.org/3/c-api/index.html). This code must call the Pyccel generated translation of the code. If that code was not generated in C then further wrappers are necessary in order to make the code compatible with C.
+
+For example in Fortran generated code there are often array arguments. These do not exist in C, therefore three or more arguments must be passed to the function instead:
+1.  The data
+2.  The size(s) of the array in each dimension
+3.  The stride(s) between elements in each dimension
+
+When defining a slice in Python it is usual to pass a fourth parameter, namely the start of the slice. Internally NumPy simplifies these three integer quantities (`start`, `end`, `step`) into just a stride and a step. This is why the start is not stored or passed as an argument.
+
+The wrapper stage takes the AST describing the generated code as an input and returns an AST which makes that code available to a target language (C or Python).
+
+It should be noted that when code is translated to Python, no wrapper is required.
+
+The entry point for the class `Wrapper` is the function `wrap`.
+
+## `_wrap`
+
+The `_wrap` function internally calls a function named `_wrap_X`, where `X` is the type of the object.
+These functions must have the form:
+```python
+def _print_ClassName(self, stmt):
+    ...
+    return Y
+```
+Each of these `_wrap_X` functions should internally call the `_wrap` function on each of the elements relevant to the wrapper to obtain AST objects which describe the same information, but in a way that is accessible from the target language.
+
+## Name Collisions
+
+While creating the wrapper, it is often necessary to create multiple variables to fully describe something which was a single object in the original code. In order to facilitate reading the code these objects should have names similar to that of the original variable. As a result a **lot** of care must be taken to avoid name collisions.
+
+In general the following rules should be respected:
+-   Any variables which we wish to be able to retrieve from the scope using the `Scope.find` function must be added to the scope with `Scope.insert_symbol`. We cannot do this with 2 variables with the same name so the names must come from the AST where collisions have already been removed.
+-   Any names saved via `Scope.insert_symbol` must be accessed using `Scope.get_expected_name` (as the name may have been changed to avoid other collisions).
+-   Any new names must be created using `Scope.get_new_name`
+
+## Fortran To C
+
+The Fortran to C wrapper wraps Fortran code to make it callable from C. This module relies heavily on [`iso_c_binding`](https://stackoverflow.com/tags/fortran-iso-c-binding/info). It is used to ensure that the functions can be called correctly and that types are compatible. We do not use `iso_fortran_binding` as it was only introduced in Fortran 2018 and we do not expect compiler support for this in the near future.
+
+### Scalar module variables
+
+Scalar module variables are already accessible from C so no extra work is done here.
+
+### Array module variables
+
+Arrays are not compatible with C. Instead the wrapper prints a function which returns a pointer to the module variable as well as its size information to make the information available from C.
+
+### Functions
+
+In the simplest case a wrapper around a function takes the same arguments, calls the function and then returns the result. The only difference between this function and the original function is the body (where the original function is called) and the `bind(c)` tag which ensures that the function name is accessible from C and is not mangled.
+
+More complex cases include functions with array arguments or results, optional variables and functions as arguments.
+
+Array arguments and results are handled similarly to array module variables, by adding additional variables for the shape and stride information.
+
+Optional arguments are passed as C pointers. An if/else block then determines whether the pointer is assigned or not. This can be quite lengthy, however it is unavoidable for compilation with intel, or nvidia. It is also unavoidable for arrays as it is important not to index an array (to access the strides) which is not present.
+
+Finally the most complex cases such as functions as arguments are simply not printed. Instead these cases raise warnings or errors to alert the user that support is missing.
diff --git a/pyccel/parser/syntax/himi.py b/obsolete/parser/syntax/himi.py
similarity index 83%
rename from pyccel/parser/syntax/himi.py
rename to obsolete/parser/syntax/himi.py
index 6f7a71b5ce..0a26f040cb 100644
--- a/pyccel/parser/syntax/himi.py
+++ b/obsolete/parser/syntax/himi.py
@@ -3,6 +3,7 @@
 # This file is part of Pyccel which is released under MIT License. See the LICENSE file or #
 # go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
 #------------------------------------------------------------------------------------------#
+#pylint: disable=pointless-string-statement
 """
 This module contains the syntax associated to the types.tx grammar
 """
@@ -13,6 +14,38 @@
 from pyccel.ast.core import Variable
 from pyccel.ast.datatypes import datatype, VariableType, FunctionType
 
+""" from pyccel.ast.datatypes
+
+class VariableType(DataType):
+    __slots__ = ('_alias','_rhs','_name')
+
+    def __init__(self, rhs, alias):
+        self._alias = alias
+        self._rhs = rhs
+        self._name = rhs._name
+
+    @property
+    def alias(self):
+        return self._alias
+
+class FunctionType(DataType):
+￼    __slots__ = ('_domain','_codomain','_domains','_name')
+￼
+￼    def __init__(self, domains):
+￼        self._domain = domains[0]
+￼        self._codomain = domains[1:]
+￼        self._domains = domains
+￼        self._name = ' -> '.join('{}'.format(V) for V in self._domains)
+￼
+￼    @property
+￼    def domain(self):
+￼        return self._domain
+￼
+￼    @property
+￼    def codomain(self):
+￼        return self._codomain
+"""
+
 
 def _construct_dtype(dtype):
     """."""
diff --git a/pyccel/ast/bind_c.py b/pyccel/ast/bind_c.py
index 7bbc976385..07f945abd5 100644
--- a/pyccel/ast/bind_c.py
+++ b/pyccel/ast/bind_c.py
@@ -3,46 +3,54 @@
 # This file is part of Pyccel which is released under MIT License. See the LICENSE file or #
 # go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
 #------------------------------------------------------------------------------------------#
+"""
+Module describing all elements of the AST needed to represent elements which appear in a Fortran-C binding
+file.
+"""
 
 from pyccel.ast.basic import Basic
-from pyccel.ast.core import CodeBlock, FunctionCall, Module
-from pyccel.ast.core import FunctionAddress
+from pyccel.ast.core import Module
 from pyccel.ast.core import FunctionDef
 from pyccel.ast.core import FunctionDefArgument, FunctionDefResult
-from pyccel.ast.core import Assign
-from pyccel.ast.core import Import
-from pyccel.ast.core import AsName
-from pyccel.ast.core import Allocate
 from pyccel.ast.datatypes import DataType, NativeInteger
-from pyccel.ast.variable import Variable
-from pyccel.parser.scope import Scope
 
 __all__ = (
     'BindCFunctionDef',
+    'BindCFunctionDefArgument',
+    'BindCFunctionDefResult',
+    'BindCModule',
     'BindCPointer',
     'CLocFunc',
-    'as_static_function',
-    'as_static_function_call',
-    'as_static_module',
-    'sanitize_arguments',
-    'wrap_array',
-    'wrap_module_array_var',
+    'C_F_Pointer',
 )
 
+
 class BindCFunctionDef(FunctionDef):
     """
-    Contains the c-compatible version of the function which is
+    Represents the definition of a C-compatible function.
+
+    Contains the C-compatible version of the function which is
     used for the wrapper.
     As compared to a normal FunctionDef, this version contains
     arguments for the shape of arrays. It should be generated by
-    calling ast.bind_c.as_static_function
+    calling `codegen.wrapper.FortranToCWrapper.wrap`.
 
     Parameters
     ----------
-    *args : See FunctionDef
+    *args : list
+        See FunctionDef.
 
     original_function : FunctionDef
-        The function from which the c-compatible version was created
+        The function from which the C-compatible version was created.
+
+    **kwargs : dict
+        See FunctionDef.
+
+    See Also
+    --------
+    pyccel.ast.core.FunctionDef
+        The class from which BindCFunctionDef inherits which contains all
+        details about the args and kwargs.
     """
     __slots__ = ('_original_function',)
     _attribute_nodes = (*FunctionDef._attribute_nodes, '_original_function')
@@ -51,333 +59,445 @@ def __init__(self, *args, original_function, **kwargs):
         self._original_function = original_function
         super().__init__(*args, **kwargs)
         assert self.name == self.name.lower()
+        assert all(isinstance(a, BindCFunctionDefArgument) for a in self._arguments)
+        assert all(isinstance(a, BindCFunctionDefResult) for a in self._results)
 
     @property
     def original_function(self):
-        """ The function which is wrapped by this BindCFunctionDef
+        """
+        The function which is wrapped by this BindCFunctionDef.
+
+        The original function which would be printed in pure Fortran which is not
+        compatible with C.
         """
         return self._original_function
 
-#=======================================================================================
-def sanitize_arguments(args):
-    """ Ensure that all arguments are of expected types (Variable or FunctionAddress)
-    """
-    _args = []
-    for a in args:
-        if isinstance(a.var, (Variable, FunctionAddress)):
-            _args.append(a.var)
+    @property
+    def bind_c_arguments(self):
+        """
+        Get the BindCFunctionDefArguments of the function.
 
-        else:
-            raise NotImplementedError('TODO for {}'.format(type(a)))
+        Return a list of all the arguments passed to the function.
+        These objects all have the type BindCFunctionDefArgument so
+        shapes and strides are hidden.
+        """
+        return self._arguments
+
+    @property
+    def bind_c_results(self):
+        """
+        Get the BindCFunctionDefResults of the function.
 
-    return _args
+        Return a list of all the results returned by the function.
+        These objects all have the type BindCFunctionDefResult so
+        shapes and strides are hidden.
+        """
+        return self._results
+
+    @property
+    def results(self):
+        """
+        List of all objects returned by the function.
 
-#=======================================================================================
-def as_static_function(func, *, mod_scope, name=None):
+        A list of all objects returned by the function including variables
+        which contain array metadata.
+        """
+        return [ai for a in self._results for ai in a.get_all_function_def_results()]
+
+    @property
+    def arguments(self):
+        """
+        List of all arguments passed to the function.
+
+        List of all arguments passed to the function including variables
+        which contain array metadata.
+        """
+        return [ai for a in self._arguments for ai in a.get_all_function_def_arguments()]
+
+# =======================================================================================
+
+
+class BindCFunctionDefArgument(FunctionDefArgument):
     """
-    Translate a FunctionDef to a BindCFunctionDef.
+    Stores all the information necessary to expose an argument to C code.
 
-    Translate a FunctionDef to a BindCFunctionDef by altering the
-    arguments to allow the function to be called from C.
-    E.g. the size of each dimension of an array is provided.
+    Arguments of a C-compatible function may need additional information
+    in order to fully construct the object. This class is mostly important
+    for array objects. These objects must pass not only the data, but also
+    meta-data. Namely the shape and strides for the array in each dimension.
+    This information is stored in this class.
 
     Parameters
     ----------
-    func : FunctionDef
-        The function to be translated.
-    mod_scope : Scope
-        The scope of the module which contains the function.
-    name : str
-        The new name of the function.
-
-    Returns
-    -------
-    BindCFunctionDef
-        The function which can be called from C.
+    var : Variable
+        The variable being passed as an argument (with a C-compatible type).
+
+    scope : pyccel.parser.scope.Scope
+        The scope in which any arguments to the function should be declared.
+        This is used to create the shape and stride variables.
+
+    original_arg_var : Variable
+        The variable which was passed to the function currently being wrapped
+        in a C-Fortran interface. This variable may have a type which is not
+        compatible with C.
+
+    **kwargs : dict
+        See FunctionDefArgument.
+
+    See Also
+    --------
+    pyccel.ast.core.FunctionDefArgument
+        The class from which BindCFunctionDefArgument inherits which
+        contains all details about the args and kwargs.
     """
+    __slots__ = ('_shape', '_strides', '_original_arg_var', '_rank')
+    _attribute_nodes = FunctionDefArgument._attribute_nodes + \
+                        ('_shape', '_strides', '_original_arg_var')
+
+    def __init__(self, var, scope, original_arg_var, **kwargs):
+        name = var.name
+        self._rank = original_arg_var.rank
+        shape   = [scope.get_temporary_variable(NativeInteger(),
+                            name=f'{name}_shape_{i+1}')
+                   for i in range(self._rank)]
+        strides = [scope.get_temporary_variable(NativeInteger(),
+                            name=f'{name}_stride_{i+1}')
+                   for i in range(self._rank)]
+        self._shape = shape
+        self._strides = strides
+        self._original_arg_var = original_arg_var
+        super().__init__(var, **kwargs)
 
-    assert(isinstance(func, FunctionDef))
-
-    args    = list(func.arguments)
-    results = [r.var for r in func.results]
-    body    = func.body.body
-    functions = func.functions
-    _results = []
-    interfaces = func.interfaces
-
-    scope = mod_scope.new_child_scope(func.name)
-
-    # Convert array results to inout arguments
-    for r in results:
-        if r.rank > 0 and r not in args:
-            #wrap the array that is returned from the original function
-            array_body, array_vars = wrap_array(r, scope, False)
-            scope.insert_variable(array_vars[-1])
-            scope.insert_variable(r)
-            body = body + tuple(array_body)
-            array_vars.pop(-1)
-            _results += array_vars
-        elif r.rank == 0:
-            _results += [r]
-
-    body = CodeBlock(body)
-
-    if name is None:
-        name = 'bind_c_{}'.format(func.name).lower()
-
-    # ...
-    results_names = [i.name for i in _results]
-    _args = []
-
-    for i_a, a in enumerate(args):
-        a = a.var
-        if not isinstance(a, (Variable, FunctionAddress)):
-            raise TypeError('Expecting a Variable or FunctionAddress type for {}'.format(a))
-        if not isinstance(a, FunctionAddress) and a.rank > 0:
-            # ...
-            additional_args = []
-            for i in range(a.rank):
-                n_name = 'n{i}_{name}'.format(name=a.name, i=i)
-                n_arg  = Variable('int', n_name)
-
-                additional_args += [n_arg]
-
-            shape_new = tuple(additional_args)
-            # ...
-
-            _args += [FunctionDefArgument(a) for a in additional_args]
-
-            a_new = Variable( a.dtype, a.name,
-                              memory_handling = a.memory_handling,
-                              is_argument = True,
-                              is_optional = a.is_optional,
-                              shape       = shape_new,
-                              rank        = a.rank,
-                              order       = a.order,
-                              precision   = a.precision)
-
-            _args.append(FunctionDefArgument(a_new))
+    @property
+    def original_function_argument_variable(self):
+        """
+        The argument which was passed to the function currently being wrapped.
 
+        The Variable which was passed to the function currently being wrapped
+        in a C-Fortran interface. This variable may have a type which is not
+        compatible with C.
+        """
+        return self._original_arg_var
+
+    @property
+    def shape(self):
+        """
+        The shape of the array argument in each dimension.
+
+        A tuple containing the variables which describe the number of
+        elements along each dimension of an array argument. These values
+        must be passed to any C-compatible function taking an array as an
+        argument.
+        """
+        return self._shape
+
+    @property
+    def strides(self):
+        """
+        The strides of the array argument in each dimension.
+
+        A tuple containing the variables which describe the strides of
+        an array argument in each dimension. These values must be passed to
+        any C-compatible function taking an array as an argument.
+        """
+        return self._strides
+
+    def get_all_function_def_arguments(self):
+        """
+        Get all argument variables which must be printed to fully describe this argument.
+
+        Get a list of all the arguments to the C-compatible function which are
+        required in order to fully describe this argument. This includes the data
+        for the object itself as well as any sizes or strides necessary to
+        define arrays.
+
+        Returns
+        -------
+        list
+            A list of FunctionDefArguments which will be arguments of a BindCFunctionDef.
+        """
+        args = [self]
+        args += [FunctionDefArgument(size) for size in self.shape]
+        args += [FunctionDefArgument(stride) for stride in self.strides]
+        return args
+
+    def __repr__(self):
+        if self.has_default:
+            argument = str(self.name)
+            value = str(self.value)
+            return f'BindCFunctionDefArgument({argument}={value}, inout={self.inout})'
+        else:
+            return f'BindCFunctionDefArgument({repr(self.name)}, inout={self.inout})'
+
+    @property
+    def inout(self):
+        """
+        Indicates whether the argument may be modified by the function.
+
+        True if the argument may be modified in the function. False if
+        the argument remains constant in the function. For array arguments
+        the inout status of the sizes and strides are also returned.
+        """
+        if self._rank:
+            return [False] + [False, False]*self._rank
         else:
-            _args.append(FunctionDefArgument(a))
-
-    args = _args
-    results = [FunctionDefResult(r) for r in _results]
-    # ...
-    return BindCFunctionDef( name, args, results, body,
-                        is_static = True,
-                        functions = functions,
-                        interfaces = interfaces,
-                        imports = func.imports,
-                        original_function = func,
-                        doc_string = func.doc_string,
-                        scope = scope
-                        )
-
-#=======================================================================================
-def as_static_module(funcs, original_module):
-    """ Create the module contained in the bind_c_mod.f90 file
-    This is the interface between the c code and the fortran code thanks
-    to iso_c_bindings
+            return super().inout
+
+# =======================================================================================
+
+
+class BindCFunctionDefResult(FunctionDefResult):
+    """
+    Stores all the information necessary to expose a result to C code.
+
+    Results of a C-compatible function may need additional information
+    in order to fully construct the object. This class is mostly important
+    for array objects. These objects must describe not only the data, but also
+    meta-data. Namely the shape for the array in each dimension.
+    This information is stored in this class.
 
     Parameters
-    ==========
-    funcs : list of FunctionDef
-            All the functions which may be exposed to c
-    original_module : Module
-            The module being wrapped
+    ----------
+    var : Variable
+        The variable being returned (with a C-compatible type).
+
+    original_res_var : Variable
+        The variable which was returned by the function currently being wrapped
+        in a C-Fortran interface. This variable may have a type which is not
+        compatible with C.
+
+    scope : pyccel.parser.scope.Scope
+        The scope in which any arguments to the function should be declared.
+        This is used to create the shape and stride variables.
+
+    **kwargs : dict
+        See FunctionDefResult.
+
+    See Also
+    --------
+    pyccel.ast.core.FunctionDefResult
+        The class from which BindCFunctionDefResult inherits which
+        contains all details about the args and kwargs.
     """
-    funcs = [f for f in funcs if not f.is_private]
-    variables = [f for f in original_module.variables if not f.is_private]
-    imports = []
-    scope = Scope(used_symbols = original_module.scope.local_used_symbols.copy())
-    bind_c_funcs = [as_static_function_call(f, original_module, scope, imports = imports) for f in funcs]
-    bind_c_arrays = [wrap_module_array_var(v, scope, original_module) for v in variables if v.rank > 0]
-    if isinstance(original_module.name, AsName):
-        name = scope.get_new_name('bind_c_{}'.format(original_module.name.target))
-    else:
-        name = scope.get_new_name('bind_c_{}'.format(original_module.name))
-    return Module(name, (), bind_c_funcs+bind_c_arrays, imports = imports, scope=scope)
-
-#=======================================================================================
-def as_static_function_call(func, mod, mod_scope, name=None, imports = None):
+    __slots__ = ('_shape', '_original_res_var')
+    _attribute_nodes = FunctionDefResult._attribute_nodes + \
+                        ('_shape', '_original_res_var')
+
+    def __init__(self, var, original_res_var, scope, **kwargs):
+        name = original_res_var.name
+        self._shape   = [scope.get_temporary_variable(NativeInteger(),
+                            name=f'{name}_shape_{i+1}')
+                         for i in range(original_res_var._rank)]
+        self._original_res_var = original_res_var
+        super().__init__(var, **kwargs)
+
+    @property
+    def original_function_result_variable(self):
+        """
+        The result returned by the function currently being wrapped.
+
+        The variable which was returned by the function currently being wrapped
+        in a C-Fortran interface. This variable may have a type which is not
+        compatible with C.
+        """
+        return self._original_res_var
+
+    @property
+    def shape(self):
+        """
+        The shape of the array result in each dimension.
+
+        A tuple containing the variables which describe the number of
+        elements along each dimension of an array result. These values
+        must be returned by any C-compatible function returning an array.
+        """
+        return self._shape
+
+    def get_all_function_def_results(self):
+        """
+        Get all result variables which must be printed to fully describe this result.
+
+        Get a list of all the results of the C-compatible function which are
+        required in order to fully describe this result. This includes the data
+        for the object itself as well as any sizes necessary to
+        define arrays.
+
+        Returns
+        -------
+        list
+            A list of FunctionDefResults which will be results of a BindCFunctionDef.
+        """
+        res = [self]
+        res += [FunctionDefResult(size) for size in self.shape]
+        return res
+
+# =======================================================================================
+
+class BindCModule(Module):
     """
-    Create a BindCFunctionDef which calls the provided FunctionDef.
+    Represents a Module which only contains functions compatible with C.
 
-    Translate a FunctionDef to a BindCFunctionDef which calls the
-    original function. A BindCFunctionDef is a FunctionDef where the
-    arguments are altered to allow the function to be called from c.
-    E.g. the size of each dimension of an array is provided.
+    Represents a Module which provides the C-Fortran interface to another module.
+    Both functions and module variables are wrapped in order to be compatible with
+    C.
 
     Parameters
     ----------
-    func : FunctionDef
-        The function to be translated.
-
-    mod : Module
-        The module which contains the function.
+    *args : tuple
+        See `pyccel.ast.core.Module`.
 
-    mod_scope : Scope
-        The scope describing the module.
+    original_module : Module
+        The Module being wrapped.
 
-    name : str
-        The new name of the function.
+    variable_wrappers : list of BindCFunctionDef
+        A list containing all the functions which expose module variables to C.
 
-    imports : list
-        An optional parameter into which any required imports
-        can be collected.
+    **kwargs : dict
+        See `pyccel.ast.core.Module`.
 
-    Returns
-    -------
-    BindCFunctionDef
-        The function which can be called from C.
+    See Also
+    --------
+    pyccel.ast.core.Module
+        The class from which BindCModule inherits which contains all details
+        about the args and kwargs.
     """
+    __slots__ = ('_orig_mod','_variable_wrappers')
+    _attribute_nodes = ('_orig_mod','_variable_wrappers')
+
+    def __init__(self, *args, original_module, variable_wrappers = (), **kwargs):
+        self._orig_mod = original_module
+        self._variable_wrappers = variable_wrappers
+        super().__init__(*args, **kwargs)
 
-    assert isinstance(func, FunctionDef)
-    assert isinstance(mod, Module)
-    mod_name = mod.scope.get_python_name(mod.name)
-
-    # from module import func
-    if imports is None:
-        local_imports = [Import(target=AsName(func, func.name), source=mod_name, mod=mod)]
-    else:
-        imports.append(Import(target=AsName(func, func.name), source=mod_name, mod=mod))
-        local_imports = ()
-
-    # function arguments
-    args = sanitize_arguments(func.arguments)
-    # function body
-    call    = FunctionCall(func, args)
-    results = [r.var for r in func.results]
-    results = results[0] if len(results) == 1 else results
-    stmt    = call if len(func.results) == 0 else Assign(results, call)
-    body    = [stmt]
-
-    # new function declaration
-    new_func = FunctionDef(func.name, func.arguments, func.results, body,
-                       functions = func.functions,
-                       interfaces = func.interfaces,
-                       imports = local_imports,
-                       doc_string = func.doc_string
-                       )
-
-    # make it compatible with c
-    static_func = as_static_function(new_func, name=name, mod_scope=mod_scope)
-
-    return static_func
-
-#=======================================================================================
+    @property
+    def original_module(self):
+        """
+        The module which was wrapped.
+
+        The original module for which this object provides the C-Fortran interface.
+        """
+        return self._orig_mod
+
+    @property
+    def variable_wrappers(self):
+        """
+        Get the wrappers which expose module variables to C.
+
+        Get a list containing all the BindCFunctionDefs which expose module variables to C.
+        """
+        return self._variable_wrappers
+
+# =======================================================================================
 
 class BindCPointer(DataType):
-    """ Datatype representing a c pointer in fortran
+    """
+    Datatype representing a C pointer in Fortran.
+
+    Datatype representing a C pointer in Fortran. This data type is defined
+    in the iso_c_binding module.
     """
     __slots__ = ()
     _name = 'bindcpointer'
 
-#=======================================================================================
+# =======================================================================================
 
 class CLocFunc(Basic):
-    """ Class representing the iso_c_binding function cloc which returns a valid
-    C pointer to the location where an object can be found
+    """
+    Creates a C-compatible pointer to the argument.
+
+    Class representing the iso_c_binding function cloc which returns a valid
+    C pointer to the location where an object can be found.
+
+    Parameters
+    ----------
+    argument : Variable
+        The object which should be pointed to.
+
+    result : Variable of dtype BindCPointer
+        The variable where the C-compatible pointer should be stored.
     """
     __slots__ = ('_arg', '_result')
     _attribute_nodes = ()
+
     def __init__(self, argument, result):
         self._arg = argument
         self._result = result
+        assert result.dtype is BindCPointer()
         super().__init__()
 
     @property
     def arg(self):
-        """ Object which will be pointed at
+        """
+        Pointer target.
+
+        Object which will be pointed at by the result pointer.
         """
         return self._arg
 
     @property
     def result(self):
-        """ The variable in which the pointer is stored
         """
-        return self._result
+        The variable where the C-compatible pointer should be stored.
 
-#=======================================================================================
+        The variable where the C-compatible pointer of dtype BindCPointer
+        should be stored.
+        """
+        return self._result
 
-def wrap_array(var, scope, persistent):
-    """ Function returning the code and local variables necessary to wrap an array
+# =======================================================================================
 
-    Parameters
-    ----------
-    var : Variable
-            The array to be wrapped
-    scope : Scope
-            The current scope (used to find valid names for variables)
-    persistent : bool
-            Indicates whether the variable is persistent in memory or
-            if it needs copying to avoid dead pointers
-
-    Results
-    -------
-    body : list
-            A list describing the lines which must be printed to wrap the array
-    variables : list
-            A list of all new variables necessary to wrap the array. The list
-            contains:
-            - The C Pointer which wraps the array
-            - Variables containing the sizes of the array
-            - The Fortran pointer which will contain a copy of the Fortran data
-              (unless the variable is persistent in memory)
+class C_F_Pointer(Basic):
     """
-    bind_var = Variable(dtype=BindCPointer(), name = scope.get_new_name('bound_'+var.name))
-    sizes = [Variable(dtype=NativeInteger(), name = scope.get_new_name()) for _ in range(var.rank)]
-    assigns = [Assign(sizes[i], var.shape[i]) for i in range(var.rank)]
-    variables = [bind_var, *sizes]
-    if not persistent:
-        ptr_var = var.clone(scope.get_new_name(var.name+'_ptr'),
-                memory_handling='alias')
-        alloc = Allocate(ptr_var, shape=var.shape, order=var.order, status='unallocated')
-        copy  = Assign(ptr_var, var)
-        c_loc = CLocFunc(ptr_var, bind_var)
-        variables.append(ptr_var)
-        body = [*assigns, alloc, copy, c_loc]
-    else:
-        c_loc = CLocFunc(var, bind_var)
-        body = [*assigns, c_loc]
-    return body, variables
-
-#=======================================================================================
-
-def wrap_module_array_var(var, scope, mod):
-    """
-    Get a function which allows a module variable to be accessed.
+    Creates a Fortran array pointer from a C pointer and size information.
 
-    Create a function which exposes a module array variable to C.
-    This allows the module variable to be accessed from the Python
-    code.
+    Represents the iso_c_binding function C_F_Pointer which takes a pointer
+    to an object in C (with dtype BindCPointer) and a list of sizes and returns
+    a Fortran array pointer.
 
     Parameters
     ----------
-    var : Variable
-            The array to be exposed.
-    scope : Scope
-            The current scope (used to find valid names for variables).
-    mod : Module
-            The module where the variable is defined.
-
-    Returns
-    -------
-    FunctionDef
-        A function which wraps an array and can be called from C.
+    c_expr : Variable of dtype BindCPointer
+        The Variable containing the C pointer.
+
+    f_expr : Variable
+        The Variable containing the resulting array.
+
+    shape : list of Variables
+        A list describing the Variables which dictate the size of the array in each dimension.
     """
-    func_name = 'bind_c_'+var.name.lower()
-    func_scope = scope.new_child_scope(func_name)
-    body, necessary_vars = wrap_array(var, func_scope, True)
-    func_scope.insert_variable(necessary_vars[0])
-    result_vars = [FunctionDefResult(v) for v in necessary_vars]
-    import_mod = Import(mod.name, AsName(var,var.name), mod=mod)
-    func = BindCFunctionDef(name = func_name,
-                  body      = body,
-                  arguments = [],
-                  results   = result_vars,
-                  imports   = [import_mod],
-                  scope = func_scope,
-                  original_function = None)
-    return func
+    __slots__ = ('_c_expr', '_f_expr', '_shape')
+    _attribute_nodes = ('_c_expr', '_f_expr', '_shape')
+
+    def __init__(self, c_expr, f_expr, shape = ()):
+        self._c_expr = c_expr
+        self._f_expr = f_expr
+        self._shape = shape
+        super().__init__()
+
+    @property
+    def c_pointer(self):
+        """
+        The Variable containing the C pointer.
+
+        The Variable of dtype BindCPointer which contains the C pointer.
+        """
+        return self._c_expr
+
+    @property
+    def f_array(self):
+        """
+        The Variable containing the resulting array.
+
+        The Variable where the array pointer will be stored.
+        """
+        return self._f_expr
+
+    @property
+    def shape(self):
+        """
+        A list of the sizes of the array in each dimension.
+
+        A list describing the Variables which are passed as arguments, in order to
+        determine the size of the array in each dimension.
+        """
+        return self._shape
diff --git a/pyccel/ast/core.py b/pyccel/ast/core.py
index d256fd1a50..04fb09c901 100644
--- a/pyccel/ast/core.py
+++ b/pyccel/ast/core.py
@@ -1232,10 +1232,12 @@ def __init__(
 
         if init_func:
             init_if = init_func.body.body[0]
-            init_cond = init_if.blocks[0].condition
-            init_var = init_cond.args[0]
-            self._variables.append(init_var)
-            self._variable_inits.append(LiteralFalse())
+            # The init function should always contain an If block unless it is part of a wrapper
+            if isinstance(init_if, If):
+                init_cond = init_if.blocks[0].condition
+                init_var = init_cond.args[0]
+                self._variables.append(init_var)
+                self._variable_inits.append(LiteralFalse())
 
         super().__init__(scope)
 
diff --git a/pyccel/ast/datatypes.py b/pyccel/ast/datatypes.py
index 6d17610a8a..b25086f384 100644
--- a/pyccel/ast/datatypes.py
+++ b/pyccel/ast/datatypes.py
@@ -21,7 +21,6 @@
 #
     'CustomDataType',
     'DataType',
-    'FunctionType',
     'NativeBool',
     'NativeComplex',
     'NativeGeneric',
@@ -34,7 +33,6 @@
     'NativeSymbol',
     'NativeVoid',
     'UnionType',
-    'VariableType',
     'DataTypeFactory',
 #
 # --------- FUNCTIONS -----------
@@ -196,36 +194,6 @@ class NativeGeneric(DataType):
     __slots__ = ()
     _name = 'Generic'
 
-
-# ...
-class VariableType(DataType):
-    __slots__ = ('_alias','_rhs','_name')
-
-    def __init__(self, rhs, alias):
-        self._alias = alias
-        self._rhs = rhs
-        self._name = rhs._name
-
-    @property
-    def alias(self):
-        return self._alias
-
-class FunctionType(DataType):
-    __slots__ = ('_domain','_codomain','_domains','_name')
-
-    def __init__(self, domains):
-        self._domain = domains[0]
-        self._codomain = domains[1:]
-        self._domains = domains
-        self._name = ' -> '.join('{}'.format(V) for V in self._domains)
-
-    @property
-    def domain(self):
-        return self._domain
-
-    @property
-    def codomain(self):
-        return self._codomain
 # ...
 
 
diff --git a/pyccel/ast/numpy_wrapper.py b/pyccel/ast/numpy_wrapper.py
index 850d8e1151..10d3f6543f 100644
--- a/pyccel/ast/numpy_wrapper.py
+++ b/pyccel/ast/numpy_wrapper.py
@@ -37,6 +37,8 @@
     #-------HELPERS ------
     'array_get_dim',
     'array_get_data',
+    'array_get_c_step',
+    'array_get_f_step',
     #-------OTHERS--------
     'get_numpy_max_acceptable_version_file'
 )
@@ -105,6 +107,18 @@ def get_numpy_max_acceptable_version_file():
                                         FunctionDefArgument(Variable(dtype=NativeInteger(), name = 'idx'))],
                            results   = [FunctionDefResult(Variable(dtype=NativeInteger(), name = 'd'))])
 
+# Return the stride of the n-th dimension : function definition in pyccel/stdlib/cwrapper/cwrapper_ndarrays.c
+array_get_c_step = FunctionDef(name    = 'nd_nstep_C',
+                           body      = [],
+                           arguments = [FunctionDefArgument(Variable(dtype=NativeVoid(), name = 'o', memory_handling='alias')),
+                                        FunctionDefArgument(Variable(dtype=NativeInteger(), name = 'idx'))],
+                           results   = [FunctionDefResult(Variable(dtype=NativeInteger(), name = 'd'))])
+array_get_f_step = FunctionDef(name    = 'nd_nstep_F',
+                           body      = [],
+                           arguments = [FunctionDefArgument(Variable(dtype=NativeVoid(), name = 'o', memory_handling='alias')),
+                                        FunctionDefArgument(Variable(dtype=NativeInteger(), name = 'idx'))],
+                           results   = [FunctionDefResult(Variable(dtype=NativeInteger(), name = 'd'))])
+
 # Return the data of ndarray : function definition in pyccel/stdlib/cwrapper/cwrapper_ndarrays.c
 array_get_data  = FunctionDef(name   = 'nd_data',
                            body      = [],
diff --git a/pyccel/codegen/printing/cwrappercode.py b/pyccel/codegen/printing/cwrappercode.py
index e5552d4dcb..a62e9871b8 100644
--- a/pyccel/codegen/printing/cwrappercode.py
+++ b/pyccel/codegen/printing/cwrappercode.py
@@ -8,7 +8,9 @@
 
 from pyccel.codegen.printing.ccode import CCodePrinter
 
-from pyccel.ast.bind_c   import as_static_function, wrap_module_array_var, BindCPointer
+from pyccel.ast.bind_c   import BindCPointer
+from pyccel.ast.bind_c   import BindCModule, BindCFunctionDef, BindCFunctionDefArgument
+from pyccel.ast.bind_c   import BindCFunctionDefResult
 
 from pyccel.ast.builtins import PythonTuple, PythonType
 
@@ -36,7 +38,7 @@
 
 from pyccel.ast.numpy_wrapper   import array_type_check
 from pyccel.ast.numpy_wrapper   import pyarray_to_ndarray
-from pyccel.ast.numpy_wrapper   import array_get_data, array_get_dim
+from pyccel.ast.numpy_wrapper   import array_get_data, array_get_dim, array_get_c_step, array_get_f_step
 
 from pyccel.ast.operators import PyccelEq, PyccelNot, PyccelOr, PyccelAssociativeParenthesis
 from pyccel.ast.operators import PyccelIsNot, PyccelLt, PyccelUnarySub
@@ -56,7 +58,7 @@
 dtype_registry = {('pyobject'     , 0) : 'PyObject',
                   ('pyarrayobject', 0) : 'PyArrayObject',
                   ('void'         , 0) : 'void',
-                  ('bind_c_ptr'   , 0) : 'void*'}
+                  ('bind_c_ptr'   , 0) : 'void'}
 
 module_imports  = [Import('numpy_version', Module('numpy_version',(),())),
             Import('numpy/arrayobject', Module('numpy/arrayobject',(),())),
@@ -124,6 +126,34 @@ def is_c_pointer(self, a):
         else:
             return CCodePrinter.is_c_pointer(self,a)
 
+    def get_python_name(self, scope, obj):
+        """
+        Get the name of object as defined in the original python code.
+
+        Get the name of the object as it was originally defined in the
+        Python code being translated. This name may have changed before
+        the printing stage in the case of name clashes or language interfaces.
+
+        Parameters
+        ----------
+        scope : pyccel.parser.scope.Scope
+            The scope where the object was defined.
+
+        obj : pyccel.ast.basic.Basic
+            The object whose name we wish to identify.
+
+        Returns
+        -------
+        str
+            The original name of the object.
+        """
+        if isinstance(obj, BindCFunctionDef):
+            return scope.get_python_name(obj.original_function.name)
+        elif isinstance(obj, BindCModule):
+            return obj.original_module.name
+        else:
+            return scope.get_python_name(obj.name)
+
     def function_signature(self, expr, print_arg_names = True):
         args = list(expr.arguments)
         if any([isinstance(a.var, FunctionAddress) for a in args]):
@@ -247,31 +277,6 @@ def get_default_assign(self, arg, func_arg, value):
         else:
             raise NotImplementedError('Default values are not implemented for this datatype : {}'.format(func_arg.dtype))
 
-    def get_static_function(self, function):
-        """
-
-        Create a static FunctionDef from the argument used for
-        C/fortran binding.
-        If target language is C return the argument function
-        If target language is Fortran, return a static function which
-        takes both the data and the shapes as arguments
-
-        Parameters
-        ----------
-        function    : FunctionDef
-            FunctionDef holding information needed to create static function
-
-        Returns   :
-        -----------
-        static_func : FunctionDef
-        """
-        if self._target_language == 'fortran':
-            static_func = as_static_function(function, mod_scope = self.scope)
-        else:
-            static_func = function
-
-        return static_func
-
     def static_function_signature(self, expr):
         """
         Get the C representation of the function signature using only basic types.
@@ -317,28 +322,53 @@ def static_function_signature(self, expr):
 
     def get_static_args(self, argument):
         """
-        Create bind_C arguments for arguments rank > 0 in fortran.
-        needed in static function call
-        func(a) ==> static_func(nd_dim(a) , nd_data(a))
-        where nd_data(a) = buffer holding data
-              nd_dim(a)  = size of array
+        Get the value(s) which should be passed for the provided argument.
+
+        Get the value(s) which should be passed to the function for the
+        specified argument. In the case of a BindCFunctionDef (when the
+        target language is Fortran) and an argument with rank > 0,
+        multiple arguments are needed:
+        - buffer holding data.
+        - shape of array in each dimension.
+        - stride for the array in each dimension.
 
         Parameters
         ----------
-        argument    : Variable
-            Variable holding information needed (rank)
+        argument : Variable
+            Variable holding information needed (rank).
 
-        Returns     : List of arguments
-            List that can contains Variables and FunctionCalls
-        -----------
+        Returns
+        -------
+        List of arguments
+            List that can contain Variables and FunctionCalls.
+
+        Examples
+        --------
+        If target language is Fortran:
+        >>> x = Variable('int', 'x', rank=2, order='c')
+        >>> self.get_static_args(x)
+        [&nd_data(x), nd_ndim(x, 0), nd_ndim(x, 1), nd_nstep_C(x, 0), nd_nstep_C(x, 1)]
+
+        If target language is C:
+        >>> x = Variable('int', 'x', rank=2, order='c')
+        >>> self.get_static_args(x)
+        [x]
         """
 
         if self._target_language == 'fortran' and argument.rank > 0:
             arg_address = ObjectAddress(argument)
-            static_args = [
+            static_args = [ObjectAddress(FunctionCall(array_get_data, [arg_address]))]
+            static_args+= [
                 FunctionCall(array_get_dim, [arg_address, i]) for i in range(argument.rank)
             ]
-            static_args.append(ObjectAddress(FunctionCall(array_get_data, [arg_address])))
+            if argument.order == 'C':
+                static_args+= [
+                    FunctionCall(array_get_c_step, [arg_address, i]) for i in range(argument.rank)
+                ]
+            else:
+                static_args+= [
+                    FunctionCall(array_get_f_step, [arg_address, i]) for i in range(argument.rank)
+                ]
         else:
             static_args = [argument]
 
@@ -367,80 +397,98 @@ def get_static_declare_type(self, variable):
         dtype = self.find_in_dtype_registry(dtype, prec)
 
         if self.is_c_pointer(variable):
-            return '{0}*'.format(dtype)
-
-        elif self._target_language == 'fortran' and variable.rank > 0:
-            return '{0}*'.format(dtype)
+            return f'{dtype}*'
 
         else:
-            return '{0}'.format(dtype)
+            return dtype
 
     def get_static_results(self, result):
         """
-        Create bind_C results for results rank > 0 in fortran.
-        needed in static function call
-        func(a) ==> static_func(a.shape[0] , &a->raw_data)
+        Get the value(s) which should be used to collect the provided result.
+
+        Get the value(s) which should be collected from the function for the
+        specified argument. In the case of a BindCFunctionDefResult (when the
+        target language is Fortran) and an argument with rank > 0,
+        multiple outputs are needed:
+        - buffer holding data.
+        - shape of array in each dimension.
+        These must then be used to create the expected result variable.
 
         Parameters
         ----------
-        result      : Variable
-            Variable holding information needed (rank)
+        result : FunctionDefResult
+            The result of the function which we want to collect.
 
         Returns
         -------
         body : list
-            Additional instructions (allocations and pointer assignments) for function body
+            Additional instructions (allocations and pointer assignments) for function body.
 
         static_results : list
-            Expanded list of function arguments corresponding to the given result
+            Expanded list of function arguments corresponding to the given result.
+
+        Examples
+        --------
+        If target language is Fortran:
+        >>> x_res = BindCFunctionDefResult(Variable(BindCPointer(), 'x_ptr', memory_handling='alias'), Variable('int', 'x', rank=2, order='c'), scope)
+        >>> self.get_static_results(x)
+        [allocate(x, [x_shape_0, x_shape_1]), x.raw_data=>x_ptr], [&nd_data(x_ptr), x_shape_0, x_shape_1]
+
+        If target language is C:
+        >>> x_res = FunctionDefResult(Variable('int', 'x', rank=2, order='c'))
+        >>> self.get_static_args(x)
+        [], [Variable(x)]
         """
 
         body = []
-
-        if self._target_language == 'fortran' and result.rank > 0:
-            sizes = [
-                     self.scope.get_temporary_variable(NativeInteger(),
-                         result.name+'_size') for _ in range(result.rank)
-                     ]
-            nd_var = self.scope.get_temporary_variable(dtype_or_var = NativeVoid(),
-                    name = result.name,
-                    memory_handling = 'alias')
-            body.append(Allocate(result, shape = sizes, order = result.order,
+        var_name = self.scope.get_expected_name(result.var.name)
+
+        if isinstance(result, BindCFunctionDefResult) and result.shape:
+            shape = [self.scope.get_temporary_variable(s) for s in result.shape]
+            orig_name = result.original_function_result_variable.name
+            orig_var = self.scope.find(self.scope.get_expected_name(orig_name), category='variables')
+            nd_var = self.scope.find(var_name, category='variables')
+            body.append(Allocate(orig_var, shape = shape, order = orig_var.order,
                 status='unallocated'))
             body.append(AliasAssign(DottedVariable(NativeVoid(), 'raw_data', memory_handling = 'alias',
-                lhs=result), nd_var))
+                lhs=orig_var), nd_var))
 
-            static_results = [ObjectAddress(nd_var), *sizes]
+            static_results = [ObjectAddress(nd_var), *shape]
 
         else:
-            static_results = [result]
+            var = self.scope.find(var_name, category='variables')
+            static_results = [var]
 
         return body, static_results
 
     def _get_static_func_call_code(self, expr, static_func_args, results):
         """
-        Get all code necessary to call the wrapped function
+        Get all code necessary to call the wrapped function.
+
+        Get the function call which calls the underlying translated function
+        being wrapped. This may involve creating new variables in order to
+        call the function in a compatible way.
 
         Parameters
         ----------
-        expr             : FunctionDef
-                           The function being wrapped
+        expr : FunctionDef
+            The function being wrapped.
+
         static_func_args : List of arguments
-                           Arguments compatible with the static function
-        results          : List of results
-                           Results of the wrapped function
+            Arguments compatible with the static function.
+
+        results : List of results
+            Results of the wrapped function.
 
         Returns
         -------
-        body             : List of Basic Nodes
-                           List of nodes describing the instructions which call the
-                           wrapped function
-
+        list of pyccel.ast.basic.Basic
+            List of nodes describing the instructions which call the
+            wrapped function.
         """
         body = []
-        static_function = self.get_static_function(expr)
         if len(results) == 0:
-            body.append(FunctionCall(static_function, static_func_args))
+            body.append(FunctionCall(expr, static_func_args))
         else:
             static_func_results = []
             for r in results:
@@ -449,7 +497,7 @@ def _get_static_func_call_code(self, expr, static_func_args, results):
                 static_func_results.extend(s)
 
             results   = static_func_results if len(static_func_results)>1 else static_func_results[0]
-            func_call = Assign(results,FunctionCall(static_function, static_func_args))
+            func_call = Assign(results,FunctionCall(expr, static_func_args))
             body.insert(0, func_call)
 
         return body
@@ -786,28 +834,33 @@ def get_PyArgParseType(self, variable):
 
         return collect_var
 
-    def get_PyBuildValue(self, variable):
+    def get_PyBuildValue(self, result):
         """
+        Get the necessary objects for calling PyBuildValue.
+
         Responsible for collecting the variable required to build the result
-        and the necessary cast function
+        using the Python function PyBuildValue. Also responsible for creating
+        the necessary cast function which creates this variable.
 
         Parameters
         ----------
-        variable : Variable
-            The variable returned by the translated function
+        result : Variable
+            The variable returned by the translated function.
 
         Returns
         -------
-        collect_var : Variable
-            The variable which will be provided to PyBuild
+        Variable
+            The variable which will be provided to PyBuild.
 
-        cast_func_stmts : functionCall
-            call to cast function responsible for the conversion of one data type into another
+        FunctionCall
+            Call to cast function responsible for the conversion of one data type into another.
         """
+        out_var = getattr(result, 'original_function_result_variable', result.var)
+        name = self.scope.get_expected_name(out_var.name)
+        variable = self.scope.find(name, category='variables')
         if variable.rank != 0:
             self.add_import(cwrapper_ndarray_import)
 
-
         cast_function = FunctionCall(C_to_Python(variable), [ObjectAddress(variable)])
 
         collect_type = PyccelPyObject()
@@ -884,25 +937,26 @@ def get_module_exec_function(self, expr, exec_func_name):
         scope.insert_variable(mod_var)
 
         # Collect module variables from translated code
-        orig_vars_to_wrap = [v for v in expr.variables if not v.is_private]
         body = []
-        if self._target_language == 'fortran':
+        if isinstance(expr, BindCModule):
+            orig_vars_to_wrap = [v for v in expr.original_module.variables if not v.is_private]
+            wrapper_funcs = {f.original_function: f for f in expr.variable_wrappers}
             # Collect python compatible module variables
             vars_to_wrap = []
             for v in orig_vars_to_wrap:
-                if v.rank > 0:
+                if v in wrapper_funcs:
                     # Get pointer to store array data
                     var = scope.get_temporary_variable(dtype_or_var = v,
                             name = v.name,
                             memory_handling = 'alias',
                             rank = 0, shape = None, order = None)
-                    # Create variables to store sizes of array
-                    sizes = [scope.get_temporary_variable(NativeInteger(),
+                    # Create variables to store the shape of the array
+                    shape = [scope.get_temporary_variable(NativeInteger(),
                             v.name+'_size') for _ in range(v.rank)]
                     # Get the bind_c function which wraps a fortran array and returns c objects
-                    var_wrapper = wrap_module_array_var(v, scope, expr)
+                    var_wrapper = wrapper_funcs[v]
                     # Call bind_c function
-                    call = Assign(PythonTuple(ObjectAddress(var), *sizes), FunctionCall(var_wrapper, ()))
+                    call = Assign(PythonTuple(ObjectAddress(var), *shape), FunctionCall(var_wrapper, ()))
                     body.append(call)
 
                     # Create ndarray to store array data
@@ -910,7 +964,7 @@ def get_module_exec_function(self, expr, exec_func_name):
                             name = v.name,
                             memory_handling = 'alias'
                             )
-                    alloc = Allocate(nd_var, shape=sizes, order=nd_var.order, status='unallocated')
+                    alloc = Allocate(nd_var, shape=shape, order=nd_var.order, status='unallocated')
                     body.append(alloc)
                     # Save raw_data into ndarray to obtain useable pointer
                     set_data = AliasAssign(DottedVariable(NativeVoid(), 'raw_data',
@@ -928,8 +982,9 @@ def get_module_exec_function(self, expr, exec_func_name):
                         w.set_fst(assign.fst)
                     vars_to_wrap.append(w)
         else:
+            orig_vars_to_wrap = [v for v in expr.variables if not v.is_private]
             vars_to_wrap = orig_vars_to_wrap
-        var_names = [str(expr.scope.get_python_name(v.name)) for v in orig_vars_to_wrap]
+        var_names = [str(self.get_python_name(expr.scope, v)) for v in orig_vars_to_wrap]
 
         # If there are any variables in the module then add them to the module object
         if vars_to_wrap:
@@ -944,8 +999,7 @@ def get_module_exec_function(self, expr, exec_func_name):
 
         if expr.init_func:
             # Call init function code
-            static_function = self.get_static_function(expr.init_func)
-            body.insert(0,FunctionCall(static_function,[],[]))
+            body.insert(0,FunctionCall(expr.init_func,[],[]))
 
         body.append(Return([LiteralInteger(0)]))
         self.exit_scope()
@@ -994,8 +1048,18 @@ def _print_Interface(self, expr):
         wrapper_results = [self.get_new_PyObject("result")]
         self.scope.insert_variable(wrapper_results[0])
 
+        if isinstance(funcs[0], BindCFunctionDef):
+            example_args = funcs[0].bind_c_arguments
+            example_arg_vars = {(a.original_function_argument_variable if isinstance(a, BindCFunctionDefArgument) else a.var): a \
+                    for a in example_args}
+            arg_names = [a.var.name for a in funcs[0].original_function.arguments]
+        else:
+            example_args = funcs[0].arguments
+            example_arg_vars = {(a.original_function_argument_variable if isinstance(a, BindCFunctionDefArgument) else a.var): a \
+                    for a in example_args}
+            arg_names = [a.name for a in example_args]
+
         # Collect argument names for PyArgParse
-        arg_names         = [a.name for a in funcs[0].arguments]
         keyword_list_name = self.scope.get_new_name('kwlist')
         keyword_list      = PyArgKeywords(keyword_list_name, arg_names)
         wrapper_body      = [keyword_list]
@@ -1008,9 +1072,9 @@ def _print_Interface(self, expr):
         default_value = {} # dict to collect all initialisation needed in the wrapper
         check_var = Variable(dtype = NativeInteger(), name = self.scope.get_new_name("check"))
         scope.insert_variable(check_var, check_var.name)
-        types_dict = OrderedDict((a.var, set()) for a in funcs[0].arguments) #dict to collect each variable possible type and the corresponding flags
+        types_dict = OrderedDict((a, set()) for a in example_arg_vars) #dict to collect each variable possible type and the corresponding flags
         # collect parse arg
-        parse_args = [self.get_PyArgParseType(a.var) for a in funcs[0].arguments]
+        parse_args = [self.get_PyArgParseType(a) for a in example_arg_vars]
 
         # Determine flags which indicate argument type
         argument_type_flags = self._determine_interface_flags(funcs)
@@ -1026,15 +1090,37 @@ def _print_Interface(self, expr):
             self.set_scope(mini_scope)
 
             # update ndarray local variables properties
-            arg_vars = {a.var: a for a in func.arguments}
-            result_vars = [r.var for r in func.results]
+            if isinstance(func, BindCFunctionDef):
+                arg_vars = {(a.original_function_argument_variable if isinstance(a, BindCFunctionDefArgument) else a.var): a \
+                        for a in func.bind_c_arguments}
+            else:
+                arg_vars = {(a.original_function_argument_variable if isinstance(a, BindCFunctionDefArgument) else a.var): a \
+                        for a in func.arguments}
+            results = func.bind_c_results if isinstance(func, BindCFunctionDef) else func.results
+            result_vars = []
+            for r in results:
+                var = r.var
+                name = var.name
+                self.scope.insert_symbol(name)
+                if var.rank > 0:
+                    v = var.clone(self.scope.get_expected_name(name), memory_handling = 'alias')
+                else:
+                    v = var.clone(self.scope.get_expected_name(name))
+                result_vars.append(v)
+                self.scope.insert_variable(v)
+                if isinstance(r, BindCFunctionDefResult) and r.shape:
+                    original_var = r.original_function_result_variable
+                    original_name = original_var.name
+                    self.scope.insert_symbol(original_name)
+                    # Declare as pointer as the array should not free its data when it goes out of scope
+                    self.scope.insert_variable(original_var.clone(self.scope.get_expected_name(original_name), memory_handling = 'alias'))
+
             local_arg_vars = {(v.clone(v.name, memory_handling='alias')
                               if isinstance(v, Variable) and v.rank > 0 or v.is_optional \
                               else v) : a for v,a in arg_vars.items()}
+
             for a in local_arg_vars:
                 mini_scope.insert_variable(a)
-            for r in result_vars:
-                mini_scope.insert_variable(r)
 
             # Loop for all args in every functions and create the corresponding condition and body
             for idx, (p_arg, (f_var, f_arg)) in enumerate(zip(parse_args, local_arg_vars.items())):
@@ -1062,11 +1148,11 @@ def _print_Interface(self, expr):
                 mini_wrapper_func_body += body
 
             # create the corresponding function call
-            mini_wrapper_func_body.extend(self._get_static_func_call_code(func, static_func_args, result_vars))
+            mini_wrapper_func_body.extend(self._get_static_func_call_code(func, static_func_args, results))
 
 
             # Loop for all res in every functions and create the corresponding body and cast
-            for r in result_vars :
+            for r in results :
                 collect_var, cast_func = self.get_PyBuildValue(r)
                 if cast_func is not None:
                     mini_wrapper_func_body.append(AliasAssign(collect_var, cast_func))
@@ -1080,6 +1166,13 @@ def _print_Interface(self, expr):
             # Call free function for C type
             mini_wrapper_func_body += [If(IfSection(PyccelIsNot(i, Nil()), [Deallocate(i)])) if self.is_c_pointer(i) \
                                         else Deallocate(i) for i in local_arg_vars if i.rank > 0]
+            if self._target_language == 'fortran':
+                dealloc_results = [self.scope.find(self.scope.get_expected_name(r.original_function_result_variable.name), category='variables')
+                                    for r in results]
+            else:
+                dealloc_results = result_vars
+            mini_wrapper_func_body += [If(IfSection(PyccelIsNot(i, Nil()), [Deallocate(i)])) if self.is_c_pointer(i) \
+                                else Deallocate(i) for i in dealloc_results if i.rank > 0]
             mini_wrapper_func_body.append(Return(wrapper_results))
             self._to_free_PyObject_list.clear()
 
@@ -1111,7 +1204,7 @@ def _print_Interface(self, expr):
 
         # Parsing Arguments
         parse_node = PyArg_ParseTupleNode(*wrapper_args[1:],
-                                          funcs[0].arguments,
+                                          example_args,
                                           parse_args, keyword_list)
 
         wrapper_body += list(default_value.values())
@@ -1177,11 +1270,12 @@ def f(a, b):
             which should be used to increment the byte flag if the argument at the same index
             has the type expected by the function.
         """
+        orig_funcs = [func.original_function if isinstance(func, BindCFunctionDef) else func for func in funcs]
         argument_type_flags = {func:[] for func in funcs}
-        nargs = len(funcs[0].arguments)
+        nargs = len(orig_funcs[0].arguments)
         step = 1
         for i in range(nargs):
-            interface_args = [func.arguments[i].var for func in funcs]
+            interface_args = [getattr(func.arguments[i], 'original_function_argument_variable', func.arguments[i].var) for func in orig_funcs]
             interface_types = [(a.dtype, a.precision) for a in interface_args]
             possible_types = list(dict.fromkeys(interface_types)) # Remove duplicates but preserve order
             for func, t in zip(funcs, interface_types):
@@ -1308,8 +1402,9 @@ def _print_FunctionDef(self, expr):
         self.set_scope(self.scope.new_child_scope(expr.name))
 
         local_arg_vars = {}
-        for a in expr.arguments:
-            v = a.var
+        args = expr.bind_c_arguments if isinstance(expr, BindCFunctionDef) else expr.arguments
+        for a in args:
+            v = (a.original_function_argument_variable if isinstance(a, BindCFunctionDefArgument) else a.var)
             if isinstance(v, Variable):
                 new_name = self.scope.get_new_name(v.name)
                 if isinstance(v, Variable) and (v.rank > 0 or v.is_optional):
@@ -1322,9 +1417,25 @@ def _print_FunctionDef(self, expr):
                 self.scope.functions[v.name] = v
                 local_arg_vars[v] = a
 
-        result_vars = [v.var.clone(self.scope.get_new_name(v.var.name)) for v in expr.results]
-        for v in result_vars:
+        results = expr.bind_c_results if isinstance(expr, BindCFunctionDef) else expr.results
+        result_vars = []
+        for r in results:
+            var = r.var
+            name = var.name
+            self.scope.insert_symbol(name)
+            if var.rank > 0:
+                v = var.clone(self.scope.get_expected_name(name), memory_handling = 'alias')
+            else:
+                v = var.clone(self.scope.get_expected_name(name))
+            result_vars.append(v)
             self.scope.insert_variable(v)
+            if isinstance(r, BindCFunctionDefResult) and r.shape:
+                original_var = r.original_function_result_variable
+                original_name = original_var.name
+                self.scope.insert_symbol(original_name)
+                # Declare as pointer as the array should not free its data when it goes out of scope
+                new_var = original_var.clone(self.scope.get_expected_name(original_name), memory_handling='alias')
+                self.scope.insert_variable(new_var)
         # update ndarray and optional local variables properties
 
         # Find a name for the wrapper function
@@ -1336,7 +1447,7 @@ def _print_FunctionDef(self, expr):
         self.scope.insert_variable(wrapper_results[0])
 
         # Collect argument names for PyArgParse
-        arg_names         = [a.var.name for a in local_arg_vars.values()]
+        arg_names         = [a.var.name for a in (expr.original_function.arguments if isinstance(expr, BindCFunctionDef) else expr.arguments)]
         keyword_list_name = self.scope.get_new_name('kwlist')
         keyword_list      = PyArgKeywords(keyword_list_name, arg_names)
 
@@ -1382,11 +1493,11 @@ def _print_FunctionDef(self, expr):
         wrapper_body.append(If(IfSection(PyccelNot(parse_node), [Return([Nil()])])))
         wrapper_body.extend(wrapper_body_translations)
 
-        wrapper_body.extend(self._get_static_func_call_code(expr, static_func_args, result_vars))
+        wrapper_body.extend(self._get_static_func_call_code(expr, static_func_args, results))
 
         # Loop over results to carry out necessary casts and collect Py_BuildValue type string
         res_args = []
-        for a in result_vars :
+        for a in results :
             collect_var, cast_func = self.get_PyBuildValue(a)
             if cast_func is not None:
                 wrapper_body.append(AliasAssign(collect_var, cast_func))
@@ -1402,6 +1513,13 @@ def _print_FunctionDef(self, expr):
         # Call free function for C type
         wrapper_body += [If(IfSection(PyccelIsNot(i, Nil()), [Deallocate(i)])) if self.is_c_pointer(i) \
                             else Deallocate(i) for i in local_arg_vars if i.rank > 0]
+        if self._target_language == 'fortran':
+            dealloc_results = [self.scope.find(self.scope.get_expected_name(r.original_function_result_variable.name), category='variables')
+                                for r in results]
+        else:
+            dealloc_results = result_vars
+        wrapper_body += [If(IfSection(PyccelIsNot(i, Nil()), [Deallocate(i)])) if self.is_c_pointer(i) \
+                            else Deallocate(i) for i in dealloc_results if i.rank > 0]
         self._to_free_PyObject_list.clear()
 
         #Return
@@ -1419,50 +1537,40 @@ def _print_FunctionDef(self, expr):
         return CCodePrinter._print_FunctionDef(self, wrapper_func)
 
     def _print_Module(self, expr):
-        scope = Scope()
+        scope = Scope(original_symbols = expr.scope.python_names.copy())
         self.set_scope(scope)
         # The initialisation and deallocation shouldn't be exposed to python
         funcs_to_wrap = [f for f in expr.funcs if f not in (expr.init_func, expr.free_func)]
 
         # Insert declared objects into scope
-        if self._target_language == 'fortran':
-            for f in expr.funcs:
-                scope.insert_symbol('bind_c_'+f.name.lower())
-            for v in expr.variables:
-                if not v.is_private:
-                    if v.rank > 0:
-                        scope.insert_symbol('bind_c_'+v.name.lower())
-                    else:
-                        scope.insert_symbol(v.name.lower())
-        else:
-            for f in expr.funcs:
-                scope.insert_symbol(f.name.lower())
-            for v in expr.variables:
-                if not v.is_private:
-                    scope.insert_symbol(v.name.lower())
-
+        variables = expr.original_module.variables if isinstance(expr, BindCModule) else expr.variables
+        for f in expr.funcs:
+            scope.insert_symbol(f.name.lower())
+        for v in variables:
+            if not v.is_private:
+                scope.insert_symbol(v.name.lower())
+
+        funcs = []
         if self._target_language == 'fortran':
             vars_to_wrap_decs = [Declare(v.dtype, v.clone(v.name.lower()), module_variable=True) \
-                                    for v in expr.variables if not v.is_private and v.rank == 0]
+                                    for v in variables if not v.is_private and v.rank == 0]
+
+            for f in expr.original_module.funcs:
+                if f.is_private:
+                    funcs.append(f)
         else:
             vars_to_wrap_decs = [Declare(v.dtype, v, module_variable=True) \
                                     for v in expr.variables if not v.is_private]
 
-        self._module_name  = expr.name
+        self._module_name  = self.get_python_name(scope, expr)
         sep = self._print(SeparatorComment(40))
 
+        function_signatures = ''.join(f'{self.static_function_signature(f)};\n' for f in expr.funcs)
         if self._target_language == 'fortran':
-            static_funcs = [self.get_static_function(f) for f in expr.funcs]
-        else:
-            static_funcs = expr.funcs
-        function_signatures = ''.join('{};\n'.format(self.static_function_signature(f)) for f in static_funcs)
-        if self._target_language == 'fortran':
-            var_wrappers = [wrap_module_array_var(v, self.scope, expr) \
-                    for v in expr.variables if not v.is_private and v.rank > 0]
-            function_signatures += ''.join('{};\n'.format(self.function_signature(v)) for v in var_wrappers)
+            function_signatures += ''.join(f'{self.static_function_signature(f)};\n' for f in expr.variable_wrappers)
 
         interface_funcs = [f.name for i in expr.interfaces for f in i.functions]
-        funcs = [*expr.interfaces, *(f for f in funcs_to_wrap if f.name not in interface_funcs)]
+        funcs += [*expr.interfaces, *(f for f in funcs_to_wrap if f.name not in interface_funcs)]
 
         self._in_header = True
         decs = ''.join('extern '+self._print(d) for d in vars_to_wrap_decs)
@@ -1477,7 +1585,7 @@ def _print_Module(self, expr):
                                      'METH_VARARGS | METH_KEYWORDS,\n'
                                      '{doc_string}\n'
                                      '}},\n').format(
-                                            name = expr.scope.get_python_name(f.name),
+                                            name = self.get_python_name(expr.scope, f),
                                             wrapper_name = self._function_wrapper_names[f.name],
                                             doc_string = self._print(LiteralString('\n'.join(f.doc_string.comments))) \
                                                         if f.doc_string else '""')
@@ -1510,7 +1618,7 @@ def _print_Module(self, expr):
                 '{method_def_name},\n'
                 '{slots_name}\n'
                 '}};\n'.format(module_def_name = module_def_name,
-                    mod_name = expr.name,
+                    mod_name = self._module_name,
                     method_def_name = method_def_name,
                     slots_name = slots_name))
 
@@ -1519,7 +1627,7 @@ def _print_Module(self, expr):
         init_func = ('PyMODINIT_FUNC PyInit_{mod_name}(void)\n{{\n'
                 'import_array();\n'
                 'return PyModuleDef_Init(&{module_def_name});\n'
-                '}}\n'.format(mod_name=expr.name,
+                '}}\n'.format(mod_name=self._module_name,
                     module_def_name = module_def_name))
 
         # Print imports last to be sure that all additional_imports have been collected
diff --git a/pyccel/codegen/printing/fcode.py b/pyccel/codegen/printing/fcode.py
index b7794c9dff..43d4d5fac5 100644
--- a/pyccel/codegen/printing/fcode.py
+++ b/pyccel/codegen/printing/fcode.py
@@ -15,7 +15,7 @@
 import functools
 
 from pyccel.ast.basic import PyccelAstNode
-from pyccel.ast.bind_c import BindCPointer
+from pyccel.ast.bind_c import BindCPointer, BindCFunctionDef, BindCFunctionDefArgument, BindCModule
 from pyccel.ast.core import get_iterable_ranges
 from pyccel.ast.core import FunctionDef, InlineFunctionDef
 from pyccel.ast.core import SeparatorComment, Comment
@@ -458,11 +458,15 @@ def _print_Module(self, expr):
         # ...
         sep = self._print(SeparatorComment(40))
         interfaces = ''
-        if expr.interfaces:
+        if expr.interfaces and not isinstance(expr, BindCModule):
             interfaces = '\n'.join(self._print(i) for i in expr.interfaces)
 
+        func_strings = []
         if expr.funcs:
-            body += '\n'.join(''.join([sep, self._print(i), sep]) for i in expr.funcs)
+            func_strings = [''.join([sep, self._print(i), sep]) for i in expr.funcs]
+        if isinstance(expr, BindCModule):
+            func_strings += [''.join([sep, self._print(i), sep]) for i in expr.variable_wrappers]
+        body += '\n'.join(func_strings)
         # ...
 
         # ...
@@ -1349,17 +1353,18 @@ def _print_Declare(self, expr):
         # ... TODO improve
         # Group the variables by intent
         var = expr.variable
+        expr_dtype = expr.dtype
         rank            = var.rank
         shape           = var.alloc_shape
-        is_target       = var.is_target
         is_const        = var.is_const
         is_optional     = var.is_optional
         is_private      = var.is_private
-        is_alias        = var.is_alias
+        is_alias        = var.is_alias and not isinstance(expr_dtype, BindCPointer)
         on_heap         = var.on_heap
         on_stack        = var.on_stack
         is_static       = expr.static
         is_external     = expr.external
+        is_target       = var.is_target and not var.is_alias
         intent          = expr.intent
         intent_in = intent and intent != 'out'
 
@@ -1369,13 +1374,11 @@ def _print_Declare(self, expr):
 
         # ... print datatype
         if isinstance(expr.dtype, CustomDataType):
-            dtype = expr.dtype
+            name   = expr_dtype.__class__.__name__
+            prefix = expr_dtype.prefix
+            alias  = expr_dtype.alias
 
-            name   = dtype.__class__.__name__
-            prefix = dtype.prefix
-            alias  = dtype.alias
-
-            if dtype.is_polymorphic or expr.passed_from_dotted:
+            if expr_dtype.is_polymorphic or expr.passed_from_dotted:
                 sig = 'class'
             else:
                 sig = 'type'
@@ -1384,9 +1387,8 @@ def _print_Declare(self, expr):
                 name = name.replace(prefix, '')
             else:
                 name = alias
-            dtype = '{0}({1})'.format(sig, name)
+            dtype = f'{sig}({name})'
         else:
-            expr_dtype = expr.dtype
             dtype = self._print(expr_dtype)
 
         # ...
@@ -1799,78 +1801,6 @@ def _print_Block(self, expr):
                  '{body}\n'
                 'end Block {name}\n').format(name=expr.name, prelude=prelude, body=body_code)
 
-    def _print_BindCFunctionDef(self, expr):
-        name = self._print(expr.name)
-        results = [r.var for r in expr.results]
-
-        self.set_scope(expr.scope)
-        self.scope.functions[expr.name] = expr
-
-        body = self._print(expr.body)
-
-        decs = OrderedDict()
-        for arg in expr.arguments:
-            if arg.inout:
-                intent='inout'
-            else:
-                intent='in'
-
-            arg_var = arg.var
-
-            if isinstance(arg_var, FunctionAddress):
-                # Functions with function addresses as arguments cannot be
-                # exposed to python so there is no need to print their signature
-                return ''
-
-            dec = Declare(arg_var.dtype, arg_var, intent=intent , static=True)
-            decs[arg_var] = dec
-
-        arguments = [a.var for a in expr.arguments]
-
-        for result in results:
-            dec = Declare(result.dtype, result, intent='out', static=True)
-            decs[result] = dec
-
-        decs.update(self._get_external_declarations())
-
-        for i in expr.local_vars:
-            dec = Declare(i.dtype, i)
-            decs[i] = dec
-        vars_to_print = self.scope.variables.values()
-        for v in vars_to_print:
-            if (v not in expr.local_vars) and (v not in expr.results) and (v not in arguments):
-                decs[v] = Declare(v.dtype,v)
-
-        func_type = 'subroutine'
-        func_end  = ''
-        if len(results) == 1 and results[0].rank == 0:
-            func_type = 'function'
-            result = results.pop()
-            func_end = 'result({0})'.format(result.name)
-            dec = Declare(result.dtype, result, static=True)
-            decs[result] = dec
-
-        # ...
-
-        interfaces = '\n'.join(self._print(i) for i in expr.interfaces)
-        arg_code  = ', '.join(self._print(i) for i in chain( arguments, results ))
-        imports   = ''.join(self._print(i) for i in expr.imports)
-        prelude   = ''.join(self._print(i) for i in decs.values())
-        body_code = body
-        doc_string = self._print(expr.doc_string) if expr.doc_string else ''
-
-        self.exit_scope()
-
-        parts = [doc_string,
-                '{0} {1}({2}) bind(c) {3}\n'.format(func_type, name, arg_code, func_end),
-                 imports,
-                'implicit none\n',
-                 prelude,
-                 interfaces,
-                 body_code,
-                 'end {} {}\n'.format(func_type, name)]
-        return '\n'.join(p for p in parts if p)
-
     def _print_FunctionAddress(self, expr):
         return expr.name
 
@@ -1930,13 +1860,22 @@ def function_signature(self, expr, name):
         for i, arg in enumerate(arguments):
             arg_var = arg.var
             if isinstance(arg_var, Variable):
-                if i == 0 and expr.cls_name:
-                    dec = Declare(arg_var.dtype, arg_var, intent='inout', passed_from_dotted = True)
-                elif arg.inout:
-                    dec = Declare(arg_var.dtype, arg_var, intent='inout')
+                if isinstance(arg, BindCFunctionDefArgument) and arg.original_function_argument_variable.rank!=0:
+                    for b_arg,inout in zip(arg.get_all_function_def_arguments(), arg.inout):
+                        v = b_arg.var
+                        if inout:
+                            dec = Declare(v.dtype, v, intent='inout')
+                        else:
+                            dec = Declare(v.dtype, v, intent='in')
+                        args_decs[v] = dec
                 else:
-                    dec = Declare(arg_var.dtype, arg_var, intent='in')
-                args_decs[arg_var] = dec
+                    if i == 0 and expr.cls_name:
+                        dec = Declare(arg_var.dtype, arg_var, intent='inout', passed_from_dotted = True)
+                    elif arg.inout:
+                        dec = Declare(arg_var.dtype, arg_var, intent='inout')
+                    else:
+                        dec = Declare(arg_var.dtype, arg_var, intent='in')
+                    args_decs[arg_var] = dec
 
         # treat case of pure function
         sig = '{0}{1} {2}'.format(rec, func_type, name)
@@ -1982,6 +1921,7 @@ def _print_FunctionDef(self, expr):
                     name = name.replace(i, _default_methods[i])
 
         sig_parts = self.function_signature(expr, name)
+        bind_c = ' bind(c)' if isinstance(expr, BindCFunctionDef) else ''
         prelude = sig_parts.pop('arg_decs')
         decs = OrderedDict()
         functions = [f for f in expr.functions if not f.is_inline]
@@ -2009,7 +1949,7 @@ def _print_FunctionDef(self, expr):
         imports = ''.join(self._print(i) for i in expr.imports)
 
         parts = [doc_string,
-                "{}({}) {}\n".format(sig_parts['sig'], sig_parts['arg_code'], sig_parts['func_end']),
+                f"{sig_parts['sig']}({sig_parts['arg_code']}){bind_c} {sig_parts['func_end']}\n",
                 imports,
                 'implicit none\n',
                 prelude,
@@ -2462,33 +2402,63 @@ def _print_Assert(self, expr):
                 '{body}\n'
                 '{epilog}\n').format(prolog=prolog, body=body, epilog=epilog)
 
+    def _handle_not_none(self, lhs, lhs_var):
+        """
+        Print code for `x is not None` statement.
+
+        Print the code which checks if x is not None. This means different
+        things depending on the type of `x`. If `x` is optional it checks
+        if it is present, if `x` is a C pointer it checks if it points at
+        anything.
+
+        Parameters
+        ----------
+        lhs : str
+            The code representing `x`.
+        lhs_var : Variable
+            The Variable `x`.
+
+        Returns
+        -------
+        str
+            The code which checks if `x is not None`.
+        """
+        if isinstance(lhs_var.dtype, BindCPointer):
+            self._constantImports.setdefault('ISO_C_Binding', set()).add('c_associated')
+            return f'c_associated({lhs})'
+        else:
+            return f'present({lhs})'
 
     def _print_PyccelIs(self, expr):
-        lhs = self._print(expr.lhs)
-        rhs = self._print(expr.rhs)
+        lhs_var = expr.lhs
+        rhs_var = expr.rhs
+        lhs = self._print(lhs_var)
+        rhs = self._print(rhs_var)
         a = expr.args[0]
         b = expr.args[1]
 
-        if isinstance(expr.rhs, Nil):
-            return '.not. present({})'.format(lhs)
+        if isinstance(rhs_var, Nil):
+            return '.not. '+ self._handle_not_none(lhs, lhs_var)
 
         if (a.dtype is NativeBool() and b.dtype is NativeBool()):
-            return '{} .eqv. {}'.format(lhs, rhs)
+            return f'{lhs} .eqv. {rhs}'
 
         errors.report(PYCCEL_RESTRICTION_IS_ISNOT,
                       symbol=expr, severity='fatal')
 
     def _print_PyccelIsNot(self, expr):
-        lhs = self._print(expr.lhs)
-        rhs = self._print(expr.rhs)
+        lhs_var = expr.lhs
+        rhs_var = expr.rhs
+        lhs = self._print(lhs_var)
+        rhs = self._print(rhs_var)
         a = expr.args[0]
         b = expr.args[1]
 
-        if isinstance(expr.rhs, Nil):
-            return 'present({})'.format(lhs)
+        if isinstance(rhs_var, Nil):
+            return self._handle_not_none(lhs, lhs_var)
 
         if a.dtype is NativeBool() and b.dtype is NativeBool():
-            return '{} .neqv. {}'.format(lhs, rhs)
+            return f'{lhs} .neqv. {rhs}'
 
         errors.report(PYCCEL_RESTRICTION_IS_ISNOT,
                       symbol=expr, severity='fatal')
@@ -3097,6 +3067,16 @@ def _print_CLocFunc(self, expr):
         self._constantImports.setdefault('ISO_C_Binding', set()).add('c_loc')
         return f'{lhs} = c_loc({rhs})\n'
 
+#=======================================================================================
+
+    def _print_C_F_Pointer(self, expr):
+        self._constantImports.setdefault('ISO_C_Binding', set()).add('C_F_Pointer')
+        shape = ','.join(self._print(s) for s in expr.shape)
+        if shape:
+            return f'call C_F_Pointer({self._print(expr.c_pointer)}, {self._print(expr.f_array)}, [{shape}])\n'
+        else:
+            return f'call C_F_Pointer({self._print(expr.c_pointer)}, {self._print(expr.f_array)})\n'
+
 #=======================================================================================
 
     def _print_PythonConjugate(self, expr):
diff --git a/pyccel/codegen/python_wrapper.py b/pyccel/codegen/python_wrapper.py
index bff17db26a..88cb8abbd8 100644
--- a/pyccel/codegen/python_wrapper.py
+++ b/pyccel/codegen/python_wrapper.py
@@ -6,10 +6,10 @@
 
 import os
 
-from pyccel.ast.bind_c                      import as_static_module
 from pyccel.ast.numpy_wrapper               import get_numpy_max_acceptable_version_file
 from pyccel.codegen.printing.fcode          import fcode
 from pyccel.codegen.printing.cwrappercode   import CWrapperCodePrinter
+from pyccel.codegen.wrapper.fortran_to_c_wrapper   import FortranToCWrapper
 from pyccel.codegen.utilities      import recompile_object
 from pyccel.codegen.utilities      import copy_internal_library
 from pyccel.codegen.utilities      import internal_libs
@@ -36,6 +36,58 @@ def create_shared_library(codegen,
                           wrapper_compiler,
                           sharedlib_modname=None,
                           verbose = False):
+    """
+    Create a shared library which can be called from Pyccel.
+
+    From a CodePrinter object describing code which has been printed
+    in a target language, create a shared library which can be
+    called from Pyccel. In order to do this the code must be wrapped.
+    First, if the code is not written in C, it must be wrapped to
+    make it callable from C. This intermediary code is printed
+    and compiled. From the C-compatible code a second (first for C)
+    wrapper is created which exposes the C code to Python. This
+    is done via the CWrapper. Finally this new code is compiled
+    to generate the required shared language.
+
+    Parameters
+    ----------
+    codegen : pyccel.codegen.printing.codeprinter.CodePrinter
+        The printer which was used to print the translated code.
+
+    main_obj : pyccel.codegen.compiling.basic.CompileObj
+        The compile object which describes the translated code.
+
+    language : str
+        The language which Pyccel translated to.
+
+    wrapper_flags : iterable
+        Any additional flags which should be used to compile the wrapper.
+
+    pyccel_dirpath : str
+        The path to the directory where the files are created and compiled.
+
+    src_compiler : pyccel.codegen.compiling.compilers.Compiler
+        The compiler which should be used to compile the library.
+
+    wrapper_compiler : pyccel.codegen.compiling.compilers.Compiler
+        The compiler which should be used to compile the wrapper.
+        Often this is the same as src_compiler but it may be different
+        when the language is not C to ensure that src_compiler can link
+        the appropriate language-specific libraries.
+
+    sharedlib_modname : str, default: None
+        The name of the shared library. The default is the name of the
+        module printed by the printer.
+
+    verbose : bool, default: False
+        Indicates if the compiling should be done with verbosity to show the
+        compiler commands.
+
+    Returns
+    -------
+    str
+        The absolute path to the shared library which was created.
+    """
 
     pyccel_stage.set_stage('cwrapper')
 
@@ -60,7 +112,8 @@ def create_shared_library(codegen,
 
     if language == 'fortran':
         # Construct static interface for passing array shapes and write it to file bind_c_MOD.f90
-        bind_c_mod = as_static_module(codegen.routines, codegen.ast)
+        wrapper = FortranToCWrapper()
+        bind_c_mod = wrapper.wrap(codegen.ast)
         bind_c_code = fcode(bind_c_mod, bind_c_mod.name)
         bind_c_filename = '{}.f90'.format(bind_c_mod.name)
 
@@ -75,6 +128,9 @@ def create_shared_library(codegen,
         src_compiler.compile_module(compile_obj=bind_c_obj,
                 output_folder=pyccel_dirpath,
                 verbose=verbose)
+        c_ast = bind_c_mod
+    else:
+        c_ast = codegen.ast
 
     #---------------------------------------
     #     Compile cwrapper from stdlib
@@ -100,7 +156,7 @@ def create_shared_library(codegen,
     codegen.ast.set_name(sharedlib_modname)
     wrapper_codegen = CWrapperCodePrinter(codegen.parser.filename, language)
     Scope.name_clash_checker = name_clash_checkers['c']
-    wrapper_code = wrapper_codegen.doprint(codegen.ast)
+    wrapper_code = wrapper_codegen.doprint(c_ast)
     if errors.has_errors():
         return
 
@@ -112,8 +168,8 @@ def create_shared_library(codegen,
     #--------------------------------------------------------
     #  Compile cwrapper_ndarrays from stdlib (if necessary)
     #--------------------------------------------------------
-    if "ndarrays" in wrapper_codegen.get_additional_imports():
-        for lib_name in ("ndarrays", "cwrapper_ndarrays"):
+    for lib_name in ("ndarrays", "cwrapper_ndarrays"):
+        if lib_name in wrapper_codegen.get_additional_imports():
             stdlib_folder, stdlib = internal_libs[lib_name]
 
             lib_dest_path = copy_internal_library(stdlib_folder, pyccel_dirpath)
diff --git a/pyccel/codegen/wrapper/__init__.py b/pyccel/codegen/wrapper/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/pyccel/codegen/wrapper/fortran_to_c_wrapper.py b/pyccel/codegen/wrapper/fortran_to_c_wrapper.py
new file mode 100644
index 0000000000..e03764933e
--- /dev/null
+++ b/pyccel/codegen/wrapper/fortran_to_c_wrapper.py
@@ -0,0 +1,451 @@
+# coding: utf-8
+#------------------------------------------------------------------------------------------#
+# This file is part of Pyccel which is released under MIT License. See the LICENSE file or #
+# go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
+#------------------------------------------------------------------------------------------#
+"""
+Module describing the code-wrapping class : FortranToCWrapper
+which creates an interface exposing Fortran code to C.
+"""
+import warnings
+from pyccel.ast.bind_c import BindCFunctionDefArgument, BindCFunctionDefResult
+from pyccel.ast.bind_c import BindCPointer, BindCFunctionDef, C_F_Pointer
+from pyccel.ast.bind_c import CLocFunc, BindCModule
+from pyccel.ast.core import Assign, FunctionCall, FunctionCallArgument
+from pyccel.ast.core import Allocate, EmptyNode, FunctionAddress
+from pyccel.ast.core import If, IfSection, Import, Interface
+from pyccel.ast.core import AsName, Module, AliasAssign
+from pyccel.ast.datatypes import NativeNumeric
+from pyccel.ast.internals import Slice
+from pyccel.ast.literals import LiteralInteger, Nil, LiteralTrue
+from pyccel.ast.operators import PyccelIsNot, PyccelMul
+from pyccel.ast.variable import Variable, IndexedElement
+from pyccel.parser.scope import Scope
+from .wrapper import Wrapper
+
+class FortranToCWrapper(Wrapper):
+    """
+    Class for creating a wrapper exposing Fortran code to C.
+
+    A class which provides all necessary functions for wrapping different AST
+    objects such that the resulting AST is C-compatible. This new AST is
+    printed as an intermediary layer.
+    """
+    def __init__(self):
+        self._additional_exprs = []
+        self._wrapper_names_dict = {}
+        super().__init__()
+
+    def _get_function_def_body(self, func, func_def_args, func_arg_to_call_arg, results, handled = ()):
+        """
+        Get the body of the bind c function definition.
+
+        Get the body of the bind c function definition by inserting if blocks
+        to check the presence of optional variables. Once we have ascertained
+        the presence of the variables the original function is called. This
+        code slices array variables to ensure the correct step.
+
+        Parameters
+        ----------
+        func : FunctionDef
+            The function which should be called.
+
+        func_def_args : list of FunctionDefArguments
+            The arguments received by the function.
+
+        func_arg_to_call_arg : dict
+            A dictionary mapping the arguments received by the function to the arguments
+            to be passed to the function call.
+
+        results : list of Variables
+            The Variables where the result of the function call will be saved.
+
+        handled : tuple
+            A list of all variables which have been handled (checked to see if they
+            are present).
+
+        Returns
+        -------
+        list
+            A list of Basic nodes describing the body of the function.
+        """
+        optional = next((a for a in func_def_args if a.original_function_argument_variable.is_optional and a not in handled), None)
+        if optional:
+            args = func_def_args.copy()
+            optional_var = optional.var
+            handled += (optional, )
+            true_section = IfSection(PyccelIsNot(optional_var, Nil()),
+                                    self._get_function_def_body(func, args, func_arg_to_call_arg, results, handled))
+            args.remove(optional)
+            false_section = IfSection(LiteralTrue(),
+                                    self._get_function_def_body(func, args, func_arg_to_call_arg, results, handled))
+            return [If(true_section, false_section)]
+        else:
+            args = [FunctionCallArgument(func_arg_to_call_arg[fa],
+                                         keyword = fa.original_function_argument_variable.name)
+                    for fa in func_def_args]
+            size = [fa.shape[::-1] if fa.original_function_argument_variable.order == 'C' else
+                    fa.shape for fa in func_def_args]
+            stride = [fa.strides[::-1] if fa.original_function_argument_variable.order == 'C' else
+                      fa.strides for fa in func_def_args]
+            orig_size = [[PyccelMul(l,s) for l,s in zip(sz, st)] for sz,st in zip(size,stride)]
+            body = [C_F_Pointer(fa.var, func_arg_to_call_arg[fa].base, s)
+                    for fa,s in zip(func_def_args, orig_size)
+                    if isinstance(func_arg_to_call_arg[fa], IndexedElement)]
+            body += [C_F_Pointer(fa.var, func_arg_to_call_arg[fa])
+                     for fa in func_def_args
+                     if not isinstance(func_arg_to_call_arg[fa], IndexedElement) \
+                        and fa.original_function_argument_variable.is_optional]
+
+            # If the function is inlined and takes an array argument create a pointer to ensure that the bounds
+            # are respected
+            if func.is_inline and any(isinstance(a.value, IndexedElement) for a in args):
+                array_args = {a: self.scope.get_temporary_variable(a.value.base, a.keyword, memory_handling = 'alias') for a in args if isinstance(a.value, IndexedElement)}
+                body += [AliasAssign(v, k.value) for k,v in array_args.items()]
+                args = [FunctionCallArgument(array_args[a], keyword=a.keyword) if a in array_args else a for a in args]
+
+            func_call = Assign(results[0], FunctionCall(func, args)) if len(results) == 1 else \
+                        Assign(results, FunctionCall(func, args))
+            return body + [func_call]
+
+    def _get_call_argument(self, bind_c_arg):
+        """
+        Get the argument which should be passed to the function call.
+
+        The FunctionDefArgument passed to the function may contain additional
+        information which should not be passed to the function being wrapped
+        (e.g. an array with strides should not pass the strides explicitly to
+        the function call, nor should it pass the entire contiguous array).
+        This function extracts the necessary information and returns the object
+        which can be passed to the function call.
+
+        Parameters
+        ----------
+        bind_c_arg : BindCFunctionDefArgument
+            The argument to the wrapped bind_c_X function.
+
+        Returns
+        -------
+        PyccelAstNode
+            An object which can be passed to a function call of the function
+            being wrapped.
+        """
+        original_arg = bind_c_arg.original_function_argument_variable
+        arg_var = self.scope.find(original_arg.name, category='variables')
+        if original_arg.is_ndarray:
+            start = LiteralInteger(1) # C_F_Pointer leads to default Fortran lbound
+            stop = None
+            indexes = [Slice(start, stop, step) for step in bind_c_arg.strides]
+            return IndexedElement(arg_var, *indexes)
+        else:
+            return arg_var
+
+    def _wrap_Module(self, expr):
+        """
+        Create a Module which is compatible with C.
+
+        Create a Module which provides an interface between C and the
+        Module described by expr. This includes wrapping functions,
+        interfaces, classes and module variables.
+
+        Parameters
+        ----------
+        expr : pyccel.ast.core.Module
+            The module to be wrapped.
+
+        Returns
+        -------
+        pyccel.ast.core.Module
+            The C-compatible module.
+        """
+        # Define scope
+        scope = expr.scope
+        mod_scope = Scope(used_symbols = scope.local_used_symbols.copy(), original_symbols = scope.python_names.copy())
+        self.scope = mod_scope
+
+        # Wrap contents
+        funcs_to_wrap = [f for f in expr.funcs if not f.is_private]
+        funcs = [self._wrap(f) for f in funcs_to_wrap]
+        if expr.init_func:
+            init_func = funcs[next(i for i,f in enumerate(funcs_to_wrap) if f == expr.init_func)]
+        else:
+            init_func = None
+        if expr.free_func:
+            free_func = funcs[next(i for i,f in enumerate(funcs_to_wrap) if f == expr.free_func)]
+        else:
+            free_func = None
+        funcs = [f for f in funcs if not isinstance(f, EmptyNode)]
+        interfaces = [self._wrap(f) for f in expr.interfaces]
+        classes = [self._wrap(f) for f in expr.classes]
+        variable_getters = [self._wrap(v) for v in expr.variables if not v.is_private]
+        variable_getters = [v for v in variable_getters if not isinstance(v, EmptyNode)]
+        imports = [Import(expr.name, target = expr, mod=expr)]
+
+        name = mod_scope.get_new_name(f'bind_c_{expr.name.target}')
+        self._wrapper_names_dict[expr.name.target] = name
+
+        self.exit_scope()
+
+        return BindCModule(name, (), funcs, variable_wrappers = variable_getters,
+                init_func = init_func, free_func = free_func,
+                interfaces = interfaces, classes = classes,
+                imports = imports, original_module = expr,
+                scope = mod_scope)
+
+    def _wrap_FunctionDef(self, expr):
+        """
+        Create a C-compatible function which executes the original function.
+
+        Create a function which can be called from C which internally calls the original
+        function. It does this by wrapping the arguments and the results and unrolling
+        the body using self._get_function_def_body to ensure optional arguments are
+        present before accessing them. With all this information a BindCFunctionDef is
+        created which is C-compatible.
+
+        Functions which cannot be wrapped raise a warning and return an EmptyNode. This
+        is the case for functions with functions as arguments.
+
+        Parameters
+        ----------
+        expr : FunctionDef
+            The function to wrap.
+
+        Returns
+        -------
+        BindCFunctionDef
+            The C-compatible function.
+        """
+        name = self.scope.get_new_name(f'bind_c_{expr.name.lower()}')
+        self._wrapper_names_dict[expr.name] = name
+
+        # Create the scope
+        func_scope = self.scope.new_child_scope(name)
+        self.scope = func_scope
+
+        self._additional_exprs = []
+
+        if any(isinstance(a.var, FunctionAddress) for a in expr.arguments):
+            warnings.warn("Functions with functions as arguments cannot be wrapped by pyccel")
+            return EmptyNode()
+
+        # Wrap the arguments and collect the expressions passed as the call argument.
+        func_arguments = [self._wrap(a) for a in expr.arguments]
+        call_arguments = [self._get_call_argument(fa) for fa in func_arguments]
+        func_to_call = {fa : ca for ca, fa in zip(call_arguments, func_arguments)}
+
+        func_results = [self._wrap_FunctionDefResult(r) for r in expr.results]
+
+        func_call_results = [r.var.clone(self.scope.get_expected_name(r.var.name)) for r in expr.results]
+
+        body = self._get_function_def_body(expr, func_arguments, func_to_call, func_call_results)
+
+        body.extend(self._additional_exprs)
+        self._additional_exprs.clear()
+
+        self.exit_scope()
+
+        func = BindCFunctionDef(name, func_arguments, func_results, body, scope=func_scope, original_function = expr,
+                doc_string = expr.doc_string)
+
+        self.scope.functions[name] = func
+
+        return func
+
+    def _wrap_Interface(self, expr):
+        """
+        Create an interface containing only C-compatible functions.
+
+        Create an interface containing only functions which can be called from C
+        from an interface which is not necessarily C-compatible.
+
+        Parameters
+        ----------
+        expr : pyccel.ast.core.Interface
+            The interface to be wrapped.
+
+        Returns
+        -------
+        pyccel.ast.core.Interface
+            The C-compatible interface.
+        """
+        functions = [self.scope.functions[self._wrapper_names_dict[f.name]] for f in expr.functions]
+        functions = [f for f in functions if not isinstance(f, EmptyNode)]
+        return Interface(expr.name, functions, expr.is_argument)
+
+    def _wrap_FunctionDefArgument(self, expr):
+        """
+        Create the equivalent BindCFunctionDefArgument for a C-compatible function.
+
+        Take a FunctionDefArgument and create a BindCFunctionDefArgument describing
+        all the information that should be passed to the C-compatible function in order
+        to be able to create the argument described by `expr`.
+
+        In the case of a scalar numerical the function simply creates a copy of the
+        variable described by the function argument in the local scope.
+
+        In the case of an array, C cannot represent the array natively. Rather it is
+        stored in a pointer. This function therefore creates a variable to represent
+        that pointer. Additionally information about the shape and strides of the array
+        are necessary, however these objects are created by the `BindCFunctionDefArgument`
+        class.
+
+        The objects which describe the argument passed to the `expr` argument of the
+        original function are also created here. However the expressions necessary to
+        collect the information from the BindCFunctionDefArgument in order to create
+        these objects are left for later. This is done to ensure that optionals are
+        handled locally to the function call. This ensures that we do not duplicate if
+        conditions.
+
+        Parameters
+        ----------
+        expr : FunctionDefArgument
+            The argument to be wrapped.
+
+        Returns
+        -------
+        BindCFunctionDefArgument
+            The C-compatible argument.
+        """
+        var = expr.var
+        name = var.name
+        self.scope.insert_symbol(name)
+        collisionless_name = self.scope.get_expected_name(var.name)
+        if var.is_ndarray or var.is_optional:
+            new_var = Variable(BindCPointer(), self.scope.get_new_name(f'bound_{name}'),
+                                is_argument = True, is_optional = False, memory_handling='alias')
+            arg_var = var.clone(collisionless_name, is_argument = False, is_optional = False,
+                                memory_handling = 'alias', allows_negative_indexes=False)
+            self.scope.insert_variable(arg_var)
+        else:
+            new_var = var.clone(collisionless_name)
+        self.scope.insert_variable(new_var)
+
+        return BindCFunctionDefArgument(new_var, value = expr.value, original_arg_var = expr.var,
+                kwonly = expr.is_kwonly, annotation = expr.annotation, scope=self.scope)
+
+    def _wrap_FunctionDefResult(self, expr):
+        """
+        Create the equivalent BindCFunctionDefResult for a C-compatible function.
+
+        Take a FunctionDefResult and create a BindCFunctionDefResult describing
+        all the information that should be returned from the C-compatible function
+        in order to fully describe the result `expr`. This function also adds any
+        expressions necessary to build the C-compatible return value to
+        `self._additional_exprs`.
+
+        In the case of a scalar numerical the function simply creates a local version
+        of the variable described by the function result and returns the
+        BindCFunctionDefResult.
+
+        In the case of an array, C cannot represent the array natively. Rather it is
+        stored in a pointer. This function therefore creates a variable to represent
+        that pointer. Additionally information about the shape and strides of the array
+        are necessary. These objects are created by the `BindCFunctionDefResult`
+        class. The assignment expressions which define the shapes and strides are
+        then stored in `self._additional_exprs` along with the allocation of the
+        pointer.
+
+        Parameters
+        ----------
+        expr : FunctionDefResult
+            The result to be wrapped.
+
+        Returns
+        -------
+        BindCFunctionDefResult
+            The C-compatible result.
+        """
+        var = expr.var
+        name = var.name
+        scope = self.scope
+        # Make name available for later
+        scope.insert_symbol(name)
+        local_var = var.clone(scope.get_expected_name(name))
+
+        if local_var.rank:
+            # Allocatable is not returned so it must appear in local scope
+            scope.insert_variable(local_var, name)
+
+            # Create the C-compatible data pointer
+            bind_var = Variable(dtype=BindCPointer(),
+                                name=scope.get_new_name('bound_'+name),
+                                is_const=True, memory_handling='alias')
+            scope.insert_variable(bind_var)
+
+            result = BindCFunctionDefResult(bind_var, var, scope)
+
+            # Save the shapes of the array
+            self._additional_exprs.extend([Assign(result.shape[i], var.shape[i]) for i in range(var.rank)])
+
+            # Create an array variable which can be passed to CLocFunc
+            ptr_var = var.clone(scope.get_new_name(name+'_ptr'),
+                                memory_handling='alias')
+            scope.insert_variable(ptr_var)
+
+            # Define the additional steps necessary to define and fill ptr_var
+            alloc = Allocate(ptr_var, shape=result.shape,
+                             order=var.order, status='unallocated')
+            copy = Assign(ptr_var, var)
+            c_loc = CLocFunc(ptr_var, bind_var)
+            self._additional_exprs.extend([alloc, copy, c_loc])
+
+            return result
+        else:
+            return BindCFunctionDefResult(local_var, var, scope)
+
+    def _wrap_Variable(self, expr):
+        """
+        Create all objects necessary to expose a module variable to C.
+
+        Create and return the objects which must be printed in the wrapping
+        module in order to expose the variable to C. In the case of scalar
+        numerical values nothing needs to be done so an EmptyNode is returned.
+        In the case of numerical arrays a C-compatible function must be created
+        which returns the array. This is necessary because built-in Fortran
+        arrays are not C-compatible. In the case of classes a C-compatible
+        function is also created which returns a pointer to the class object.
+
+        Parameters
+        ----------
+        expr : pyccel.ast.variables.Variable
+            The module variable.
+
+        Returns
+        -------
+        pyccel.ast.basic.Basic
+            The AST object describing the code which must be printed in
+            the wrapping module to expose the variable.
+        """
+        if expr.rank == 0 and expr.dtype in NativeNumeric:
+            return EmptyNode()
+        else:
+            scope = self.scope
+            func_name = scope.get_new_name('bind_c_'+expr.name.lower())
+            func_scope = scope.new_child_scope(func_name)
+            mod = expr.get_user_nodes(Module)[0]
+            import_mod = Import(mod.name, AsName(expr,expr.name), mod=mod)
+            func_scope.imports['variables'][expr.name] = expr
+
+            # Create the data pointer
+            bind_var = Variable(dtype=BindCPointer(),
+                                name=scope.get_new_name('bound_'+expr.name),
+                                is_const=True, memory_handling='alias')
+            func_scope.insert_variable(bind_var)
+
+            result = BindCFunctionDefResult(bind_var, expr, func_scope)
+            if expr.rank == 0:
+                #assigns = []
+                #c_loc = CLocFunc(expr, bind_var)
+                raise NotImplementedError("Classes cannot be wrapped")
+            else:
+                assigns = [Assign(result.shape[i], expr.shape[i]) for i in range(expr.rank)]
+                c_loc = CLocFunc(expr, bind_var)
+            body = [*assigns, c_loc]
+            func = BindCFunctionDef(name = func_name,
+                          body      = body,
+                          arguments = [],
+                          results   = [result],
+                          imports   = [import_mod],
+                          scope = func_scope,
+                          original_function = expr)
+            return func
diff --git a/pyccel/codegen/wrapper/wrapper.py b/pyccel/codegen/wrapper/wrapper.py
new file mode 100644
index 0000000000..a260e2b12b
--- /dev/null
+++ b/pyccel/codegen/wrapper/wrapper.py
@@ -0,0 +1,113 @@
+# coding: utf-8
+#------------------------------------------------------------------------------------------#
+# This file is part of Pyccel which is released under MIT License. See the LICENSE file or #
+# go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
+#------------------------------------------------------------------------------------------#
+"""
+Module describing the base code-wrapping class : Wrapper.
+"""
+
+from pyccel.parser.scope      import Scope
+from pyccel.errors.errors     import Errors
+from pyccel.errors.messages   import PYCCEL_RESTRICTION_TODO
+
+__all__ = ["Wrapper"]
+
+errors = Errors()
+
+
+class Wrapper:
+    """
+    The base class for code-wrapping subclasses.
+
+    The base class for any classes designed to create a wrapper around code.
+    Such wrappers are necessary to create an interface between two different
+    languages.
+    """
+    start_language = None
+    target_language = None
+
+    def __init__(self):
+        self._scope = None
+
+    @property
+    def scope(self):
+        """
+        Get the current scope.
+
+        Get the scope for the current context.
+
+        See Also
+        --------
+        pyccel.parser.scope.Scope
+            The type of the returned object.
+        """
+        return self._scope
+
+    @scope.setter
+    def scope(self, scope):
+        assert isinstance(scope, Scope)
+        self._scope = scope
+
+    def exit_scope(self):
+        """
+        Exit the current scope and return to the enclosing scope.
+
+        Exit the current scope and set the scope back to the value
+        of the enclosing scope.
+        """
+        self._scope = self._scope.parent_scope
+
+    def wrap(self, expr):
+        """
+        Get the wrapped version of the AST object.
+
+        Return the AST object which allows the object `expr` printed
+        in the start language to be accessed from the target language.
+
+        Parameters
+        ----------
+        expr : pyccel.ast.basic.Basic
+            The expression that should be wrapped.
+
+        Returns
+        -------
+        pyccel.ast.basic.Basic
+            The AST which describes the object that lets you
+            access the expression.
+        """
+        return self._wrap(expr)
+
+    def _wrap(self, expr):
+        """
+        Get the wrapped version of the AST object.
+
+        Private function returning the AST object which is used to access
+        the object `expr` from the target language.
+
+        Parameters
+        ----------
+        expr : pyccel.ast.basic.Basic
+            The expression that should be wrapped.
+
+        Returns
+        -------
+        pyccel.ast.basic.Basic
+            The AST which describes the object that lets you
+            access the expression.
+        """
+        classes = type(expr).mro()
+        for cls in classes:
+            wrap_method = '_wrap_' + cls.__name__
+            if hasattr(self, wrap_method):
+                return getattr(self, wrap_method)(expr)
+
+        return self._wrap_not_supported(expr)
+
+    def _wrap_not_supported(self, expr):
+        """ Print an error message if the wrap function for the type
+        is not implemented """
+        msg = f'_wrap_{type(expr).__name__} is not yet implemented for wrapper : {type(self)}\n'
+        errors.report(msg+PYCCEL_RESTRICTION_TODO, symbol = expr,
+                severity='fatal')
+
diff --git a/pyccel/compilers/default_compilers.py b/pyccel/compilers/default_compilers.py
index c04566a22b..50ac871e49 100644
--- a/pyccel/compilers/default_compilers.py
+++ b/pyccel/compilers/default_compilers.py
@@ -14,7 +14,7 @@
               'mpi_exec' : 'mpif90',
               'language': 'fortran',
               'module_output_flag': '-J',
-              'debug_flags': ("-fcheck=bounds",),
+              'debug_flags': ("-fcheck=bounds","-g","-O0"),
               'release_flags': ("-O3","-funroll-loops",),
               'general_flags' : ('-fPIC',),
               'standard_flags' : ('-std=f2003',),
@@ -41,7 +41,7 @@
               'mpi_exec' : 'mpiifort',
               'language': 'fortran',
               'module_output_flag': '-module',
-              'debug_flags': ("-check=bounds",),
+              'debug_flags': ("-check=bounds","-g","-O0"),
               'release_flags': ("-O3","-funroll-loops",),
               'general_flags' : ('-fPIC',),
               'standard_flags' : ('-std=f2003',),
@@ -60,7 +60,7 @@
               'mpi_exec' : 'pgfortran',
               'language': 'fortran',
               'module_output_flag': '-module',
-              'debug_flags': ("-Mbounds",),
+              'debug_flags': ("-Mbounds","-g","-O0"),
               'release_flags': ("-O3","-Munroll",),
               'general_flags' : ('-fPIC',),
               'standard_flags' : ('-Mstandard',),
@@ -78,7 +78,7 @@
               'mpi_exec' : 'nvfort',
               'language': 'fortran',
               'module_output_flag': '-module',
-              'debug_flags': ("-Mbounds",),
+              'debug_flags': ("-Mbounds","-g","-O0"),
               'release_flags': ("-O3","-Munroll",),
               'general_flags' : ('-fPIC',),
               'standard_flags' : ('-Mstandard',),
@@ -95,7 +95,7 @@
 gcc_info = {'exec' : 'gcc',
             'mpi_exec' : 'mpicc',
             'language': 'c',
-            'debug_flags': ("-g",),
+            'debug_flags': ("-g","-O0"),
             'release_flags': ("-O3","-funroll-loops",),
             'general_flags' : ('-fPIC',),
             'standard_flags' : ('-std=c99',),
@@ -126,7 +126,7 @@
 icc_info = {'exec' : 'icc',
             'mpi_exec' : 'mpiicc',
             'language': 'c',
-            'debug_flags': ("-g",),
+            'debug_flags': ("-g","-O0"),
             'release_flags': ("-O3","-funroll-loops",),
             'general_flags' : ('-fPIC',),
             'standard_flags' : ('-std=c99',),
@@ -143,7 +143,7 @@
 pgcc_info = {'exec' : 'pgcc',
             'mpi_exec' : 'pgcc',
             'language': 'c',
-            'debug_flags': ("-g",),
+            'debug_flags': ("-g","-O0"),
             'release_flags': ("-O3","-Munroll",),
             'general_flags' : ('-fPIC',),
             'standard_flags' : ('-std=c99',),
@@ -160,7 +160,7 @@
 nvc_info = {'exec' : 'nvc',
             'mpi_exec' : 'nvc',
             'language': 'c',
-            'debug_flags': ("-g",),
+            'debug_flags': ("-g","-O0"),
             'release_flags': ("-O3","-Munroll",),
             'general_flags' : ('-fPIC',),
             'standard_flags' : ('-std=c99',),
diff --git a/pyccel/naming/cnameclashchecker.py b/pyccel/naming/cnameclashchecker.py
index 7bc31965bb..19710c7ff7 100644
--- a/pyccel/naming/cnameclashchecker.py
+++ b/pyccel/naming/cnameclashchecker.py
@@ -6,11 +6,15 @@
 """
 Handles name clash problems in C
 """
-from pyccel.utilities.metaclasses import Singleton
-from pyccel.utilities.strings import create_incremented_string
+from .languagenameclashchecker import LanguageNameClashChecker
 
-class CNameClashChecker(metaclass = Singleton):
-    """ Class containing functions to help avoid problematic names in C
+class CNameClashChecker(LanguageNameClashChecker):
+    """
+    Class containing functions to help avoid problematic names in C.
+
+    A class which provides functionalities to check or propose variable names and
+    verify that they do not cause name clashes. Name clashes may be due to
+    new variables, or due to the use of reserved keywords.
     """
     # Keywords as mentioned on https://en.cppreference.com/w/c/keyword
     keywords = set(['isign', 'fsign', 'csign', 'auto', 'break', 'case', 'char', 'const',
@@ -42,20 +46,29 @@ def has_clash(self, name, symbols):
                any(name == s for s in symbols)
 
     def get_collisionless_name(self, name, symbols):
-        """ Get the name that will be used in the fortran code
+        """
+        Get a valid name which doesn't collision with symbols or C keywords.
+
+        Find a new name based on the suggested name which will not cause
+        conflicts with C keywords, does not appear in the provided symbols,
+        and is a valid name in C code.
+
+        Parameters
+        ----------
+        name : str
+            The suggested name.
+        symbols : set
+            Symbols which should be considered as collisions.
+
+        Returns
+        -------
+        str
+            A new name which is collision free.
         """
         if len(name)>4 and all(name[i] == '_' for i in (0,1,-1,-2)):
             # Ignore magic methods
             return name
         if name[0] == '_':
             name = 'private'+name
-        prefix = name
-        coll_symbols = self.keywords.copy()
-        coll_symbols.update(s.lower() for s in symbols)
-        if prefix in coll_symbols:
-            counter = 1
-            new_name, counter = create_incremented_string(coll_symbols,
-                    prefix = prefix, counter = counter)
-            name = name+new_name[-5:]
-        return name
+        return self._get_collisionless_name(name, symbols)
 
diff --git a/pyccel/naming/fortrannameclashchecker.py b/pyccel/naming/fortrannameclashchecker.py
index f6c041cd06..62a518e5d0 100644
--- a/pyccel/naming/fortrannameclashchecker.py
+++ b/pyccel/naming/fortrannameclashchecker.py
@@ -4,15 +4,20 @@
 # go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
 #------------------------------------------------------------------------------------------#
 """
-Handles name clash problems in Fortran
+Handles name clash problems in Fortran.
 """
 import warnings
 
-from pyccel.utilities.metaclasses import Singleton
-from pyccel.utilities.strings import create_incremented_string
+from .languagenameclashchecker import LanguageNameClashChecker
 
-class FortranNameClashChecker(metaclass = Singleton):
-    """ Class containing functions to help avoid problematic names in C
+class FortranNameClashChecker(LanguageNameClashChecker):
+    """
+    Class containing functions to help avoid problematic names in Fortran.
+
+    A class which provides functionalities to check or propose variable names and
+    verify that they do not cause name clashes. Name clashes may be due to
+    capitalisation (as Fortran is not case-sensitive), or due to the use of reserved
+    keywords.
     """
     # Keywords as mentioned on https://fortranwiki.org/fortran/show/Keywords
     # Intrinsic functions as mentioned on https://pages.mtu.edu/~shene/COURSES/cs201/NOTES/chap02/funct.html
@@ -40,7 +45,7 @@ class FortranNameClashChecker(metaclass = Singleton):
             'unlock', 'test', 'abs', 'sqrt', 'sin', 'cos', 'tan',
             'asin', 'acos', 'atan', 'exp', 'log', 'int', 'nint',
             'floor', 'fraction', 'real', 'max', 'mod', 'count',
-            'pack', 'numpy_sign'])
+            'pack', 'numpy_sign', 'c_associated', 'c_loc', 'c_f_pointer', 'c_ptr'])
 
     def has_clash(self, name, symbols):
         """ Indicate whether the proposed name causes any clashes
@@ -50,21 +55,31 @@ def has_clash(self, name, symbols):
                any(name == s.lower() for s in symbols)
 
     def get_collisionless_name(self, name, symbols):
-        """ Get the name that will be used in the fortran code
+        """
+        Get a valid name which doesn't collide with symbols or Fortran keywords.
+
+        Find a new name based on the suggested name which will not cause
+        conflicts with Fortran keywords, does not conflict with the provided symbols,
+        and is a valid name in Fortran code.
+
+        Parameters
+        ----------
+        name : str
+            The suggested name.
+        symbols : set
+            Symbols which should be considered as collisions.
+
+        Returns
+        -------
+        str
+            A new name which is collision free.
         """
         if len(name)>4 and all(name[i] == '_' for i in (0,1,-1,-2)):
             # Ignore magic methods
             return name
         if name[0] == '_':
             name = 'private'+name
-        prefix = name.lower()
-        coll_symbols = self.keywords.copy()
-        coll_symbols.update(s.lower() for s in symbols)
-        if prefix in coll_symbols:
-            counter = 1
-            new_name, counter = create_incremented_string(coll_symbols,
-                    prefix = prefix, counter = counter)
-            name = name+new_name[-5:]
+        name = self._get_collisionless_name(name, symbols)
         if len(name) > 96:
             warnings.warn("Name {} is too long for Fortran. This may cause compiler errors".format(name))
         return name
diff --git a/pyccel/naming/languagenameclashchecker.py b/pyccel/naming/languagenameclashchecker.py
new file mode 100644
index 0000000000..a606db003d
--- /dev/null
+++ b/pyccel/naming/languagenameclashchecker.py
@@ -0,0 +1,47 @@
+# coding: utf-8
+#------------------------------------------------------------------------------------------#
+# This file is part of Pyccel which is released under MIT License. See the LICENSE file or #
+# go to https://github.com/pyccel/pyccel/blob/master/LICENSE for full license details.     #
+#------------------------------------------------------------------------------------------#
+"""
+Superclass for handling name clash problems.
+"""
+from pyccel.utilities.metaclasses import Singleton
+from pyccel.utilities.strings import create_incremented_string
+
+class LanguageNameClashChecker(metaclass = Singleton):
+    """
+    Class containing functions to help avoid problematic names in a target language.
+
+    A super class which provides functionalities to check or propose variable names and
+    verify that they do not cause name clashes. Name clashes may be due to
+    a variety of reasons which vary from language to language.
+    """
+    keywords = None
+
+    def _get_collisionless_name(self, name, symbols):
+        """
+        Get a name which doesn't collision with keywords or symbols.
+
+        Find a new name based on the suggested name which does not collision
+        with the language keywords or the provided symbols.
+
+        Parameters
+        ----------
+        name : str
+            The suggested name.
+        symbols : set
+            Symbols which should be considered as collisions.
+
+        Returns
+        -------
+        str
+            A new name which is collision free.
+        """
+        coll_symbols = self.keywords.copy()
+        coll_symbols.update(s.lower() for s in symbols)
+        if self.has_clash(name, coll_symbols): #pylint: disable=no-member
+            counter = 1
+            name, counter = create_incremented_string(coll_symbols,
+                    prefix = name, counter = counter, name_clash_checker = self)
+        return name
diff --git a/pyccel/naming/pythonnameclashchecker.py b/pyccel/naming/pythonnameclashchecker.py
index cf97b8aa90..4c44189bed 100644
--- a/pyccel/naming/pythonnameclashchecker.py
+++ b/pyccel/naming/pythonnameclashchecker.py
@@ -6,11 +6,15 @@
 """
 Handles name clash problems in Python
 """
-from pyccel.utilities.metaclasses import Singleton
-from pyccel.utilities.strings import create_incremented_string
+from .languagenameclashchecker import LanguageNameClashChecker
 
-class PythonNameClashChecker(metaclass = Singleton):
-    """ Class containing functions to help avoid problematic names in Python
+class PythonNameClashChecker(LanguageNameClashChecker):
+    """
+    Class containing functions to help avoid problematic names in Python.
+
+    A class which provides functionalities to check or propose variable names and
+    verify that they do not cause name clashes. Name clashes may arise when
+    generating names for new variables.
     """
     keywords = set()
 
@@ -21,15 +25,25 @@ def has_clash(self, name, symbols):
                any(name == s for s in symbols)
 
     def get_collisionless_name(self, name, symbols):
-        """ Get the name that will be used in the fortran code
         """
-        prefix = name
-        coll_symbols = self.keywords.copy()
-        coll_symbols.update(s.lower() for s in symbols)
-        if prefix in coll_symbols:
-            counter = 1
-            new_name, counter = create_incremented_string(coll_symbols,
-                    prefix = prefix, counter = counter)
-            name = name+new_name[-5:]
-        return name
+        Get a valid name which doesn't collision with symbols.
+
+        Find a new name based on the suggested name which does not
+        appear in the provided symbols. It is not necessary to exclude
+        keywords for names which were either originally valid Python
+        names, or internally generated names.
+
+        Parameters
+        ----------
+        name : str
+            The suggested name.
+        symbols : set
+            Symbols which should be considered as collisions.
+
+        Returns
+        -------
+        str
+            A new name which is collision free.
+        """
+        return self._get_collisionless_name(name, symbols)
 
diff --git a/pyccel/parser/scope.py b/pyccel/parser/scope.py
index 669ce59b96..74d79df50d 100644
--- a/pyccel/parser/scope.py
+++ b/pyccel/parser/scope.py
@@ -22,14 +22,33 @@
 errors = Errors()
 
 class Scope(object):
-    """ Class representing all objects defined within a given
-    scope
+    """
+    Class representing all objects defined within a given scope.
+
+    This class provides all necessary functionalities for creating new object
+    names without causing name clashes. It also stores all objects defined
+    within the scope. This allows us to search for variables only in relevant
+    scopes.
 
     Parameters
     ----------
-    decorators : dict
-                 A dictionary of any decorators which operate on
-                 objects in this scope
+    decorators : dict, default: ()
+        A dictionary of any decorators which operate on objects in this scope.
+
+    is_loop : bool, default: False
+        Indicates if the scope represents a loop (in Python variables declared
+        in loops are not scoped to the loop).
+
+    parent_scope : Scope, default: None
+        The enclosing scope.
+
+    used_symbols : set, default: None
+        A set of all the names which we know will appear in the scope and which
+        we therefore want to avoid when creating new names.
+
+    original_symbols : dict, default: None
+        A dictionary which maps names used in the code to the original name used
+        in the Python code.
     """
     allow_loop_scoping = False
     name_clash_checker = PythonNameClashChecker()
@@ -42,7 +61,7 @@ class Scope(object):
             'macros','templates','headers','decorators',
             'cls_constructs')
 
-    def __init__(self, *, decorators=None, is_loop = False,
+    def __init__(self, *, decorators = (), is_loop = False,
                     parent_scope = None, used_symbols = None,
                     original_symbols = None):
 
@@ -60,8 +79,7 @@ def __init__(self, *, decorators=None, is_loop = False,
 
         self._dummy_counter = 0
 
-        if decorators:
-            self._locals['decorators'].update(decorators)
+        self._locals['decorators'].update(decorators)
 
         # TODO use another name for headers
         #      => reserved keyword, or use __
@@ -407,18 +425,27 @@ def local_used_symbols(self):
 
     def get_new_incremented_symbol(self, prefix, counter):
         """
-        Creates a new name by adding a numbered suffix to the provided prefix.
+        Create a new name by adding a numbered suffix to the provided prefix.
+
+        Create a new name which does not clash with any existing names by
+        adding a numbered suffix to the provided prefix.
+
+        Parameters
+        ----------
+        prefix : str
+            The prefix from which the new name will be created.
 
-          Parameters
-          ----------
-          prefix : str
+        counter : int
+            The starting point for the incrementation.
 
-          Returns
-          -------
-          new_name     : str
+        Returns
+        -------
+        pyccel.ast.internals.PyccelSymbol
+            The newly created name.
         """
 
-        new_name, counter = create_incremented_string(self.local_used_symbols.values(), prefix = prefix)
+        new_name, counter = create_incremented_string(self.local_used_symbols.values(),
+                                    prefix = prefix, counter = counter, name_clash_checker = self.name_clash_checker)
 
         new_symbol = PyccelSymbol(new_name, is_temp=True)
 
@@ -428,20 +455,24 @@ def get_new_incremented_symbol(self, prefix, counter):
 
     def get_new_name(self, current_name = None):
         """
+        Get a new name which does not clash with any names in the current context.
+
         Creates a new name. A current_name can be provided indicating the name the
         user would like to use if possible. If this name is not available then it
         will be used as a prefix for the new name.
         If no current_name is provided, then the standard prefix is used, and the
         dummy counter is used and updated to facilitate finding the next value of
-        this common case
+        this common case.
 
-          Parameters
-          ----------
-          current_name : str
+        Parameters
+        ----------
+        current_name : str, default: None
+            The name the user would like to use if possible.
 
-          Returns
-          -------
-          new_name     : PyccelSymbol
+        Returns
+        -------
+        PyccelSymbol
+            The new name which will be printed in the code.
         """
         if current_name is not None and not self.name_clash_checker.has_clash(current_name, self.all_used_symbols):
             new_name = PyccelSymbol(current_name)
@@ -452,13 +483,13 @@ def get_new_name(self, current_name = None):
             # Avoid confusing names by also searching in parent scopes
             new_name, self._dummy_counter = create_incremented_string(self.all_used_symbols,
                                                 prefix = current_name,
-                                                counter = self._dummy_counter)
+                                                counter = self._dummy_counter,
+                                                name_clash_checker = self.name_clash_checker)
         else:
             # When a name is suggested, try to stick to it
             new_name,_ = create_incremented_string(self.all_used_symbols, prefix = current_name)
 
         new_name = PyccelSymbol(new_name, is_temp = True)
-
         self.insert_symbol(new_name)
 
         return new_name
diff --git a/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.c b/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.c
index 112bffbc07..9d998ba6a1 100644
--- a/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.c
+++ b/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.c
@@ -463,7 +463,6 @@ int     nd_ndim(t_ndarray *a, int n)
  *
  * 	Parameters	:
  *		a 	  : python array object
- *      index : dimension index
  * 	Returns		:
  *		return NULL if object is NULL or the data of the array
  * reference of the used c/numpy api function
@@ -477,3 +476,48 @@ void    *nd_data(t_ndarray *a)
 
 	return a->raw_data;
 }
+
+/*
+ * Function: nd_step_C
+ * --------------------
+ * Return the step in the nth dimension for a C-ordered array (decreasing strides)
+ *
+ * 	Parameters	:
+ *		a 	  : python array object
+ * 		index	  : dimension index
+ * 	Returns		:
+ *		return 1 if object is NULL or the step along the indexed dimension
+ */
+int     nd_nstep_C(t_ndarray *a, int n)
+{
+	if (a == NULL || a->length == 0)
+		return 1;
+
+	int step = a->strides[n];
+	for (int i = n+1; i<a->nd; ++i) {
+		step /= a->shape[i];
+	}
+	return step > 0 ? step : 1;
+}
+/*
+ * Function: nd_step_F
+ * --------------------
+ * Return the step in the nth dimension for a F-ordered array (increasing strides)
+ *
+ * 	Parameters	:
+ *		a 	  : python array object
+ * 		index	  : dimension index
+ * 	Returns		:
+ *		return 1 if object is NULL or the step along the indexed dimension
+ */
+int     nd_nstep_F(t_ndarray *a, int n)
+{
+	if (a == NULL || a->length == 0)
+		return 1;
+
+	int step = a->strides[n];
+	for (int i = 0; i<n; ++i) {
+		step /= a->shape[i];
+	}
+	return step > 0 ? step : 1;
+}
diff --git a/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.h b/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.h
index 3404a3301a..c92c005a0c 100644
--- a/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.h
+++ b/pyccel/stdlib/cwrapper_ndarrays/cwrapper_ndarrays.h
@@ -49,5 +49,7 @@ bool	is_numpy_array(PyObject *o, int dtype, int rank, int flag);
 
 void    *nd_data(t_ndarray *a);
 int     nd_ndim(t_ndarray *a, int n);
+int     nd_nstep_C(t_ndarray *a, int n);
+int     nd_nstep_F(t_ndarray *a, int n);
 
 #endif
diff --git a/pyccel/utilities/metaclasses.py b/pyccel/utilities/metaclasses.py
index 7158429747..c23b631d6e 100644
--- a/pyccel/utilities/metaclasses.py
+++ b/pyccel/utilities/metaclasses.py
@@ -12,12 +12,17 @@
 )
 
 class Singleton(type):
-    """ Indicates that there is only one instance of the class
+    """
+    Metaclass indicating that there is only one instance of the class.
+
+    A metaclass which ensures that only one instance of the class is ever
+    created. Trying to create a second instance will result in accessing
+    the first.
     """
     _instances = {}
-    def __call__(cls, *args, **kwargs):
+    def __call__(cls):
         if cls not in cls._instances:
-            cls._instances[cls] = super().__call__(*args, **kwargs)
+            cls._instances[cls] = super().__call__()
         return cls._instances[cls]
 
 class ArgumentSingleton(type):
diff --git a/pyccel/utilities/strings.py b/pyccel/utilities/strings.py
index 5d9964756f..13deff4ae3 100644
--- a/pyccel/utilities/strings.py
+++ b/pyccel/utilities/strings.py
@@ -6,31 +6,36 @@
 """ Module containing helper functions for managing strings
 """
 
-def create_incremented_string(forbidden_exprs, prefix = 'Dummy', counter = 1):
-    """This function takes a prefix and a counter and uses them to construct
+def create_incremented_string(forbidden_exprs, prefix = 'Dummy', counter = 1, name_clash_checker = None):
+    """
+    Create a new unique string by incrementing a prefix.
+
+    This function takes a prefix and a counter and uses them to construct
     a new name of the form:
             prefix_counter
     Where counter is formatted to fill 4 characters
     The new name is checked against a list of forbidden expressions. If the
     constructed name is forbidden then the counter is incremented until a valid
-    name is found
-
-      Parameters
-      ----------
-      forbidden_exprs : Set
-                        A set of all the values which are not valid solutions to this problem
-      prefix          : str
-                        The prefix used to begin the string
-      counter         : int
-                        The expected value of the next name
+    name is found.
 
-      Returns
-      ----------
-      name            : str
-                        The incremented string name
-      counter         : int
-                        The expected value of the next name
+    Parameters
+    ----------
+    forbidden_exprs : set
+        A set of all the values which are not valid solutions to this problem.
+    prefix : str
+        The prefix used to begin the string.
+    counter : int
+        The expected value of the next name.
+    name_clash_checker : pyccel.naming.languagenameclashchecker.LanguageNameClashChecker
+        A class instance providing access to a `has_clash` function which determines
+        if names clash in a given language.
 
+    Returns
+    -------
+    name : str
+        The incremented string name.
+    counter : int
+        The expected value of the next name.
     """
     nDigits = 4
 
@@ -40,8 +45,13 @@ def create_incremented_string(forbidden_exprs, prefix = 'Dummy', counter = 1):
     name_format = "{prefix}_{counter:0="+str(nDigits)+"d}"
     name = name_format.format(prefix=prefix, counter = counter)
     counter += 1
-    while name in forbidden_exprs:
-        name = name_format.format(prefix=prefix, counter = counter)
-        counter += 1
+    if name_clash_checker:
+        while name_clash_checker.has_clash(name, forbidden_exprs):
+            name = name_format.format(prefix=prefix, counter = counter)
+            counter += 1
+    else:
+        while name in forbidden_exprs:
+            name = name_format.format(prefix=prefix, counter = counter)
+            counter += 1
 
     return name, counter


Test patch:
diff --git a/tests/parser/test_himi.py b/obsolete/tests/test_himi.py
similarity index 100%
rename from tests/parser/test_himi.py
rename to obsolete/tests/test_himi.py
diff --git a/tests/epyccel/test_arrays.py b/tests/epyccel/test_arrays.py
index ad417f0977..519a5040ad 100644
--- a/tests/epyccel/test_arrays.py
+++ b/tests/epyccel/test_arrays.py
@@ -85,6 +85,20 @@ def test_array_int32_1d_scalar_add(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_1d_scalar_add_stride(language):
+
+    f1 = arrays.array_int32_1d_scalar_add
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [1,2,3,4,5,6,7,8], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[::3], a)
+    f2(x2[::3], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_1d_scalar_sub(language):
 
     f1 = arrays.array_int32_1d_scalar_sub
@@ -99,6 +113,20 @@ def test_array_int32_1d_scalar_sub(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_1d_scalar_sub_stride(language):
+
+    f1 = arrays.array_int32_1d_scalar_sub
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [1,2,3,4,5,6,7,8,9], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[::2], a)
+    f2(x2[::2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_1d_scalar_mul(language):
 
     f1 = arrays.array_int32_1d_scalar_mul
@@ -113,6 +141,20 @@ def test_array_int32_1d_scalar_mul(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_1d_scalar_mul_stride(language):
+
+    f1 = arrays.array_int32_1d_scalar_mul
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [1,2,3,4,5,6,7,8,9], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[3:7:2], a)
+    f2(x2[3:7:2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_1d_scalar_div(language):
 
     f1 = arrays.array_int32_1d_scalar_div
@@ -141,6 +183,20 @@ def test_array_int32_1d_scalar_idiv(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_1d_scalar_idiv_stride(language):
+
+    f1 = arrays.array_int32_1d_scalar_idiv
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [1,2,3,4,5,6,7,8,9], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = 1, high = 1e9, dtype = np.int32)
+
+    f1(x1[:3:2], a)
+    f2(x2[:3:2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_1d_add(language):
 
     f1 = arrays.array_int32_1d_add
@@ -264,6 +320,20 @@ def test_array_int32_2d_C_scalar_add(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_2d_C_scalar_add_stride(language):
+
+    f1 = arrays.array_int32_2d_C_scalar_add
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3], [4,5,6]], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[::2,:], a)
+    f2(x2[::2,:], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_C_scalar_sub(language):
 
     f1 = arrays.array_int32_2d_C_scalar_sub
@@ -278,6 +348,30 @@ def test_array_int32_2d_C_scalar_sub(language):
 
     assert np.array_equal( x1, x2 )
 
+@pytest.mark.parametrize( 'language', [
+        pytest.param("c", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.c]),
+        pytest.param("fortran", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.fortran]),
+        pytest.param("python", marks = pytest.mark.python)
+    ]
+)
+def test_array_int32_2d_C_scalar_sub_stride(language):
+
+    f1 = arrays.array_int32_2d_C_scalar_sub
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3,7], [4,5,6,8]], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[:,::2], a)
+    f2(x2[:,::2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_C_scalar_mul(language):
 
     f1 = arrays.array_int32_2d_C_scalar_mul
@@ -292,6 +386,20 @@ def test_array_int32_2d_C_scalar_mul(language):
 
     assert np.array_equal( x1, x2 )
 
+def test_array_int32_2d_C_scalar_mul_stride(language):
+
+    f1 = arrays.array_int32_2d_C_scalar_mul
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3], [4,5,6]], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[1:,:], a)
+    f2(x2[1:,:], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_C_scalar_idiv(language):
 
     f1 = arrays.array_int32_2d_C_scalar_idiv
@@ -306,6 +414,30 @@ def test_array_int32_2d_C_scalar_idiv(language):
 
     assert np.array_equal( x1, x2 )
 
+@pytest.mark.parametrize( 'language', [
+        pytest.param("c", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.c]),
+        pytest.param("fortran", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.fortran]),
+        pytest.param("python", marks = pytest.mark.python)
+    ]
+)
+def test_array_int32_2d_C_scalar_idiv_stride(language):
+
+    f1 = arrays.array_int32_2d_C_scalar_idiv
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3], [4,5,6]], dtype=np.int32 )
+    x2 = np.copy(x1)
+    a = randint(low = 1, high = 1e9, dtype = np.int32)
+
+    f1(x1[:,1:], a)
+    f2(x2[:,1:], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_C_add(language):
 
     f1 = arrays.array_int32_2d_C_add
@@ -380,6 +512,30 @@ def test_array_int32_2d_F_scalar_add(language):
 
     assert np.array_equal( x1, x2 )
 
+@pytest.mark.parametrize( 'language', [
+        pytest.param("c", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.c]),
+        pytest.param("fortran", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.fortran]),
+        pytest.param("python", marks = pytest.mark.python)
+    ]
+)
+def test_array_int32_2d_F_scalar_add_stride(language):
+
+    f1 = arrays.array_int32_2d_F_scalar_add
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3], [4,5,6], [7,8,9]], dtype=np.int32, order='F' )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[1::2,::2], a)
+    f2(x2[1::2,::2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_F_scalar_sub(language):
 
     f1 = arrays.array_int32_2d_F_scalar_sub
@@ -394,6 +550,30 @@ def test_array_int32_2d_F_scalar_sub(language):
 
     assert np.array_equal( x1, x2 )
 
+@pytest.mark.parametrize( 'language', [
+        pytest.param("c", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.c]),
+        pytest.param("fortran", marks = [
+            pytest.mark.xfail(reason="Ordering is unknown on non-contiguous array"),
+            pytest.mark.fortran]),
+        pytest.param("python", marks = pytest.mark.python)
+    ]
+)
+def test_array_int32_2d_F_scalar_sub_stride(language):
+
+    f1 = arrays.array_int32_2d_F_scalar_sub
+    f2 = epyccel( f1 , language = language)
+
+    x1 = np.array( [[1,2,3], [4,5,6], [7,8,9]], dtype=np.int32, order='F' )
+    x2 = np.copy(x1)
+    a = randint(low = -1e9, high = 1e9, dtype = np.int32)
+
+    f1(x1[::2,1::2], a)
+    f2(x2[::2,1::2], a)
+
+    assert np.array_equal( x1, x2 )
+
 def test_array_int32_2d_F_scalar_mul(language):
 
     f1 = arrays.array_int32_2d_F_scalar_mul