diff --git a/CMakeLists.txt b/CMakeLists.txt
index aec33e15b..995cd9e95 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -35,6 +35,7 @@ ecbuild_add_option( FEATURE SINGLE_PRECISION
                     DEFAULT ON
                     DESCRIPTION "Support for Single Precision" )
 
+
 if( HAVE_SINGLE_PRECISION )
   set( single "single" )
 endif()
@@ -64,6 +65,13 @@ ecbuild_add_option( FEATURE ETRANS
                     DEFAULT OFF
                     DESCRIPTION "Include Limited-Area-Model Transforms" )
 
+	   
+ecbuild_add_option( FEATURE ECTRANS4PY
+                    DEFAULT OFF
+		    CONDITION HAVE_ETRANS
+                    DESCRIPTION "Compile ectrans4py interface routines for python binding w/ ctypesForFortran" )
+
+
 ectrans_find_lapack()
 
 ### Add sources and tests
diff --git a/MANIFEST.in b/MANIFEST.in
new file mode 100644
index 000000000..1e1af65a5
--- /dev/null
+++ b/MANIFEST.in
@@ -0,0 +1,8 @@
+include AUTHORS
+include CMakeLists.txt
+include LICENSE
+include README.md
+include VERSION
+recursive-include src *
+recursive-include cmake *
+recursive-include tests *
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 000000000..a003b0c80
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,4 @@
+[build-system]
+requires = ["setuptools", "wheel", "scikit-build"]
+build-backend = "setuptools.build_meta"
+
diff --git a/setup.py b/setup.py
new file mode 100644
index 000000000..80a8788eb
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,32 @@
+import os
+import re
+import ast
+from skbuild import setup
+
+def get_version():   # remove this part 
+    version_file = os.path.join("src", "ectrans4py", "__init__.py")
+    with open(version_file, "r", encoding="utf-8") as f:
+        content = f.read()
+        version_match = re.search(r"^__version__\s*=\s*['\"]([^'\"]*)['\"]", content, re.M)
+        if version_match:
+            return version_match.group(1)
+        raise RuntimeError("Unable to find version string.")
+
+version=get_version()
+# ectrans4py package : 
+setup(
+    name="ectrans4py",
+    version=version, 
+    packages=['ectrans4py'],
+    cmake_minimum_required_version="3.13",
+    cmake_args=[
+        '-DENABLE_ETRANS=ON',  
+        '-DENABLE_ECTRANS4PY=ON',
+        '-DENABLE_SINGLE_PRECISION=OFF',
+        '-DENABLE_OMP=OFF',
+    ],
+    package_dir={"": "src"},
+    cmake_install_dir="src/ectrans4py",
+    setup_requires=["scikit-build", "setuptools"],
+    install_requires=["numpy", "ctypesforfortran==1.1.3"],
+)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 81883e1d2..706806cd8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -13,4 +13,7 @@ if( HAVE_TRANSI )
 endif()
 if( HAVE_ETRANS )
   add_subdirectory(etrans)
-endif()
\ No newline at end of file
+endif()
+if(HAVE_ECTRANS4PY)
+  add_subdirectory(ectrans4py)
+endif()
diff --git a/src/ectrans4py/CMakeLists.txt b/src/ectrans4py/CMakeLists.txt
new file mode 100644
index 000000000..857d7d609
--- /dev/null
+++ b/src/ectrans4py/CMakeLists.txt
@@ -0,0 +1,20 @@
+if(HAVE_ETRANS)
+    # (using CMAKE_CURRENT_SOURCE_DIR is necessary because sources are in a different directory than the target library (trans_${prec}))
+    ecbuild_list_add_pattern(
+        LIST ectrans4py_src
+        GLOB ${CMAKE_CURRENT_SOURCE_DIR}/*.F90
+        QUIET
+    )
+
+    set(HAVE_dp ${HAVE_DOUBLE_PRECISION})
+    set(prec dp)
+
+    if(HAVE_${prec})
+        # Add sources
+        target_sources(trans_${prec} PRIVATE ${ectrans4py_src})
+    endif()
+
+else()
+    ecbuild_critical("To activate the ectrans Python interface, you must enable the ETRANS option.")
+endif()
+
diff --git a/src/ectrans4py/__init__.py b/src/ectrans4py/__init__.py
new file mode 100644
index 000000000..5591f3465
--- /dev/null
+++ b/src/ectrans4py/__init__.py
@@ -0,0 +1,389 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Copyright (c) Météo France (2014-)
+# This software is governed by the CeCILL-C license under French law.
+# http://www.cecill.info
+"""
+ialsptrans4py:
+
+Contains the interface to spectral transforms from the IAL/ecTrans.
+Note that this is temporary between the former package arpifs4py and a direct python interface to ecTrans.
+
+Actual .so library should be in one of the preinstalled paths or in a directory specified via LD_LIBRARY_PATH
+"""
+
+from __future__ import print_function, absolute_import, unicode_literals, division
+
+import os
+import numpy as np
+import ctypesForFortran
+from ctypesForFortran import addReturnCode, treatReturnCode, IN, OUT
+
+
+
+__version__ = "2.0.1"
+
+# Shared objects library
+########################
+shared_objects_library = os.environ.get('IALSPTRANS4PY_SO', None)
+if shared_objects_library is None or not os.path.exists(shared_objects_library):
+    # not specified or path does not exist : find in known locations
+    so_basename = "libtrans_dp.so"  # local name in the directory
+    LD_LIBRARY_PATH = [p for p in os.environ.get('LD_LIBRARY_PATH', '').split(':') if p != '']
+    potential_locations = LD_LIBRARY_PATH + [
+            os.path.join(os.path.dirname(os.path.realpath(__file__)), 'lib'),   # FIXEME : but requiere changes in CMakeLists.txt
+        os.path.join(os.path.dirname(os.path.realpath(__file__)), 'lib64'),     # force one libdir directory name ! 
+#        "/home/common/epygram/public/EPyGrAM/libs4py",  # CNRM
+#        "/home/gmap/mrpe/mary/public/EPyGrAM/libs4py",  # belenos/taranis
+#        "/home/acrd/public/EPyGrAM/libs4py",  # ECMWF's Atos aa-ad
+        ]
+    for _libs4py_dir in potential_locations:
+        shared_objects_library = os.path.join(_libs4py_dir, so_basename)
+        if os.path.exists(shared_objects_library):
+            break
+        else:
+            shared_objects_library = None
+    if shared_objects_library is None:
+        msg = ' '.join(["'{}' was not found in any of potential locations: {}.",
+                        "You can specify a different location using env var LD_LIBRARY_PATH",
+                        "or specify a precise full path using env var IALSPTRANS4PY_SO."]).format(
+                                so_filename, str(potential_locations))
+        raise FileNotFoundError(msg)
+else:
+    so_basename = os.path.basename(shared_objects_library)
+ctypesFF, handle = ctypesForFortran.ctypesForFortranFactory(shared_objects_library)
+
+# Initialization
+################
+
+def init_env(omp_num_threads=None,
+             no_mpi=False):
+    """
+    Set adequate environment for the inner libraries.
+
+    :param int omp_num_threads: sets OMP_NUM_THREADS
+    :param bool no_mpi: environment variable DR_HOOK_NOT_MPI set to 1
+    """
+    # because arpifs library is compiled with MPI & openMP
+    if omp_num_threads is not None:
+        os.environ['OMP_NUM_THREADS'] = str(omp_num_threads)
+    if no_mpi:
+        os.environ['DR_HOOK_NOT_MPI'] = '1'
+
+# Transforms interfaces
+#######################
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def etrans_inq4py(KSIZEI, KSIZEJ,
+                 KPHYSICALSIZEI, KPHYSICALSIZEJ,
+                 KTRUNCX, KTRUNCY,
+                 KNUMMAXRESOL,
+                 PDELATX, PDELATY):
+    """
+    Simplified wrapper to ETRANS_INQ.
+
+    Args:\n
+    1,2) KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+    3,4) KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+    5,6) KTRUNCX, KTRUNCY: troncatures
+    7) KNUMMAXRESOL: maximum number of troncatures handled
+    8,9) PDELTAX, PDELTAY: resolution along x,y axis
+
+    Returns:\n
+    1) KGPTOT: number of gridpoints
+    2) KSPEC: number of spectral coefficients
+    """
+    return ([KSIZEI, KSIZEJ,
+             KPHYSICALSIZEI, KPHYSICALSIZEJ,
+             KTRUNCX, KTRUNCY,
+             KNUMMAXRESOL,
+             PDELATX, PDELATY],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.float64, None, IN),
+             (np.float64, None, IN),
+             (np.int64, None, OUT),
+             (np.int64, None, OUT)],
+            None)
+
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def trans_inq4py(KSIZEJ, KTRUNC, KSLOEN, KLOEN, KNUMMAXRESOL):
+    """
+    Simplified wrapper to TRANS_INQ.
+
+    Args:\n
+    1) KSIZEJ: number of latitudes in grid-point space
+    2) KTRUNC: troncature
+    3) KSLOEN: Size of KLOEN
+    4) KLOEN: number of points on each latitude row
+    5) KNUMMAXRESOL: maximum number of troncatures handled
+
+    Returns:\n
+    1) KGPTOT: number of gridpoints
+    2) KSPEC: number of spectral coefficients
+    3) KNMENG: cut-off zonal wavenumber
+    """
+    return ([KSIZEJ, KTRUNC, KSLOEN, KLOEN, KNUMMAXRESOL],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, (KSLOEN,), IN),
+             (np.int64, None, IN),
+             (np.int64, None, OUT),
+             (np.int64, None, OUT),
+             (np.int64, (KSLOEN,), OUT)],
+            None)
+
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def sp2gp_lam4py(KSIZEI, KSIZEJ,
+                   KPHYSICALSIZEI, KPHYSICALSIZEJ,
+                   KTRUNCX, KTRUNCY,
+                   KNUMMAXRESOL,
+                   KSIZE,
+                   LGRADIENT,
+                   LREORDER,
+                   PDELTAX, PDELTAY,
+                   PSPEC):
+    """
+    Transform spectral coefficients into grid-point values.
+
+    Args:\n
+    1,2) KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+    3,4) KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+    5,6) KTRUNCX, KTRUNCY: troncatures
+    7) KNUMMAXRESOL: maximum number of troncatures handled
+    8) KSIZE: size of PSPEC
+    9) LGRADIENT: gradient computation
+    10) LREORDER: reorder spectral coefficients or not
+    11,12) PDELTAX,PDELTAY: resolution along x,y axis
+    13) PSPEC: spectral coefficient array
+
+    Returns:\n
+    1) PGPT: grid-point field
+    2) PGPTM: N-S derivative if LGRADIENT
+    3) PGPTL: E-W derivative if LGRADIENT
+    """
+    return ([KSIZEI, KSIZEJ,
+             KPHYSICALSIZEI, KPHYSICALSIZEJ,
+             KTRUNCX, KTRUNCY,
+             KNUMMAXRESOL,
+             KSIZE,
+             LGRADIENT,
+             LREORDER,
+             PDELTAX, PDELTAY,
+             PSPEC],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (bool, None, IN),
+             (bool, None, IN),
+             (np.float64, None, IN),
+             (np.float64, None, IN),
+             (np.float64, (KSIZE,), IN),
+             (np.float64, (KSIZEI * KSIZEJ,), OUT),
+             (np.float64, (KSIZEI * KSIZEJ,), OUT),
+             (np.float64, (KSIZEI * KSIZEJ,), OUT)],
+            None)
+
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def gp2sp_lam4py(KSIZE,
+                   KSIZEI, KSIZEJ,
+                   KPHYSICALSIZEI, KPHYSICALSIZEJ,
+                   KTRUNCX, KTRUNCY,
+                   KNUMMAXRESOL,
+                   PDELTAX, PDELTAY,
+                   LREORDER,
+                   PGPT):
+    """
+    Transform grid point values into spectral coefficients.
+
+    Args:\n
+    1) KSIZE: size of spectral field
+    2,3) KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+    4,5) KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+    6,7) KTRUNCX, KTRUNCY: troncatures
+    8) KNUMMAXRESOL: maximum number of troncatures handled
+    9,10) PDELTAX, PDELTAY: resolution along x,y axis
+    11) LREORDER: reorder spectral coefficients or not
+    12) PGPT: grid-point field
+
+    Returns:\n
+    1) PSPEC: spectral coefficient array
+    """
+    return ([KSIZE,
+             KSIZEI, KSIZEJ,
+             KPHYSICALSIZEI, KPHYSICALSIZEJ,
+             KTRUNCX, KTRUNCY,
+             KNUMMAXRESOL,
+             PDELTAX, PDELTAY,
+             LREORDER,
+             PGPT],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.float64, None, IN),
+             (np.float64, None, IN),
+             (bool, None, IN),
+             (np.float64, (KSIZEI * KSIZEJ,), IN),
+             (np.float64, (KSIZE,), OUT)],
+            None)
+
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def sp2gp_gauss4py(KSIZEJ,
+                     KTRUNC,
+                     KNUMMAXRESOL,
+                     KGPTOT,
+                     KSLOEN,
+                     KLOEN,
+                     KSIZE,
+                     LGRADIENT,
+                     LREORDER,
+                     PSPEC):
+    """
+    Transform spectral coefficients into grid-point values.
+
+    Args:\n
+    1) KSIZEJ: Number of latitudes
+    2) KTRUNC: troncature
+    3) KNUMMAXRESOL: maximum number of troncatures handled
+    4) KGPTOT: number of grid-points
+    5) KSLOEN: Size of KLOEN
+    6) KLOEN:
+    7) KSIZE: Size of PSPEC
+    8) LGRADIENT: compute derivatives
+    9) LREORDER: reorder spectral coefficients or not
+    10) PSPEC: spectral coefficient array
+
+    Returns:\n
+    1) PGPT: grid-point field
+    2) PGPTM: N-S derivative if LGRADIENT
+    3) PGPTL: E-W derivative if LGRADIENT
+    """
+    return ([KSIZEJ,
+             KTRUNC,
+             KNUMMAXRESOL,
+             KGPTOT,
+             KSLOEN,
+             KLOEN,
+             KSIZE,
+             LGRADIENT,
+             LREORDER,
+             PSPEC],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, (KSLOEN,), IN),
+             (np.int64, None, IN),
+             (bool, None, IN),
+             (bool, None, IN),
+             (np.float64, (KSIZE,), IN),
+             (np.float64, (KGPTOT,), OUT),
+             (np.float64, (KGPTOT,), OUT),
+             (np.float64, (KGPTOT,), OUT)],
+            None)
+
+
+@treatReturnCode
+@ctypesFF()
+@addReturnCode
+def gp2sp_gauss4py(KSPEC,
+                     KSIZEJ,
+                     KTRUNC,
+                     KNUMMAXRESOL,
+                     KSLOEN,
+                     KLOEN,
+                     KSIZE,
+                     LREORDER,
+                     PGPT):
+    """
+    Transform grid-point values into spectral coefficients.
+
+    Args:\n
+    1) KSPEC: size of spectral coefficients array
+    2) KSIZEJ: Number of latitudes
+    3) KTRUNC: troncature
+    4) KNUMMAXRESOL: maximum number of troncatures handled
+    5) KSLOEN: Size of KLOEN
+    6) KLOEN
+    7) KSIZE: Size of PGPT
+    8) LREORDER: reorder spectral coefficients or not
+    9) PGPT: grid-point field
+
+    Returns:\n
+    1) PSPEC: spectral coefficient array
+    """
+    return ([KSPEC,
+             KSIZEJ,
+             KTRUNC,
+             KNUMMAXRESOL,
+             KSLOEN,
+             KLOEN,
+             KSIZE,
+             LREORDER,
+             PGPT],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.int64, (KSLOEN,), IN),
+             (np.int64, None, IN),
+             (bool, None, IN),
+             (np.float64, (KSIZE,), IN),
+             (np.float64, (KSPEC,), OUT)],
+            None)
+
+
+@ctypesFF()
+def sp2gp_fft1d4py(KSIZES, KTRUNC, PSPEC, KSIZEG):
+    """
+    Transform spectral coefficients into grid-point values,
+    for a 1D array (vertical section academic model)
+
+    Args:\n
+    1) KSIZES size of PSPEC
+    2) KTRUNC: troncature
+    3) PSPEC: spectral coefficient array
+    4) KSIZEG: size of grid-point field (with extension zone)
+
+    Returns:\n
+    1) PGPT: grid-point field
+    """
+    return ([KSIZES, KTRUNC, PSPEC, KSIZEG],
+            [(np.int64, None, IN),
+             (np.int64, None, IN),
+             (np.float64, (KSIZES,), IN),
+             (np.int64, None, IN),
+             (np.float64, (KSIZEG,), OUT)],
+            None)
diff --git a/src/ectrans4py/etrans_inq4py.F90 b/src/ectrans4py/etrans_inq4py.F90
new file mode 100644
index 000000000..7f2113fba
--- /dev/null
+++ b/src/ectrans4py/etrans_inq4py.F90
@@ -0,0 +1,66 @@
+SUBROUTINE ETRANS_INQ4PY(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
+                       &KTRUNCX, KTRUNCY, KNUMMAXRESOL, PDELTAX, PDELTAY, &
+                       &KGPTOT, KSPEC)
+! ** PURPOSE
+!    Simplified wrapper to ETRANS_INQ
+!
+! ** DUMMY ARGUMENTS
+!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+!    KTRUNCX, KTRUNCY: troncatures
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    PDELTAX: x resolution
+!    PDELTAY: y resolution
+!    KGPTOT: number of gridpoints
+!    KSPEC: number of spectral coefficients
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan., S. Riette: PDELTAX and PDELTAY added
+!
+! I. Dummy arguments declaration
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSIZEI, KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+REAL(KIND=8), INTENT(IN) :: PDELTAX
+REAL(KIND=8), INTENT(IN) :: PDELTAY
+INTEGER(KIND=8), INTENT(OUT) :: KGPTOT
+INTEGER(KIND=8), INTENT(OUT) :: KSPEC
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+LOGICAL :: LLSTOP
+INTEGER :: IIDENTRESOL
+INTEGER, DIMENSION(1) :: ILOEN
+INTEGER :: IGPTOT, ISPEC
+
+#include "etrans_inq.h"
+
+ISIZEI=KSIZEI
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=KPHYSICALSIZEI
+IPHYSICALSIZEJ=KPHYSICALSIZEJ
+ITRUNCX=KTRUNCX
+ITRUNCY=KTRUNCY
+INUMMAXRESOL=KNUMMAXRESOL
+
+! III. Setup
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
+                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
+IF (.NOT. LLSTOP) THEN
+  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
+  KGPTOT=IGPTOT
+  KSPEC=ISPEC
+ENDIF
+!
+END SUBROUTINE ETRANS_INQ4PY
diff --git a/src/ectrans4py/gp2sp_gauss4py.F90 b/src/ectrans4py/gp2sp_gauss4py.F90
new file mode 100644
index 000000000..76fff02c8
--- /dev/null
+++ b/src/ectrans4py/gp2sp_gauss4py.F90
@@ -0,0 +1,113 @@
+SUBROUTINE GP2SP_GAUSS4PY(KRETURNCODE, KSPEC, KSIZEJ, KTRUNC, KNUMMAXRESOL, KSLOEN, KLOEN, KSIZE, LREORDER, PGPT, PSPEC)
+! ** PURPOSE
+!    Transform spectral coefficients into grid-point values
+!
+! ** DUMMY ARGUMENTS
+!    KRETURNCODE: error code
+!    KSPEC: size of spectral coefficients array
+!    KSIZEJ: Number of latitudes
+!    KTRUNC: troncature
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    KSLOEN: Size ok KLOEN
+!    KLOEN
+!    KSIZE: Size of PGPT
+!    LREORDER: switch to reorder spectral coefficients or not
+!    PGPT: grid-point field
+!    PSPEC: spectral coefficient array
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan. 2016, S. Riette: w_spec_setup interface modified
+!    March, 2016, A.Mary: LREORDER
+!
+! I. Dummy arguments declaration
+USE PARKIND1, ONLY : JPRB
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSPEC
+INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNC
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+INTEGER(KIND=8), INTENT(IN) :: KSLOEN
+INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
+INTEGER(KIND=8), INTENT(IN) :: KSIZE
+LOGICAL, INTENT(IN) :: LREORDER
+REAL(KIND=8), DIMENSION(KSIZE), INTENT(IN) :: PGPT
+REAL(KIND=8), DIMENSION(KSPEC), INTENT(OUT) :: PSPEC
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+LOGICAL :: LLSTOP
+INTEGER :: IIDENTRESOL
+INTEGER :: JI, JM, JN
+INTEGER, DIMENSION(0:KTRUNC) :: NASM0
+REAL(KIND=JPRB), DIMENSION(1, SIZE(PGPT)) :: ZSPBUF !size over-evaluated
+REAL(KIND=JPRB), DIMENSION(SIZE(PGPT), 1, 1) :: ZGPBUF
+REAL(KIND=8) :: ZDELTAX, ZDELTAY
+
+#include "trans_inq.h"
+#include "dir_trans.h"
+KRETURNCODE=0
+ILOEN(:)=KLOEN(:)
+ISIZEI=0
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=0
+IPHYSICALSIZEJ=0
+ITRUNCX=KTRUNC
+ITRUNCY=0
+INUMMAXRESOL=KNUMMAXRESOL
+!
+! III. Setup
+ZDELTAX=0.
+ZDELTAY=0.
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
+                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
+!
+! IV. Transformation
+! IV.a Shape of coefficient array
+IF (.NOT. LLSTOP) THEN
+  JI=1
+  DO JN=0, KTRUNC
+    NASM0(JN)=JI
+    JI=JI+1+JN+(JN+1)
+  ENDDO
+ENDIF
+
+! IV.b Direct transform
+IF (.NOT. LLSTOP) THEN
+  ZGPBUF(:,1,1)=REAL(PGPT(:),KIND=JPRB)
+  CALL DIR_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
+ENDIF
+
+! IV.c Reordering
+IF (LREORDER) THEN
+  IF(.NOT. LLSTOP) THEN
+    PSPEC(:)=0.
+    JI=1
+    DO JM=0, KTRUNC
+      DO JN=JM, KTRUNC
+        PSPEC(NASM0(JN)+JM)=REAL(ZSPBUF(1,JI),KIND=8)
+        JI=JI+1
+        IF(JM/=0) THEN
+          PSPEC(NASM0(JN)-JM)=REAL(ZSPBUF(1,JI),KIND=8)
+        ENDIF
+        JI=JI+1
+      ENDDO
+    ENDDO
+    IF(JI-1/=KSPEC) THEN
+      PRINT*, "Internal error in GP2SP_GAUSS4PY (spectral reordering)"
+      KRETURNCODE=-999
+    ENDIF
+  ENDIF
+ELSE
+  PSPEC(1:KSPEC) = REAL(ZSPBUF(1,1:KSPEC),KIND=8)
+ENDIF
+
+END SUBROUTINE GP2SP_GAUSS4PY
diff --git a/src/ectrans4py/gp2sp_lam4py.f90 b/src/ectrans4py/gp2sp_lam4py.f90
new file mode 100644
index 000000000..036a4674b
--- /dev/null
+++ b/src/ectrans4py/gp2sp_lam4py.f90
@@ -0,0 +1,121 @@
+SUBROUTINE GP2SP_LAM4PY(KRETURNCODE, KSIZE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
+                         &KTRUNCX, KTRUNCY, KNUMMAXRESOL, PDELTAX, PDELTAY, LREORDER, PGPT, PSPEC)
+! ** PURPOSE
+!    Transform grid point values into spectral coefficients
+!
+! ** DUMMY ARGUMENTS
+!    KRETURNCODE: error code
+!    KSIZE: size of spectral field
+!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+!    KTRUNCX, KTRUNCY: troncatures
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    PDELTAX: x resolution
+!    PDELTAY: y resolution
+!    LREORDER: switch to reorder spectral coefficients or not
+!    PGPT: grid-point field
+!    PSPEC: spectral coefficient array
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan., S. Riette: PDELTAX and PDELTAY added
+!    March, 2016, A.Mary: LREORDER
+!
+! I. Dummy arguments declaration
+USE PARKIND1, ONLY : JPRB
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSIZE, KSIZEI, KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+REAL(KIND=8), INTENT(IN) :: PDELTAX
+REAL(KIND=8), INTENT(IN) :: PDELTAY
+LOGICAL, INTENT(IN) :: LREORDER
+REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(IN) :: PGPT
+REAL(KIND=8), DIMENSION(KSIZE), INTENT(OUT) :: PSPEC
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
+INTEGER :: IGPTOT, ISPEC
+INTEGER, DIMENSION(0:KTRUNCY) :: ISPECINI, ISPECEND
+REAL(KIND=JPRB), DIMENSION(1, KSIZEI*KSIZEJ) :: ZSPBUF !size over-evaluated
+REAL(KIND=JPRB), DIMENSION(KSIZEI*KSIZEJ, 1, 1) :: ZGPBUF
+INTEGER :: JI, JM, JN, IIDENTRESOL
+LOGICAL :: LLSTOP
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+INTEGER, DIMENSION(1) :: ILOEN
+
+#include "edir_trans.h"
+#include "etrans_inq.h"
+
+KRETURNCODE=0
+LLSTOP=.FALSE.
+ISIZEI=KSIZEI
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=KPHYSICALSIZEI
+IPHYSICALSIZEJ=KPHYSICALSIZEJ
+ITRUNCX=KTRUNCX
+ITRUNCY=KTRUNCY
+INUMMAXRESOL=KNUMMAXRESOL
+
+! III. Setup
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
+                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
+
+! IV. Transformation
+
+! IV.a Shape of coefficient array
+!IGPTOT is the total number of points in grid-point space
+!ISPEC is the number of spectral coefficients
+!IESM0(m) is the index of spectral coefficient (m,0) in model
+!ISPECINI(n) is the index of the first of the 4 spectral coefficient (0,n) in FA file
+!ISPECEND(n) is the index of the last of the last 4 spectral coefficients (:,n) in FA file
+IF (.NOT. LLSTOP) THEN
+  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
+  JI=1
+  DO JN=0, ITRUNCY
+    ISPECINI(JN)=(JI-1)*4+1
+    JI=JI+COUNT(IESM0(1:ITRUNCX)-IESM0(0:ITRUNCX-1)>JN*4)
+    IF (ISPEC-IESM0(ITRUNCX)>JN*4) JI=JI+1
+    ISPECEND(JN)=(JI-1)*4
+  ENDDO
+ENDIF
+
+! III.b transform
+IF (.NOT. LLSTOP) THEN
+  ZGPBUF(:,1,1)=REAL(PGPT(:),KIND=JPRB)
+  CALL EDIR_TRANS(PSPSCALAR=ZSPBUF(:,:), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
+ENDIF
+
+! III.c Reordering
+! reorder Aladin :  file ordering = coeffs per blocks of m, 4 reals per coeff
+!           Aladin array ordering = coeffs per blocks of n, 4 reals per coeff
+IF (LREORDER) THEN
+  IF (.NOT. LLSTOP) THEN
+    JI=1
+    PSPEC(:)=0.
+    DO JM=0,ITRUNCX*4+4,4
+      DO JN=0,ITRUNCY
+        IF (ISPECINI(JN)+JM+3<=ISPECEND(JN)) THEN
+          PSPEC(ISPECINI(JN)+JM:ISPECINI(JN)+JM+3) = REAL(ZSPBUF(1,JI:JI+3),KIND=8)
+          JI=JI+4
+        ENDIF
+      ENDDO
+    ENDDO
+    IF(JI/=ISPEC+1) THEN
+      PRINT*, "Internal error in GP2SP_LAM4PY (spectral reordering)"
+      KRETURNCODE=-999
+    ENDIF
+  ENDIF
+ELSE
+  PSPEC(1:KSIZE) = REAL(ZSPBUF(1,1:KSIZE),KIND=8)
+ENDIF
+
+END SUBROUTINE GP2SP_LAM4PY
diff --git a/src/ectrans4py/sp2gp_fft1d4py.F90 b/src/ectrans4py/sp2gp_fft1d4py.F90
new file mode 100644
index 000000000..060f14f4d
--- /dev/null
+++ b/src/ectrans4py/sp2gp_fft1d4py.F90
@@ -0,0 +1,114 @@
+SUBROUTINE SP2GP_FFT1D4PY(KSIZES, KTRUNC, PSPEC, KSIZEG, PGPT)
+! ** PURPOSE
+!    Transform spectral coefficients into grid-point values,
+!    for a 1D array (vertical section academic model)
+!
+! ** DUMMY ARGUMENTS
+!    KSIZES size of PSPEC
+!    KTRUNC: troncature
+!    PSPEC: spectral coefficient array
+!    KSIZEG: size of grid-point field (with extension zone)
+!    PGPT: grid-point field
+!
+! ** AUTHOR
+!    26 March 2015, A. Mary, from utilities/pinuts/module/fa_datas_mod.F90
+!
+! ** MODIFICATIONS
+!
+! I. Dummy arguments declaration
+IMPLICIT NONE
+
+INTEGER(KIND=8), INTENT(IN) :: KSIZES
+INTEGER(KIND=8), INTENT(IN) :: KTRUNC
+REAL(KIND=8), DIMENSION(KSIZES), INTENT(IN) :: PSPEC
+INTEGER(KIND=8), INTENT(IN) :: KSIZEG
+REAL(KIND=8), DIMENSION(KSIZEG), INTENT(OUT) :: PGPT
+
+INTEGER(KIND=8) :: NSEFRE, NFTM, NDGLSUR
+REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: SP2, TRIGSE
+INTEGER(KIND=8), DIMENSION(:), ALLOCATABLE  :: NFAXE
+INTEGER(KIND=8), PARAMETER :: NZERO=0
+
+NDGLSUR = KSIZEG+MOD(KSIZEG,2)+2
+NFTM    = 2*(KTRUNC+1)
+ALLOCATE(SP2(NDGLSUR*NFTM))
+SP2     = 0.0
+SP2     = CONVRT2FFT(PSPEC,NZERO,KTRUNC,NDGLSUR)
+ALLOCATE(NFAXE(1:10))
+ALLOCATE(TRIGSE(1:KSIZEG)) 
+CALL SET99(TRIGSE,NFAXE,KSIZEG)
+CALL FFT992(SP2(1:KSIZEG), TRIGSE, NFAXE, 1, NDGLSUR, KSIZEG, 1, 1)
+DEALLOCATE(TRIGSE)
+DEALLOCATE(NFAXE)
+PGPT(:) = SP2(1:KSIZEG)
+
+CONTAINS
+
+! from utilities/pinuts/module/fa_datas_mod.F90
+! and utilities/pinuts/module/array_lib_mod.F90
+
+FUNCTION CONVRT2FFT(IN,X,Y,N) RESULT(OU)
+REAL(KIND=8),DIMENSION(:),INTENT(IN)      :: IN
+INTEGER(KIND=8),INTENT(IN)                       :: X, Y, N
+REAL(KIND=8),DIMENSION(N*2*(X+1))         :: OU
+
+INTEGER(KIND=8),DIMENSION(2*(X+1),(N/2))         :: MINQ 
+INTEGER(KIND=8),DIMENSION((N/2),2*(X+1))         :: TMINQ 
+REAL(KIND=8),DIMENSION(2*(X+1),(N/2))     :: OMINQ, EMINQ
+REAL(KIND=8),DIMENSION((N/2),2*(X+1))     :: TOMINQ, TEMINQ   
+REAL(KIND=8),DIMENSION(N*(X+1))           :: OINI, EINI
+REAL(KIND=8), PARAMETER                   :: ZZERO=0.0
+
+CALL SPLIT_ODEV(IN,OINI,EINI)
+MINQ   = MASQ(X,Y,N)
+OMINQ  = UNPACK(OINI,MINQ == 1,ZZERO)
+TOMINQ = TRANSPOSE(OMINQ)
+EMINQ  = UNPACK(EINI,MINQ == 1,ZZERO)
+TEMINQ = TRANSPOSE(EMINQ)
+TMINQ  = 1
+OINI   = PACK(TOMINQ,TMINQ > 0)
+EINI   = PACK(TEMINQ,TMINQ > 0)
+OU     = MIX_ODEV(OINI,EINI)
+END FUNCTION CONVRT2FFT
+
+FUNCTION MASQ(X,Y,N) RESULT(T)
+INTEGER(KIND=8),INTENT(IN)                       :: X, Y, N
+INTEGER(KIND=8),DIMENSION(1:2*(X+1),1:(N/2))     :: T
+
+INTEGER(KIND=8)                                  :: I, J
+INTEGER(KIND=8),DIMENSION(0:X)                   :: KM
+INTEGER(KIND=8),DIMENSION(0:Y)                   :: KN
+CALL ELLIPS64(X,Y,KN,KM) 
+T = 0
+DO I=0,Y
+  DO J=0,2*KN(I)+1
+    T(J+1,I+1)=1
+  END DO
+END DO
+END FUNCTION MASQ
+
+FUNCTION MIX_ODEV(TO,TE) RESULT(T)
+REAL(KIND=8),DIMENSION(:),INTENT(IN)        :: TO,TE
+REAL(KIND=8),DIMENSION(SIZE(TO)+SIZE(TE))   :: T
+
+INTEGER(KIND=8) :: I
+
+DO I=1,(SIZE(TO)+SIZE(TE))/2
+  T((2*I)-1)=TE(I)
+  T(2*I)=TO(I)
+END DO
+END FUNCTION MIX_ODEV
+
+SUBROUTINE SPLIT_ODEV(T,TO,TE)
+REAL(KIND=8),DIMENSION(:),INTENT(IN)          :: T
+REAL(KIND=8),DIMENSION(SIZE(T)/2),INTENT(OUT) :: TO,TE
+
+INTEGER(KIND=8) :: I
+
+DO I=1,SIZE(T)/2
+  TO(I)=T(2*I)
+  TE(I)=T((2*I)-1)
+END DO
+END SUBROUTINE SPLIT_ODEV
+
+END SUBROUTINE SP2GP_FFT1D4PY
\ No newline at end of file
diff --git a/src/ectrans4py/sp2gp_gauss4py.F90 b/src/ectrans4py/sp2gp_gauss4py.F90
new file mode 100644
index 000000000..61186f53f
--- /dev/null
+++ b/src/ectrans4py/sp2gp_gauss4py.F90
@@ -0,0 +1,123 @@
+SUBROUTINE SP2GP_GAUSS4PY(KRETURNCODE, KSIZEJ, KTRUNC, KNUMMAXRESOL, KGPTOT, KSLOEN, KLOEN, KSIZE, &
+                          & LGRADIENT, LREORDER, PSPEC, PGPT, PGPTM, PGPTL)
+! ** PURPOSE
+!    Transform spectral coefficients into grid-point values
+!
+! ** DUMMY ARGUMENTS
+!    KSIZEJ: Number of latitudes
+!    KTRUNC: troncature
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    KGPTOT: number of grid-points
+!    KSLOEN: Size of KLOEN
+!    KLOEN:
+!    KSIZE: Size of PSPEC
+!    LREORDER: switch to reorder spectral coefficients or not
+!    LGRADIENT: switch to compute or not gradient
+!    PSPEC: spectral coefficient array
+!    PGPT: grid-point field
+!    PGPTM: N-S derivative if LGRADIENT
+!    PGPTL: E-W derivative if LGRADIENT
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan., S. Riette: w_spec_setup interface modified
+!    March, 2016, A.Mary: LREORDER
+!    Sept., 2016, A.Mary: LGRADIENT
+!
+! I. Dummy arguments declaration
+USE PARKIND1, ONLY : JPRB
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNC
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+INTEGER(KIND=8), INTENT(IN) :: KGPTOT
+INTEGER(KIND=8), INTENT(IN) :: KSLOEN
+INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
+INTEGER(KIND=8), INTENT(IN) :: KSIZE
+LOGICAL, INTENT(IN) :: LGRADIENT
+LOGICAL, INTENT(IN) :: LREORDER
+REAL(KIND=8), DIMENSION(KSIZE),  INTENT(IN)  :: PSPEC
+REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPT
+REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPTM
+REAL(KIND=8), DIMENSION(KGPTOT), INTENT(OUT) :: PGPTL
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+LOGICAL :: LLSTOP
+INTEGER :: IIDENTRESOL
+INTEGER :: JI, JM, JN
+INTEGER, DIMENSION(0:KTRUNC) :: NASM0
+REAL(KIND=8), DIMENSION(1, KSIZE) :: ZSPBUF
+REAL(KIND=JPRB), DIMENSION(KGPTOT, 3, 1) :: ZGPBUF
+REAL(KIND=8) :: ZDELTAX, ZDELTAY
+#include "trans_inq.h"
+#include "inv_trans.h"
+
+ILOEN(:)=KLOEN(:)
+ISIZEI=0
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=0
+IPHYSICALSIZEJ=0
+ITRUNCX=KTRUNC
+ITRUNCY=0
+INUMMAXRESOL=KNUMMAXRESOL
+!
+! III. Setup
+ZDELTAX=0.
+ZDELTAY=0.
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
+                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
+!
+! IV. Transformation
+IF (LREORDER) THEN
+  ! IV.a Shape of coefficient array
+  IF (.NOT. LLSTOP) THEN
+    JI=1
+    DO JN=0, KTRUNC
+      NASM0(JN)=JI
+      JI=JI+1+JN+(JN+1)
+    ENDDO
+  ENDIF
+
+  ! IV.b Reordering
+  IF(.NOT. LLSTOP) THEN
+    ZSPBUF(1,:)=0.
+    JI=1
+    DO JM=0, KTRUNC
+      DO JN=JM, KTRUNC
+        ZSPBUF(1,JI)=PSPEC(NASM0(JN)+JM)
+        JI=JI+1
+        IF(JM==0) THEN
+          ZSPBUF(1,JI)=0
+        ELSE
+          ZSPBUF(1,JI)=PSPEC(NASM0(JN)-JM)
+        ENDIF
+        JI=JI+1
+      ENDDO
+    ENDDO
+  ENDIF
+ELSE
+  ZSPBUF(1,:) = PSPEC(:)
+ENDIF
+
+! IV.c Inverse transform
+IF (.NOT. LLSTOP) THEN
+  IF (.NOT. LGRADIENT) THEN
+    CALL INV_TRANS(PSPSCALAR=REAL(ZSPBUF(:,:),KIND=JPRB), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
+    PGPT(:)=REAL(ZGPBUF(:,1,1),KIND=8)
+  ELSE
+    CALL INV_TRANS(PSPSCALAR=REAL(ZSPBUF(:,:),KIND=JPRB), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL, LDSCDERS=.TRUE.)
+    PGPT(:)=REAL(ZGPBUF(:,1,1),KIND=8)
+    PGPTM(:)=REAL(ZGPBUF(:,2,1),KIND=8)
+    PGPTL(:)=REAL(ZGPBUF(:,3,1),KIND=8)
+  ENDIF
+ENDIF
+END SUBROUTINE SP2GP_GAUSS4PY
diff --git a/src/ectrans4py/sp2gp_lam4py.F90 b/src/ectrans4py/sp2gp_lam4py.F90
new file mode 100644
index 000000000..17657966f
--- /dev/null
+++ b/src/ectrans4py/sp2gp_lam4py.F90
@@ -0,0 +1,140 @@
+SUBROUTINE SP2GP_LAM4PY(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
+                         &KTRUNCX, KTRUNCY, KNUMMAXRESOL, KSIZE, LGRADIENT, LREORDER, PDELTAX, PDELTAY, &
+                         &PSPEC, PGPT, PGPTM, PGPTL)
+! ** PURPOSE
+!    Transform spectral coefficients into grid-point values
+!
+! ** DUMMY ARGUMENTS
+!    KRETURNCODE: error code
+!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone)
+!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field
+!    KTRUNCX, KTRUNCY: troncatures
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    KSIZE: size of PSPEC
+!    LREORDER: switch to reorder spectral coefficients or not
+!    LGRADIENT: switch to compute or not gradient
+!    PDELTAX: x resolution
+!    PDELTAY: y resolution
+!    PSPEC: spectral coefficient array
+!    PGPT: grid-point field
+!    PGPTM: N-S derivative if LGRADIENT
+!    PGPTL: E-W derivative if LGRADIENT
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    5 Jan., S. Riette: PDELTAX, PDELTAY, LGRADIENT, PGPTM and PGPTL added
+!    March, 2016, A.Mary: LREORDER
+!
+! I. Dummy arguments declaration
+USE PARKIND1, ONLY : JPRB
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSIZEI, KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNCX, KTRUNCY
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+INTEGER(KIND=8), INTENT(IN) :: KSIZE
+LOGICAL, INTENT(IN) :: LGRADIENT
+LOGICAL, INTENT(IN) :: LREORDER
+REAL(KIND=8), INTENT(IN) :: PDELTAX
+REAL(KIND=8), INTENT(IN) :: PDELTAY
+REAL(KIND=8), DIMENSION(KSIZE), INTENT(IN) :: PSPEC
+REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPT
+REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPTM
+REAL(KIND=8), DIMENSION(KSIZEI*KSIZEJ), INTENT(OUT) :: PGPTL
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(0:KTRUNCX) :: IESM0
+INTEGER :: IGPTOT, ISPEC
+INTEGER, DIMENSION(0:KTRUNCY) :: ISPECINI, ISPECEND
+REAL(KIND=8), DIMENSION(1, KSIZE) :: ZSPBUF
+REAL(KIND=JPRB), DIMENSION(KSIZEI*KSIZEJ, 3, 1) :: ZGPBUF
+INTEGER :: JI, JM, JN, IINDEX, IIDENTRESOL
+LOGICAL :: LLSTOP
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+INTEGER, DIMENSION(1) :: ILOEN
+
+#include "einv_trans.h"
+#include "etrans_inq.h"
+
+KRETURNCODE=0
+LLSTOP=.FALSE.
+ISIZEI=KSIZEI
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=KPHYSICALSIZEI
+IPHYSICALSIZEJ=KPHYSICALSIZEJ
+ITRUNCX=KTRUNCX
+ITRUNCY=KTRUNCY
+INUMMAXRESOL=KNUMMAXRESOL
+ILOEN(:)=0
+
+! III. Setup
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .TRUE., 1, &
+                  &PDELTAX, PDELTAY, IIDENTRESOL, LLSTOP)
+
+! IV. Transformation
+
+! IV.a Shape of coefficient array
+!IGPTOT is the total number of points in grid-point space
+!ISPEC is the number of spectral coefficients
+!IESM0(m) is the index of spectral coefficient (m,0) in model
+!ISPECINI(n) is the index of the first of the 4 spectral coefficient (0,n) in FA file
+!ISPECEND(n) is the index of the last of the last 4 spectral coefficients (:,n) in FA file
+IF (.NOT. LLSTOP) THEN
+  CALL ETRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KESM0=IESM0)
+  JI=1
+  DO JN=0, ITRUNCY
+    ISPECINI(JN)=(JI-1)*4+1
+    JI=JI+COUNT(IESM0(1:ITRUNCX)-IESM0(0:ITRUNCX-1)>JN*4)
+    IF (ISPEC-IESM0(ITRUNCX)>JN*4) JI=JI+1
+    ISPECEND(JN)=(JI-1)*4
+  ENDDO
+ENDIF
+
+! III.b Reordering
+! reorder Aladin :  file ordering = coeffs per blocks of m, 4 reals per coeff
+!           Aladin array ordering = coeffs per blocks of n, 4 reals per coeff
+IF (LREORDER) THEN
+  IF (.NOT. LLSTOP) THEN
+    ZSPBUF(:,:)=0.
+    JI=1
+    DO JM=0,ITRUNCX+1
+      DO JN=0,ITRUNCY
+        IF (ISPECINI(JN)+JM*4+3<=ISPECEND(JN)) THEN
+          DO IINDEX=ISPECINI(JN)+JM*4, ISPECINI(JN)+JM*4+3
+            ZSPBUF(1,JI)=PSPEC(IINDEX)
+            JI=JI+1
+          ENDDO
+        ENDIF
+      ENDDO
+    ENDDO
+    IF (JI/=ISPEC+1) THEN
+      PRINT*, "Internal error in SP2GP_LAM4PY (spectral reordering)"
+      KRETURNCODE=-999
+      LLSTOP=.TRUE.
+    ENDIF
+  ENDIF
+ELSE
+  ZSPBUF(1,:) = PSPEC(:)
+ENDIF
+
+! III.c Inverse transform
+IF (.NOT. LLSTOP) THEN
+  IF (.NOT. LGRADIENT) THEN
+      CALL EINV_TRANS(PSPSCALAR=REAL(ZSPBUF(:,:),KIND=JPRB), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL)
+      PGPT(:)=REAL(ZGPBUF(:,1,1),KIND=8)
+  ELSE
+      CALL EINV_TRANS(PSPSCALAR=REAL(ZSPBUF(:,:),KIND=JPRB), PGP=ZGPBUF(:,:,:), KRESOL=IIDENTRESOL, LDSCDERS=.TRUE.)
+      PGPT(:)=REAL(ZGPBUF(:,1,1),KIND=8)
+      PGPTM(:)=REAL(ZGPBUF(:,2,1),KIND=8)
+      PGPTL(:)=REAL(ZGPBUF(:,3,1),KIND=8)
+  ENDIF
+ENDIF
+
+END SUBROUTINE SP2GP_LAM4PY
diff --git a/src/ectrans4py/spec_setup4py.F90 b/src/ectrans4py/spec_setup4py.F90
new file mode 100644
index 000000000..644962e3a
--- /dev/null
+++ b/src/ectrans4py/spec_setup4py.F90
@@ -0,0 +1,160 @@
+SUBROUTINE SPEC_SETUP4PY(KRETURNCODE, KSIZEI, KSIZEJ, KPHYSICALSIZEI, KPHYSICALSIZEJ, &
+                       &KTRUNCX, KTRUNCY, KNUMMAXRESOL, KLOEN, LDLAM, &
+                       &KSIZEKLOEN, PDELTAX, PDELTAY, &
+                       &KIDENTRESOL, LDSTOP)
+! ** PURPOSE
+!    Setup spectral transform for LAM and global
+!
+! ** DUMMY ARGUMENTS
+!    KRETURNCODE: error code
+!    KSIZEI, KSIZEJ: size of grid-point field (with extension zone for LAM), put max size for KSIZEI in global
+!    KPHYSICALSIZEI, KPHYSICALSIZEJ: size of physical part of grid-point field for LAM (put 0 for global)
+!    KTRUNCX, KTRUNCY: troncatures for LAM (only KTRUNCX is used for global, put 0 for KTRUNCY)
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    KLOEN: number of points on each latitude row
+!    KSIZEKLOEN: size of KLOEN array
+!    PDELTAX: x resolution
+!    PDELTAY: y resolution
+!    LDLAM: LAM (.TRUE.) or global (.FALSE.)
+!    KIDENTRESOL: identification of resolution
+!    LDSTOP: exception raised?
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan 2016, S. Riette: PDELTAX and PDELTAY added
+!    31 Jan 2019 R. El Khatib fix for single precision compilation
+!
+! I. Dummy arguments declaration
+USE PARKIND1, ONLY : JPRB
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER, INTENT(IN) :: KSIZEI, KSIZEJ
+INTEGER, INTENT(IN) :: KPHYSICALSIZEI, KPHYSICALSIZEJ
+INTEGER, INTENT(IN) :: KTRUNCX, KTRUNCY
+INTEGER, INTENT(IN) :: KNUMMAXRESOL
+INTEGER, DIMENSION(KSIZEKLOEN), INTENT(IN) :: KLOEN
+LOGICAL, INTENT(IN) :: LDLAM
+INTEGER, INTENT(IN) :: KSIZEKLOEN
+REAL(KIND=8), INTENT(IN) :: PDELTAX
+REAL(KIND=8), INTENT(IN) :: PDELTAY
+INTEGER, INTENT(OUT) :: KIDENTRESOL
+LOGICAL, INTENT(OUT) :: LDSTOP
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(2*KSIZEJ) :: ILOEN
+INTEGER :: JI
+LOGICAL, SAVE :: LLFIRSTCALL=.TRUE.
+LOGICAL :: LLNEWRESOL
+INTEGER, SAVE :: INBRESOL=0
+INTEGER(KIND=8) :: ICODEILOEN
+INTEGER, SAVE :: INUMMAXRESOLSAVE=-1
+INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ITRUNCXSAVE, ITRUNCYSAVE, &
+                                            IPHYSICALSIZEISAVE, &
+                                            IPHYSICALSIZEJSAVE, &
+                                            ISIZEISAVE, ISIZEJSAVE, &
+                                            IIDENTRESOLSAVE
+INTEGER(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: ILOENSAVE
+REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: ZDELTAXSAVE, &
+                                                        ZDELTAYSAVE
+REAL(KIND=8) :: ZEXWN, ZEYWN
+
+#include "setup_trans0.h"
+#include "esetup_trans.h"
+#include "setup_trans.h"
+
+KRETURNCODE=0
+LDSTOP=.FALSE.
+! III. Setup
+
+! III.a Setup LAM and global spectral transform - all resolutions
+! Maximum number of resolution is set now and cannot be change anymore
+IF (LLFIRSTCALL) THEN
+  !This code is called only once, whatever is the number of resolutions
+  CALL SETUP_TRANS0(KPRINTLEV=0, LDMPOFF=.TRUE., KMAX_RESOL=KNUMMAXRESOL)
+  ALLOCATE(ITRUNCXSAVE(KNUMMAXRESOL))
+  ALLOCATE(ITRUNCYSAVE(KNUMMAXRESOL))
+  ALLOCATE(IPHYSICALSIZEISAVE(KNUMMAXRESOL))
+  ALLOCATE(IPHYSICALSIZEJSAVE(KNUMMAXRESOL))
+  ALLOCATE(ISIZEJSAVE(KNUMMAXRESOL))
+  ALLOCATE(ISIZEISAVE(KNUMMAXRESOL))
+  ALLOCATE(ILOENSAVE(KNUMMAXRESOL))
+  ALLOCATE(IIDENTRESOLSAVE(KNUMMAXRESOL))
+  ALLOCATE(ZDELTAXSAVE(KNUMMAXRESOL))
+  ALLOCATE(ZDELTAYSAVE(KNUMMAXRESOL))
+  ITRUNCXSAVE=-1
+  ITRUNCYSAVE=-1
+  IPHYSICALSIZEISAVE=-1
+  IPHYSICALSIZEJSAVE=-1
+  ISIZEJSAVE=-1
+  ISIZEISAVE=-1
+  ILOENSAVE=-1
+  IIDENTRESOLSAVE=-1
+  ZDELTAXSAVE=-1.
+  ZDELTAXSAVE=-1.
+  LLFIRSTCALL=.FALSE.
+  INUMMAXRESOLSAVE=KNUMMAXRESOL
+ENDIF
+!
+! III.b Is-it a new resolution?
+LLNEWRESOL=.TRUE.
+IF(LDLAM) THEN
+  ILOEN(:)=KSIZEI
+ELSE
+  ILOEN(:)=0
+  ILOEN(1:MIN(SIZE(ILOEN),SIZE(KLOEN)))=KLOEN(1:MIN(SIZE(ILOEN),SIZE(KLOEN)))
+ENDIF
+ICODEILOEN=0
+DO JI=1, SIZE(ILOEN)
+  ICODEILOEN=ICODEILOEN+ILOEN(JI)*JI**4
+ENDDO
+DO JI=1, INBRESOL
+  IF (KTRUNCX==ITRUNCXSAVE(JI) .AND. KTRUNCY==ITRUNCYSAVE(JI) .AND. &
+   &KPHYSICALSIZEI==IPHYSICALSIZEISAVE(JI) .AND. &
+   &KPHYSICALSIZEJ==IPHYSICALSIZEJSAVE(JI) .AND. &
+   &KSIZEJ==ISIZEJSAVE(JI) .AND. KSIZEI==ISIZEISAVE(JI) .AND. &
+   &ICODEILOEN==ILOENSAVE(JI) .AND. &
+   &PDELTAX==ZDELTAXSAVE(JI) .AND. PDELTAY==ZDELTAYSAVE(JI)) THEN
+    KIDENTRESOL=IIDENTRESOLSAVE(JI)
+    LLNEWRESOL=.FALSE.
+  ENDIF
+ENDDO
+IF(LLNEWRESOL) THEN
+  INBRESOL=INBRESOL+1
+  IF(INBRESOL>INUMMAXRESOLSAVE) THEN
+    PRINT*, "Error in SPEC_SETUP4PY : Maximum number of resolution is exceeded."
+    KRETURNCODE=-999
+    LDSTOP=.TRUE.
+  ENDIF
+ENDIF
+!
+! III.c Setup LAM or global spectral transform - once by resolution
+IF(LLNEWRESOL .AND. .NOT. LDSTOP) THEN
+  ! The following code is exectuded once for each resolution
+  ITRUNCXSAVE(INBRESOL)=KTRUNCX
+  ITRUNCYSAVE(INBRESOL)=KTRUNCY
+  IPHYSICALSIZEISAVE(INBRESOL)=KPHYSICALSIZEI
+  IPHYSICALSIZEJSAVE(INBRESOL)=KPHYSICALSIZEJ
+  ISIZEISAVE(INBRESOL)=KSIZEI
+  ISIZEJSAVE(INBRESOL)=KSIZEJ
+  ILOENSAVE(INBRESOL)=ICODEILOEN
+  ZDELTAXSAVE(INBRESOL)=PDELTAX
+  ZDELTAYSAVE(INBRESOL)=PDELTAY
+  IF(LDLAM) THEN
+    ZEXWN=2*3.141592653589797/(KSIZEI*PDELTAX)
+    ZEYWN=2*3.141592653589797/(KSIZEJ*PDELTAY)
+    CALL ESETUP_TRANS(KMSMAX=ITRUNCXSAVE(INBRESOL), KSMAX=ITRUNCYSAVE(INBRESOL), &
+                     &KDGUX=IPHYSICALSIZEJSAVE(INBRESOL), &
+                     &KDGL=ISIZEJSAVE(INBRESOL), KLOEN=ILOEN(:), KRESOL=IIDENTRESOLSAVE(INBRESOL), &
+                     &PEXWN=REAL(ZEXWN,KIND=JPRB), PEYWN=REAL(ZEYWN,KIND=JPRB))
+  ELSE
+    PRINT*, "Setup spectral transform"
+    CALL SETUP_TRANS(KSMAX=ITRUNCXSAVE(INBRESOL), KDGL=ISIZEJSAVE(INBRESOL), &
+                    &KLOEN=ILOEN(1:ISIZEJSAVE(INBRESOL)), KRESOL=IIDENTRESOLSAVE(INBRESOL))
+    PRINT*, "End Setup spectral transform"
+  ENDIF
+  KIDENTRESOL=IIDENTRESOLSAVE(INBRESOL)
+ENDIF
+END SUBROUTINE SPEC_SETUP4PY
+
diff --git a/src/ectrans4py/trans_inq4py.F90 b/src/ectrans4py/trans_inq4py.F90
new file mode 100644
index 000000000..f989ef175
--- /dev/null
+++ b/src/ectrans4py/trans_inq4py.F90
@@ -0,0 +1,70 @@
+SUBROUTINE TRANS_INQ4PY(KRETURNCODE, KSIZEJ, KTRUNC, KSLOEN, KLOEN, KNUMMAXRESOL, &
+                      &KGPTOT, KSPEC, KNMENG)
+! ** PURPOSE
+!    Simplified wrapper to TRANS_INQ
+!
+! ** DUMMY ARGUMENTS
+!    KSIZEJ: number of latitudes in grid-point space
+!    KTRUNC: troncature
+!    KSLOEN: Size of KLOEN
+!    KLOEN: number of points on each latitude row
+!    KNUMMAXRESOL: maximum number of troncatures handled
+!    KGPTOT: number of gridpoints
+!    KSPEC: number of spectral coefficients
+!    KNMENG: cut-off zonal wavenumber
+!
+! ** AUTHOR
+!    9 April 2014, S. Riette
+!
+! ** MODIFICATIONS
+!    6 Jan., S. Riette: w_spec_setup interfaced modified
+!
+! I. Dummy arguments declaration
+IMPLICIT NONE
+INTEGER(KIND=8), INTENT(OUT) :: KRETURNCODE
+INTEGER(KIND=8), INTENT(IN) :: KSIZEJ
+INTEGER(KIND=8), INTENT(IN) :: KTRUNC
+INTEGER(KIND=8), INTENT(IN) :: KSLOEN
+INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(IN) :: KLOEN
+INTEGER(KIND=8), INTENT(IN) :: KNUMMAXRESOL
+INTEGER(KIND=8), INTENT(OUT) :: KGPTOT
+INTEGER(KIND=8), INTENT(OUT) :: KSPEC
+INTEGER(KIND=8), DIMENSION(KSLOEN), INTENT(OUT) :: KNMENG
+!
+! II. Local variables declaration
+INTEGER, DIMENSION(SIZE(KLOEN)) :: ILOEN
+INTEGER :: ISIZEI, ISIZEJ, &
+         & IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+         & ITRUNCX, ITRUNCY, &
+         & INUMMAXRESOL
+LOGICAL :: LLSTOP
+INTEGER :: IIDENTRESOL
+INTEGER :: IGPTOT, ISPEC
+INTEGER, DIMENSION(SIZE(KLOEN)) :: INMENG
+REAL(KIND=8) :: ZDELTAX, ZDELTAY
+#include "trans_inq.h"
+
+ILOEN(:)=KLOEN(:)
+ISIZEI=0
+ISIZEJ=KSIZEJ
+IPHYSICALSIZEI=0
+IPHYSICALSIZEJ=0
+ITRUNCX=KTRUNC
+ITRUNCY=0
+INUMMAXRESOL=KNUMMAXRESOL
+INMENG(:)=KNMENG(:)
+!
+! III. Setup
+ZDELTAX=0.
+ZDELTAY=0.
+CALL SPEC_SETUP4PY(KRETURNCODE, ISIZEI, ISIZEJ, IPHYSICALSIZEI, IPHYSICALSIZEJ, &
+                  &ITRUNCX, ITRUNCY, INUMMAXRESOL, ILOEN, .FALSE., SIZE(ILOEN), &
+                  &ZDELTAX, ZDELTAY, IIDENTRESOL, LLSTOP)
+IF (.NOT. LLSTOP) THEN
+  CALL TRANS_INQ(KRESOL=IIDENTRESOL, KGPTOT=IGPTOT, KSPEC=ISPEC, KNMENG=INMENG)
+  KGPTOT=IGPTOT
+  KSPEC=ISPEC
+  KNMENG=INMENG
+ENDIF
+!
+END SUBROUTINE TRANS_INQ4PY