diff --git a/.gitignore b/.gitignore
index bcecafd7..066ab9c0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,8 @@
 *.o
 *.a
+*.so
 *~
+.*.swp
 *_out.nc
 config_*.nam
 *_sza.nc
@@ -11,3 +13,14 @@ mod
 /bin/ecrad
 /bin/ecrad_ifs
 /bin/ecrad_ifs_blocked
+
+# Python binding
+/pyecrad/VERSION
+/pyecrad/data
+/tmp
+/dist
+/wheelhouse
+__pycache__
+pyecrad.egg-info
+/test/pyecrad/control.nc
+/test/pyecrad/experiment.nc
diff --git a/Makefile b/Makefile
index 2afecc0f..2ce418fb 100644
--- a/Makefile
+++ b/Makefile
@@ -73,6 +73,7 @@ export FCFLAGS = $(WARNFLAGS) $(BASICFLAGS) $(CPPFLAGS) -I../include \
 	$(OPTFLAGS) $(DEBUGFLAGS) $(NETCDF_INCLUDE) $(OMPFLAG)
 export LIBS    = $(LDFLAGS) -L../lib -lradiation -lutilities \
 	-lifsrrtm -lifsaux $(FCLIBS) $(NETCDF_LIB) $(OMPFLAG)
+export SHAREDLIBFLAGS
 
 # Do we include Dr Hook from ECMWF's fiat library?
 ifdef FIATDIR
@@ -158,6 +159,11 @@ ifsdriver: libifsaux libifsrrtm libutilities libradiation libifs
 test_programs: driver
 	cd driver && $(MAKE) test_programs
 
+python: build
+	cd driver && $(MAKE) python
+	cp -f VERSION pyecrad/
+	rm -rf pyecrad/data; cp -r data pyecrad/
+
 symlinks: clean-symlinks
 	cd practical && ln -s ../bin/ecrad
 	cd practical && ln -s ../data
@@ -197,6 +203,9 @@ clean-mods:
 clean-symlinks:
 	rm -f practical/ecrad practical/data
 
+clean-python:
+	rm -rf tmp dist wheelhouse
+
 clean-autosaves:
 	rm -f *~ .gitignore~ */*~ */*/*~
 
diff --git a/Makefile_include.gfortran b/Makefile_include.gfortran
index 5eaa3fd4..c52b82f3 100644
--- a/Makefile_include.gfortran
+++ b/Makefile_include.gfortran
@@ -8,7 +8,9 @@ CPPFLAGS = -cpp
 # flag "-fconvert=big-endian" here because the RRTM input files are
 # big endian Fortran unformatted files, but now the file ordering has
 # been specified at the OPEN command so no compiler flags are needed.
-BASICFLAGS = -J../mod -fno-range-check
+BASICFLAGS = -J../mod -fno-range-check -fPIC
+
+SHAREDLIBFLAGS = --shared
 
 # OpenMP flag; type "make OMPFLAG=-DNO_OPENMP" to compile with OpenMP
 # disabled
diff --git a/README.md b/README.md
index 591018f1..89490792 100644
--- a/README.md
+++ b/README.md
@@ -16,7 +16,8 @@ optics, cloud optics and solver are completely separated (see
 `radiation/radiation_interface.F90` where they are called in sequence),
 thereby facilitating future changes where different gas models or
 solvers may be switched in and out independently. The offline code is
-parallelized using OpenMP.
+parallelized using OpenMP. In addition, the offline code can be called
+from the pyecrad python package.
 
 Five solvers are currently available:
 
@@ -94,6 +95,8 @@ The subdirectories are as follows:
 
 - `practical` - exercises to get started with ecRad
 
+- `pyecrad` - python source code for the pyecrad package
+
 
 ## Compilation
 
@@ -138,6 +141,15 @@ Fortran compiler.
    `$FIATDIR/module/fiat/yomhook.mod` can be found at build time.
 
 
+## Python package
+
+The python package can be installed by entering the command `pip install .`
+in the root directory of ecrad.
+
+Wheels are built and/or pushed on pypi using the `pyecrad_wheel.sh` script.
+Help on this command can be obtained with `pyecrad_wheel.sh -h`.
+
+
 ## Testing
 
 The offline driver is run via
diff --git a/driver/Makefile b/driver/Makefile
index 9670382a..6095cde4 100644
--- a/driver/Makefile
+++ b/driver/Makefile
@@ -10,8 +10,9 @@ OBJECTS := $(SOURCES:.F90=.o)
 EXECUTABLE = ../bin/ecrad
 IFS_EXECUTABLE = ../bin/ecrad_ifs
 IFS_BLOCKED_EXECUTABLE = ../bin/ecrad_ifs_blocked
+PYSO = ../pyecrad/libecrad4py.so
 
-all: driver ifs_driver test_programs
+all: driver ifs_driver test_programs python
 
 driver: $(EXECUTABLE)
 
@@ -19,6 +20,8 @@ ifs_driver: $(IFS_EXECUTABLE) $(IFS_BLOCKED_EXECUTABLE)
 
 test_programs: $(TEST_PROGRAMS)
 
+python: $(PYSO)
+
 # Link ecrad executable; add "-lifs" if you want to use the "satur"
 # routine in ecrad_driver.F90
 $(EXECUTABLE): $(OBJECTS) ../lib/*.a ecrad_driver.o
@@ -36,6 +39,9 @@ test_%: test_%.F90 ../lib/*.a
 #$(TEST): $(TEST).F90 ../lib/*.a
 #	$(FC) $(FCFLAGS) $(TEST).F90 $(LIBS) -o $(TEST)
 
+$(PYSO): ../lib/*.a ecrad4py.o
+	$(FC) $(SHAREDLIBFLAGS) $(FCFLAGS) ecrad4py.o $(OBJECTS) -lifs $(LIBS) -o $(PYSO)
+
 # Note that the dependence on mod files can mean that rerunning "make"
 # recreates the executable
 %.o: %.F90 ../lib/*.a
@@ -43,11 +49,11 @@ test_%: test_%.F90 ../lib/*.a
 
 clean:
 	rm -f *.o $(EXECUTABLE) $(IFS_EXECUTABLE) $(IFS_BLOCKED_EXECUTABLE) \
-		$(TEST_PROGRAMS)
+		$(TEST_PROGRAMS) $(PYSO)
 
 ecrad_driver.o: ecrad_driver_config.o ecrad_driver_read_input.o
 ecrad_ifs_driver.o: ecrad_driver_config.o ecrad_driver_read_input.o
 ecrad_ifs_driver_blocked.o: ecrad_driver_config.o ecrad_driver_read_input.o ifs_blocking.o
 ecrad_driver_read_input.o ifs_blocking.o: ecrad_driver_config.o
 
-.PHONY: driver ifs_driver test_programs all
+.PHONY: driver ifs_driver test_programs python all
diff --git a/driver/ecrad4py.F90 b/driver/ecrad4py.F90
new file mode 100644
index 00000000..43439c1b
--- /dev/null
+++ b/driver/ecrad4py.F90
@@ -0,0 +1,401 @@
+! ecrad4py.F90 - Driver for python ECRAD radiation scheme
+!
+! (C) Copyright 2014- ECMWF.
+!
+! This software is licensed under the terms of the Apache Licence Version 2.0
+! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+!
+! In applying this licence, ECMWF does not waive the privileges and immunities
+! granted to it by virtue of its status as an intergovernmental organisation
+! nor does it submit to any jurisdiction.
+!
+! Author:  Sébastien Riette
+! Email:   sebastien.riette@meteo.fr
+!
+! This module contains two subroutine:
+! 1) setup with namelist filereading (to be called once)
+! 2) run to run the scheme on input profile
+!
+! The strucuture of the code is taken from the offline driver
+
+module ecrad4py
+
+  use radiation_config,         only : config_type
+
+  implicit none
+
+  type(config_type), dimension(:), allocatable :: config
+
+  contains
+
+    subroutine setup(namelist_file_name, directory_name, iconfig) bind(C, name='setup_')
+      use radiation_interface,      only : setup_radiation
+      use, intrinsic :: iso_c_binding, only: c_char, c_int64_t
+      character(kind=c_char), dimension(512), intent(in) :: namelist_file_name
+      character(kind=c_char), dimension(511), intent(in) :: directory_name
+      integer(kind=c_int64_t), intent(out) :: iconfig
+      character(len=size(namelist_file_name)) :: string_namelist_file_name
+      character(len=size(directory_name)) :: string_directory_name
+      integer :: jk
+      type(config_type), dimension(:), allocatable :: config_tmp
+
+      do jk=1, size(namelist_file_name)
+        string_namelist_file_name(jk:jk)=namelist_file_name(jk)
+      enddo
+      do jk=1, size(directory_name)
+        string_directory_name(jk:jk)=directory_name(jk)
+      enddo
+
+      if (.not. allocated(config)) then
+          iconfig = 1
+          allocate(config(iconfig))
+      else
+          allocate(config_tmp(size(config)))
+          config_tmp(:) = config(:)
+          iconfig = size(config) + 1
+          deallocate(config)
+          allocate(config(iconfig))
+          config(:iconfig - 1) = config_tmp(:)
+          deallocate(config_tmp)
+      endif
+
+      ! Read "radiation" namelist into radiation configuration type
+      call config(iconfig)%read(file_name=string_namelist_file_name)
+      config(iconfig)%directory_name = string_directory_name
+
+      ! Setup the radiation scheme: load the coefficients for gas and
+      ! cloud optics, currently from RRTMG
+      call setup_radiation(config(iconfig))
+    end subroutine setup
+
+    function get_IVolumeMixingRatio() bind(C, name='get_IVolumeMixingRatio_')
+      use radiation_gas,               only : IVolumeMixingRatio
+      use, intrinsic :: iso_c_binding, only: c_int64_t
+      integer(kind=c_int64_t) :: get_IVolumeMixingRatio
+      get_IVolumeMixingRatio=int(IVolumeMixingRatio, kind(get_IVolumeMixingRatio))
+    end function get_IVolumeMixingRatio
+
+    function get_IMassMixingRatio() bind(C, name='get_IMassMixingRatio_')
+      use radiation_gas,               only : IMassMixingRatio
+      use, intrinsic :: iso_c_binding, only: c_int64_t
+      integer(kind=c_int64_t) :: get_IMassMixingRatio
+      get_IMassMixingRatio=int(IMassMixingRatio, kind(get_IMassMixingRatio))
+    end function get_IMassMixingRatio
+
+    subroutine run(iconfig, ncol, nlev, nblocksize, pressure_hl, temperature_hl, solar_irradiance, &
+                  &spectral_solar_cycle_multiplier, &
+                  &cos_solar_zenith_angle, cloud_fraction, fractional_std, &
+                  &q_liquid, re_liquid, q_ice, re_ice, iseed, overlap_param, &
+                  &skin_temperature, nalbedobands, sw_albedo, sw_albedo_direct, &
+                  &nemissivitygpoints, lw_emissivity, &
+                  &q_unit, q, co2_unit, co2, o3_unit, o3, n2o_unit, n2o, &
+                  &co_unit, co, ch4_unit, ch4, o2_unit, o2, cfc11_unit, cfc11, &
+                  &cfc12_unit, cfc12, hcfc22_unit, hcfc22, ccl4_unit, ccl4, no2_unit, no2, &
+                  &naerosols, aerosols, inv_cloud_effective_size, inv_inhom_effective_size, &
+                  &lw_up, lw_dn, lw_up_clear, lw_dn_clear, cloud_cover_lw, &
+                  &sw_up, sw_dn, sw_up_clear, sw_dn_clear, cloud_cover_sw) bind(C, name='run_')
+
+      use, intrinsic :: iso_c_binding, only: c_int64_t, c_double
+      use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+      use parkind1,                 only : jprb, jprd ! Working/double precision
+
+      use radiation_interface,      only : radiation, set_gas_units
+      use radiation_thermodynamics, only : thermodynamics_type
+      use radiation_single_level,   only : single_level_type
+      use radiation_cloud,          only : cloud_type
+      use radiation_config,         only : ISolverSPARTACUS, IGasModelMonochromatic
+      use radiation_aerosol,        only : aerosol_type
+      use radiation_gas,            only : gas_type, &
+                                         & IVolumeMixingRatio, IMassMixingRatio, &
+                                         & IH2O, ICO2, IO3, IN2O, ICO, ICH4, IO2, ICFC11, ICFC12, &
+                                         & IHCFC22, ICCl4, INO2, GasName, GasLowerCaseName, NMaxGases
+      use radiation_flux,           only : flux_type
+      use radiation_io,             only : radiation_abort
+
+      integer(kind=c_int64_t), intent(in) :: iconfig, ncol, nlev, nblocksize
+      real(kind=c_double), dimension(ncol, nlev+1), intent(in) :: pressure_hl ! pressure (Pa) on half-levels
+      real(kind=c_double), dimension(ncol, nlev+1), intent(in) :: temperature_hl ! temperature (K) on half-levels
+      real(kind=c_double), intent(in) :: solar_irradiance ! solar irradiance (W m-2)
+      real(kind=c_double), intent(in) :: spectral_solar_cycle_multiplier ! +1.0 solar max, -1.0 min, 0.0 mean solar spectrum
+      real(kind=c_double), dimension(ncol), intent(in) :: cos_solar_zenith_angle ! cosine of the solar zenith angle
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: cloud_fraction ! cloud fraction
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: fractional_std ! fractional standard deviation of in-cloud water content
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: q_liquid ! liquid specific content (kg/kg)
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: re_liquid ! liquid effective radius (m)
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: q_ice ! ice specific content (kg/kg)
+      real(kind=c_double), dimension(ncol, nlev), intent(in) :: re_ice ! ice effective radius (m)
+      integer(kind=c_int64_t), dimension(ncol), intent(in) :: iseed ! Seed for random number generator in McICA
+      real(kind=c_double), dimension(ncol, nlev-1) :: overlap_param ! overlap of cloud boundaries
+      real(kind=c_double), dimension(ncol), intent(in) :: skin_temperature ! Ts (K)
+      integer(kind=c_int64_t), intent(in) :: nalbedobands
+      real(kind=c_double), dimension(ncol, nalbedobands), intent(in) :: sw_albedo
+      real(kind=c_double), dimension(ncol, nalbedobands), optional, intent(in) :: sw_albedo_direct
+      integer(kind=c_int64_t), intent(in) :: nemissivitygpoints
+      real(kind=c_double), dimension(ncol, nemissivitygpoints), intent(in) :: lw_emissivity
+      integer(kind=c_int64_t), optional, intent(in) :: q_unit ! vapour unit: mass or volume mixing ratio
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: q ! vapour mass/volume mixing ratio value
+      integer(kind=c_int64_t), optional, intent(in) :: co2_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) ::co2
+      integer(kind=c_int64_t), optional, intent(in) :: o3_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: o3
+      integer(kind=c_int64_t), optional, intent(in) :: n2o_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: n2o
+      integer(kind=c_int64_t), optional, intent(in) :: co_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: co
+      integer(kind=c_int64_t), optional, intent(in) :: ch4_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: ch4
+      integer(kind=c_int64_t), optional, intent(in) :: o2_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: o2
+      integer(kind=c_int64_t), optional, intent(in) :: cfc11_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: cfc11
+      integer(kind=c_int64_t), optional, intent(in) :: cfc12_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: cfc12
+      integer(kind=c_int64_t), optional, intent(in) :: hcfc22_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: hcfc22
+      integer(kind=c_int64_t), optional, intent(in) :: ccl4_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: ccl4
+      integer(kind=c_int64_t), optional, intent(in) :: no2_unit
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: no2
+      integer(kind=c_int64_t), intent(in) :: naerosols
+      real(kind=c_double), dimension(ncol, nlev, naerosols), optional, intent(in) :: aerosols
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: inv_cloud_effective_size
+      real(kind=c_double), dimension(ncol, nlev), optional, intent(in) :: inv_inhom_effective_size
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: lw_up
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: lw_dn
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: lw_up_clear
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: lw_dn_clear
+      real(kind=c_double), dimension(ncol), intent(out) :: cloud_cover_lw
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: sw_up
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: sw_dn
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: sw_up_clear
+      real(kind=c_double), dimension(ncol, nlev+1), intent(out) :: sw_dn_clear
+      real(kind=c_double), dimension(ncol), intent(out) :: cloud_cover_sw
+
+      type(single_level_type) :: single_level
+      type(thermodynamics_type) :: thermodynamics
+      type(cloud_type), target  :: cloud
+      type(gas_type)            :: gas
+      type(aerosol_type)        :: aerosol
+      type(flux_type)           :: flux
+
+      character(40)             :: gas_var_name
+      integer                   :: jgas
+      integer                   :: jblock, nblock, istartcol, iendcol
+
+      ! Pressure and temperature (SI units) are on half-levels, i.e. of length (ncol, nlev+1)
+      thermodynamics%pressure_hl = pressure_hl
+      thermodynamics%temperature_hl = temperature_hl
+
+      ! --------------------------------------------------------
+      ! Sun
+      ! --------------------------------------------------------
+
+      single_level%solar_irradiance = solar_irradiance
+      single_level%spectral_solar_cycle_multiplier = spectral_solar_cycle_multiplier
+      single_level%cos_sza = cos_solar_zenith_angle
+
+      ! --------------------------------------------------------
+      ! Cloud 
+      ! --------------------------------------------------------
+
+      if (config(iconfig)%do_clouds) then
+        cloud%fraction = cloud_fraction
+        cloud%fractional_std = fractional_std
+        cloud%ntype=2
+        allocate(cloud%mixing_ratio(int(ncol), nlev, cloud%ntype))
+        allocate(cloud%effective_radius(int(ncol), nlev, cloud%ntype))
+
+        cloud%mixing_ratio(:,:,1) = q_liquid
+        cloud%mixing_ratio(:,:,2) = q_ice
+        cloud%effective_radius(:,:,1) = re_liquid
+        cloud%effective_radius(:,:,2) = re_ice
+        cloud%q_liq  => cloud%mixing_ratio(:,:,1)
+        cloud%q_ice  => cloud%mixing_ratio(:,:,2)
+        cloud%re_liq => cloud%effective_radius(:,:,1)
+        cloud%re_ice => cloud%effective_radius(:,:,2)
+
+        call single_level%init_seed_simple(1, int(ncol)) ! Simple initialization of the seeds for the Monte Carlo scheme
+        single_level%iseed = int(iseed)
+
+        cloud%overlap_param = overlap_param
+
+        ! --------------------------------------------------------
+        ! Cloud properties needed by SPARTACUS
+        ! --------------------------------------------------------
+        if (config(iconfig)%i_solver_sw == ISolverSPARTACUS .or. \
+            config(iconfig)%i_solver_lw == ISolverSPARTACUS) then
+          if (present(inv_cloud_effective_size)) then
+            cloud%inv_cloud_effective_size = inv_cloud_effective_size
+            if (present(inv_inhom_effective_size)) then
+              cloud%inv_inhom_effective_size = inv_inhom_effective_size
+            endif
+          else
+            call radiation_abort('inv_cloud_effective_size array absent with SPARTACUS solver')
+          endif
+        endif ! Spartacus
+      endif ! Cloud
+
+      ! --------------------------------------------------------
+      ! Surface properties
+      ! --------------------------------------------------------
+
+      single_level%is_simple_surface = .true.
+
+      single_level%skin_temperature = skin_temperature
+      single_level%sw_albedo = sw_albedo
+      if(present(sw_albedo_direct)) then
+        single_level%sw_albedo_direct = sw_albedo_direct
+      endif
+      single_level%lw_emissivity = lw_emissivity
+
+      ! --------------------------------------------------------
+      ! Aerosol and gas concentrations
+      ! --------------------------------------------------------
+
+      if (config(iconfig)%use_aerosols) then
+        if (present(aerosols)) then
+          aerosol%mixing_ratio = aerosols
+        else
+          call radiation_abort('aerosols array absent with config%use_aerosols==.true.')
+        endif
+      endif
+
+      ! Water vapour and ozone are always in terms of mass mixing ratio
+      ! (kg/kg) and always 2D arrays with dimensions (ncol, nlev), unlike
+      ! other gases (see below)
+
+      call gas%allocate(int(ncol), int(nlev))
+
+      ! Loop through all radiatively important gases
+      do jgas = 1, NMaxGases
+        if (jgas == IH2O .and. present(q_unit) .and. present(q)) then
+          call gas%put(IH2O, int(q_unit), q)
+        else if (jgas == ICO2 .and. present(co2_unit) .and. present(co2)) then
+          call gas%put(ICO2, int(co2_unit), co2)
+        else if (jgas == IO3 .and. present(o3_unit) .and. present(o3)) then
+          call gas%put(IO3, int(o3_unit), o3)
+        else if (jgas == IN2O .and. present(n2o_unit) .and. present(n2o)) then
+          call gas%put(IN2O, int(n2o_unit), n2o)
+        else if (jgas == ICO .and. present(co_unit) .and. present(co)) then
+          call gas%put(ICO, int(co_unit), co)
+        else if (jgas == ICH4 .and. present(ch4_unit) .and. present(ch4)) then
+          call gas%put(ICH4, int(ch4_unit), ch4)
+        else if (jgas == IO2 .and. present(o2_unit) .and. present(o2)) then
+          call gas%put(IO2, int(o2_unit), o2)
+        else if (jgas == ICFC11 .and. present(cfc11_unit) .and. present(cfc11)) then
+          call gas%put(ICFC11, int(cfc11_unit), cfc11)
+        else if (jgas == ICFC12 .and. present(cfc12_unit) .and. present(cfc12)) then
+          call gas%put(ICFC12, int(cfc12_unit), cfc12)
+        else if (jgas == IHCFC22 .and. present(hcfc22_unit) .and. present(hcfc22)) then
+          call gas%put(IHCFC22, int(hcfc22_unit), hcfc22)
+        else if (jgas == ICCl4 .and. present(ccl4_unit) .and. present(ccl4)) then
+          call gas%put(ICCl4, int(ccl4_unit), ccl4)
+        else if (jgas == INO2 .and. present(no2_unit) .and. present(no2)) then
+          call gas%put(INO2, int(no2_unit), no2)
+        else
+          gas_var_name = trim(GasLowerCaseName(jgas))
+          !irank = file%get_rank(trim(gas_var_name))
+          !if (irank == 0) then
+          !  ! Store this as a well-mixed gas
+          !  call gas%put_well_mixed(jgas, IVolumeMixingRatio, well_mixed_gas_vmr)
+          !else if (irank == 2) then
+          !  call gas%put(jgas, IVolumeMixingRatio, gas_mr)
+          !else if (irank > 0) then
+          !  write(nulout,'(a,a,a)')  '***  Error: ', trim(gas_var_name), ' does not have 0 or 2 dimensions'
+          !  stop
+          !end if
+        end if
+      end do
+
+      ! --------------------------------------------------------
+      ! Call radiation scheme
+      ! --------------------------------------------------------
+    
+      ! Ensure the units of the gas mixing ratios are what is required
+      ! by the gas absorption model
+      call set_gas_units(config(iconfig), gas)
+    
+      ! Compute saturation with respect to liquid (needed for aerosol hydration) call...
+      call thermodynamics%calc_saturation_wrt_liquid(1, int(ncol))
+    
+      ! Allocate memory for the flux profiles, which may include arrays
+      ! of dimension n_bands_sw/n_bands_lw, so must be called after
+      ! setup_radiation
+      call flux%allocate(config(iconfig), 1, int(ncol), int(nlev))
+        
+      ! Call the ECRAD radiation scheme
+      if (nblocksize <= 0) then
+        call radiation(int(ncol), int(nlev), 1, int(ncol), &                              
+             &  config(iconfig), single_level, thermodynamics, gas, cloud, aerosol, flux)
+      else
+        nblock = (int(ncol) - 1 + int(nblocksize)) / int(nblocksize)
+        !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME)
+        do jblock = 1, nblock
+          istartcol = (jblock - 1) * int(nblocksize) + 1
+          iendcol = min(istartcol + int(nblocksize) - 1, int(ncol))
+          call radiation(int(ncol), int(nlev), istartcol, iendcol, &
+               &  config(iconfig), single_level, thermodynamics, gas, cloud, aerosol, flux)
+        enddo
+      endif
+    
+      ! --------------------------------------------------------
+      ! Output
+      ! --------------------------------------------------------
+      if (config(iconfig)%i_gas_model_lw == IGasModelMonochromatic .and. &
+         &config(iconfig)%mono_lw_wavelength > 0.0_jprb) then
+          call radiation_abort('lw flux unit is W m-3')
+      endif
+    
+      if (config(iconfig)%do_lw) then
+        lw_up = flux%lw_up
+        lw_dn = flux%lw_dn
+        if (config(iconfig)%do_clear) then
+          lw_up_clear = flux%lw_up_clear
+          lw_dn_clear = flux%lw_dn_clear
+        else
+          lw_up_clear = ieee_value(lw_up_clear, ieee_quiet_nan)
+          lw_dn_clear = ieee_value(lw_dn_clear, ieee_quiet_nan)
+        endif
+        if (config(iconfig)%do_clouds) then
+          cloud_cover_lw = flux%cloud_cover_lw
+        else
+          cloud_cover_lw = ieee_value(cloud_cover_lw, ieee_quiet_nan)
+        endif
+      else
+        lw_up = ieee_value(lw_up, ieee_quiet_nan)
+        lw_dn = ieee_value(lw_dn, ieee_quiet_nan)
+        lw_up_clear = ieee_value(lw_up_clear, ieee_quiet_nan)
+        lw_dn_clear = ieee_value(lw_dn_clear, ieee_quiet_nan)
+        cloud_cover_lw = ieee_value(cloud_cover_lw, ieee_quiet_nan)
+      endif
+      if (config(iconfig)%do_sw) then
+         sw_up = flux%sw_up
+         sw_dn = flux%sw_dn
+         if (config(iconfig)%do_clear) then
+           sw_up_clear = flux%sw_up_clear
+           sw_dn_clear = flux%sw_dn_clear
+         else
+           sw_up_clear = ieee_value(sw_up_clear, ieee_quiet_nan)
+           sw_dn_clear = ieee_value(sw_dn_clear, ieee_quiet_nan)
+         endif
+         if (config(iconfig)%do_clouds) then
+          cloud_cover_sw = flux%cloud_cover_sw
+        else
+          cloud_cover_sw = ieee_value(cloud_cover_sw, ieee_quiet_nan)
+        endif
+      else
+        sw_up = ieee_value(sw_up, ieee_quiet_nan)
+        sw_dn = ieee_value(sw_dn, ieee_quiet_nan)
+        sw_up_clear = ieee_value(sw_up_clear, ieee_quiet_nan)
+        sw_dn_clear = ieee_value(sw_dn_clear, ieee_quiet_nan)
+        cloud_cover_sw = ieee_value(cloud_cover_sw, ieee_quiet_nan)
+      endif
+      if (config(iconfig)%do_clouds) then
+        deallocate(cloud%mixing_ratio)
+        deallocate(cloud%effective_radius)
+      endif
+      call gas%deallocate()
+      call flux%deallocate()
+    
+    end subroutine run
+end module ecrad4py
diff --git a/pyecrad/README b/pyecrad/README
new file mode 100644
index 00000000..d62ef235
--- /dev/null
+++ b/pyecrad/README
@@ -0,0 +1,116 @@
+This directory holds the python code needed for the pyecrad package.
+
+The class that allows you to run ecrad is Ecrad and can be used as
+shown in the following examples.
+
+Example 1: namelist file, input and output as netcdf files
+```
+import pyecrad
+
+ecrad = pyecrad.Ecrad('namelist.nam')
+ecrad.driver('input.nc', 'output.nc')
+```
+
+Example 2: equivalent to example 1 but using explicitely read and
+           write methods
+```
+import pyecrad
+
+ecrad = pyecrad.Ecrad('namelist.nam')
+input = ecrad.read('input.nc')
+result = ecrad.exec(**input)
+ecrad.write('output.nc', result)
+```
+
+Example 3: all in memory example
+```
+import matplotlib.pyplot as plt
+import numpy
+import pyecrad
+
+# Build namelist as a dict
+namelist = {'radiation': {
+    # GENERAL
+    #
+    'iverbose'            : 1, 
+    'iverbosesetup'       : 1,
+    'do_surface_sw_spectral_flux' : False,  # Save surface fluxes in each band?
+    #
+    # CLOUDS
+    #
+    'use_general_cloud_optics' : False,
+    'ice_model_name'      : "Fu-IFS",       # Can be "Fu-IFS" or "Yi"
+    'sw_solver_name'      : "Tripleclouds", # "Tripleclouds", "McICA" or "SPARTACUS"
+    'lw_solver_name'      : "Tripleclouds", # "Tripleclouds", "McICA" or "SPARTACUS"
+    'overlap_scheme_name' : "Exp-Ran",      # McICA also accepts Max-Ran or Exp-Exp
+    'do_lw_cloud_scattering' : True,        # Clouds scatter in the longwave?
+    'gas_model_name'      : "RRTMG-IFS",    # "RRTMG-IFS" or "ECCKD"
+    #
+    # AEROSOLS
+    #
+    'use_aerosols'             : True,      # Radiation sees aerosols?
+    'use_general_aerosol_optics':False,
+    'do_lw_aerosol_scattering' : False,     # Aerosols scatter in the longwave?
+    #
+    # 11 IFS aerosol mixing ratios are stored in the ecRad input file: 1-3
+    # Sea salt, 4-6 mineral dust, 7 hydrophilic organics, 8 hydrophobic
+    # organics, 9 hydrophilic black carbon, 10 hydrophobic black carbon, 11
+    # ammonium sulfate
+    'n_aerosol_types'  : 11,   # Number of aerosol types in input file
+    #
+    # The aerosol optical properties are in this file:
+    'aerosol_optics_override_file_name' : 'aerosol_ifs_rrtm_46R1_with_NI_AM.nc',
+    #
+    # For each of the 11 mixing ratios in the input file, we need to map to
+    # one of the optical properties, where negative numbers index
+    # hydrophilic aerosol types and positive numbers index hydrophobic
+    # aerosol types, e.g. 11=black carbon, -5=sulphate.
+    'i_aerosol_type_map' : [-1, -2, -3, 1, 2, 3, -4, 10, 11, 11, -5],
+}}
+
+# Dimensions
+nlev, ncol = 5, 1
+
+# Input profiles
+input = {
+    'pressure_hl': numpy.array([[10000., 20000., 30000., 50000., 70000., 100000.]]).T,  # nlev + 1
+    'temperature_hl': numpy.array([[215., 215., 230., 250., 270., 290.]]).T,  # nlev + 1
+    'solar_irradiance': 1366.0,
+    'spectral_solar_cycle_multiplier': 0.,
+    'cos_solar_zenith_angle': numpy.array([0.9]),
+    'cloud_fraction': numpy.array([[0., 0.5, 1., 0.5, 0.]]).T,
+    'fractional_std': numpy.array([[1., 1., 1., 1., 1.]]).T,
+    'q_liquid': numpy.array([[0., 2.E-4, 5.E-4, 2.E-4, 0.]]).T,
+    're_liquid': numpy.array([[0., 1.E-5, 1.E-5, 1.E-5, 0.]]).T,
+    'q_ice': numpy.array([[0., 0., 0., 0., 0.]]).T,
+    're_ice': numpy.array([[0., 0., 0., 0., 0.]]).T,
+    'overlap_param': numpy.array([[0.5, 0.5, 0.5, 0.5]]).T,  # nlev - 1
+    'skin_temperature': numpy.array([320.]),
+    'sw_albedo': numpy.array([[0.5]]),
+    'lw_emissivity': numpy.array([[0.9]]),
+    'q_unit': pyecrad.IMassMixingRatio,
+    'q': numpy.array([[8.E-5, 1.E-4, 5.E-4, 2.E-3, 6.E-3]]).T,
+    'aerosols': numpy.zeros((namelist['radiation']['n_aerosol_types'], nlev, ncol))
+}
+for gas in ('co2', 'o3', 'n2o', 'co', 'ch4', 'o2', 'cfc11',                                                                                                                                       
+            'cfc12', 'hcfc22', 'ccl4', 'no2'):
+    input[gas + '_unit'] = pyecrad.IVolumeMixingRatio
+    input[gas] = numpy.array([[0., 0., 0., 0., 0.]]).T
+
+# ecrad execution
+ecrad = pyecrad.Ecrad(namelist)
+result = ecrad.exec(**input)
+
+# Plot results
+print(f'cloud_cover_sw: {result['cloud_cover_sw']}')
+print(f'cloud_cover_lw: {result['cloud_cover_lw']}')
+fig, ax = plt.subplots()
+for k, v in result.items():
+    if k not in ('cloud_cover_sw', 'cloud_cover_lw', 'pressure_hl'):
+        ax.plot(v, result['pressure_hl'], label=k)
+
+ax.set_xlabel('Radiation fluxes (W)')
+ax.set_ylabel('Pressure (Pa)')
+ax.legend()
+plt.show()
+```
diff --git a/pyecrad/__init__.py b/pyecrad/__init__.py
new file mode 100755
index 00000000..c67cb1b1
--- /dev/null
+++ b/pyecrad/__init__.py
@@ -0,0 +1,25 @@
+#!/usr/bin/env python3
+
+"""
+Python wrapper for ecrad
+
+(C) Copyright 2026- ECMWF.
+
+This software is licensed under the terms of the Apache Licence Version 2.0
+which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+In applying this licence, ECMWF does not waive the privileges and immunities
+granted to it by virtue of its status as an intergovernmental organisation
+nor does it submit to any jurisdiction.
+
+Author:  Sébastien Riette
+Email:   sebastien.riette@meteo.fr
+License: see the COPYING file for details
+"""
+
+from .ecrad4py import get_IVolumeMixingRatio, get_IMassMixingRatio
+IVolumeMixingRatio = get_IVolumeMixingRatio()
+IMassMixingRatio = get_IMassMixingRatio()
+del get_IVolumeMixingRatio, get_IMassMixingRatio
+
+from .ecrad import Ecrad
diff --git a/pyecrad/ecrad.py b/pyecrad/ecrad.py
new file mode 100644
index 00000000..cff61278
--- /dev/null
+++ b/pyecrad/ecrad.py
@@ -0,0 +1,219 @@
+"""
+Ecrad class
+
+(C) Copyright 2026- ECMWF.
+
+This software is licensed under the terms of the Apache Licence Version 2.0
+which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+In applying this licence, ECMWF does not waive the privileges and immunities
+granted to it by virtue of its status as an intergovernmental organisation
+nor does it submit to any jurisdiction.
+
+Author:  Sébastien Riette
+Email:   sebastien.riette@meteo.fr
+License: see the COPYING file for details
+"""
+
+import tempfile
+import f90nml
+import netCDF4
+import numpy
+from .ecrad4py import setup, run
+from . import IVolumeMixingRatio, IMassMixingRatio
+
+class Ecrad():
+    """
+    Ecrad class
+    """
+    def __init__(self, namelist, data_directory=None):
+        """
+        Ecrad class
+        :param namelist: namelist as a dictionary or a file name
+        :param data_directory: (optional) ecrad data directory
+        """
+        if isinstance(namelist, dict):
+            with tempfile.NamedTemporaryFile() as namel:
+                nml = f90nml.read(namel.name)
+                for k, v in namelist.items():
+                    nml[k] = v
+                nml.write(namel.name, force=True)
+                self._iconfig = setup(namel.name, data_directory)
+        else:
+            self._iconfig = setup(namelist, data_directory)
+        self.IVolumeMixingRatio = IVolumeMixingRatio
+        self.IMassMixingRatio = IMassMixingRatio
+
+    def exec(self,
+             pressure_hl, temperature_hl, solar_irradiance,
+             spectral_solar_cycle_multiplier,
+             cos_solar_zenith_angle, cloud_fraction, fractional_std,
+             q_liquid, re_liquid, q_ice, re_ice, overlap_param,
+             skin_temperature, sw_albedo, lw_emissivity, sw_albedo_direct=None,
+             q_unit=None, q=None, co2_unit=None, co2=None,
+             o3_unit=None, o3=None, n2o_unit=None, n2o=None,
+             co_unit=None, co=None, ch4_unit=None, ch4=None,
+             o2_unit=None, o2=None, cfc11_unit=None, cfc11=None,
+             cfc12_unit=None, cfc12=None, hcfc22_unit=None, hcfc22=None,
+             ccl4_unit=None, ccl4=None, no2_unit=None, no2=None,
+             aerosols=None,
+             inv_cloud_effective_size=None, inv_inhom_effective_size=None ,
+             nblocksize=32, iseed=None):
+        """
+        Main method executing ecrad
+        :return: dictionary hoding the different fields computed by ecrad
+        """
+        # Guess sizes from fields
+        nlev = pressure_hl.shape[0] - 1
+        ncol = pressure_hl.shape[1]
+        if aerosols is None:
+            naerosols = 0
+        else:
+            naerosols = aerosols.shape[0]
+        nemissivitygpoints = lw_emissivity.shape[0]
+        nalbedobands = sw_albedo.shape[0]
+
+        # Default value for iseed
+        if iseed is None:
+            iseed = numpy.ones((ncol, ))
+
+        result = run(self._iconfig, ncol, nlev, nblocksize,
+                     pressure_hl, temperature_hl, solar_irradiance,
+                     spectral_solar_cycle_multiplier,
+                     cos_solar_zenith_angle, cloud_fraction, fractional_std,
+                     q_liquid, re_liquid, q_ice, re_ice, iseed, overlap_param,
+                     skin_temperature, nalbedobands, sw_albedo, sw_albedo_direct,
+                     nemissivitygpoints, lw_emissivity,
+                     q_unit, q, co2_unit, co2,
+                     o3_unit, o3, n2o_unit, n2o,
+                     co_unit, co, ch4_unit, ch4,
+                     o2_unit, o2, cfc11_unit, cfc11,
+                     cfc12_unit, cfc12, hcfc22_unit, hcfc22,
+                     ccl4_unit, ccl4, no2_unit, no2,
+                     naerosols, aerosols,
+                     inv_cloud_effective_size, inv_inhom_effective_size)
+        # Converts the result into a dictionary
+        result = {name: result[i] for i, name in enumerate(
+            ('lw_up', 'lw_dn', 'lw_up_clear', 'lw_dn_clear', 'cloud_cover_lw',
+             'sw_up', 'sw_dn', 'sw_up_clear', 'sw_dn_clear', 'cloud_cover_sw'))}
+        # Add the input pressure field to the output
+        result['pressure_hl'] = pressure_hl
+
+        return result
+
+    def driver(self, input_file, output_file, iseed=None, nblocksize=32):
+        """
+        Read input netCDF file, run ecrad and write output netCDF file
+        :param input_file: input netCDF file
+        :param output_file: netCDF output file
+        :param iseed: (optional) seed
+        :param nblocksize: (optional) block size
+        """
+        return self.write(output_file, self.exec(iseed=iseed, nblocksize=nblocksize,
+                                                 **self.read(input_file)))
+
+    def read(self, input_file):
+        """
+        Read a netCDF file
+        :param input_file: netCDF input file
+        :return:  kkwargs arguments for the exec method
+        """
+        with netCDF4.Dataset(input_file, 'r') as nci:
+            if len(nci['sw_albedo'].shape) == 2:
+                sw_albedo = nci['sw_albedo'][...]
+            else:
+                sw_albedo = nci['sw_albedo'][...][numpy.newaxis, ...]
+            if len(nci['lw_emissivity'].shape) == 2:
+                lw_emissivity = nci['lw_emissivity']
+            else:
+                lw_emissivity = nci['lw_emissivity'][...][numpy.newaxis, ...]
+            fractional_std = numpy.ones(nci['cloud_fraction'].shape)
+            solar_irradiance = 1366.0  # default value set in ecrad_driver_read_input.F90
+            spectral_solar_cycle_multiplier = 0.  # default value set in ecrad_driver_read_input.F90
+
+            shape = (nci.dimensions['level'].size, nci.dimensions['column'].size)
+            gases = {}
+            for gas in ('co2', 'o3', 'n2o', 'co', 'ch4', 'o2', 'cfc11',
+                      'cfc12', 'hcfc22', 'ccl4', 'no2'):
+                if gas + '_vmr' in nci.variables:
+                    gases[gas + '_unit'] = IVolumeMixingRatio
+                    if len(nci[gas + '_vmr'].shape) == 0:
+                        gases[gas] = numpy.ndarray(shape)
+                        gases[gas][...] = nci[gas + '_vmr'][...]
+                    else:
+                        gases[gas] = nci[gas + '_vmr'][...].T
+                elif gas + '_mmr' in nci.variables:
+                    gases[gas + '_unit'] = IMassMixingRatio
+                    gases[gas] = nci[gas + '_mmr'][...].T
+
+            aerosols = numpy.moveaxis(nci['aerosol_mmr'][...], [0, 1, 2], [0, 2, 1])
+
+            result = {'pressure_hl': nci['pressure_hl'][...].T,
+                      'temperature_hl': nci['temperature_hl'][...].T,
+                      'solar_irradiance': solar_irradiance,
+                      'spectral_solar_cycle_multiplier':spectral_solar_cycle_multiplier,
+                      'cos_solar_zenith_angle': nci['cos_solar_zenith_angle'][...],
+                      'cloud_fraction': nci['cloud_fraction'][...].T,
+                      'fractional_std': fractional_std.T,
+                      'q_liquid': nci['q_liquid'][...].T,
+                      're_liquid': nci['re_liquid'][...].T,
+                      'q_ice': nci['q_ice'][...].T,
+                      're_ice': nci['re_ice'][...].T,
+                      'overlap_param': nci['overlap_param'][...].T,
+                      'skin_temperature': nci['skin_temperature'][...],
+                      'sw_albedo': sw_albedo,
+                      'lw_emissivity': lw_emissivity,
+                      'q_unit': IMassMixingRatio,
+                      'q': nci['q'][...].T,
+                      'aerosols': aerosols.T,
+                     }
+            result.update(gases)
+        return result
+
+    def write(self, output_file, result):
+        """
+        Write the output to a netCDF file
+        :param output_file: output netCDF file name
+        :param result: output dictionary from the exec method
+        """
+        with netCDF4.Dataset(output_file, 'w') as nco:
+            # Add dimensions
+            nco.createDimension('column', result['pressure_hl'].shape[1])
+            nco.createDimension('half_level', result['pressure_hl'].shape[0])
+
+            # Save outputs
+            for name, long_name, units, standard_name, vdim, result_name in \
+                [('flux_up_lw', 'Upwelling longwave flux', 'W m-2',
+                  'upwelling_longwave_flux_in_air', 'half_level', 'lw_up'),
+                 ('flux_dn_lw', 'Downwelling longwave flux', 'W m-2',
+                  'downwelling_longwave_flux_in_air', 'half_level', 'lw_dn'),
+                 ('flux_up_lw_clear', 'Upwelling clear-sky longwave flux', 'W m-2',
+                  None, 'half_level', 'lw_up_clear'),
+                 ('flux_dn_lw_clear', 'Downwelling clear-sky longwave flux', 'W m-2',
+                  None, 'half_level', 'lw_dn_clear'),
+                 ('cloud_cover_lw', 'Total cloud cover diagnosed by longwave solver', '1',
+                  'cloud_area_fraction', None, 'cloud_cover_lw'),
+
+                 ('flux_up_sw', 'Upwelling shortwave flux', 'W m-2',
+                  'upwelling_shortwave_flux_in_air', 'half_level', 'sw_up'),
+                 ('flux_dn_sw', 'Downwelling shortwave flux', 'W m-2',
+                  'downwelling_shortwave_flux_in_air', 'half_level', 'sw_dn'),
+                 ('flux_up_sw_clear', 'Upwelling clear-sky shortwave flux', 'W m-2',
+                  None, 'half_level', 'sw_up_clear'),
+                 ('flux_dn_sw_clear', 'Downwelling clear-sky shortwave flux', 'W m-2',
+                  None, 'half_level', 'sw_dn_clear'),
+                 ('cloud_cover_sw', 'Total cloud cover diagnosed by shortwave solver', '1',
+                  'cloud_area_fraction', None, 'cloud_cover_sw'),
+                 ('pressure_hl', 'Pressure at half-levels', 'Pa',
+                  None, 'half_level', 'pressure_hl'),
+                ]:
+                if vdim is None:
+                    dimension = (nco.dimensions['column'], )
+                else:
+                    dimension = (nco.dimensions['column'], nco.dimensions[vdim])
+                flux = nco.createVariable(name, result[result_name].dtype, dimension)
+                flux[...] = result[result_name].T
+                flux.long_name = long_name
+                flux.units = units
+                if standard_name is not None:
+                    flux.standard_name = standard_name
diff --git a/pyecrad/ecrad4py.py b/pyecrad/ecrad4py.py
new file mode 100644
index 00000000..37d1674c
--- /dev/null
+++ b/pyecrad/ecrad4py.py
@@ -0,0 +1,175 @@
+"""
+Low level API to ecrad
+
+(C) Copyright 2026- ECMWF.
+
+This software is licensed under the terms of the Apache Licence Version 2.0
+which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+In applying this licence, ECMWF does not waive the privileges and immunities
+granted to it by virtue of its status as an intergovernmental organisation
+nor does it submit to any jurisdiction.
+
+Author:  Sébastien Riette
+Email:   sebastien.riette@meteo.fr
+License: see the COPYING file for details
+"""
+
+import os
+import numpy
+import ctypesForFortran
+from ctypesForFortran import (IN, OUT, INOUT, MISSING, MANDATORY_AFTER_OPTIONAL as MAO,
+                              string2array)
+IN = ctypesForFortran.IN
+OUT = ctypesForFortran.OUT
+INOUT = ctypesForFortran.INOUT
+MISSING = ctypesForFortran.MISSING
+MAO = ctypesForFortran.MANDATORY_AFTER_OPTIONAL
+
+ctypesFF, handle = ctypesForFortran.ctypesForFortranFactory(
+    os.path.join(os.path.dirname(os.path.abspath(__file__)), "./libecrad4py.so"))
+
+version_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), "VERSION")
+if os.path.exists(version_file):
+    with open(version_file, 'r', encoding='utf-8') as f:
+        __version__ = f.read().strip()
+else:
+    __version__ = 'unknown'
+
+
+def close():
+    """Close the shared lib"""
+    ctypesForFortran.dlclose(handle)
+
+
+@ctypesFF()
+def setup(namelist_file_name, directory_name=None):
+    """
+    Ecrad setup
+    :param namelist_file_name: namelist file name
+    :param directory_name: path to the data directory
+    """
+    if directory_name is None:
+        directory_name = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data")
+    return ([string2array(namelist_file_name, 512),
+             string2array(directory_name, 511)],
+            [(str, (1, 512), IN),
+             (str, (1, 511), IN),
+             (numpy.int64, None, OUT)],
+            None)
+
+
+@ctypesFF()
+def get_IVolumeMixingRatio():
+    """
+    Ecrad value for volume mixing ratio
+    """
+    return ([], [], (numpy.int64, None))
+
+
+@ctypesFF()
+def get_IMassMixingRatio():
+    """
+    Ecrad value for mass mixing ratio
+    """
+    return ([], [], (numpy.int64, None))
+
+
+@ctypesFF(castInput=True, indexing='C')
+def run(iconfig, ncol, nlev, nblocksize, pressure_hl, temperature_hl, solar_irradiance,
+        spectral_solar_cycle_multiplier,
+        cos_solar_zenith_angle, cloud_fraction, fractional_std,
+        q_liquid, re_liquid, q_ice, re_ice, iseed, overlap_param,
+        skin_temperature, nalbedobands, sw_albedo, sw_albedo_direct=None,
+        nemissivitygpoints=MAO, lw_emissivity=MAO,
+        q_unit=None, q=None, co2_unit=None, co2=None,
+        o3_unit=None, o3=None, n2o_unit=None, n2o=None,
+        co_unit=None, co=None, ch4_unit=None, ch4=None,
+        o2_unit=None, o2=None, cfc11_unit=None, cfc11=None,
+        cfc12_unit=None, cfc12=None, hcfc22_unit=None, hcfc22=None,
+        ccl4_unit=None, ccl4=None, no2_unit=None, no2=None,
+        naerosols=0, aerosols=None,
+        inv_cloud_effective_size=None, inv_inhom_effective_size=None):
+
+    """
+    Ecrad simulation
+    """
+    def n2m(x):
+        """Converts None to MISSING"""
+        return MISSING if x is None else x
+    return ([iconfig, ncol, nlev, nblocksize, pressure_hl, temperature_hl, solar_irradiance,
+             spectral_solar_cycle_multiplier,
+             cos_solar_zenith_angle, cloud_fraction, fractional_std,
+             q_liquid, re_liquid, q_ice, re_ice, numpy.array(iseed), overlap_param,
+             skin_temperature, nalbedobands, sw_albedo, n2m(sw_albedo_direct),
+             nemissivitygpoints, lw_emissivity,
+             n2m(q_unit), n2m(q), n2m(co2_unit), n2m(co2), n2m(o3_unit), n2m(o3),
+             n2m(n2o_unit), n2m(n2o), n2m(co_unit), n2m(co), n2m(ch4_unit), n2m(ch4),
+             n2m(o2_unit), n2m(o2), n2m(cfc11_unit), n2m(cfc11), n2m(cfc12_unit), n2m(cfc12),
+             n2m(hcfc22_unit), n2m(hcfc22), n2m(ccl4_unit), n2m(ccl4), n2m(no2_unit), n2m(no2),
+             naerosols, n2m(aerosols),
+             n2m(inv_cloud_effective_size), n2m(inv_inhom_effective_size)],
+            [(numpy.int64, None, IN),  # iconfig
+             (numpy.int64, None, IN),  # ncol
+             (numpy.int64, None, IN),  # nlev
+             (numpy.int64, None, IN),  # nblocksize
+             (numpy.float64, (nlev + 1, ncol), IN),  # pressure (Pa) on half-levels
+             (numpy.float64, (nlev + 1, ncol), IN),  # temperature (K) on half-levels
+             (numpy.float64, None, IN),  # solar irradiance (W m-2)
+             (numpy.float64, None, IN),  # spectral_solar_cycle_multiplier ! +1.0 solar max, -1.0 min, 0.0 mean solar spectrum
+             (numpy.float64, (ncol, ), IN),  # cosine of the solar zenith angle
+             (numpy.float64, (nlev, ncol), IN),  # cloud fraction
+             (numpy.float64, (nlev, ncol), IN),  # fractional standard deviation of in-cloud water content
+             (numpy.float64, (nlev, ncol), IN),  # liquid specific content (kg/kg)
+             (numpy.float64, (nlev, ncol), IN),  # liquid effective radius (m)
+             (numpy.float64, (nlev, ncol), IN),  # ice specific content (kg/kg)
+             (numpy.float64, (nlev, ncol), IN),  # ice effective radius (m)
+             (numpy.int64, (ncol, ), IN),  # Seed for random number generator in McICA
+             (numpy.float64, (nlev - 1, ncol), IN),  # overlap_param ! overlap of cloud boundaries
+             (numpy.float64, (ncol, ), IN),  # skin_temperature ! Ts (K)
+             (numpy.int64, None, IN),  # nalbedobands
+             (numpy.float64, (nalbedobands, ncol), IN),  # sw_albedo
+             (numpy.float64, (nalbedobands, ncol), IN),  # sw_albedo_direct
+             (numpy.int64, None, IN),  # nemissivitygpoints
+             (numpy.float64, (nemissivitygpoints, ncol), IN),  # lw_emissivity
+             (numpy.int64, None, IN),  # vapour unit
+             (numpy.float64, (nlev, ncol), IN),  # vapour mixing ratio
+             (numpy.int64, None, IN),  # co2 unit
+             (numpy.float64, (nlev, ncol), IN),  # co2 mixing ratio
+             (numpy.int64, None, IN),  # o3 unit
+             (numpy.float64, (nlev, ncol), IN),  # o3 mixing ratio
+             (numpy.int64, None, IN),  # n2o_unit
+             (numpy.float64, (nlev, ncol), IN),  # n2o mixing ratio
+             (numpy.int64, None, IN),  # co_unit
+             (numpy.float64, (nlev, ncol), IN),  # co mixing ratio
+             (numpy.int64, None, IN),  # ch4_unit
+             (numpy.float64, (nlev, ncol), IN),  # ch4 mixing ratio
+             (numpy.int64, None, IN),  # o2_unit
+             (numpy.float64, (nlev, ncol), IN),  # o2 mixing ratio
+             (numpy.int64, None, IN),  # cfc11_unit
+             (numpy.float64, (nlev, ncol), IN),  # cfc11 mixing ratio
+             (numpy.int64, None, IN),  # cfc12_unit
+             (numpy.float64, (nlev, ncol), IN),  # cfc12 mixing ratio
+             (numpy.int64, None, IN),  # hcfc22_unit
+             (numpy.float64, (nlev, ncol), IN),  # hcfc22 mixing ratio
+             (numpy.int64, None, IN),  # ccl4_unit
+             (numpy.float64, (nlev, ncol), IN),  # ccl4 mixing ratio
+             (numpy.int64, None, IN),  # no2_unit
+             (numpy.float64, (nlev, ncol), IN),  # no2 mixing ratio
+             (numpy.int64, None, IN),  # naerosols
+             (numpy.float64, (naerosols, nlev, ncol), IN),  # aerosols mixing ratio
+             (numpy.float64, (nlev, ncol), IN),  # inv_cloud_effective_size
+             (numpy.float64, (nlev, ncol), IN),  # inv_inhom_effective_size
+
+             (numpy.float64, (nlev + 1, ncol), OUT),  # lw_up
+             (numpy.float64, (nlev + 1, ncol), OUT),  # lw_dn
+             (numpy.float64, (nlev + 1, ncol), OUT),  # lw_up_clear
+             (numpy.float64, (nlev + 1, ncol), OUT),  # lw_dn_clear
+             (numpy.float64, (ncol, ), OUT),  # cloud_cover_lw
+             (numpy.float64, (nlev + 1, ncol), OUT),  # sw_up
+             (numpy.float64, (nlev + 1, ncol), OUT),  # sw_dn
+             (numpy.float64, (nlev + 1, ncol), OUT),  # sw_up_clear
+             (numpy.float64, (nlev + 1, ncol), OUT),  # sw_dn_clear
+             (numpy.float64, (ncol, ), OUT),  # cloud_cover_sw
+
+            ], None)
diff --git a/pyecrad_wheel.sh b/pyecrad_wheel.sh
new file mode 100755
index 00000000..3fba9a2e
--- /dev/null
+++ b/pyecrad_wheel.sh
@@ -0,0 +1,208 @@
+#!/usr/bin/env bash
+
+# Wheel build script
+#
+# (C) Copyright 2026- ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# In applying this licence, ECMWF does not waive the privileges and immunities
+# granted to it by virtue of its status as an intergovernmental organisation
+# nor does it submit to any jurisdiction.
+#
+# Author:  Sébastien Riette
+# Email:   sebastien.riette@meteo.fr
+# License: see the COPYING file for details
+
+
+function usage {
+  echo "
+$0 target [pypi]
+$0 -h
+
+This script builds the python wheel using a container with low version of GLIBC
+allowing execution on a wide range of linux systems.
+
+Script can be called with up to two arguments: target and pypi
+
+target can be:
+- get: to get the container and dependencies source code
+- deps: to build the dependencies
+- lib: to build the library
+- wheel: to build the wheels and the source distribution
+- pypi: to push the source distributions and the wheels on PyPI
+        in this case, the pypi argument must be provided
+- all: to perform all these actions
+
+pypi is the repository on which wheels are pushed ('pypi' and 'all' targets)."
+}
+
+# versions of dependencies
+netcdf_version="v4.9.2"
+netcdf_uri=http://github.com/Unidata/netcdf-c
+netcdffortran_version="v4.6.2"
+netcdffortran_uri=https://github.com/Unidata/netcdf-fortran
+hdf5=https://resources.unidata.ucar.edu/netcdf/netcdf-4/hdf5-1.8.15.tar.bz2
+curl=https://curl.haxx.se/download/curl-8.5.0.tar.gz
+
+# container image to be used
+container_uri=docker://quay.io/pypa/manylinux2014_x86_64
+
+# local certificate (for pip)
+local_certificate=/etc/ssl/certs/ca-certificates.crt
+
+# python versions for which to build a wheel
+#export python_versions="cp310-cp310 cp311-cp311 cp312-cp312 cp313-cp313 cp313-cp313t"
+# This module does not depend on the exact version used, we pick one
+export python_versions="cp312-cp312"
+
+# end of user customisation
+# =============================================================================
+
+# working directories
+export ECRAD_ROOT=$( cd -- "$( dirname -- "${BASH_SOURCE[0]}" )" &> /dev/null && pwd )
+export TMPWORKDIR=$ECRAD_ROOT/tmp
+export SRC=$TMPWORKDIR/src
+mkdir -p $TMPWORKDIR
+mkdir -p $SRC
+
+# container stuff
+export CONTAINER_SIF_PATH="$TMPWORKDIR/$(echo $container_uri | awk -F '/' '{print $NF}').sif"
+export CONTAINER_ROOT=/work  # bind path from within the container
+export SINGULARITY_BINDPATH=$ECRAD_ROOT:$CONTAINER_ROOT  # binding out:in
+export SINGULARITY_TMPDIR=$TMPWORKDIR/singularity
+export REQUESTS_CA_BUNDLE=$TMPWORKDIR/ca-certificates.crt  # accessible copy within the container
+mkdir -p $SINGULARITY_TMPDIR
+
+# cmake/make stuff
+export BUILD=$CONTAINER_ROOT/tmp/build
+export INSTALL_DIR=$CONTAINER_ROOT/tmp/install  # path within the container for the sake of wheel link edition
+export CMAKE_INSTALL_PREFIX=$INSTALL_DIR
+
+# end of definitions
+# =============================================================================
+
+# target action
+if [ "$1" == '-h' ]; then
+  usage
+  exit
+else
+  target=$1
+  shift
+  if [ "$target" != "all" ] && [ "$target" != "get" ] && [ "$target" != "deps" ] && \
+     [ "$target" != "lib" ] && [ "$target" != "wheel" ]  && [ "$target" != "pypi" ]; then
+    usage
+    echo
+    echo "First argument must be one of 'all', 'get', 'deps', 'lib', 'wheel', 'pypi'"
+    exit 1
+  fi
+  if [ "$target" == "pypi" ] || [ "$target" == "all" ]; then
+    pypi=$1
+    shift
+    if [ "$pypi" == "" ]; then
+      usage
+      echo
+      echo "The pypi argument is needed to push the wheel"
+      exit 2
+    fi
+  fi
+fi
+
+# ecrad version
+version=$(cat VERSION)
+
+# get container and dependencies
+if [ "$target" == "get" ] || [ "$target" == "all" ]; then
+  # clone dependencies packages (out of container)
+  [ -d $SRC/netcdf ] && rm -rf $SRC/netcdf
+  git clone $netcdf_uri -b $netcdf_version $SRC/netcdf
+
+  [ -d $SRC/netcdffortran ] && rm -rf $SRC/netcdffortran
+  git clone $netcdffortran_uri -b $netcdffortran_version $SRC/netcdffortran
+
+  rm -rf $SRC/hdf5-* $SRC/hdf5.tgz
+  wget $hdf5 -O $SRC/hdf5.tgz
+  (cd $SRC; tar xf hdf5.tgz)
+
+  rm -rf $SRC/curl-* $SRC/curl.tgz
+  wget $curl -O $SRC/curl.tgz
+  (cd $SRC; tar xf curl.tgz; rm -f curl.tgz)
+
+  # get container
+  singularity pull $CONTAINER_SIF_PATH $container_uri
+
+  # copy certificate for pip install from within the container
+  cp $local_certificate $REQUESTS_CA_BUNDLE
+fi
+
+# build dependencies
+if [ "$target" == "deps" ] || [ "$target" == "all" ]; then
+  cat - <<..EOF > ${TMPWORKDIR}/build_deps.sh
+  set -e
+
+  cd $SRC/curl-*
+  ./configure --prefix=$INSTALL_DIR --without-ssl
+  make -j 4
+  make install
+
+  cd $SRC/hdf5-*
+  ./configure --prefix=$INSTALL_DIR
+  make -j 4
+  make install
+
+  rm -rf $BUILD
+  mkdir -p $BUILD
+  cd $BUILD
+  cmake -S $SRC/netcdf -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR
+  cmake --build $BUILD
+  cmake --install $BUILD
+
+  rm -rf $BUILD
+  mkdir -p $BUILD
+  cd $BUILD
+  cmake -S $SRC/netcdffortran -DCMAKE_INSTALL_PREFIX=$INSTALL_DIR
+  cmake --build $BUILD
+  cmake --install $BUILD
+..EOF
+  chmod +x ${TMPWORKDIR}/build_deps.sh
+  singularity run -i $CONTAINER_SIF_PATH ${TMPWORKDIR}/build_deps.sh
+fi
+
+# build ecrad lib
+if [ "$target" == "lib" ] || [ "$target" == "all" ]; then
+  cat - <<..EOF > ${TMPWORKDIR}/build_lib.sh
+  export PATH=\$PATH:$INSTALL_DIR/bin # to find nf-config
+  make python
+..EOF
+  chmod +x ${TMPWORKDIR}/build_lib.sh
+  singularity run -i $CONTAINER_SIF_PATH ${TMPWORKDIR}/build_lib.sh
+fi
+
+# build wheels
+if [ "$target" == "wheel" ] || [ "$target" == "all" ]; then
+  cat - <<..EOF > ${TMPWORKDIR}/build_wheel.sh
+  export PATH=\$PATH:$INSTALL_DIR/bin # to find nf-config
+  export LD_LIBRARY_PATH=$INSTALL_DIR/lib:$INSTALL_DIR/lib64 # for auditwheel repair
+  for p in $python_versions; do
+    /opt/python/\$p/bin/pip install build
+    /opt/python/\$p/bin/pip install auditwheel
+    /opt/python/\$p/bin/python -m build --wheel
+    wheel=dist/pyecrad-${version}-\$p-linux_x86_64.whl
+    if [ -f \$wheel ]; then
+      /opt/python/\$p/bin/python -m auditwheel repair \$wheel
+    fi
+  done
+  wheel=dist/pyecrad-${version}-py3-none-any.whl
+  if [ -f \$wheel ]; then
+    /opt/python/\$p/bin/python -m auditwheel repair \$wheel
+  fi
+..EOF
+  chmod +x ${TMPWORKDIR}/build_wheel.sh
+  singularity run -i $CONTAINER_SIF_PATH ${TMPWORKDIR}/build_wheel.sh
+fi
+
+# distribute wheels
+if [ "$target" == "pypi" ] || [ "$target" == "all" ]; then
+  python3 -m twine upload --repository $pypi wheelhouse/pyecrad-${version}-py3-none-manylinux2014_x86_64.manylinux_2_17_x86_64.whl
+fi
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..15e40201
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,31 @@
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "pyecrad"
+dynamic = ["version"]
+dependencies = ["ctypesForFortran >=1.3.0, !=2.0.*, !=2.1.*",
+                "numpy",
+                "f90nml",
+                "netCDF4",
+               ]
+authors = [
+  {name = "Sébastien Riette", email = "sebastien.riette@meteo.fr" },
+]
+description = "ecRad python binding"
+readme = "README.md"
+requires-python = ">= 3.10"
+license = "Apache-2.0"
+
+[tool.setuptools.dynamic]
+version = {file = "VERSION"}
+
+[tool.setuptools.packages.find]
+include = ["pyecrad*"]
+
+[tool.setuptools.package-data]
+pyecrad = ['libecrad4py.so',
+           'VERSION',
+           'data/*',
+          ]
diff --git a/setup.py b/setup.py
new file mode 100644
index 00000000..f118c946
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,21 @@
+"""
+Build extension for pyecrad
+"""
+
+import subprocess
+from setuptools import setup
+from setuptools.command.build import build
+
+
+class PyecradBuild(build):
+    """
+    Custom class to invoke make python
+    """
+    def run(self):
+        """
+        Method actually doing the build
+        """
+        subprocess.run(['make', 'python'], check=True)
+
+
+setup(cmdclass={"build": PyecradBuild})
diff --git a/test/pyecrad/Makefile b/test/pyecrad/Makefile
new file mode 100644
index 00000000..48e17ef4
--- /dev/null
+++ b/test/pyecrad/Makefile
@@ -0,0 +1,19 @@
+all: compare
+
+compare: control experiment
+	python3 ../common/nccmp.py control.nc experiment.nc
+
+clean:
+	rm -f control.nc experiment.nc
+
+control:
+	ecrad config.nam era5slice.nc control.nc
+
+experiment:
+	# Run ecrad
+	driver.py config.nam era5slice.nc experiment.nc
+	# Add the history attribute for the nccmp.py script
+	python3 -c "import netCDF4; \
+	            nc = netCDF4.Dataset('experiment.nc', 'a'); \
+	            nc.history='pyecrad test'; \
+	            nc.close()"
diff --git a/test/pyecrad/config.nam b/test/pyecrad/config.nam
new file mode 100644
index 00000000..7c8317cd
--- /dev/null
+++ b/test/pyecrad/config.nam
@@ -0,0 +1,84 @@
+! Configuration namelists for ecRad radiation scheme
+!
+! The "radiation_driver" namelist controls the behaviour of the fortran
+! driver routine, including parallelization options and overriding numbers
+! read from the NetCDF input file. Because this namelist is used to
+! run the reference fortran implementation (which uses this namelist)
+! to evaluate the python implementation (which do not use it), every
+! modification done here must be transferred to the python script.
+!
+! The "radiation" namelist controls the behaviour of the radiative
+! transfer algorithm itself, it is used by both versions. Any line
+! prefixed by "!" is ignored. If a namelist parameter is missing then
+! ecRad will use a default value.  For most parameters you can see
+! what ecRad has used from the output printed to the terminal when it
+! runs.
+!
+! The "iverbose*" parameters specify the verbosity level: 0=none,
+! 1=warning, 2=info, 3=progress, 4=detailed, 5=debug
+
+&radiation_driver
+!
+! GENERAL
+!
+iverbose    	= 2, 
+do_parallel     = true,      ! Use OpenMP parallelization, if possible?
+nblocksize      = 32,        ! Number of columns to process per thread
+experiment_name = "Control", ! Written to output file, used in plot legends
+!
+! SCALE OR OVERRIDE ECRAD INPUTS
+!
+fractional_std  = 1, ! Fractional standard dev. of in-cloud water content
+!
+! The following settings configure the SPARTACUS solver
+cloud_separation_scale_toa     = 14000.0,
+cloud_separation_scale_surface = 2500.0,
+cloud_separation_scale_power   = 3.5,
+cloud_inhom_separation_factor  = 0.75,
+!
+
+! Writing fluxes in double precision removes noise from differences in
+! mesospheric heating rates
+do_write_double_precision = true,
+/
+
+&radiation
+!
+! GENERAL
+!
+iverbose            = 1, 
+iverbosesetup       = 1,
+directory_name      = "data",         ! Location of configuration files
+do_surface_sw_spectral_flux = false,  ! Save surface fluxes in each band?
+!
+! CLOUDS
+!
+use_general_cloud_optics = false,
+ice_model_name      = "Fu-IFS",       ! Can be "Fu-IFS" or "Yi"
+sw_solver_name      = "Tripleclouds", ! "Tripleclouds", "McICA" or "SPARTACUS"
+lw_solver_name      = "Tripleclouds", ! "Tripleclouds", "McICA" or "SPARTACUS"
+overlap_scheme_name = "Exp-Ran",      ! McICA also accepts Max-Ran or Exp-Exp
+do_lw_cloud_scattering = true,        ! Clouds scatter in the longwave?
+gas_model_name      = "RRTMG-IFS",    ! "RRTMG-IFS" or "ECCKD"
+!
+! AEROSOLS
+!
+use_aerosols             = true,      ! Radiation sees aerosols?
+use_general_aerosol_optics=false,
+do_lw_aerosol_scattering = false,     ! Aerosols scatter in the longwave?
+!
+! 11 IFS aerosol mixing ratios are stored in the ecRad input file: 1-3
+! Sea salt, 4-6 mineral dust, 7 hydrophilic organics, 8 hydrophobic
+! organics, 9 hydrophilic black carbon, 10 hydrophobic black carbon, 11
+! ammonium sulfate
+n_aerosol_types  = 11,   ! Number of aerosol types in input file
+!
+! The aerosol optical properties are in this file:
+aerosol_optics_override_file_name = 'aerosol_ifs_rrtm_46R1_with_NI_AM.nc'
+!
+! For each of the 11 mixing ratios in the input file, we need to map to
+! one of the optical properties, where negative numbers index
+! hydrophilic aerosol types and positive numbers index hydrophobic
+! aerosol types, e.g. 11=black carbon, -5=sulphate.
+i_aerosol_type_map = -1, -2, -3, 1, 2, 3, -4, 10, 11, 11, -5,
+/
diff --git a/test/pyecrad/data b/test/pyecrad/data
new file mode 120000
index 00000000..e67b4559
--- /dev/null
+++ b/test/pyecrad/data
@@ -0,0 +1 @@
+../../data
\ No newline at end of file
diff --git a/test/pyecrad/driver.py b/test/pyecrad/driver.py
new file mode 100755
index 00000000..65d22e5f
--- /dev/null
+++ b/test/pyecrad/driver.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+"""
+This script mimics the behaviour of the offline driver in order
+to test the python implementation.
+
+(C) Copyright 2026- ECMWF.
+
+This software is licensed under the terms of the Apache Licence Version 2.0
+which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+
+In applying this licence, ECMWF does not waive the privileges and immunities
+granted to it by virtue of its status as an intergovernmental organisation
+nor does it submit to any jurisdiction.
+
+Author:  Sébastien Riette
+Email:   sebastien.riette@meteo.fr
+License: see the COPYING file for details
+"""
+
+import importlib.util
+import sys
+
+# Import pyecrad module
+spec = importlib.util.spec_from_file_location('pyecrad', '../../pyecrad/__init__.py')
+pyecrad = importlib.util.module_from_spec(spec)
+sys.modules['pyecrad'] = pyecrad
+spec.loader.exec_module(pyecrad)
+
+if __name__ == '__main__':
+    import argparse
+    parser = argparse.ArgumentParser(description='Offline driver using pyecrad')
+    parser.add_argument('namelist', help='Namelist file name')
+    parser.add_argument('input', help='Input netCDF file')
+    parser.add_argument('output', help='Output netCDF file')
+    parser.add_argument('--nblocksize', type=int, default=32,
+                        help='Block size')
+    args = parser.parse_args()
+
+    ecrad = pyecrad.Ecrad(args.namelist)
+    ecrad.driver(args.input, args.output, nblocksize=args.nblocksize)
diff --git a/test/pyecrad/ecrad b/test/pyecrad/ecrad
new file mode 120000
index 00000000..68e0999c
--- /dev/null
+++ b/test/pyecrad/ecrad
@@ -0,0 +1 @@
+../../bin/ecrad
\ No newline at end of file
diff --git a/test/pyecrad/era5slice.nc b/test/pyecrad/era5slice.nc
new file mode 120000
index 00000000..1777a735
--- /dev/null
+++ b/test/pyecrad/era5slice.nc
@@ -0,0 +1 @@
+../../practical/era5slice.nc
\ No newline at end of file
diff --git a/test/pyecrad/test_pyecrad.py b/test/pyecrad/test_pyecrad.py
new file mode 100644
index 00000000..b24d6245
--- /dev/null
+++ b/test/pyecrad/test_pyecrad.py
@@ -0,0 +1,29 @@
+"""
+pytest script
+"""
+
+import os
+import subprocess
+import netCDF4
+
+
+def test_pyecrad():
+    """
+    Builds the reference and test file and compare them
+    """
+    cwd = os.path.dirname(os.path.abspath(__file__))
+
+    # reference file
+    subprocess.run(['ecrad', 'config.nam', 'era5slice.nc', 'control.nc'],
+                   cwd=cwd, check=True)
+
+    # test file
+    subprocess.run(['driver.py', 'config.nam', 'era5slice.nc',
+                    'experiment.nc'], cwd=cwd, check=True)
+    with netCDF4.Dataset(os.path.join(cwd, 'experiment.nc'), 'a') as nc:
+        nc.history = 'pyecrad test'
+
+    # comparison
+    assert subprocess.run(['python3', '../common/nccmp.py',
+                           'control.nc', 'experiment.nc'],
+                          cwd=cwd, check=False).returncode == 0