Repository: jouvetg/igm (language: python)

Issue:
Added custom modules, added param_file CLI argument
- tested custom modules for working directory (folder and file)
- tested error messages for custom modules
TODO
- allow custom modules from various paths in CLI as an argument maybe?
- check and fix dependencies
- look into override_json function
- look into parsing help arguments (parsing multiple times and general sub-parser structure)
- implement unit testing


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/install_cuda_ubuntu.sh b/.github/workflows/install_cuda_ubuntu.sh
new file mode 100755
index 0000000..8e9a7b1
--- /dev/null
+++ b/.github/workflows/install_cuda_ubuntu.sh
@@ -0,0 +1,13 @@
+#!/bin/bash
+# Update the system
+sudo apt-get update
+sudo apt-get upgrade -y
+
+# Install dependencies
+sudo apt-get install -y build-essential dkms
+
+# Install CUDA Toolkit
+sudo apt-get install -y nvidia-cuda-toolkit
+
+# Verify installation
+nvidia-smi
diff --git a/.github/workflows/python-publish.yml b/.github/workflows/python-publish.yml
index 37d7cb0..ba334a7 100644
--- a/.github/workflows/python-publish.yml
+++ b/.github/workflows/python-publish.yml
@@ -11,32 +11,53 @@ name: Upload Python Package
 on:
   push:
     branches:
-      - main
+      - reorganize
 
 jobs:
   deploy:
 
-    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        os: [ubuntu-latest, windows-latest, macos-latest]
 
+    runs-on: ${{ matrix.os }}
     steps:
     - uses: actions/checkout@v3
     - name: Set up Python
       uses: actions/setup-python@v3
       with:
             python-version: '<3.11'
+    # Should integration tests be run here or just assume the user has a working environment?
+    # - name: Install CUDA
+    #   run: |
+    #     ./.github/workflows/install_cuda_ubuntu.sh
+    #     nvidia-smi
+    # - uses: Jimver/cuda-toolkit@v0.2.11
+    #   id: cuda-toolkit
+    #   with:
+    #     cuda: '12.1.0'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
         pip install build
     - name: Build package
       run: python -m build
+    - name: Install package
+      run: pip install .
 #    - name: Publish package to testpypi
 #      uses: pypa/gh-action-pypi-publish@release/v1
 #      with:
 #        password: ${{ secrets.PYPI_API_TOKEN }}
 #        repository-url: https://test.pypi.org/legacy/
-    - name: Publish package to pypi
-      uses: pypa/gh-action-pypi-publish@release/v1
-      with:
-        password: ${{ secrets.IGM_TOKEN }}
-        repository-url: https://upload.pypi.org/legacy/
+    - name: Run tests
+      run: |
+        pip install pytest
+        python -c "import os; print(os.getcwd())"
+        cd tests
+        python -m pytest -W ignore::DeprecationWarning -p no:cacheprovider
+  # Not deploying pypi at the moment as I am testing some CI/CD with unit-testing (and its a feature branch anyways...)
+    # - name: Publish package to pypi
+    #   uses: pypa/gh-action-pypi-publish@release/v1
+    #   with:
+    #     password: ${{ secrets.IGM_TOKEN }}
+    #     repository-url: https://upload.pypi.org/legacy/
diff --git a/.gitignore b/.gitignore
index 7e8600c..252ec97 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,7 +1,33 @@
 __pycache__
-igm_model.egg-info
+igm_model.egg-info/
 dist
 sol_analytic
 paper.tex
-fig
+fig/
 biblio.bib
+pix2pixhd-texture-model/
+texture-aif-test
+checkpoints/
+logdir/
+# texture/
+smb_accpdd/
+clim_glacialindex/
+
+
+# hide output files and data
+input_saved.nc
+RGI60-11.01450/
+output.nc
+computational-statistics.txt
+memory-statistics.txt
+computational-pie.png
+memory-pie.png
+examples/params.json
+clean.sh
+params_saved.json
+
+smb_accpdd/
+clim_glacialindex/
+input_saved.nc
+RGI60-11.01450/
+.vscode/
\ No newline at end of file
diff --git a/CHANGES.txt b/CHANGES.txt
deleted file mode 100644
index 139597f..0000000
--- a/CHANGES.txt
+++ /dev/null
@@ -1,2 +0,0 @@
-
-
diff --git a/README.md b/README.md
index 11fe953..1a73841 100644
--- a/README.md
+++ b/README.md
@@ -28,5 +28,4 @@ IGM has a [discord channel](https://discord.com/) for quick support, getting lat
 
 # Contact
 
-Feel free to drop me an email for any questions, bug reports, or ideas of model extension: guillaume.jouvet at unil.ch
-
+Feel free to drop me an email for any questions, bug reports, or ideas of model extension: guillaume.jouvet at unil.ch
\ No newline at end of file
diff --git a/doc/anim_mayavi.md b/doc/anim_mayavi.md
index f605097..8a70abf 100644
--- a/doc/anim_mayavi.md
+++ b/doc/anim_mayavi.md
@@ -18,3 +18,10 @@ pip install mayavi pyqt5
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/anim_plotly.md b/doc/anim_plotly.md
index c900899..d3b4119 100644
--- a/doc/anim_plotly.md
+++ b/doc/anim_plotly.md
@@ -28,3 +28,10 @@ This module was implemented by [Oskar Herrmann](https://github.com/ho11laqe).
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/anim_video.md b/doc/anim_video.md
index 6c4046b..b50da41 100644
--- a/doc/anim_video.md
+++ b/doc/anim_video.md
@@ -13,3 +13,10 @@ This module depends on `xarray` library.
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/avalanche.md b/doc/avalanche.md
index 420b329..c073937 100644
--- a/doc/avalanche.md
+++ b/doc/avalanche.md
@@ -15,5 +15,12 @@ program and implemented in IGM by Jürgen Mey with support from Guillaume Jouvet
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--avalanche_update_freq`|`1`|Update avalanche each X years (1)|
 ||`--avalanche_angleOfRepose`|`30`|Angle of repose (30°)|
diff --git a/doc/clim_oggm.md b/doc/clim_oggm.md
index 403b6d3..1176afa 100644
--- a/doc/clim_oggm.md
+++ b/doc/clim_oggm.md
@@ -32,6 +32,13 @@ If parameter `clim_oggm_clim_trend_array` is set to empty list `[]`, then it wil
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--clim_oggm_update_freq`|`1`|Update the climate each X years|
 ||`--smb_oggm_file`|`smb_oggm_param.txt`|Name of the imput file for the climate outide the given datatime frame (time, delta_temp, prec_scali)|
 ||`--clim_oggm_clim_trend_array`|`[['time', 'delta_temp', 'prec_scal'], [1900, 0.0, 1.0], [2020, 0.0, 1.0]]`|Define climate trend outside available time window|
diff --git a/doc/enthalpy.md b/doc/enthalpy.md
index ea4bbec..be0cd10 100644
--- a/doc/enthalpy.md
+++ b/doc/enthalpy.md
@@ -186,6 +186,13 @@ Here one updates the enthalpy as many times as the ice flow, we assume that the
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--enth_water_density`|`1000`|Constant of the Water density [kg m-3]|
 ||`--enth_spy`|`31556926`|Number of seconds per years [s y-1]|
 ||`--enth_ki`|`2.1`|Conductivity of cold ice [W m-1 K-1] (Aschwanden and al, JOG, 2012)|
diff --git a/doc/flow_dt_thk.md b/doc/flow_dt_thk.md
index 4373c79..cc9c97b 100644
--- a/doc/flow_dt_thk.md
+++ b/doc/flow_dt_thk.md
@@ -11,6 +11,13 @@ This IGM module wraps up module iceflow, time and thk (check for each modules fo
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--iflo_type`|`emulated`|Type of iceflow: it can emulated (default), solved, or in diagnostic mode to investigate the fidelity of the emulator towads the solver|
 ||`--iflo_pretrained_emulator`||Do we take a pretrained emulator or start from scratch?|
 ||`--iflo_emulator`|``|Directory path of the deep-learning pretrained ice flow model, take from the library if empty string|
diff --git a/doc/generate_all_params_doc1.py b/doc/generate_all_params_doc1.py
index a89543e..7068eab 100644
--- a/doc/generate_all_params_doc1.py
+++ b/doc/generate_all_params_doc1.py
@@ -5,10 +5,11 @@
 
 for m in ['preproc','process','postproc']:
     
-    L = glob.glob('../igm/modules/'+m+'/*.py')
+    L = glob.glob('../igm/modules/'+m+'/*/*.py')
     L = [l.split('/')[-1].split('.')[0] for l in L]
-    L.remove('__init__')
-    
+    while '__init__' in L:
+        L.remove('__init__')
+     
     for l in L:
 
         print(m,l)
@@ -17,19 +18,20 @@
 
             print("import sys", file=f)
             print("sys.path.append('../')", file=f)
-            print("import argparse", file=f)
-            print("import igm", file=f)
-            print("parser = argparse.ArgumentParser()", file=f)
-            print("igm.params_"+l+"(parser)", file=f)
+            print("import igm", file=f)            
+            print("parser = igm.params_core()", file=f)
+            print("params, __ = parser.parse_known_args()", file=f)
+            print('params.modules_'+m+'=["'+l+'"]', file=f)
+            print("imported_modules = igm.load_modules(params)", file=f)
+            print("igm.modules."+m+"."+l+".params(parser)", file=f)
             print("params = parser.parse_args()", file=f)
-            
-         
+
         os.system("argmark -f igm-run.py")
         os.system("echo ' ' > argmark1.md")
         os.system("echo '# Parameters: ' >> argmark1.md")
         os.system("sed -z 's/.*Arguments//' argmark.md >> argmark1.md")
 
-        os.system('cat ../igm/modules/'+m+'/'+l+'.md argmark1.md > ' + l + '.md')
+        os.system('cat ../igm/modules/'+m+'/'+l+'/'+l+'.md argmark1.md > ' + l + '.md')
 
 os.system('rm igm-run.py argmark.md argmark1.md')
  
diff --git a/doc/gflex.md b/doc/gflex.md
index 7238261..47b6a7c 100644
--- a/doc/gflex.md
+++ b/doc/gflex.md
@@ -23,5 +23,12 @@ On the other hand, we do not expect a frequent update, therefore, this module sh
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--gflex_update_freq`|`100.0`|Update gflex each X years (1)|
 ||`--gflex_default_Te`|`50000`|Default value for Te (Elastic thickness [m]) if not given as ncdf file|
diff --git a/doc/glerosion.md b/doc/glerosion.md
index d2e440c..f1b5fa4 100644
--- a/doc/glerosion.md
+++ b/doc/glerosion.md
@@ -16,6 +16,13 @@ By default, we use the parameters from
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--glerosion_cst`|`2.7e-07`|Erosion multiplicative factor, here taken from Herman, F. et al.               Erosion by an Alpine glacier. Science 350, 193–195 (2015)|
 ||`--glerosion_exp`|`2`|Erosion exponent factor, here taken from Herman, F. et al.                Erosion by an Alpine glacier. Science 350, 193–195 (2015)|
 ||`--glerosion_update_freq`|`1`|Update the erosion only each X years (Default: 100)|
diff --git a/doc/iceflow.md b/doc/iceflow.md
index e485ee0..10de11e 100644
--- a/doc/iceflow.md
+++ b/doc/iceflow.md
@@ -61,6 +61,13 @@ When treating ery large arrays, retraining must be done sequentially patch-wise
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--iflo_type`|`emulated`|Type of iceflow: it can emulated (default), solved, or in diagnostic mode to investigate the fidelity of the emulator towads the solver|
 ||`--iflo_pretrained_emulator`||Do we take a pretrained emulator or start from scratch?|
 ||`--iflo_emulator`|``|Directory path of the deep-learning pretrained ice flow model, take from the library if empty string|
diff --git a/doc/iceflow_v1.md b/doc/iceflow_v1.md
index 9375eb3..87538b4 100644
--- a/doc/iceflow_v1.md
+++ b/doc/iceflow_v1.md
@@ -37,6 +37,13 @@ as output: (ubar,vbar, uvelsurf, vvelsurf, uvelbase, vvelbase)
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--init_strflowctrl`|`78`|Initial strflowctrl (default 78)|
 ||`--iflo_emulator`|`f15_cfsflow_GJ_22_a`|Directory path of the deep-learning ice flow model,               create a new if empty string|
 ||`--iflo_init_slidingco`|`0`|Initial sliding coeeficient slidingco (default: 0)|
diff --git a/doc/include_icemask.md b/doc/include_icemask.md
index c13d3e3..d85df4a 100644
--- a/doc/include_icemask.md
+++ b/doc/include_icemask.md
@@ -32,5 +32,12 @@ Author: Andreas Henz, andreas.henz@geo.uzh.ch  (06.09.2023)
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--mask_shapefile`|`icemask.shp`|Icemask input file (default: icemask.shp)|
 ||`--mask_invert`||Invert ice mask if the mask is where the ice should be (default: False)|
diff --git a/doc/load_ncdf.md b/doc/load_ncdf.md
index 6c0149c..73892de 100644
--- a/doc/load_ncdf.md
+++ b/doc/load_ncdf.md
@@ -14,6 +14,13 @@ This module depends on `netCDF4`.
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--lncd_input_file`|`input.nc`|NetCDF input data file|
 ||`--lncd_coarsen`|`1`|Coarsen the data from NetCDF file by a certain (integer) number: 2 would be twice coarser ignore data each 2 grid points|
 ||`--lncd_crop`|`False`|Crop the data from NetCDF file with given top/down/left/right bounds|
diff --git a/doc/load_tif.md b/doc/load_tif.md
index f81be58..5fecb9c 100644
--- a/doc/load_tif.md
+++ b/doc/load_tif.md
@@ -12,6 +12,13 @@ This module depends on `rasterio`.
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--ltif_coarsen`|`1`|coarsen the data to a coarser resolution (default: 1), e.g. 2 would be twice coarser ignore data each 2 grid points|
 ||`--ltif_crop`|`False`|Crop the data with xmin, xmax, ymin, ymax (default: False)|
 ||`--ltif_xmin`|`-100000000000000000000`|crop_xmin|
diff --git a/doc/oggm_shop.md b/doc/oggm_shop.md
index 29963be..4d66ec5 100644
--- a/doc/oggm_shop.md
+++ b/doc/oggm_shop.md
@@ -21,6 +21,13 @@ The module depends (of course) of `oggm` library. Unfortunatly the module works
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--oggm_RGI_ID`|`RGI60-11.01450`|RGI ID|
 ||`--oggm_preprocess`||Use preprocessing|
 ||`--oggm_dx`|`100`|Spatial resolution (need preprocess false to change it)|
diff --git a/doc/optimize.md b/doc/optimize.md
index ea6aa01..ccd7b06 100644
--- a/doc/optimize.md
+++ b/doc/optimize.md
@@ -139,6 +139,13 @@ You may monitor the data assimilation during the inverse modelling in several wa
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--opti_vars_to_save`|`['usurf', 'thk', 'slidingco', 'velsurf_mag', 'velsurfobs_mag', 'divflux', 'icemask']`|List of variables to be recorded in the ncdef file|
 ||`--opti_init_zero_thk`|`False`|Initialize the optimization with zero ice thickness|
 ||`--opti_regu_param_thk`|`10.0`|Regularization weight for the ice thickness in the optimization|
diff --git a/doc/optimize_v1.md b/doc/optimize_v1.md
index 00d4236..32bdb67 100644
--- a/doc/optimize_v1.md
+++ b/doc/optimize_v1.md
@@ -127,6 +127,13 @@ You may monitor the data assimilation during the inverse modelling in several wa
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--opti_vars_to_save`|`['usurf', 'thk', 'strflowctrl', 'arrhenius', 'slidingco', 'velsurf_mag', 'velsurfobs_mag', 'divflux']`|List of variables to be recorded in the ncdef file|
 ||`--opti_thr_strflowctrl`|`78.0`|threshold value for strflowctrl|
 ||`--opti_init_zero_thk`|`False`|Initialize the optimization with zero ice thickness|
diff --git a/doc/particles.md b/doc/particles.md
index a2d0749..1d26e21 100644
--- a/doc/particles.md
+++ b/doc/particles.md
@@ -25,6 +25,13 @@ The module needs horizontal velocities (state.U), as well as vertical speeds (st
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--part_tracking_method`|`simple`|Method for tracking particles (simple or 3d)|
 ||`--part_frequency_seeding`|`50`|Frequency of seeding (unit : year)|
 ||`--part_density_seeding`|`0.2`|Density of seeding (1 means we seed all pixels, 0.2 means we seed each 5 grid cell, ect.)|
diff --git a/doc/particles_v1.md b/doc/particles_v1.md
index abd493c..8eb095c 100644
--- a/doc/particles_v1.md
+++ b/doc/particles_v1.md
@@ -12,6 +12,13 @@ associated with iceflow_v1 (check at the doc of particles).
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--part_tracking_method`|`3d`|Method for tracking particles (3d or simple)|
 ||`--part_frequency_seeding`|`10`|Frequency of seeding (default: 10)|
 ||`--part_density_seeding`|`0.2`|Density of seeding (default: 0.2)|
diff --git a/doc/plot2d.md b/doc/plot2d.md
index d9d31d9..023d0bf 100644
--- a/doc/plot2d.md
+++ b/doc/plot2d.md
@@ -14,6 +14,13 @@ If the `particles` module is activated, one may plot particles on the top settin
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--plt2d_editor`|`vs`|Optimized for VS code (vs) or spyder (sp) for live plot|
 ||`--plt2d_live`||Display plots live the results during computation instead of making png|
 ||`--plt2d_particles`||Display particles is True, does not display if False|
diff --git a/doc/pretraining.md b/doc/pretraining.md
index 6146a15..ff125ce 100644
--- a/doc/pretraining.md
+++ b/doc/pretraining.md
@@ -46,6 +46,13 @@ The path (or name of the data folder) must be pass in parameter `data_dir`.
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--data_dir`|`surflib3d_shape_100`|Directory of the data of the glacier catalogu|
 ||`--batch_size`|`1`|Batch size|
 ||`--freq_test`|`20`|Frequence of the test|
diff --git a/doc/print_comp.md b/doc/print_comp.md
index 0492da8..43b51da 100644
--- a/doc/print_comp.md
+++ b/doc/print_comp.md
@@ -12,3 +12,10 @@ Note: These numbers must be interepreted with care: Leaks of computational times
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/print_info.md b/doc/print_info.md
index 0c9981c..032e49d 100644
--- a/doc/print_info.md
+++ b/doc/print_info.md
@@ -10,3 +10,10 @@ This IGM modules prints in the terminal basic live-time modelling informations s
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/rockflow.md b/doc/rockflow.md
index b92bf14..dde88a3 100644
--- a/doc/rockflow.md
+++ b/doc/rockflow.md
@@ -11,4 +11,11 @@ This module extends the ice flow outside the glaciated area, by giving a constan
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--rock_flow_speed`|`1`|Speed of rock flow along the slope in m/y|
diff --git a/doc/smb_oggm.md b/doc/smb_oggm.md
index d57e728..c9117c9 100644
--- a/doc/smb_oggm.md
+++ b/doc/smb_oggm.md
@@ -13,6 +13,13 @@ where $P_i^{sol}$ is the is the monthly solid precipitation, $T_i$ is the monthl
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--smb_oggm_update_freq`|`1`|Update the mass balance each X years |
 ||`--smb_oggm_ice_density`|`910.0`|Density of ice for conversion of SMB into ice equivalent|
 ||`--smb_oggm_wat_density`|`1000.0`|Density of water|
diff --git a/doc/smb_simple.md b/doc/smb_simple.md
index 22dc202..b267931 100644
--- a/doc/smb_simple.md
+++ b/doc/smb_simple.md
@@ -39,6 +39,13 @@ If one has provided in input an "icemask" field, then this module will compute n
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--smb_simple_update_freq`|`1`|Update the mass balance each X years (1)|
 ||`--smb_simple_file`|`smb_simple_param.txt`|Name of the imput file for the simple mass balance model (time, gradabl, gradacc, ela, accmax)|
 ||`--smb_simple_array`|`[]`|Time dependent parameters for simple mass balance model (time, gradabl, gradacc, ela, accmax)|
diff --git a/doc/thk.md b/doc/thk.md
index 61099e2..a65a8bc 100644
--- a/doc/thk.md
+++ b/doc/thk.md
@@ -28,3 +28,10 @@ This IGM module solves the mass conservation of ice to update the thickness from
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/time.md b/doc/time.md
index 7799f09..795baf1 100644
--- a/doc/time.md
+++ b/doc/time.md
@@ -33,6 +33,13 @@ A bit more details on the time step stability conditionsis given in the followin
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--time_start`|`2000.0`|Start modelling time|
 ||`--time_end`|`2100.0`|End modelling time|
 ||`--time_save`|`10`|Save result frequency for many modules (in year)|
diff --git a/doc/vert_flow.md b/doc/vert_flow.md
index 1b52dd0..95315dc 100644
--- a/doc/vert_flow.md
+++ b/doc/vert_flow.md
@@ -11,3 +11,10 @@ This IGM module computes the vertical component (providing state.W) of the veloc
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
diff --git a/doc/write_ncdf.md b/doc/write_ncdf.md
index 21f8669..9334ba4 100644
--- a/doc/write_ncdf.md
+++ b/doc/write_ncdf.md
@@ -11,5 +11,12 @@ This IGM module writes 2D field variables defined in the paramer list `wncd_vars
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--wncd_output_file`|`output.nc`|Output ncdf data file|
 ||`--wncd_vars_to_save`|`['topg', 'usurf', 'thk', 'smb', 'velbar_mag', 'velsurf_mag', 'uvelsurf', 'vvelsurf', 'wvelsurf']`|List of variables to be recorded in the ncdf file|
diff --git a/doc/write_particles.md b/doc/write_particles.md
index 8f0895b..ebaba26 100644
--- a/doc/write_particles.md
+++ b/doc/write_particles.md
@@ -21,4 +21,11 @@ providing in turn the particle ID, x,y,z positions, the relative height within t
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--wpar_add_topography`||Add topg|
diff --git a/doc/write_tif.md b/doc/write_tif.md
index 735275c..c0603ed 100644
--- a/doc/write_tif.md
+++ b/doc/write_tif.md
@@ -15,4 +15,11 @@ This module depends on the `rasterio` library.
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--wtif_vars_to_save`|`['usurf', 'thk']`|List of variables to be recorded in the NetCDF file|
diff --git a/doc/write_ts.md b/doc/write_ts.md
index 3a09049..d9dd23e 100644
--- a/doc/write_ts.md
+++ b/doc/write_ts.md
@@ -10,4 +10,11 @@ This IGM module writes time serie variables (ice glaciated area and volume) into
 |short|long|default|help|
 | :--- | :--- | :--- | :--- |
 |`-h`|`--help`||show this help message and exit|
+||`--working_dir`|``|Working directory (default empty string)|
+||`--modules_preproc`|`['oggm_shop']`|List of pre-processing modules|
+||`--modules_process`|`['iceflow', 'time', 'thk']`|List of processing modules|
+||`--modules_postproc`|`['write_ncdf', 'plot2d', 'print_info']`|List of post-processing modules|
+||`--logging`||Activate the looging|
+||`--logging_file`|``|Logging file name, if empty it prints in the screen|
+||`--print_params`||Print definitive parameters in a file for record|
 ||`--wts_output_file`|`output_ts.nc`|Output ncdf data file (time serie)|
diff --git a/examples/aletsch-1880-2100/clim_aletsch.py b/examples/aletsch-1880-2100/clim_aletsch.py
index 0d05211..b252362 100755
--- a/examples/aletsch-1880-2100/clim_aletsch.py
+++ b/examples/aletsch-1880-2100/clim_aletsch.py
@@ -9,7 +9,7 @@
 import tensorflow as tf
 import time
 
-def params_clim_aletsch(parser):
+def params(parser):
 
     # CLIMATE PARAMETERS
     parser.add_argument(
@@ -25,7 +25,7 @@ def params_clim_aletsch(parser):
         help="Give the resolution the climate forcing should be (monthly=12, daily=365)",
     ) 
 
-def  initialize_clim_aletsch(params,state):
+def  initialize(params,state):
     """
         load climate data to run the Aletsch Glacier simulation
     """
@@ -35,7 +35,7 @@ def  initialize_clim_aletsch(params,state):
 
     # read temperature and precipitation data from temp_prec.dat
     temp_prec = np.loadtxt(
-        os.path.join(params.working_dir, "temp_prec.dat"),
+        "temp_prec.dat",
         dtype=np.float32,
         skiprows=2,
     )
@@ -79,7 +79,7 @@ def  initialize_clim_aletsch(params,state):
     state.tcomp_clim_aletsch = []
 
 
-def update_clim_aletsch(params,state):
+def update(params,state):
 
     if ((state.t - state.tlast_clim_aletsch) >= params.clim_update_freq):
 
@@ -124,6 +124,6 @@ def update_clim_aletsch(params,state):
         state.tcomp_clim_aletsch[-1] *= -1
 
 
-def  finalize_clim_aletsch(params,state):
+def  finalize(params,state):
     pass
  
diff --git a/examples/aletsch-1880-2100/smb_accmelt.py b/examples/aletsch-1880-2100/smb_accmelt.py
index 4de5226..5972675 100755
--- a/examples/aletsch-1880-2100/smb_accmelt.py
+++ b/examples/aletsch-1880-2100/smb_accmelt.py
@@ -29,7 +29,7 @@
 import time
 
  ## add custumized smb function
-def params_smb_accmelt(parser):  
+def params(parser):  
 
     parser.add_argument(
         "--smb_simple_update_freq",
@@ -83,6 +83,7 @@ def params_smb_accmelt(parser):
         type=float,
         default=1.0, 
     )
+    
     parser.add_argument(
         "--smb_accmelt_ice_density",
         type=float,
@@ -95,15 +96,14 @@ def params_smb_accmelt(parser):
         default=1000.0,
         help="Density of water",
     )
-    
-def initialize_smb_accmelt(params,state):
+
+
+def initialize(params,state):
     """
         load smb data to run the Aletsch Glacier simulation 
     """
 
-    nc = Dataset(
-        os.path.join(params.working_dir, 'massbalance.nc'), "r"
-    )
+    nc = Dataset('massbalance.nc', "r")
     x = np.squeeze(nc.variables["x"]).astype("float32")
     y = np.squeeze(nc.variables["y"]).astype("float32")
 
@@ -115,7 +115,7 @@ def initialize_smb_accmelt(params,state):
     nc.close()
     
     
-    nc = Dataset( os.path.join(params.working_dir, 'bassin.nc'), "r" )
+    nc = Dataset('bassin.nc', "r" )
     x = np.squeeze(nc.variables["x"]).astype("float32")
     y = np.squeeze(nc.variables["y"]).astype("float32")
 
@@ -144,7 +144,7 @@ def initialize_smb_accmelt(params,state):
     # read mass balance parameters
     state.mb_parameters = tf.Variable(
         np.loadtxt(
-            os.path.join(params.working_dir, "mbparameter.dat"),
+            "mbparameter.dat",
             skiprows=2,
             dtype=np.float32,
         )
@@ -176,7 +176,7 @@ def initialize_smb_accmelt(params,state):
 # Note that tf.function works best with TensorFlow ops; NumPy and Python calls are converted to constants.
 # Therefore: you must make sure any variables are TensorFlow Tensor (and not Numpy)
 # @tf.function()
-def update_smb_accmelt(params,state):
+def update(params,state):
 
     if ((state.t - state.tlast_smb) >= params.smb_simple_update_freq):
 
@@ -257,5 +257,5 @@ def update_smb_accmelt(params,state):
         state.tcomp_smb[-1] *= -1
 
 
-def finalize_smb_accmelt(params,state):
+def finalize(params,state):
     pass
diff --git a/examples/aletsch-1880-2100/track_usurf_obs.py b/examples/aletsch-1880-2100/track_usurf_obs.py
index c1626b6..303a874 100755
--- a/examples/aletsch-1880-2100/track_usurf_obs.py
+++ b/examples/aletsch-1880-2100/track_usurf_obs.py
@@ -22,12 +22,12 @@
 from netCDF4 import Dataset
 
 
-def params_track_usurf_obs(parser):
+def params(parser):
     pass
 
-def initialize_track_usurf_obs(params,state):
+def initialize(params,state):
     
-    nc = Dataset( os.path.join(params.working_dir, 'past_surf.nc'), "r" )
+    nc = Dataset('past_surf.nc', "r" )
 #    for y in [1880,1926,1957,1980,1999,2009,2017]:
 #        vars(state)['surf_'+str(y)] = np.squeeze( nc.variables['surf_'+str(y)] ).astype("float32") 
     for v in nc.variables:
@@ -38,9 +38,9 @@ def initialize_track_usurf_obs(params,state):
  
     # load the surface toporgaphy available at given year
     state.usurf = vars(state)['surf_'+str(int(params.time_start))]
-    state.thk   = state.usurf-state.topg
+    state.thk   = state.usurf -state.topg
 
-def update_track_usurf_obs(params,state):
+def update(params,state):
 
     if state.t in [1880,1926,1957,1980,1999,2009,2017]:
 
@@ -53,5 +53,5 @@ def update_track_usurf_obs(params,state):
                 % (state.t, mean, std, vol) )
 
 
-def  finalize_track_usurf_obs(params,state):
+def finalize(params,state):
     pass
diff --git a/examples/paleo-alps/smb_signal.py b/examples/paleo-alps/smb_signal.py
index 36808e9..a921648 100755
--- a/examples/paleo-alps/smb_signal.py
+++ b/examples/paleo-alps/smb_signal.py
@@ -7,7 +7,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-def params_smb_signal(parser):
+def params(parser):
 
     parser.add_argument(
         "--pdela",
@@ -40,7 +40,7 @@ def params_smb_signal(parser):
         help="Maximum accumulation",
     )
 
-def initialize_smb_signal(params,state):
+def initialize(params,state):
     """
         Retrieve the Temperature difference from the EPICA signal
     """
@@ -50,7 +50,7 @@ def initialize_smb_signal(params,state):
     dT   = ss[:,3]          # extract the dT, i.e. global temp. difference
     state.dT =  interp1d(time,dT, fill_value=(dT[0], dT[-1]),bounds_error=False)
 
-def update_smb_signal(params,state):
+def update(params,state):
     """
         mass balance 'signal'
     """
@@ -65,5 +65,5 @@ def update_smb_signal(params,state):
     state.smb  *= tf.where(tf.less(state.smb , 0), params.gradabl, params.gradacc)
     state.smb  = tf.clip_by_value(state.smb , -100, params.maxacc)
 
-def finalize_smb_signal(params, state):
+def finalize(params, state):
     pass
diff --git a/examples/quick-demo-mysmb/mysmb.py b/examples/quick-demo-mysmb/mysmb.py
index e690f6a..8681b83 100755
--- a/examples/quick-demo-mysmb/mysmb.py
+++ b/examples/quick-demo-mysmb/mysmb.py
@@ -5,13 +5,13 @@
 import math
 
  ## add custumized smb function
-def params_mysmb(parser):  
+def params(parser):  
     parser.add_argument("--meanela", type=float, default=3000 )
 
-def initialize_mysmb(params,state):
+def initialize(params,state):
     params.meanela = np.quantile(state.usurf[state.thk>10],0.2)
 
-def update_mysmb(params,state):
+def update(params,state):
     # perturabe the ELA with sinusional signal 
     ELA = ( params.meanela - 750*math.sin((state.t/100)*math.pi) )
     # compute smb linear with elevation with 2 acc & abl gradients
@@ -23,5 +23,5 @@ def update_mysmb(params,state):
     if hasattr(state, "icemask"):
         state.smb  = tf.where((state.smb<0)|(state.icemask>0.5),state.smb,-10)
     
-def finalize_mysmb(params,state):
+def finalize(params,state):
     pass
diff --git a/examples/synthetic/make_synthetic.py b/examples/synthetic/make_synthetic.py
index d128802..25fa4ac 100755
--- a/examples/synthetic/make_synthetic.py
+++ b/examples/synthetic/make_synthetic.py
@@ -10,11 +10,11 @@
 
 from igm.modules.utils import complete_data
 
-def params_make_synthetic(parser):
+def params(parser):
     pass
 
 
-def initialize_make_synthetic(params, state):
+def initialize(params, state):
 
     x = np.arange(0, 100) * 100  # make x-axis, lenght 10 km, resolution 100 m
     y = np.arange(0, 200) * 100  # make y-axis, lenght 20 km, resolution 100 m
@@ -33,10 +33,10 @@ def initialize_make_synthetic(params, state):
     complete_data(state)
 
 
-def update_make_synthetic(params, state):
+def update(params, state):
     pass
 
 
-def finalize_make_synthetic(params, state):
+def finalize(params, state):
     pass
     
diff --git a/igm/__init__.py b/igm/__init__.py
index 1e91dba..dcb8f00 100755
--- a/igm/__init__.py
+++ b/igm/__init__.py
@@ -3,44 +3,32 @@
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
 
-from igm.common import *
+import sys
+import os
+sys.path.append(os.getcwd()) # i guess by default sys.path does not have cwd. not sure if this is the convention though...
+    
+from . import modules
+# from .modules import (
+# 	preproc,
+# 	process,
+# 	postproc,
+    
+# )
 
-from igm.modules.utils import *
+from .common import (
+    State,
+    params_core,
+    load_modules,
+    add_logger,
+    print_params,
+    load_dependent_modules,
+    get_modules_list,
+    load_user_defined_params,
+    run_intializers,
+    run_processes,
+    run_finalizers,
+    setup_igm,
+    add_logger
+)
 
-from igm.modules.preproc.oggm_shop import *
-from igm.modules.preproc.load_ncdf import *
-from igm.modules.preproc.load_tif import *
-from igm.modules.preproc.optimize_v1 import *
-from igm.modules.preproc.optimize import *
-from igm.modules.preproc.include_icemask import *
-from igm.modules.preproc.pretraining import *
-from igm.modules.preproc.infersmb import *
-
-from igm.modules.process.iceflow_v1 import *
-from igm.modules.process.iceflow import *
-from igm.modules.process.vert_flow import *
-from igm.modules.process.smb_simple import *
-from igm.modules.process.thk import *
-from igm.modules.process.time import *
-from igm.modules.process.particles_v1 import *
-from igm.modules.process.particles import *
-from igm.modules.process.glerosion import *
-from igm.modules.process.rockflow import *
-from igm.modules.process.flow_dt_thk import *
-from igm.modules.process.avalanche import *
-from igm.modules.process.gflex import *
-from igm.modules.process.enthalpy import *
-from igm.modules.process.clim_oggm import *
-from igm.modules.process.smb_oggm import *
-from igm.modules.process.read_output import *
-
-from igm.modules.postproc.print_info import *
-from igm.modules.postproc.print_comp import *
-from igm.modules.postproc.write_tif import *
-from igm.modules.postproc.write_ncdf import *
-from igm.modules.postproc.write_ts import *
-from igm.modules.postproc.plot2d  import *
-from igm.modules.postproc.write_particles import *
-from igm.modules.postproc.anim_mayavi import *
-from igm.modules.postproc.anim_video import *
-from igm.modules.postproc.anim_plotly import *
\ No newline at end of file
+from . import emulators
diff --git a/igm/common.py b/igm/common.py
index 90ebe1a..594e31d 100755
--- a/igm/common.py
+++ b/igm/common.py
@@ -5,28 +5,35 @@
 Published under the GNU GPL (Version 3), check at the LICENSE file
 """
 
-import os, glob, json, sys, inspect, re
+import os, json
+from json import JSONDecodeError
 import importlib
-import argparse
-from igm.modules.utils import *
+from argparse import ArgumentParser, Namespace
 import igm
+import tensorflow as tf
+from pathlib import Path
+from functools import partial
+from typing import List, Any, Dict, Tuple
+from types import ModuleType
+import logging
+import warnings
 
 
-# The (empty) class state serves in fact to use a dictionnary, but it looks like a class
 class State:
-    def __init__(self):
-        self.__dict__ = dict()
+    pass
 
 
 # this create core parameters for any IGM run
 def params_core():
-    parser = argparse.ArgumentParser(description="IGM")
+    parser = ArgumentParser(
+        description="IGM", conflict_handler="resolve"
+    )  # automatically overrides repeated/older parameters! Vaid solution?
 
     parser.add_argument(
-        "--working_dir",
+        "--param_file",
         type=str,
-        default="",
-        help="Working directory (default empty string)",
+        default="params.json",
+        help="Path for the JSON parameter file.",
     )
     parser.add_argument(
         "--modules_preproc",
@@ -48,9 +55,15 @@ def params_core():
     )
     parser.add_argument(
         "--logging",
-        type=str2bool,
+        action="store_true",
         default=False,
-        help="Activate the looging",
+        help="Activates the logging (default value is False)",
+    )
+    parser.add_argument(
+        "--logging_level",
+        default=30,
+        type=int,
+        help="Determine logging level used for logger",
     )
     parser.add_argument(
         "--logging_file",
@@ -60,132 +73,289 @@ def params_core():
     )
     parser.add_argument(
         "--print_params",
-        type=str,
+        action="store_false",
         default=True,
         help="Print definitive parameters in a file for record",
     )
     parser.add_argument(
-        "--gpu",
+        "--gpu_info",
+        action="store_true",
+        default=False,
+        help="Print CUDA and GPU information to the screen",
+    )
+    parser.add_argument(
+        "--gpu_id",
         type=int,
         default=0,
-        help="Id of the GPU to use (default 0)",
-        
+        help="Id of the GPU to use (default is 0)",
+    )
+    parser.add_argument(
+        "--saved_params_filename",
+        type=str,
+        default="params_saved",
+        help="Name of the file to store the parameters used in the run (current working directory)",
     )
-    return parser
-
-
-# Function to remove comments from a JSON string
-def remove_comments(json_str):
-    lines = json_str.split("\n")
-    cleaned_lines = [line for line in lines if not line.strip().startswith("#")]
-    return "\n".join(cleaned_lines)
-
-
-def overide_from_json_file(parser, check_if_params_exist=True):
-    # get the path of the json file
-    param_file = os.path.join(parser.parse_args(args=[]).working_dir, "params.json")
-
-    # load the given parameters from the json file
-    with open(param_file, "r") as json_file:
-        json_text = json_file.read()
-
-    # Remove comments from the JSON string
-    json_without_comments = remove_comments(json_text)
 
-    # Parse the modified JSON string
-    try:
-        dic_params = json.loads(json_without_comments)
-    except json.JSONDecodeError as e:
-        print(f"Error decoding JSON: {e}")
+    return parser
 
-    # list only the parameters registered so far
-    LIST = list(vars(parser.parse_args(args=[])).keys())
 
-    if "time_step" in dic_params["modules_process"]:
-        import sys
+def setup_igm(
+    state: State, parser: ArgumentParser
+) -> Tuple[List[ModuleType], Namespace, State]:
+    params, _ = parser.parse_known_args()
 
-        print(
-            " ------- CHECK THE WIKI TO CHANGE YOUR PARAM FILE AND USER MODULES -------------"
-        )
-        sys.exit(
-            "WARNING: the parameter time_step is deprecated, UPDATE you params.json file"
-        )
+    if params.gpu_info:
+        print_gpu_info()
 
-    if check_if_params_exist:
-        for key in dic_params.keys():
-            if not key in LIST:
-                print(
-                    "WARNING: the following parameters of the json file do not exist in igm: ",
-                    key,
-                )
+    if params.logging:
+        add_logger(params=params, state=state)
+        tf.get_logger().setLevel(params.logging_level)
+    # else:
+    #     os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3" # this is a little dangerous as the default value (for everyone) is to hide info, warnings, etc! Commenting it out for now
 
-    # keep only the parameters to overide hat were registerd so far
-    filtered_dict = {key: value for key, value in dic_params.items() if key in LIST}
+    modules_dict = get_modules_list(params.param_file)
+    imported_modules = load_modules(modules_dict)
+    imported_modules = load_dependent_modules(imported_modules)
 
-    parser.set_defaults(**filtered_dict)
+    for module in imported_modules:
+        module.params(parser)
 
+    core_and_module_params = parser.parse_args()
+    params = load_user_defined_params(
+        param_file=core_and_module_params.param_file,
+        params_dict=vars(core_and_module_params),
+    )
 
-# this add a logger to the state
-def load_custom_module(params, module):
-    if os.path.exists(os.path.join(params.working_dir, module + ".py")):
-        sys.path.append(params.working_dir)
-        custmod = importlib.import_module(module)
+    parser.set_defaults(**params)
+    params = parser.parse_args()
 
-        def is_user_defined_function(obj):
-            return inspect.isfunction(obj) and inspect.getmodule(obj) == custmod
+    return imported_modules, params, state
 
-        LIST = [
-            name
-            for name, obj in inspect.getmembers(custmod)
-            if is_user_defined_function(obj)
-        ]
 
-        for st in LIST:
-            vars(igm)[st] = vars(custmod)[st]
+def run_intializers(modules: List, params: Any, state: State) -> None:
+    for module in modules:
+        module.initialize(params, state)
 
 
-# this function checks if there dependency modules that are not in the module list already
-# dependent modules are listed by function dependency_MODULENAME, and exist only if this function exists.
-def find_dependent_modules(modules):
-    dm = sum(
-        [
-            getattr(igm, "dependency_" + m)()
-            for m in modules
-            if hasattr(igm, "dependency_" + m)
-        ],
-        [],
-    )
+def run_processes(modules: List, params: Any, state: State) -> None:
+    if hasattr(state, "t"):
+        while state.t < params.time_end:
+            for module in modules:
+                module.update(params, state)
 
-    return [m for m in dm if m not in modules]
 
+def run_finalizers(modules: List, params: Any, state: State) -> None:
+    for module in modules:
+        module.finalize(params, state)
 
-# this add a logger to the state
-def add_logger(params, state, logging_level="INFO"):
-    import logging
 
+def add_logger(params, state) -> None:
     if params.logging_file == "":
         pathf = ""
     else:
-        pathf = os.path.join(params.working_dir, params.logging_file)
+        pathf = params.logging_file
 
     logging.basicConfig(
         filename=pathf,
         encoding="utf-8",
         filemode="w",
-        level=getattr(logging, logging_level),
+        level=params.logging_level,
         format="%(asctime)s - %(levelname)s - %(message)s",
     )
+    logging.root.setLevel(params.logging_level)
 
-    state.logger = logging.getLogger("my_logger")
-
+    state.logger = logging.getLogger("igm_logger")
     if not pathf == "":
         os.system("echo rm " + pathf + " >> clean.sh")
 
 
-# Print parameters in screen and a dedicated file
-def print_params(params):
-    param_file = os.path.join(params.working_dir, "params_saved.json")
+def remove_comments(json_str) -> str:
+    lines = json_str.split("\n")
+    cleaned_lines = [
+        line for line in lines if not line.strip().startswith(("//", "#"))
+    ]  # ! TODO: Add blocks comments...
+    cleaned_text = "\n".join(cleaned_lines)
+    return cleaned_text
+
 
+def load_json_file(
+    param_file: str, remove_param_comments: bool = True
+) -> Dict[str, Any]:
+    with open(param_file, "r") as json_file:
+        json_text = json_file.read()
+
+    if remove_param_comments:
+        json_text = remove_comments(json_text)
+
+    dic_params = json.loads(json_text)
+    return dic_params
+
+
+def get_modules_list(params_path: str) -> dict[str, List[str]]:
+    try:
+        with open(params_path) as f:
+            # params_dict = json.load(f) #re-instate if you want to enforce no comments in the json file(remove all the lines that load the clean json in this case...)
+            json_text = f.read()
+            json_cleaned = remove_comments(json_text)
+            params_dict = json.loads(json_cleaned)
+            module_dict = {
+                "modules_preproc": params_dict["modules_preproc"],
+                "modules_process": params_dict["modules_process"],
+                "modules_postproc": params_dict["modules_postproc"],
+            }
+
+            return module_dict
+    except JSONDecodeError as e:
+        raise JSONDecodeError(
+            msg="For the following line, please check the 1) JSON file structure and/or 2) make sure there are no comments (//, #, etc.)",
+            doc=e.doc,
+            pos=e.pos,
+        )
+
+
+def load_user_defined_params(param_file: str, params_dict: dict[str, Any]):
+    try:
+        json_defined_params = load_json_file(param_file=param_file)
+    except JSONDecodeError as e:
+        raise JSONDecodeError(
+            msg=f"Error decoding JSON: Please make sure your file path is correct or your json is properly formatted.",
+            doc=e.doc,
+            pos=e.pos,
+        )
+
+    unrecognized_json_arguments = {
+        k: v for k, v in json_defined_params.items() if k not in params_dict
+    }
+    params_dict.update(json_defined_params)
+
+    for key in unrecognized_json_arguments.keys():
+        warnings.warn(
+            f"The following argument specified in the JSON file does not exist among the core arguments nor the modules you have chosen. Ignoring: {key}"
+        )
+
+    return params_dict
+
+
+def load_modules(modules_dict: Dict) -> List[ModuleType]:
+    """Returns a list of actionable modules to then apply the update, initialize, finalize functions on for IGM."""
+
+    imported_preproc_modules = load_modules_from_directory(
+        modules_list=modules_dict["modules_preproc"], module_folder="preproc"
+    )
+    imported_process_modules = load_modules_from_directory(
+        modules_list=modules_dict["modules_process"], module_folder="process"
+    )
+    imported_postproc_modules = load_modules_from_directory(
+        modules_list=modules_dict["modules_postproc"], module_folder="postproc"
+    )
+    # ? Should we have custom modules in a seperate folder?
+    # imported_custom_modules = load_modules_from_directory(
+    #     modules_list=modules_dict.modules_custom, module_folder=params.modules_custom_folder
+    # )
+
+    return (
+        imported_preproc_modules
+        + imported_process_modules
+        + imported_postproc_modules  # + imported_custom_modules
+    )
+
+
+def validate_module(module) -> None:
+    """Validates that a module has the required functions to be used in IGM."""
+    required_functions = ["params", "initialize", "finalize", "update"]
+    for function in required_functions:
+        if not hasattr(module, function):
+            raise AttributeError(
+                f"Module {module} is missing the required function ({function}). If it is a custom python package, make sure to include it in the __init__.py file.",
+                f"Please see https://github.com/jouvetg/igm/wiki/5.-Custom-modules-(coding) for more information on how to construct custom modules.",
+            )
+
+
+def load_modules_from_directory(
+    modules_list: List[str], module_folder: str
+) -> List[ModuleType]:
+    imported_modules = []
+    for module_name in modules_list:
+        module_path = f"igm.modules.{module_folder}.{module_name}"
+        try:
+            module = importlib.import_module(module_path)
+        except ModuleNotFoundError:
+            logging.info(
+                f"Error importing module: {module_path}, checking for custom package in current working directory."
+            )
+            try:
+                logging.info(
+                    f"Trying to import custom module from current working directory (folder or .py): {module_name}"
+                )
+                module = importlib.import_module(module_name)
+            except ModuleNotFoundError:
+                raise ModuleNotFoundError(
+                    f"Can not find module {module_name}. Make sure it is either in the 1) {Path(igm.__file__).parent}/modules/{module_folder} directory or 2) in your current working directory."
+                )
+
+        validate_module(module)
+        imported_modules.append(module)
+
+    return imported_modules
+
+
+def has_dependencies(module: Any) -> bool:
+    if hasattr(module, "dependencies"):
+        return True
+    return False
+
+
+def load_dependent_modules(imported_modules: List) -> List[ModuleType]:
+    imported_dependencies = set()
+    for module in imported_modules:
+        if has_dependencies(module):
+            module_dependencies = module.dependencies
+            directories_to_search = [
+                "preproc",
+                "process",
+                "postproc",
+            ]  # will also check current working directory if any of them fail
+            load_modules_partial = partial(
+                load_modules_from_directory, module_dependencies
+            )
+            for directory in directories_to_search:
+                try:
+                    dependent_module = load_modules_partial(module_folder=directory)
+                    imported_dependencies.add(
+                        dependent_module[0]
+                    )  # [0] because it returns a singleton list
+                    logging.info(
+                        f"Found dependencies in directory {directory} or current working directory. Checking next module for dependencies."
+                    )
+                    break
+                except ModuleNotFoundError:
+                    logging.info(
+                        f"Could not find dependencies in directory {directory} or current working directory. Checking next IGM directory."
+                    )
+
+    return imported_modules + list(imported_dependencies)
+
+
+def print_gpu_info() -> None:
+    gpus = tf.config.experimental.list_physical_devices("GPU")
+    print(f"{'CUDA Enviroment':-^150}")
+    tf.sysconfig.get_build_info().pop("cuda_compute_capabilities", None)
+    print(f"{json.dumps(tf.sysconfig.get_build_info(), indent=2, default=str)}")
+    print(f"{'Available GPU Devices':-^150}")
+    for gpu in gpus:
+        gpu_info = {"gpu_id": gpu.name, "device_type": gpu.device_type}
+        device_details = tf.config.experimental.get_device_details(gpu)
+        gpu_info.update(device_details)
+
+        print(f"{json.dumps(gpu_info, indent=2, default=str)}")
+    print(f"{'':-^150}")
+
+
+# Print parameters in screen and a dedicated file
+def print_params(params: Namespace) -> None:
+    param_file = params.saved_params_filename
+    param_file = param_file + ".json"
+    
     # load the given parameters
     with open(param_file, "w") as json_file:
         json.dump(params.__dict__, json_file, indent=2)
diff --git a/igm/igm_run.py b/igm/igm_run.py
index 30c41d8..7653af3 100755
--- a/igm/igm_run.py
+++ b/igm/igm_run.py
@@ -3,69 +3,31 @@
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
 
-import numpy as np
-import matplotlib.pyplot as plt
 import tensorflow as tf
-import os
-import igm
+from igm import (
+    State,
+    params_core,
+    print_params,
+    run_intializers,
+    run_processes,
+    run_finalizers,
+    setup_igm,
+)
+
+
+def main() -> None:
+    state = State()  # class acting as a dictionary
+    parser = params_core()
+    imported_modules, params, state = setup_igm(state=state, parser=parser)
 
-
-def main():
-    
-    print("-----------------------------------------------------------------")
-    print("Num GPUs Available: ", len(tf.config.list_physical_devices("GPU")))
-    print("-----------------------------------------------------------------")
-    
-    # Collect defaults, overide from json file, and parse all core parameters
-    parser = igm.params_core()
-    igm.overide_from_json_file(parser, check_if_params_exist=False)
-    params, __ = parser.parse_known_args()  # args=[] add this for jupyter notebook
-
-    # get the list of all modules in order
-    modules = params.modules_preproc + params.modules_process + params.modules_postproc
-
-    # load custom modules from file (must be called my_module_name.py) to igm
-    for module in modules:
-        igm.load_custom_module(params, module)
-
-    # get the list of all dependent modules, which parameters must be called too
-    dependent_modules = igm.find_dependent_modules(modules)
-
-    # Collect defaults, overide from json file, and parse all specific module parameters
-    for module in modules + dependent_modules:
-        getattr(igm, "params_" + module)(parser)
-    igm.overide_from_json_file(parser, check_if_params_exist=True)
-    params = parser.parse_args()  # args=[] add this for jupyter notebook
-
-    # print definive parameters in a file for record
     if params.print_params:
-        igm.print_params(params)
-
-    # Define a state class/dictionnary that contains all the data
-    state = igm.State()
-
-    # if logging is activated, add a logger to the state
-    if params.logging:
-        igm.add_logger(params, state)
-    else:
-        os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
+        print_params(params=params)
 
     # Place the computation on your device GPU ('/GPU:0') or CPU ('/CPU:0')
-    with tf.device("/GPU:"+str(params.gpu)):
-        # Initialize all the model components in turn
-        for module in modules:
-            getattr(igm, "initialize_" + module)(params, state)
-
-        # Time loop, perform the simulation until reaching the defined end time
-        if hasattr(state, "t"):
-            while state.t < params.time_end:
-                # Update each model components in turn
-                for module in modules:
-                    getattr(igm, "update_" + module)(params, state)
-
-        # Finalize each module in turn
-        for module in modules:
-            getattr(igm, "finalize_" + module)(params, state)
+    with tf.device(f"/GPU:{params.gpu_id}"):  # type: ignore for linting checks
+        run_intializers(imported_modules, params, state)
+        run_processes(imported_modules, params, state)
+        run_finalizers(imported_modules, params, state)
 
 
 if __name__ == "__main__":
diff --git a/igm/instructed_oggm.py b/igm/instructed_oggm.py
index a7c79e2..f730b10 100755
--- a/igm/instructed_oggm.py
+++ b/igm/instructed_oggm.py
@@ -76,8 +76,16 @@ def __init__(
         
         """
 
-        parser = argparse.ArgumentParser(description="IGM")
-        igm.params_iceflow(parser)
+        # parser = argparse.ArgumentParser(description="IGM")
+        parser = igm.params_core()
+        
+        params, __ = parser.parse_known_args()  # args=[] add this for jupyter notebook
+             
+        imported_modules = igm.load_modules(params)
+
+        for module in imported_modules:
+            module.params(parser)
+        
         self.params = parser.parse_args(args=[])
 
         self.state = igm.State()
@@ -111,7 +119,7 @@ def __init__(
 
         self.icemask = mb_filter
 
-        igm.initialize_iceflow(self.params, self.state)
+        igm.modules.process.iceflow.initialize(self.params, self.state)
 
     def step(self, dt):
         # recast glacier variables into igm-like variables
@@ -120,7 +128,7 @@ def step(self, dt):
         self.state.usurf.assign(self.surface_h)
 
         # compute ubar and vbar
-        igm.update_iceflow(self.params, self.state)
+        igm.modules.process.iceflow.update(self.params, self.state)
 
         # retrurn the divergence of the flux using upwind fluxes
         divflux = (
diff --git a/igm/modules/__init__.py b/igm/modules/__init__.py
index cb01f05..1c50675 100755
--- a/igm/modules/__init__.py
+++ b/igm/modules/__init__.py
@@ -2,3 +2,9 @@
 
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
+
+from . import (
+    preproc,
+    process,
+    postproc,
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/anim_mayavi/__init__.py b/igm/modules/postproc/anim_mayavi/__init__.py
new file mode 100644
index 0000000..a2ab92c
--- /dev/null
+++ b/igm/modules/postproc/anim_mayavi/__init__.py
@@ -0,0 +1,6 @@
+from .anim_mayavi import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/anim_mayavi.md b/igm/modules/postproc/anim_mayavi/anim_mayavi.md
similarity index 100%
rename from igm/modules/postproc/anim_mayavi.md
rename to igm/modules/postproc/anim_mayavi/anim_mayavi.md
diff --git a/igm/modules/postproc/anim_mayavi.py b/igm/modules/postproc/anim_mayavi/anim_mayavi.py
similarity index 81%
rename from igm/modules/postproc/anim_mayavi.py
rename to igm/modules/postproc/anim_mayavi/anim_mayavi.py
index 33c1e57..ff99dc6 100755
--- a/igm/modules/postproc/anim_mayavi.py
+++ b/igm/modules/postproc/anim_mayavi/anim_mayavi.py
@@ -3,26 +3,23 @@
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
 
-import numpy as np
-import os, sys, shutil
-import time
 import tensorflow as tf
 import matplotlib.pyplot as plt
 
 
-def params_anim_mayavi(parser):
+def params(parser):
     pass
 
 
-def initialize_anim_mayavi(params, state):
+def initialize(params, state):
     pass
 
 
-def update_anim_mayavi(params, state):
+def update(params, state):
     pass
 
 
-def finalize_anim_mayavi(params, state):
+def finalize(params, state):
     from mayavi import mlab
     import xarray as xr
     import numpy as np
@@ -30,9 +27,7 @@ def finalize_anim_mayavi(params, state):
 
     plt.close("all")
 
-    ds = xr.open_dataset(
-        os.path.join(params.working_dir, "output.nc"), engine="netcdf4"
-    )
+    ds = xr.open_dataset("output.nc", engine="netcdf4")
 
     X, Y = np.meshgrid(ds.x, ds.y)
 
@@ -57,7 +52,7 @@ def finalize_anim_mayavi(params, state):
 
     VEC = hasattr(ds, "uvelsurf") & hasattr(ds, "vvelsurf") & hasattr(ds, "wvelsurf")
 
-    PAR = os.path.isdir(os.path.join(params.working_dir, "trajectories"))
+    PAR = os.path.isdir("trajectories")
 
     if VEC:
         quiv = mlab.quiver3d(
@@ -70,7 +65,7 @@ def finalize_anim_mayavi(params, state):
         )
 
     # if PAR:
-    #     f = os.path.join(params.working_dir, "trajectories",'traj-'+"{:06.0f}".format(TIME[0])+".csv")
+    #     f = os.path.join("trajectories",'traj-'+"{:06.0f}".format(TIME[0])+".csv")
     #     XYZ = np.loadtxt(f,skiprows=1,delimiter=',')
     #     pt3d = mlab.points3d(XYZ[:,1], XYZ[:,2], XYZ[:,3], colormap="RdBu",mode='point')
     #     pt3d.actor.property.point_size = 5
@@ -87,7 +82,7 @@ def anim(VEC, PAR):
                 quiv.mlab_source.v = np.array(tf.expand_dims(ds.vvelsurf[i], axis=0))
                 quiv.mlab_source.w = np.array(tf.expand_dims(ds.wvelsurf[i], axis=0))
             # if PAR:
-            #     f = os.path.join(params.working_dir, "trajectories",'traj-'+"{:06.0f}".format(TIME[i])+".csv")
+            #     f = os.path.join("trajectories",'traj-'+"{:06.0f}".format(TIME[i])+".csv")
             #     XYZ = np.loadtxt(f,skiprows=1,delimiter=',')
             #     pt3d.mlab_source.x = XYZ[:,1]
             #     pt3d.mlab_source.y = XYZ[:,2]
diff --git a/igm/modules/postproc/anim_plotly/__init__.py b/igm/modules/postproc/anim_plotly/__init__.py
new file mode 100644
index 0000000..9fa0075
--- /dev/null
+++ b/igm/modules/postproc/anim_plotly/__init__.py
@@ -0,0 +1,6 @@
+from .anim_plotly import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/anim_plotly.md b/igm/modules/postproc/anim_plotly/anim_plotly.md
similarity index 100%
rename from igm/modules/postproc/anim_plotly.md
rename to igm/modules/postproc/anim_plotly/anim_plotly.md
diff --git a/igm/modules/postproc/anim_plotly.py b/igm/modules/postproc/anim_plotly/anim_plotly.py
similarity index 97%
rename from igm/modules/postproc/anim_plotly.py
rename to igm/modules/postproc/anim_plotly/anim_plotly.py
index 7e771d0..3cc4aeb 100755
--- a/igm/modules/postproc/anim_plotly.py
+++ b/igm/modules/postproc/anim_plotly/anim_plotly.py
@@ -2,25 +2,25 @@
 import copy
 
 import numpy as np
-import os, sys, shutil
+import os
 import json
 from types import SimpleNamespace
 import xarray as xr
 
 
-def params_anim_plotly(parser):
+def params(parser):
     pass
 
 
-def initialize_anim_plotly(params, state):
+def initialize(params, state):
     pass
 
 
-def update_anim_plotly(params, state):
+def update(params, state):
     pass
 
 
-def finalize_anim_plotly(params, state):
+def finalize(params, state):
     import plotly.graph_objects as go
     import math
     from dash import Dash, dcc, html, Input, Output
@@ -107,9 +107,7 @@ def updata_graph(property, camera_angle, camera_height):
         params = json.loads(json_text, object_hook=lambda d: SimpleNamespace(**d))
 
         # read output.nc
-        ds = xr.open_dataset(
-            os.path.join(params.working_dir, "output.nc"), engine="netcdf4"
-        )
+        ds = xr.open_dataset("output.nc", engine="netcdf4")
 
         # get attributes from ds
         bedrock = np.array(ds.topg[0])
diff --git a/igm/modules/postproc/anim_video/__init__.py b/igm/modules/postproc/anim_video/__init__.py
new file mode 100644
index 0000000..5638e9c
--- /dev/null
+++ b/igm/modules/postproc/anim_video/__init__.py
@@ -0,0 +1,6 @@
+from .anim_video import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/anim_video.md b/igm/modules/postproc/anim_video/anim_video.md
similarity index 100%
rename from igm/modules/postproc/anim_video.md
rename to igm/modules/postproc/anim_video/anim_video.md
diff --git a/igm/modules/postproc/anim_video.py b/igm/modules/postproc/anim_video/anim_video.py
similarity index 73%
rename from igm/modules/postproc/anim_video.py
rename to igm/modules/postproc/anim_video/anim_video.py
index 5656a56..348c004 100755
--- a/igm/modules/postproc/anim_video.py
+++ b/igm/modules/postproc/anim_video/anim_video.py
@@ -4,31 +4,27 @@
 # Published under the GNU GPL (Version 3), check at the LICENSE file
 
 import numpy as np
-import os, sys, shutil
-import time
-import tensorflow as tf
+import os
 import matplotlib.pyplot as plt
 
 
-def params_anim_video(parser):
+def params(parser):
     pass
 
 
-def initialize_anim_video(params, state):
+def initialize(params, state):
     pass
 
 
-def update_anim_video(params, state):
+def update(params, state):
     pass
 
 
-def finalize_anim_video(params, state):
+def finalize(params, state):
     import xarray as xr
     from matplotlib import animation
 
-    ds = xr.open_dataset(
-        os.path.join(params.working_dir, "output.nc"), engine="netcdf4"
-    )
+    ds = xr.open_dataset("output.nc", engine="netcdf4")
 
     tas = ds.thk
 
@@ -65,8 +61,6 @@ def animate(frame):
         interval=500,  # ms between frames
     )
 
-    ani.save(os.path.join(params.working_dir, "animation.mp4"))
+    ani.save("animation.mp4")
 
-    os.system(
-        "echo rm " + os.path.join(params.working_dir, "animation.mp4") + " >> clean.sh"
-    )
+    os.system("echo rm " + "animation.mp4" + " >> clean.sh")
diff --git a/igm/modules/postproc/plot2d/__init__.py b/igm/modules/postproc/plot2d/__init__.py
new file mode 100644
index 0000000..4be0376
--- /dev/null
+++ b/igm/modules/postproc/plot2d/__init__.py
@@ -0,0 +1,6 @@
+from .plot2d import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/plot2d.md b/igm/modules/postproc/plot2d/plot2d.md
similarity index 100%
rename from igm/modules/postproc/plot2d.md
rename to igm/modules/postproc/plot2d/plot2d.md
diff --git a/igm/modules/postproc/plot2d.py b/igm/modules/postproc/plot2d/plot2d.py
similarity index 88%
rename from igm/modules/postproc/plot2d.py
rename to igm/modules/postproc/plot2d/plot2d.py
index b9e1e7f..eed9a28 100755
--- a/igm/modules/postproc/plot2d.py
+++ b/igm/modules/postproc/plot2d/plot2d.py
@@ -13,7 +13,7 @@
 from igm.modules.utils import *
 
 
-def params_plot2d(parser):
+def params(parser):
     parser.add_argument(
         "--plt2d_editor",
         type=str,
@@ -46,7 +46,7 @@ def params_plot2d(parser):
     )
 
 
-def initialize_plot2d(params, state):
+def initialize(params, state):
     state.extent = [np.min(state.x), np.max(state.x), np.min(state.y), np.max(state.y)]
 
     if params.plt2d_editor == "vs":
@@ -59,14 +59,10 @@ def initialize_plot2d(params, state):
     state.ax.axis("off")
     state.ax.set_aspect("equal")
 
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, params.plt2d_var + "*.png")
-        + " >> clean.sh"
-    )
+    os.system("echo rm " + params.plt2d_var + "*.png" + " >> clean.sh")
 
 
-def update_plot2d(params, state):
+def update(params, state):
     if state.saveresult:
         state.tcomp_plot2d.append(time.time())
 
@@ -120,10 +116,7 @@ def update_plot2d(params, state):
                 display(state.fig)
         else:
             plt.savefig(
-                os.path.join(
-                    params.working_dir,
-                    params.plt2d_var + "-" + str(state.t.numpy()).zfill(4) + ".png",
-                ),
+                params.plt2d_var + "-" + str(state.t.numpy()).zfill(4) + ".png",
                 bbox_inches="tight",
                 pad_inches=0.2,
             )
@@ -132,5 +125,5 @@ def update_plot2d(params, state):
         state.tcomp_plot2d[-1] *= -1
 
 
-def finalize_plot2d(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/postproc/print_comp/__init__.py b/igm/modules/postproc/print_comp/__init__.py
new file mode 100644
index 0000000..b1646f5
--- /dev/null
+++ b/igm/modules/postproc/print_comp/__init__.py
@@ -0,0 +1,6 @@
+from .print_comp import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/print_comp.md b/igm/modules/postproc/print_comp/print_comp.md
similarity index 100%
rename from igm/modules/postproc/print_comp.md
rename to igm/modules/postproc/print_comp/print_comp.md
diff --git a/igm/modules/postproc/print_comp.py b/igm/modules/postproc/print_comp/print_comp.py
similarity index 80%
rename from igm/modules/postproc/print_comp.py
rename to igm/modules/postproc/print_comp/print_comp.py
index 36c5ab7..cd859d7 100755
--- a/igm/modules/postproc/print_comp.py
+++ b/igm/modules/postproc/print_comp/print_comp.py
@@ -10,19 +10,19 @@
 import matplotlib.pyplot as plt
 
 
-def params_print_comp(parser):
+def params(parser):
     pass
 
 
-def initialize_print_comp(params, state):
+def initialize(params, state):
     pass
 
 
-def update_print_comp(params, state):
+def update(params, state):
     pass
 
 
-def finalize_print_comp(params, state):
+def finalize(params, state):
     ################################################################
 
     size_of_tensor = {}
@@ -41,16 +41,12 @@ def finalize_print_comp(params, state):
     )
 
     print("Memory statistics report:")
-    with open(os.path.join(params.working_dir, "memory-statistics.txt"), "w") as f:
+    with open("memory-statistics.txt", "w") as f:
         for key, value in size_of_tensor.items():
             print("     %24s  |  size : %8.4f Gb " % (key, value), file=f)
             print("     %24s  |  size : %8.4f Gb  " % (key, value))
 
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, "memory-statistics.txt")
-        + " >> clean.sh"
-    )
+    os.system("echo rm " + "memory-statistics.txt" + " >> clean.sh" )
 
     _plot_memory_pie(params, state)
 
@@ -61,9 +57,7 @@ def finalize_print_comp(params, state):
     state.tcomp_all = [np.sum([np.sum(getattr(state, m)) for m in modules])]
 
     print("Computational statistics report:")
-    with open(
-        os.path.join(params.working_dir, "computational-statistics.txt"), "w"
-    ) as f:
+    with open("computational-statistics.txt", "w") as f:
         for m in modules:
             CELA = (
                 m[6:],
@@ -81,11 +75,7 @@ def finalize_print_comp(params, state):
                 % CELA
             )
 
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, "computational-statistics.txt")
-        + " >> clean.sh"
-    )
+    os.system( "echo rm " + "computational-statistics.txt" + " >> clean.sh" )
 
     _plot_computational_pie(params, state)
 
@@ -129,15 +119,10 @@ def my_autopct(pct):
     plt.setp(autotexts, size=8, weight="bold")
     #    ax.set_title("Matplotlib bakery: A pie")
     plt.tight_layout()
-    plt.savefig(os.path.join(params.working_dir, "computational-pie.png"), pad_inches=0)
+    plt.savefig("computational-pie.png", pad_inches=0)
     plt.close("all")
 
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, "computational-pie.png")
-        + " >> clean.sh"
-    )
-
+    os.system("echo rm computational-pie.png >> clean.sh")
 
 def _plot_memory_pie(params, state):
     """
@@ -185,9 +170,7 @@ def my_autopct(pct):
     plt.setp(autotexts, size=8, weight="bold")
     #    ax.set_title("Matplotlib bakery: A pie")
     plt.tight_layout()
-    plt.savefig(os.path.join(params.working_dir, "memory-pie.png"), pad_inches=0)
+    plt.savefig("memory-pie.png", pad_inches=0)
     plt.close("all")
 
-    os.system(
-        "echo rm " + os.path.join(params.working_dir, "memory-pie.png") + " >> clean.sh"
-    )
+    os.system("echo rm memory-pie.png >> clean.sh")
diff --git a/igm/modules/postproc/print_info/__init__.py b/igm/modules/postproc/print_info/__init__.py
new file mode 100644
index 0000000..1768657
--- /dev/null
+++ b/igm/modules/postproc/print_info/__init__.py
@@ -0,0 +1,6 @@
+from .print_info import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/print_info.md b/igm/modules/postproc/print_info/print_info.md
similarity index 100%
rename from igm/modules/postproc/print_info.md
rename to igm/modules/postproc/print_info/print_info.md
diff --git a/igm/modules/postproc/print_info.py b/igm/modules/postproc/print_info/print_info.py
similarity index 84%
rename from igm/modules/postproc/print_info.py
rename to igm/modules/postproc/print_info/print_info.py
index 7751b4a..fe55b8d 100755
--- a/igm/modules/postproc/print_info.py
+++ b/igm/modules/postproc/print_info/print_info.py
@@ -9,17 +9,17 @@
 import matplotlib.pyplot as plt
 
 
-def params_print_info(parser):
+def params(parser):
     pass
 
 
-def initialize_print_info(params, state):
+def initialize(params, state):
     print(
         "IGM %s :         Iterations   |         Time (y)     |     Time Step (y)   |   Ice Volume (km^3) "
     )
 
 
-def update_print_info(params, state):
+def update(params, state):
     """
     This serves to print key info on the fly during computation
     """
@@ -36,5 +36,5 @@ def update_print_info(params, state):
         )
 
 
-def finalize_print_info(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/postproc/texture/__init__.py b/igm/modules/postproc/texture/__init__.py
new file mode 100644
index 0000000..fd3a582
--- /dev/null
+++ b/igm/modules/postproc/texture/__init__.py
@@ -0,0 +1,1 @@
+from .texture import params, initialize, update, finalize
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/constants.py b/igm/modules/postproc/texture/constants.py
new file mode 100644
index 0000000..0b8857e
--- /dev/null
+++ b/igm/modules/postproc/texture/constants.py
@@ -0,0 +1,20 @@
+from typing import Tuple
+from dataclasses import dataclass
+
+@dataclass(frozen=True)
+class FeatureConstants:
+    topg: Tuple[float, float] = (4151.99267578125, 0.0)
+    water: Tuple[float, float] = (1.0, 0.0)
+    vx: Tuple[float, float] = (3.5998013782501204, -4.1134811401367175)
+    thk: Tuple[float, float] = (125.9672691345204, 0.0)
+    prec: Tuple[float, float] = (463.0692443847656, 2.575910019874573)
+    ndvi: Tuple[float, float] = (221.4811553955078, -0.6113666892051697)
+    vy: Tuple[float, float] = (6.441854915618933, -1.314163446426393)
+    temp: Tuple[float, float] = (22.149993896484375, -8.147315979003906)
+
+
+@dataclass(frozen=True)
+class ImageConstants:
+    red: Tuple[float, float] = (255.0, 0.0)
+    blue: Tuple[float, float] = (255.0, 0.0)
+    green: Tuple[float, float] = (255.0, 0.0)
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/emulator.py b/igm/modules/postproc/texture/emulator.py
new file mode 100644
index 0000000..a8b7bb2
--- /dev/null
+++ b/igm/modules/postproc/texture/emulator.py
@@ -0,0 +1,19 @@
+import tensorflow as tf
+from abc import ABC, abstractmethod
+
+
+class Emulator(ABC):
+    @abstractmethod
+    def apply(self, array):
+        pass
+
+
+class LinearRegressor(Emulator):
+    def __init__(self, coefficients, b) -> None:
+        self.coefficients = coefficients
+        self.b = b
+
+    def apply(self, array: tf.Tensor) -> tf.Tensor:
+        y = tf.matmul(self.coefficients, array) + self.b
+
+        return y
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/exporter.py b/igm/modules/postproc/texture/exporter.py
new file mode 100644
index 0000000..d551099
--- /dev/null
+++ b/igm/modules/postproc/texture/exporter.py
@@ -0,0 +1,51 @@
+import tensorflow as tf
+import matplotlib.pyplot as plt
+from abc import ABC, abstractmethod
+
+from .normalizer import Normalizer
+
+class Exporter(ABC):
+    @abstractmethod
+    def export(self, image: tf.Tensor, filename: str) -> None:
+        pass
+
+
+class PngExporter(Exporter):
+    def convert_8bit(self, image: tf.Tensor, image_normalizer: Normalizer) -> tf.Tensor:
+        out = image_normalizer.unnormalize(image)
+        out = tf.cast(out, tf.uint8)
+        out = tf.squeeze(out).numpy()
+
+        return out
+
+    def unpad(self, image: tf.Tensor, padding_parameters) -> tf.Tensor:
+        n_dimensions = len(padding_parameters)
+        for i, dimension in enumerate(padding_parameters):
+            if dimension != (0, 0):
+                begin = [0] * n_dimensions
+                size = [-1] * n_dimensions
+                begin[i] = dimension[0]
+                size[i] = image.shape[i] - dimension[0] - dimension[1]
+                unpadded_image = tf.slice(image, begin=begin, size=size)
+
+        return unpadded_image
+
+    def export(self, image: tf.Tensor, filename: str) -> None:
+        fig = plt.figure()
+        plt.imshow(image, origin="lower")
+        plt.axis("off")
+        plt.savefig(
+            filename,
+            bbox_inches="tight",
+            # dpi=300,
+            pad_inches=0,
+        )
+        plt.close(fig)
+
+
+class TiffExporter(Exporter):
+    def export(self, image: tf.Tensor, filename: str):
+        raise NotImplementedError(
+            "Tiff/tif exporter not implemented yet. Choose (png) please."
+        )
+        # TODO: export image as tiff
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/image_data.py b/igm/modules/postproc/texture/image_data.py
new file mode 100644
index 0000000..c3cdbf0
--- /dev/null
+++ b/igm/modules/postproc/texture/image_data.py
@@ -0,0 +1,90 @@
+import tensorflow as tf
+
+from typing import Tuple
+from dataclasses import dataclass
+
+from .utils import resize_image
+from igm import State
+
+@dataclass
+class ImageData:
+    values: tf.Tensor
+    height: int
+    width: int
+    state: State
+    upsampled: bool = False
+    square: bool = False
+ 
+    def upsample(self, height: int, width: int) -> None:
+        self.values = resize_image(self.values, height, width)
+        self.upsampled = True
+
+    def compute_shape(self, resolution: int) -> Tuple[int, int]:
+        self.upsample_resolution = resolution
+        if self.height == self.width:
+            self.height = resolution
+            self.width = resolution
+            self.square = True
+
+            return resolution, resolution
+        elif self.height < self.width:
+            ratio = resolution / self.width
+            new_width = int(tf.math.ceil(ratio * self.width))
+            new_height = int(tf.math.ceil(ratio * self.height))
+
+            if new_height % 2 != 0:
+                new_height += 1
+
+            self.height = new_height
+            self.width = new_width
+
+            return new_width, new_height
+        else:
+            ratio = resolution / self.height
+            new_width = int(tf.math.ceil(ratio * self.width))
+            new_height = int(tf.math.ceil(ratio * self.height))
+
+            if new_width % 2 != 0:
+                new_width += 1
+
+            self.height = new_height
+            self.width = new_width
+
+            return new_width, new_height
+
+    def make_square(self, method: str = "reflect"):
+        if not self.upsampled:
+            self.upsample(self.height, self.width)
+            padding_parameters = [
+                (0, 0),
+                (0, 0),
+                (0, 0),
+                (0, 0),
+            ]  # [bs, h, w, c]
+        if self.height > self.width:
+            pad_width = (self.height - self.width) // 2
+            padding_parameters = [
+                (0, 0),
+                (0, 0),
+                (pad_width, pad_width),
+                (0, 0),
+            ]  # [bs, h, w + padding, c]
+        else:
+            pad_width = (self.width - self.height) // 2
+            padding_parameters = [
+                (0, 0),
+                (pad_width, pad_width),
+                (0, 0),
+                (0, 0),
+            ]  # [bs, h + padding, w, c]
+
+        self.values = tf.pad(self.values, padding_parameters, mode=method)
+
+        return padding_parameters
+
+    def run(self) -> None:
+        new_width, new_height = self.compute_shape(resolution=self.upsample_resolution)
+        self.upsample(width=new_width, height=new_height)
+        if not self.square:
+            padding_parameters = self.make_square()
+            self.state.texture_exporter.padding_parameters = padding_parameters
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/normalizer.py b/igm/modules/postproc/texture/normalizer.py
new file mode 100644
index 0000000..4050961
--- /dev/null
+++ b/igm/modules/postproc/texture/normalizer.py
@@ -0,0 +1,53 @@
+import tensorflow as tf
+from abc import ABC, abstractmethod
+
+from .constants import FeatureConstants, ImageConstants
+
+class Normalizer(ABC):
+    @abstractmethod
+    def normalize(self, image: tf.Tensor) -> tf.Tensor:
+        pass
+
+    @abstractmethod
+    def unnormalize(self, image: tf.Tensor) -> tf.Tensor:
+        pass
+
+
+
+
+from dataclasses import fields
+class FeatureNormalizer(Normalizer):
+    def __init__(self, constants: FeatureConstants):
+        self.constants = constants
+        self.constant_names = [field.name for field in fields(constants)]
+
+    def normalize(self, image: tf.Tensor, variable_name: str) -> tf.Tensor:
+        array_max, array_min = getattr(self.constants, variable_name)
+        normalized_image = (
+            2 * ((image - array_min) / (array_max - array_min)) - 1
+        )  # normalize to [-1, 1]
+
+        return normalized_image
+
+    def unnormalize(self, image: tf.Tensor) -> tf.Tensor:
+        return image
+
+    def normalize_all(self, image: tf.Tensor) -> tf.Tensor:
+        numpy_image = image.numpy()
+        for i, feature in enumerate(self.constant_names):
+            numpy_image[..., i] = self.normalize(numpy_image[..., i], feature)
+            # TODO: make it not cast to numpy and maybe define a custom tf.map or tf.function...
+
+        return tf.convert_to_tensor(numpy_image, dtype=tf.float32)
+
+class ImageNormalizer(Normalizer):
+    def __init__(self, constants: ImageConstants):
+        self.constants = constants
+
+    def normalize(self, image: tf.Tensor) -> tf.Tensor:
+        return image
+
+    def unnormalize(self, image: tf.TensorArray) -> tf.Tensor:
+        image = 255 * (image * 0.5 + 0.5)
+
+        return image
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd.py b/igm/modules/postproc/texture/pix2pixhd.py
new file mode 100644
index 0000000..f23ac68
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd.py
@@ -0,0 +1,49 @@
+from .normalizer import Normalizer
+from .image_data import ImageData
+from .exporter import Exporter
+from typing import Any
+import logging
+
+class Pix2PixHDPipeline:
+    def __init__(
+        self,
+        feature_normalizer: Normalizer,
+        image_normalizer: Normalizer,
+        image: ImageData,
+        exporter: Exporter,
+        model: Any,
+        state: Any,
+        params: Any,
+
+    ):
+        self.feature_normalizer: Normalizer = feature_normalizer
+        self.image_normalizer: Normalizer = image_normalizer
+        self.image: ImageData = image
+        self.exporter: Exporter = exporter
+        self.model: Any = model
+        self.state: Any = state
+        self.params: Any = params
+
+    def run(self) -> None:
+        year = int(round(self.state.t.numpy()))
+
+        self.image.values = self.feature_normalizer.normalize_all(self.image.values)
+
+        image = self.model(self.image.values, training=False)
+
+        logging.info(f"Output Image shape from model: {image.shape}")
+        if not self.image.square:
+            image = self.exporter.unpad(
+                image=image, padding_parameters=self.exporter.padding_parameters
+            )
+
+        logging.info(f"Output Image shape after unpadding (if needed): {image.shape}")
+        # Converts raw model output [-1,1] to [0,255] for imshow
+        image = self.exporter.convert_8bit(
+            image=image, image_normalizer=self.image_normalizer
+        )
+        logging.info(f"Image type after converting to 8bit (if needed): {image.dtype}")
+        logging.info(f"First column from Image after converting to 8bit (if needed): {image[:,0]}")
+        self.exporter.export(
+            image=image, filename=f"texture_{year}.{self.params.texture_format}"
+        )
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/generator.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/generator.py
new file mode 100644
index 0000000..c713040
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/generator.py
@@ -0,0 +1,178 @@
+import tensorflow as tf
+import tensorflow.keras as K
+from tensorflow_addons.layers import InstanceNormalization
+from tensorflow.keras import layers
+# from tensorflow.keras.layers import BatchNormalization
+# from tensorflow.keras.layers import Input
+
+from .layer_functions.reflection_pad2d import ReflectionPad2d
+from .layer_functions.resnet_block import ResnetBlock
+from .layer_functions.tanh import Tanh
+from .layer_functions.bilinear_upsampling import BilinearUpSampling
+
+weight_init = {}
+weight_init['conv'] = tf.random_normal_initializer(0.0, 0.02)
+weight_init['bn_gamma'] = tf.random_normal_initializer(1.0, 0.02)
+weight_init['bn_beta'] = tf.zeros_initializer()
+
+
+class LocalEnhancer(K.Model):
+    def __init__(self, input_nc, output_nc, ngf=32, n_downsample_global=3, n_blocks_global=9,
+                 n_local_enhancers=1, n_blocks_local=3, norm_layer=InstanceNormalization, padding_type='REFLECT'):
+        super(LocalEnhancer, self).__init__()
+        
+        # Setup variables
+        self.n_local_enhancers = n_local_enhancers
+        paddings = (3, 3)
+        activation = layers.LeakyReLU(0.2)
+        # self.split_model = split_model
+        
+        # Setup Global Generator to use
+        ngf_global = ngf * (2**n_local_enhancers)
+
+        model_global = GlobalGenerator(input_nc, output_nc, ngf=ngf_global, n_downsampling=n_downsample_global, n_blocks=n_blocks_global, norm_layer=norm_layer).layers[0].layers
+
+        model_global = model_global[:-3] # get rid of final convolution layers
+        self.model = K.Sequential(model_global)
+
+        # # Single Local Enhancer
+        # Downsampling layers
+        ngf_global = ngf * (2**(n_local_enhancers-1))
+        
+        model_downsample = K.Sequential(name="DownSampler")
+        model_downsample.add(ReflectionPad2d(paddings))
+        model_downsample.add(layers.Conv2D(ngf_global, 7, kernel_initializer=weight_init['conv']))
+        model_downsample.add(norm_layer()) # can't seem to figure out how to alter RandomNormal mean/std
+        model_downsample.add(activation)
+        model_downsample.add(layers.Conv2D(ngf_global * 2, 3, strides=2, padding='same', kernel_initializer=weight_init['conv']))
+        model_downsample.add(norm_layer()) # can't seem to figure out how to alter RandomNormal mean/std
+        model_downsample.add(activation)
+
+        # Residual blocks
+        model_upsample = K.Sequential(name="ResnetAndUpsampler")
+        for i in range(n_blocks_local):
+            model_upsample.add(ResnetBlock(ngf_global * 2, padding_type=padding_type, norm_layer=norm_layer))
+
+        # Upsampling layers
+        # model_upsample.add(layers.Conv2DTranspose(ngf_global, 3, strides=2, padding='same', output_padding=1, kernel_initializer=weight_init['conv'])) # change this to bilinear?
+        model_upsample.add(BilinearUpSampling(ngf=ngf_global, mult=2, weight_init=weight_init, name=f"UpConv2D_0")) # FIX THIS TO WORK WITH LOCAL
+        model_upsample.add(norm_layer()) # can't seem to figure out how to alter RandomNormal mean/std
+        model_upsample.add(activation)
+
+        # Final Convolutional
+        model_upsample.add(ReflectionPad2d(paddings))
+        model_upsample.add(layers.Conv2D(output_nc, 7, kernel_initializer=weight_init['conv'])) # change this to bilinear?
+        model_upsample.add(Tanh(name='LocalOutput'))
+
+        # Create model attributes for calling function
+        self.model_downsampler = model_downsample
+        self.model_upsampler = model_upsample
+
+        # Layer to reduce resolution by half...
+        self.downsample = layers.AveragePooling2D(3, strides=2, padding='same')
+
+    def call(self, x):
+
+        # Initial input (full resolution image)
+        input = [x]
+        
+        # Creating low resolution input -> input is now [high_res, low_res]
+        low_res_input = self.downsample(input[-1])
+        input.append(low_res_input)
+
+        # Pass low resolution input into global model to get feature map
+        global_feature_map = self.model(input[-1])
+
+        # Get downsampling and upsampling parts of the local enhancer
+        model_downsample = getattr(self, 'model_downsampler')
+        model_upsample = getattr(self, 'model_upsampler')
+        
+        # Get high resolution input and pass it through the local enhancers downsampler
+        output = model_upsample(model_downsample(input[0]) + global_feature_map)
+
+        return output
+    # def call(self, x):
+    #     # Initial input (full resolution image)
+    #     input = [x]
+    #     #print("am I spliting?", self.split_model)
+    #     # Creating low resolution input -> input is now [high_res, low_res]
+    #     low_res_input = self.downsample(input[-1])
+    #     input.append(low_res_input)
+
+    #     # Pass low resolution input into global model to get feature map
+    #     if self.split_model:
+    #         with tf.device("/gpu:0"):
+    #             global_feature_map = self.model(input[-1])
+    #     else:
+    #         global_feature_map = self.model(input[-1])
+
+    #     # Get downsampling and upsampling parts of the local enhancer
+    #     model_downsample = getattr(self, 'model_downsampler')
+    #     model_upsample = getattr(self, 'model_upsampler')
+        
+    #     # Get high resolution input and pass it through the local enhancers downsampler
+    #     if self.split_model:
+    #      #   print("yep, im splitting")
+    #         with tf.device("/gpu:1"):
+    #             output = model_upsample(model_downsample(input[0]) + global_feature_map)
+    #     else:
+    #         output = model_upsample(model_downsample(input[0]) + global_feature_map)
+
+    #     return output
+
+class GlobalGenerator(K.Model):
+    def __init__(self, input_nc, output_nc, ngf=64, n_downsampling=3, n_blocks=9,
+                 norm_layer=InstanceNormalization, # instance norm?
+                 padding_type='REFLECT'):
+        assert(n_blocks >= 0)
+        super(GlobalGenerator, self).__init__()
+        
+        # Setup initial variables
+        self.padding_type = padding_type
+        self.output_nc = output_nc
+        paddings = (3, 3)
+
+        # Setup Global Model
+        model = K.Sequential(name="GlobalGenerator")
+        model.add(layers.InputLayer([None, None, input_nc], name=f"Input"))
+        
+        # Initial Block
+        model.add(ReflectionPad2d(paddings, name=f"InReflectPad"))
+        model.add(layers.Conv2D(ngf, 7, kernel_initializer=weight_init['conv'], name=f"InConv2D"))
+        model.add(norm_layer(name=f"InNorm"))
+        model.add(layers.LeakyReLU(0.2, name=f"InActivation"))
+
+        # Downsampling blocks
+        for i in range(n_downsampling):
+            mult = 2**i
+            # model.add(ReflectionPad2d(paddings, name=f"ReflectPad_{i}"))
+            model.add(layers.Conv2D(ngf * mult * 2, 3, strides=2, padding='same', kernel_initializer=weight_init['conv'], name=f"Conv2D_{i}"))
+            model.add(norm_layer(name=f"DownNorm_{i}"))
+            model.add(layers.LeakyReLU(0.2, name=f"DownActivation_{i}"))
+
+        # Resnet blocks
+        mult = 2**n_downsampling
+        for i in range(n_blocks):
+            model.add(ResnetBlock(dim=ngf * mult, padding_type=padding_type, norm_layer=norm_layer, activation=layers.ReLU(), name=f"Resnet_{i}"))
+
+        # Upsampling blocks
+        for i in range(n_downsampling):
+            mult = 2**(n_downsampling - i)
+            # model.add(layers.Conv2DTranspose(int(ngf * mult / 2), 3, strides=2, padding='same', output_padding=1, kernel_initializer=weight_init['conv'], name=f"TransConv2D_{i}"))
+            model.add(BilinearUpSampling(ngf=ngf, mult=mult, weight_init=weight_init))
+            model.add(norm_layer(name=f"UpNorm_{i}"))
+            model.add(layers.LeakyReLU(0.2, name=f"UpActivation_{i}"))
+            
+        # Final layers
+        model.add(ReflectionPad2d(paddings, name=f"OutReflectPad"))
+        model.add(layers.Conv2D(self.output_nc, 7, kernel_initializer=weight_init['conv'], name=f"OutConv2D"))
+        model.add(Tanh(name='GlobalOutput'))
+
+        # Save global model
+        self.model = model
+
+    # @tf.function
+    def call(self, x):
+        x = self.model(x)
+
+        return x
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/bilinear_upsampling.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/bilinear_upsampling.py
new file mode 100644
index 0000000..19baa31
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/bilinear_upsampling.py
@@ -0,0 +1,18 @@
+from tensorflow.keras import layers
+
+class BilinearUpSampling(layers.Layer):
+    def __init__(self, ngf, mult, weight_init, **kwargs):
+        super(BilinearUpSampling, self).__init__(**kwargs)
+        self.ngf = ngf
+        self.mult = mult
+        self.weight_init = weight_init
+        # self.name = name
+
+    def build(self, input_shape):
+        self.upsampling = layers.UpSampling2D(size=2, interpolation='nearest')
+        self.conv = layers.Conv2D(filters=int(self.ngf * self.mult / 2), kernel_size=3, padding='same', kernel_initializer=self.weight_init['conv'])
+
+    def call(self, inputs):
+        x = self.upsampling(inputs)
+        x = self.conv(x)
+        return x
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/reflection_pad2d.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/reflection_pad2d.py
new file mode 100644
index 0000000..21d4850
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/reflection_pad2d.py
@@ -0,0 +1,17 @@
+import tensorflow as tf
+from tensorflow.keras import layers
+
+
+class ReflectionPad2d(layers.Layer):
+    def __init__(self, paddings=(1, 1), **kwargs):
+        self.paddings = tuple(paddings)
+        self.input_spec = [layers.InputSpec(ndim=4)]
+        super(ReflectionPad2d, self).__init__(**kwargs)
+
+    def compute_output_shape(self, s):
+        return (s[0], s[1]+self.paddings[0], s[2]+self.paddings[1], s[3])
+
+    def call(self, x):
+        w_pad, h_pad = self.paddings
+        x = tf.pad(x, [[0, 0], [h_pad, h_pad], [w_pad, w_pad], [0, 0]], 'REFLECT')
+        return x
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/resnet_block.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/resnet_block.py
new file mode 100644
index 0000000..88bdd1e
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/resnet_block.py
@@ -0,0 +1,32 @@
+import tensorflow as tf
+import tensorflow.keras as K
+from tensorflow.keras import layers
+
+weight_init = {}
+weight_init['conv'] = tf.random_normal_initializer(0.0, 0.02)
+weight_init['bn_gamma'] = tf.random_normal_initializer(1.0, 0.02)
+weight_init['bn_beta'] = tf.zeros_initializer()
+
+class ResnetBlock(layers.Layer):
+    def __init__(self, dim, padding_type, norm_layer,
+                 activation=layers.ReLU(), use_dropout=False, **kwargs):
+        super(ResnetBlock, self).__init__(**kwargs)
+        # CONSTANT REFLECT SYMMETRIC
+        model = K.Sequential((layers.ZeroPadding2D(padding=(1, 1)),
+                              layers.Conv2D(dim, 3, kernel_initializer=weight_init['conv']),
+                              norm_layer(),
+                              activation
+                              ))
+        if use_dropout:
+            model.add(layers.Dropout(0.5))
+        model.add(layers.ZeroPadding2D([1, 1]))
+        model.add(layers.Conv2D(dim, 3, kernel_initializer=weight_init['conv']))
+        model.add(norm_layer())
+        self.model = model
+
+    def call(self, x):
+        identity = x
+        x = self.model(x)
+        out = x + identity
+
+        return out
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/tanh.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/tanh.py
new file mode 100644
index 0000000..a46cee4
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/layer_functions/tanh.py
@@ -0,0 +1,10 @@
+import tensorflow.keras as K
+from tensorflow.keras import layers
+
+
+class Tanh(layers.Layer):
+    def __init__(self, **kwargs):
+        super(Tanh, self).__init__(**kwargs)
+
+    def call(self, x):
+        return K.activations.tanh(x)
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/pix2pixhd_model_assets/loading.py b/igm/modules/postproc/texture/pix2pixhd_model_assets/loading.py
new file mode 100644
index 0000000..4b09bf7
--- /dev/null
+++ b/igm/modules/postproc/texture/pix2pixhd_model_assets/loading.py
@@ -0,0 +1,21 @@
+import tensorflow as tf
+import os
+
+def load_model_test(checkpoint, checkpoint_dir):
+    
+    # if model.opt.load_pretrain_dir == '':
+    #     print("Training from scratch.")
+    #     return 0
+    # else:
+    #     if model.opt.which_epoch == 'latest':
+    model_path = tf.train.latest_checkpoint(checkpoint_dir)
+    #     else:
+    # model_path = os.path.join(model_path)
+
+    try:
+        checkpoint.restore(model_path).expect_partial()
+        print("Successfully restored model.")
+    except AssertionError as error:
+        print(error)
+    except:
+        raise ImportError("Could not find model to load. Check pretrain model path.")
diff --git a/igm/modules/postproc/texture/preparer.py b/igm/modules/postproc/texture/preparer.py
new file mode 100644
index 0000000..edbf567
--- /dev/null
+++ b/igm/modules/postproc/texture/preparer.py
@@ -0,0 +1,79 @@
+import tensorflow as tf
+
+from typing import Tuple, Any
+from .emulator import Emulator
+
+class Pix2PixHDImagePreparer:
+    def __init__(self, params: Any, state: Any, ndvi_emulator: Emulator = None) -> None:
+        if ndvi_emulator is None and (not hasattr(state, "ndvi")):
+            raise ValueError(
+                "Emulator is required for NDVI if no NDVI is provided."
+            )  # ? getting a type error with None here... could be a bug to look into
+        self.state = state
+        self.params = params
+        self.ndvi_emulator = ndvi_emulator
+
+    def get_features(self) -> None:
+        topg = self.state.topg
+        thk = self.state.thk
+        vx = self.state.U[-1]
+        vy = self.state.V[-1]
+
+        try:
+            air_temp = self.state.air_temp
+        except AttributeError:
+            raise AttributeError("Air temperature is required for the texture model. Please include it in an input file with the variable name 'air_temp'. \
+                                Also, note that air_temp will be averaged to obtain the annually averaged air_temp for the spatial grid.")
+        
+        try:
+            precipitation = self.state.precipitation
+        except AttributeError:
+            raise AttributeError("Precipitation is required for the texture model. Please include it in an input file with the variable name 'precipitation'. \
+                                Also, note that precipitation will be averaged to obtain the annually averaged precipitation for the spatial grid.")
+            
+        # if the user supplies monthly, quarerly, or seasonal data, average to get annual data
+        if air_temp.shape[0] > 1:
+            temp = tf.math.reduce_mean(air_temp, axis=0)
+        else:
+            temp = air_temp
+        if precipitation.shape[0] > 1:
+            prec = tf.math.reduce_mean(precipitation, axis=0)
+        else:
+            prec = precipitation
+
+        prec = self.prec_units(prec)
+
+        if hasattr(self.state, "ndvi"):
+            ndvi = self.state.ndvi
+        else:
+            ndvi = self.compute_ndvi(temp, prec)
+
+        if hasattr(self.state, "water"):
+            water = self.state.water
+        else:
+            water = self.compute_water(topg.shape)
+
+        self.image = tf.stack([topg, water, vx, thk, prec, ndvi, vy, temp], axis=0)
+
+    def prec_units(self, prec: tf.Tensor) -> tf.Tensor:
+        return (prec / self.params.divide_by_density) * (1000 / 12)  # (m / yr. -> mm / month)
+
+    def compute_ndvi(self, temp: tf.Tensor, prec: tf.Tensor) -> tf.Tensor:
+        temp_flat = tf.reshape(temp, [-1])
+        prec_flat = tf.reshape(prec, [-1])
+        climate = tf.stack((temp_flat, prec_flat))
+
+        ndvi = self.ndvi_emulator.apply(climate)
+        ndvi = tf.reshape(ndvi, [temp.shape[0], temp.shape[1]])
+        ndvi = tf.clip_by_value(ndvi, 0, 255)
+
+        return ndvi
+
+    def compute_water(self, tensor_shape: Tuple[int, int]):
+        return tf.zeros(shape=tensor_shape)  # assumes no water for now...
+
+    def prepare_batch(self) -> tf.Tensor:
+        image = tf.transpose(self.image, [1, 2, 0])
+        image = tf.expand_dims(image, axis=0)
+
+        return image
\ No newline at end of file
diff --git a/igm/modules/postproc/texture/texture.py b/igm/modules/postproc/texture/texture.py
new file mode 100644
index 0000000..eef8c0c
--- /dev/null
+++ b/igm/modules/postproc/texture/texture.py
@@ -0,0 +1,149 @@
+from tensorflow.keras.models import load_model
+import tensorflow as tf
+
+from typing import Any
+import time
+import os
+import igm
+import logging
+
+from .normalizer import FeatureNormalizer, ImageNormalizer
+from .utils import TextureModelNotFoundError, resize_image
+from .emulator import LinearRegressor
+from .exporter import PngExporter, TiffExporter
+from .constants import FeatureConstants, ImageConstants
+from .preparer import Pix2PixHDImagePreparer
+from .image_data import ImageData
+from .pix2pixhd import Pix2PixHDPipeline
+from .pix2pixhd_model_assets.generator import LocalEnhancer
+from .pix2pixhd_model_assets.loading import load_model_test
+
+TEXTURE_DEFAULT_DIR = igm.__path__[0] + "/modules/postproc/texture/"
+TEXTURE_DEFAULT_PATH = os.path.join(TEXTURE_DEFAULT_DIR, "pix2pixhd-texture-model")
+TEXTURE_CKPT_DIR = os.path.join(TEXTURE_DEFAULT_DIR, "checkpoints")
+
+def params(parser: Any) -> None:
+    parser.add_argument(
+        "--texture_format",
+        type=str,
+        default="png",
+        help="Format of the texture image (png, tif, or tiff)",
+    )
+    parser.add_argument(
+        "--texture_model_path",
+        type=str,
+        default=TEXTURE_DEFAULT_PATH,
+        help="Name of the folder for the texture model (tf format)",
+    )
+    parser.add_argument(
+        "--texture_verbosity",
+        type=int,
+        default=30,
+        help="Python Logger verbosity level (10=DEBUG, 20=INFO, 30=WARNING, 40=ERROR, 50=CRITICAL)",
+    )
+    
+    parser.add_argument(
+        "--divide_by_density",
+        type=float,
+        default=1,
+        help="This parameter solves an imcompatibility (this option will be removed in the future)",
+    )
+    # TODO: Add an option where they can overwrite certain features with their own data (e.g. ndvi, water, etc.), so one can you high resolution topg for example...
+    # TODO: Add a logger for verbosity
+
+def initialize(params: Any, state: Any) -> None:
+    state.tcomp_texture = []
+    logging.basicConfig(level=params.texture_verbosity)
+    if not os.path.exists(params.texture_model_path):
+        model_url = "https://drive.google.com/drive/folders/1UP761XQpD4UvqNtbKNO20EXrxmAd4uoA?usp=sharing"
+
+        # TODO: Only available to download folder if you use external packages (i do not know, but for now, I will let the user manually download from google drive)
+
+        raise TextureModelNotFoundError(
+            f"Model not found.\n\nPlease download the model\n{model_url})\nand place the downloaded folder in the following directory:\n{TEXTURE_DEFAULT_DIR}"
+        )
+
+    # state.texture_model = load_model(params.texture_model_path, compile=False)
+    state.texture_model = LocalEnhancer(input_nc=8, output_nc=3, ngf=32, n_downsample_global=4, n_blocks_global=9, n_local_enhancers=1, n_blocks_local=3)
+    checkpoint_dict = {"generator": state.texture_model}
+    checkpoint = tf.train.Checkpoint(**checkpoint_dict)
+
+    __ = load_model_test(checkpoint, TEXTURE_CKPT_DIR)
+
+    feature_constants = FeatureConstants()
+    image_constants = ImageConstants()
+    state.feature_normalizer = FeatureNormalizer(constants=feature_constants)
+    state.image_normalizer = ImageNormalizer(constants=image_constants)
+
+    b = 15.335304
+    coefficients = tf.constant([2.0448825, 0.62188774], shape=(1, 2))
+    state.ndvi_emulator = LinearRegressor(coefficients=coefficients, b=b)
+
+    if params.texture_format == "png":
+        state.texture_exporter = PngExporter()
+    elif params.texture_format == ("tif" or "tiff"):
+        state.texture_exporter = TiffExporter()
+    else:
+        raise NotImplementedError(
+            "Texture format not implemented. Please choose one of the following: (png, tif, or tiff)"
+        )
+def is_power_of_two(number):
+    """Checks if a number is a power of two"""
+    return (number & (number - 1)) == 0
+
+def nearest_power_of_two(number):
+    """Finds nearest power of two"""
+    import math
+    return 2 ** round(math.log2(number))
+
+
+def update(params: Any, state: Any) -> None:
+    if state.saveresult:
+        state.tcomp_texture.append(time.time())
+
+        preparer = Pix2PixHDImagePreparer(
+            state=state, params=params, ndvi_emulator=state.ndvi_emulator
+        )
+        preparer.get_features()
+        image = preparer.prepare_batch()
+
+        logging.info(f"Input Image shape (before resizing): {image.shape}")
+        data = ImageData(values=image, height=image.shape[1], width=image.shape[2], state=state)
+        
+        # TODO: Check if the image is a power of 2...
+        if data.height > data.width:
+            if not is_power_of_two(data.height):
+                resolution = nearest_power_of_two(data.height)
+        elif data.height < data.width:
+            if not is_power_of_two(data.width):
+                resolution = nearest_power_of_two(data.width)
+        else:
+            resolution = data.height # or data.width but it assumes its a square
+        logging.info(f"Long side resolution for resizing: {resolution}")
+        new_width, new_height = data.compute_shape(
+            resolution=resolution
+        )
+        data.upsample(width=new_width, height=new_height)
+        logging.info(f"Input Image shape (after resizing): {data.values.shape}")
+        if not data.square:
+            padding_parameters = data.make_square()
+            state.texture_exporter.padding_parameters = padding_parameters
+
+        pipeline = Pix2PixHDPipeline(
+            feature_normalizer=state.feature_normalizer,
+            image_normalizer=state.image_normalizer,
+            image=data,
+            model=state.texture_model,
+            exporter=state.texture_exporter,
+            state=state,
+            params=params,
+        )
+
+        pipeline.run()
+
+        state.tcomp_texture[-1] -= time.time()
+        state.tcomp_texture[-1] *= -1
+
+
+def finalize(params: Any, state: Any) -> None:
+    pass
diff --git a/igm/modules/postproc/texture/utils.py b/igm/modules/postproc/texture/utils.py
new file mode 100644
index 0000000..1dbc95d
--- /dev/null
+++ b/igm/modules/postproc/texture/utils.py
@@ -0,0 +1,16 @@
+import tensorflow as tf
+
+class TextureModelNotFoundError(Exception):
+    pass
+
+def resize_image(
+    input_image: tf.Tensor,
+    height: int,
+    width: int,
+    upsampling_method=tf.image.ResizeMethod.BICUBIC,
+):
+    input_image = tf.image.resize(
+        input_image, [height, width], method=upsampling_method
+    )
+
+    return input_image
\ No newline at end of file
diff --git a/igm/modules/postproc/write_ncdf/__init__.py b/igm/modules/postproc/write_ncdf/__init__.py
new file mode 100644
index 0000000..33838c2
--- /dev/null
+++ b/igm/modules/postproc/write_ncdf/__init__.py
@@ -0,0 +1,6 @@
+from .write_ncdf import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/write_ncdf.md b/igm/modules/postproc/write_ncdf/write_ncdf.md
similarity index 100%
rename from igm/modules/postproc/write_ncdf.md
rename to igm/modules/postproc/write_ncdf/write_ncdf.md
diff --git a/igm/modules/postproc/write_ncdf.py b/igm/modules/postproc/write_ncdf/write_ncdf.py
similarity index 91%
rename from igm/modules/postproc/write_ncdf.py
rename to igm/modules/postproc/write_ncdf/write_ncdf.py
index b2251c2..5953ad6 100755
--- a/igm/modules/postproc/write_ncdf.py
+++ b/igm/modules/postproc/write_ncdf/write_ncdf.py
@@ -14,7 +14,7 @@
 from igm.modules.utils import getmag
 
 
-def params_write_ncdf(parser):
+def params(parser):
     parser.add_argument(
         "--wncd_output_file",
         type=str,
@@ -24,7 +24,7 @@ def params_write_ncdf(parser):
 
     parser.add_argument(
         "--wncd_vars_to_save",
-        type=list,
+        nargs='+',
         default=[
             "topg",
             "usurf",
@@ -40,14 +40,10 @@ def params_write_ncdf(parser):
     )
 
 
-def initialize_write_ncdf(params, state):
+def initialize(params, state):
     state.tcomp_write_ncdf = []
 
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, params.wncd_output_file)
-        + " >> clean.sh"
-    )
+    os.system("echo rm " + params.wncd_output_file + " >> clean.sh")
 
     # give information on variables for output ncdf, TODO: IMPROVE
     state.var_info_ncdf_ex = {}
@@ -87,7 +83,7 @@ def initialize_write_ncdf(params, state):
     state.var_info_ncdf_ex["weight_particles"] = ["weight_particles", "no"]
 
 
-def update_write_ncdf(params, state):
+def update(params, state):
     if state.saveresult:
         state.tcomp_write_ncdf.append(time.time())
 
@@ -112,11 +108,7 @@ def update_write_ncdf(params, state):
             if hasattr(state, "logger"):
                 state.logger.info("Initialize NCDF ex output Files")
 
-            nc = Dataset(
-                os.path.join(params.working_dir, params.wncd_output_file),
-                "w",
-                format="NETCDF4",
-            )
+            nc = Dataset(params.wncd_output_file, "w", format="NETCDF4")
 
             nc.createDimension("time", None)
             E = nc.createVariable("time", np.dtype("float32").char, ("time",))
@@ -172,11 +164,7 @@ def update_write_ncdf(params, state):
                     "Write NCDF ex file at time : " + str(state.t.numpy())
                 )
 
-            nc = Dataset(
-                os.path.join(params.working_dir, params.wncd_output_file),
-                "a",
-                format="NETCDF4",
-            )
+            nc = Dataset( params.wncd_output_file, "a", format="NETCDF4" )
 
             d = nc.variables["time"][:].shape[0]
             nc.variables["time"][d] = state.t.numpy()
@@ -194,5 +182,5 @@ def update_write_ncdf(params, state):
         state.tcomp_write_ncdf[-1] *= -1
 
 
-def finalize_write_ncdf(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/postproc/write_particles/__init__.py b/igm/modules/postproc/write_particles/__init__.py
new file mode 100644
index 0000000..a09a54a
--- /dev/null
+++ b/igm/modules/postproc/write_particles/__init__.py
@@ -0,0 +1,6 @@
+from .write_particles import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/write_particles.md b/igm/modules/postproc/write_particles/write_particles.md
similarity index 100%
rename from igm/modules/postproc/write_particles.md
rename to igm/modules/postproc/write_particles/write_particles.md
diff --git a/igm/modules/postproc/write_particles.py b/igm/modules/postproc/write_particles/write_particles.py
similarity index 78%
rename from igm/modules/postproc/write_particles.py
rename to igm/modules/postproc/write_particles/write_particles.py
index 6703851..d393279 100755
--- a/igm/modules/postproc/write_particles.py
+++ b/igm/modules/postproc/write_particles/write_particles.py
@@ -5,13 +5,13 @@
 
 import numpy as np
 import os
-import datetime, time
+import time
 import tensorflow as tf
 
 from igm.modules.utils import *
 
 
-def params_write_particles(parser):
+def params(parser):
     parser.add_argument(
         "--wpar_add_topography",
         type=str2bool,
@@ -20,22 +20,18 @@ def params_write_particles(parser):
     )
 
 
-def initialize_write_particles(params, state):
+def initialize(params, state):
     state.tcomp_write_particles = []
 
-    directory = os.path.join(params.working_dir, "trajectories")
+    directory = "trajectories"
     if os.path.exists(directory):
         shutil.rmtree(directory)
     os.mkdir(directory)
 
-    os.system(
-        "echo rm -r "
-        + os.path.join(params.working_dir, "trajectories")
-        + " >> clean.sh"
-    )
+    os.system( "echo rm -r " + "trajectories" + " >> clean.sh" )
 
     if params.wpar_add_topography:
-        ftt = os.path.join(params.working_dir, "trajectories", "topg.csv")
+        ftt = os.path.join("trajectories", "topg.csv")
         array = tf.transpose(
             tf.stack(
                 [state.X[state.X > 0], state.Y[state.X > 0], state.topg[state.X > 0]]
@@ -44,12 +40,11 @@ def initialize_write_particles(params, state):
         np.savetxt(ftt, array, delimiter=",", fmt="%.2f", header="x,y,z")
 
 
-def update_write_particles(params, state):
+def update(params, state):
     if state.saveresult:
         state.tcomp_write_particles.append(time.time())
 
         f = os.path.join(
-            params.working_dir,
             "trajectories",
             "traj-" + "{:06d}".format(int(state.t.numpy())) + ".csv",
         )
@@ -71,13 +66,12 @@ def update_write_particles(params, state):
         )
         np.savetxt(f, array, delimiter=",", fmt="%.2f", header="Id,x,y,z,rh,t,englt")
 
-        ft = os.path.join(params.working_dir, "trajectories", "time.dat")
+        ft = os.path.join("trajectories", "time.dat")
         with open(ft, "a") as f:
             print(state.t.numpy(), file=f)
 
         if params.wpar_add_topography:
             ftt = os.path.join(
-                params.working_dir,
                 "trajectories",
                 "usurf-" + "{:06d}".format(int(state.t.numpy())) + ".csv",
             )
@@ -96,5 +90,5 @@ def update_write_particles(params, state):
         state.tcomp_write_particles[-1] *= -1
 
 
-def finalize_write_particles(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/postproc/write_tif/__init__.py b/igm/modules/postproc/write_tif/__init__.py
new file mode 100644
index 0000000..537d26b
--- /dev/null
+++ b/igm/modules/postproc/write_tif/__init__.py
@@ -0,0 +1,6 @@
+from .write_tif import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/write_tif.md b/igm/modules/postproc/write_tif/write_tif.md
similarity index 100%
rename from igm/modules/postproc/write_tif.md
rename to igm/modules/postproc/write_tif/write_tif.md
diff --git a/igm/modules/postproc/write_tif.py b/igm/modules/postproc/write_tif/write_tif.py
similarity index 85%
rename from igm/modules/postproc/write_tif.py
rename to igm/modules/postproc/write_tif/write_tif.py
index 101dd78..92b2ac3 100755
--- a/igm/modules/postproc/write_tif.py
+++ b/igm/modules/postproc/write_tif/write_tif.py
@@ -8,7 +8,7 @@
 import tensorflow as tf
 
 
-def params_write_tif(parser):
+def params(parser):
     parser.add_argument(
         "--wtif_vars_to_save",
         type=list,
@@ -20,18 +20,16 @@ def params_write_tif(parser):
     )
 
 
-def initialize_write_tif(params, state):
+def initialize(params, state):
     pass
 
 
-def update_write_tif(params, state):
+def update(params, state):
     import rasterio
 
     if state.saveresult:
         for var in params.wtif_vars_to_save:
-            file = os.path.join(
-                params.working_dir, var + "-" + str(int(state.t)).zfill(6) + ".tif"
-            )
+            file = var + "-" + str(int(state.t)).zfill(6) + ".tif"
 
             os.system("echo rm " + file + " >> clean.sh")
 
@@ -60,5 +58,5 @@ def update_write_tif(params, state):
             del src
 
 
-def finalize_write_tif(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/postproc/write_ts/__init__.py b/igm/modules/postproc/write_ts/__init__.py
new file mode 100644
index 0000000..f4d7900
--- /dev/null
+++ b/igm/modules/postproc/write_ts/__init__.py
@@ -0,0 +1,6 @@
+from .write_ts import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/postproc/write_ts.md b/igm/modules/postproc/write_ts/write_ts.md
similarity index 100%
rename from igm/modules/postproc/write_ts.md
rename to igm/modules/postproc/write_ts/write_ts.md
diff --git a/igm/modules/postproc/write_ts.py b/igm/modules/postproc/write_ts/write_ts.py
similarity index 73%
rename from igm/modules/postproc/write_ts.py
rename to igm/modules/postproc/write_ts/write_ts.py
index 1c62637..77fc6fc 100755
--- a/igm/modules/postproc/write_ts.py
+++ b/igm/modules/postproc/write_ts/write_ts.py
@@ -4,13 +4,11 @@
 # Published under the GNU GPL (Version 3), check at the LICENSE file
 
 import numpy as np
-import os, sys, shutil
-import matplotlib.pyplot as plt
-import tensorflow as tf
+import os
 from netCDF4 import Dataset
 
 
-def params_write_ts(parser):
+def params(parser):
     parser.add_argument(
         "--wts_output_file",
         type=str,
@@ -19,19 +17,15 @@ def params_write_ts(parser):
     )
 
 
-def initialize_write_ts(params, state):
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, params.wts_output_file)
-        + " >> clean.sh"
-    )
+def initialize(params, state):
+    os.system( "echo rm " + params.wts_output_file + " >> clean.sh" )
 
     state.var_info_ncdf_ts = {}
     state.var_info_ncdf_ts["vol"] = ["Ice volume", "km^3"]
     state.var_info_ncdf_ts["area"] = ["Glaciated area", "km^2"]
 
 
-def update_write_ts(params, state):
+def update(params, state):
     if state.saveresult:
         vol = np.sum(state.thk) * (state.dx**2) / 10**9
         area = np.sum(state.thk > 1) * (state.dx**2) / 10**6
@@ -42,11 +36,7 @@ def update_write_ts(params, state):
             if hasattr(state, "logger"):
                 state.logger.info("Initialize NCDF ts output Files")
 
-            nc = Dataset(
-                os.path.join(params.working_dir, params.wts_output_file),
-                "w",
-                format="NETCDF4",
-            )
+            nc = Dataset( params.wts_output_file,"w", format="NETCDF4" )
 
             nc.createDimension("time", None)
             E = nc.createVariable("time", np.dtype("float32").char, ("time",))
@@ -68,11 +58,7 @@ def update_write_ts(params, state):
                     "Write NCDF ts file at time : " + str(state.t.numpy())
                 )
 
-            nc = Dataset(
-                os.path.join(params.working_dir, params.wts_output_file),
-                "a",
-                format="NETCDF4",
-            )
+            nc = Dataset( params.wts_output_file, "a", format="NETCDF4" )
             d = nc.variables["time"][:].shape[0]
 
             nc.variables["time"][d] = state.t.numpy()
@@ -81,5 +67,5 @@ def update_write_ts(params, state):
             nc.close()
 
 
-def finalize_write_ts(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/preproc/__init__.py b/igm/modules/preproc/__init__.py
index cb01f05..ad781bc 100755
--- a/igm/modules/preproc/__init__.py
+++ b/igm/modules/preproc/__init__.py
@@ -2,3 +2,23 @@
 
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
+
+# from . import (
+# 	include_icemask,
+# 	load_ncdf,
+# 	load_tif,
+# 	optimize,
+# 	optimize_v1,
+# 	oggm_shop,
+# 	pretraining,
+# )
+
+# __all__ = ['load_ncdf',
+#            'load_tif',
+#            'include_icemask',
+#            'optimize',
+#            'optimize_v1',
+#            'oggm_shop',
+#            'pretraining']
+
+# this above would cause conflicts with finalize, initialize, params, update functions... (i think)
\ No newline at end of file
diff --git a/igm/modules/preproc/include_icemask/__init__.py b/igm/modules/preproc/include_icemask/__init__.py
new file mode 100644
index 0000000..1b6015d
--- /dev/null
+++ b/igm/modules/preproc/include_icemask/__init__.py
@@ -0,0 +1,6 @@
+from .include_icemask import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/include_icemask.md b/igm/modules/preproc/include_icemask/include_icemask.md
similarity index 100%
rename from igm/modules/preproc/include_icemask.md
rename to igm/modules/preproc/include_icemask/include_icemask.md
diff --git a/igm/modules/preproc/include_icemask.py b/igm/modules/preproc/include_icemask/include_icemask.py
similarity index 95%
rename from igm/modules/preproc/include_icemask.py
rename to igm/modules/preproc/include_icemask/include_icemask.py
index 238450b..f7dd227 100755
--- a/igm/modules/preproc/include_icemask.py
+++ b/igm/modules/preproc/include_icemask/include_icemask.py
@@ -30,7 +30,7 @@
 from shapely.geometry import Point
 
 
-def params_include_icemask(parser):
+def params(parser):
     parser.add_argument(
         "--mask_shapefile",
         type=str,
@@ -46,7 +46,7 @@ def params_include_icemask(parser):
     )
 
 
-def initialize_include_icemask(params, state):
+def initialize(params, state):
     # read_shapefile
     gdf = read_shapefile(params.mask_shapefile)
 
@@ -83,11 +83,11 @@ def initialize_include_icemask(params, state):
     state.icemask = tf.constant(mask_values)
 
 
-def update_include_icemask(params, state):
+def update(params, state):
     pass
 
 
-def finalize_include_icemask(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/preproc/load_ncdf/__init__.py b/igm/modules/preproc/load_ncdf/__init__.py
new file mode 100644
index 0000000..cc5b0b1
--- /dev/null
+++ b/igm/modules/preproc/load_ncdf/__init__.py
@@ -0,0 +1,6 @@
+from .load_ncdf import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/load_ncdf.md b/igm/modules/preproc/load_ncdf/load_ncdf.md
similarity index 100%
rename from igm/modules/preproc/load_ncdf.md
rename to igm/modules/preproc/load_ncdf/load_ncdf.md
diff --git a/igm/modules/preproc/load_ncdf.py b/igm/modules/preproc/load_ncdf/load_ncdf.py
similarity index 94%
rename from igm/modules/preproc/load_ncdf.py
rename to igm/modules/preproc/load_ncdf/load_ncdf.py
index bad17fa..b43acde 100755
--- a/igm/modules/preproc/load_ncdf.py
+++ b/igm/modules/preproc/load_ncdf/load_ncdf.py
@@ -14,7 +14,7 @@
 from igm.modules.utils import *
 
 
-def params_load_ncdf(parser):
+def params(parser):
     parser.add_argument(
         "--lncd_input_file",
         type=str,
@@ -59,11 +59,11 @@ def params_load_ncdf(parser):
     )
 
 
-def initialize_load_ncdf(params, state):
+def initialize(params, state):
     if hasattr(state, "logger"):
         state.logger.info("LOAD NCDF file")
 
-    nc = Dataset(os.path.join(params.working_dir, params.lncd_input_file), "r")
+    nc = Dataset(params.lncd_input_file, "r")
 
     x = np.squeeze(nc.variables["x"]).astype("float32")
     y = np.squeeze(nc.variables["y"]).astype("float32")
@@ -131,9 +131,9 @@ def initialize_load_ncdf(params, state):
     complete_data(state)
 
 
-def update_load_ncdf(params, state):
+def update(params, state):
     pass
 
 
-def finalize_load_ncdf(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/preproc/load_tif/__init__.py b/igm/modules/preproc/load_tif/__init__.py
new file mode 100644
index 0000000..e94f53d
--- /dev/null
+++ b/igm/modules/preproc/load_tif/__init__.py
@@ -0,0 +1,6 @@
+from .load_tif import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/load_tif.md b/igm/modules/preproc/load_tif/load_tif.md
similarity index 100%
rename from igm/modules/preproc/load_tif.md
rename to igm/modules/preproc/load_tif/load_tif.md
diff --git a/igm/modules/preproc/load_tif.py b/igm/modules/preproc/load_tif/load_tif.py
similarity index 94%
rename from igm/modules/preproc/load_tif.py
rename to igm/modules/preproc/load_tif/load_tif.py
index dcf43a3..4e4ff2f 100755
--- a/igm/modules/preproc/load_tif.py
+++ b/igm/modules/preproc/load_tif/load_tif.py
@@ -10,7 +10,7 @@
 from igm.modules.utils import *
 
 
-def params_load_tif(parser):
+def params(parser):
     parser.add_argument(
         "--ltif_coarsen",
         type=int,
@@ -49,10 +49,10 @@ def params_load_tif(parser):
     )
 
 
-def initialize_load_tif(params, state):
+def initialize(params, state):
     import rasterio
 
-    files = glob.glob(os.path.join(params.working_dir, "*.tif"))
+    files = glob.glob("*.tif")
 
     for file in files:
         var = os.path.split(file)[-1].split(".")[0]
@@ -110,9 +110,9 @@ def initialize_load_tif(params, state):
     complete_data(state)
 
 
-def update_load_tif(params, state):
+def update(params, state):
     pass
 
 
-def finalize_load_tif(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/preproc/oggm_shop/__init__.py b/igm/modules/preproc/oggm_shop/__init__.py
new file mode 100644
index 0000000..774e79f
--- /dev/null
+++ b/igm/modules/preproc/oggm_shop/__init__.py
@@ -0,0 +1,6 @@
+from .oggm_shop import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/oggm_shop.md b/igm/modules/preproc/oggm_shop/oggm_shop.md
similarity index 100%
rename from igm/modules/preproc/oggm_shop.md
rename to igm/modules/preproc/oggm_shop/oggm_shop.md
diff --git a/igm/modules/preproc/oggm_shop.py b/igm/modules/preproc/oggm_shop/oggm_shop.py
similarity index 96%
rename from igm/modules/preproc/oggm_shop.py
rename to igm/modules/preproc/oggm_shop/oggm_shop.py
index 1f4a330..304f072 100755
--- a/igm/modules/preproc/oggm_shop.py
+++ b/igm/modules/preproc/oggm_shop/oggm_shop.py
@@ -13,7 +13,7 @@
 from igm.modules.utils import complete_data
 
 
-def params_oggm_shop(parser):
+def params(parser):
     # aletsch  RGI60-11.01450
     # malspina RGI60-01.13696
     # brady    RGI60-01.20796
@@ -73,7 +73,7 @@ def params_oggm_shop(parser):
     )
 
 
-def initialize_oggm_shop(params, state):
+def initialize(params, state):
     import json
 
     _oggm_util([params.oggm_RGI_ID], params)
@@ -81,7 +81,7 @@ def initialize_oggm_shop(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Prepare data using oggm and glathida")
 
-    nc = Dataset(os.path.join(params.working_dir, params.oggm_RGI_ID, "gridded_data.nc"), "r+")
+    nc = Dataset(os.path.join(params.oggm_RGI_ID, "gridded_data.nc"), "r+")
 
     x = np.squeeze(nc.variables["x"]).astype("float32")
     y = np.flip(np.squeeze(nc.variables["y"]).astype("float32"))
@@ -193,9 +193,7 @@ def initialize_oggm_shop(params, state):
         var_info["icemask"] = ["Ice mask", "no unit"]
         var_info["dhdt"] = ["Ice thickness change", "m/y"]
 
-        nc = Dataset(
-            os.path.join(params.working_dir, "input_saved.nc"), "w", format="NETCDF4"
-        )
+        nc = Dataset( "input_saved.nc", "w", format="NETCDF4" )
 
         nc.createDimension("y", len(y))
         yn = nc.createVariable("y", np.dtype("float32").char, ("y",))
@@ -227,11 +225,11 @@ def initialize_oggm_shop(params, state):
 #            state.logger.setLevel(params.logging_level)
 
 
-def update_oggm_shop(params, state):
+def update(params, state):
     pass
 
 
-def finalize_oggm_shop(params, state):
+def finalize(params, state):
     pass
 
 
@@ -377,17 +375,13 @@ def _oggm_util(RGIs, params):
                                                 gdirs, informed_threestep=True)
 
     source_folder = gdirs[0].get_filepath("gridded_data").split("gridded_data.nc")[0]
-    destination_folder = os.path.join(params.working_dir, params.oggm_RGI_ID)
+    destination_folder = params.oggm_RGI_ID
 
     if os.path.exists(destination_folder):
         shutil.rmtree(destination_folder)
     shutil.copytree(source_folder, destination_folder)
 
-    os.system(
-        "echo rm -r "
-        + os.path.join(params.working_dir, params.oggm_RGI_ID)
-        + " >> clean.sh"
-    )
+    os.system( "echo rm -r " + params.oggm_RGI_ID + " >> clean.sh" )
 
 
 def _read_glathida(x, y, usurf, proj, path_glathida, state):
diff --git a/igm/modules/preproc/optimize/__init__.py b/igm/modules/preproc/optimize/__init__.py
new file mode 100644
index 0000000..4ac343b
--- /dev/null
+++ b/igm/modules/preproc/optimize/__init__.py
@@ -0,0 +1,8 @@
+from .optimize import (
+	params,
+	initialize,
+	finalize,
+	update,
+)
+
+# dependencies = ["iceflow"] Can NOT use if you also inherit params from another module in this modules params (conflict will occur with same name)
\ No newline at end of file
diff --git a/igm/modules/preproc/optimize.md b/igm/modules/preproc/optimize/optimize.md
similarity index 100%
rename from igm/modules/preproc/optimize.md
rename to igm/modules/preproc/optimize/optimize.md
diff --git a/igm/modules/preproc/optimize.py b/igm/modules/preproc/optimize/optimize.py
similarity index 95%
rename from igm/modules/preproc/optimize.py
rename to igm/modules/preproc/optimize/optimize.py
index 7a626ef..2c3ebbc 100755
--- a/igm/modules/preproc/optimize.py
+++ b/igm/modules/preproc/optimize/optimize.py
@@ -14,14 +14,21 @@
 from netCDF4 import Dataset
 
 from igm.modules.utils import *
-from igm.modules.process.iceflow import *
+from igm.modules.process.iceflow import initialize as initialize_iceflow
+from igm.modules.process.iceflow import params as params_iceflow
+from igm.modules.process.iceflow.iceflow import (
+    fieldin_to_X,
+    update_2d_iceflow_variables,
+    iceflow_energy_XY,
+    Y_to_UV,
+    save_iceflow_model
+)
 
 
-def dependency_optimize():
-    return ["iceflow"]
+def params(parser):
+    # dependency on iceflow parameters...
+    params_iceflow(parser)
 
-
-def params_optimize(parser):
     parser.add_argument(
         "--opti_vars_to_save",
         type=list,
@@ -201,7 +208,7 @@ def params_optimize(parser):
         help="uniformize the density of thkobs",
     )
 
-def initialize_optimize(params, state):
+def initialize(params, state):
     """
     This function does the data assimilation (inverse modelling) to optimize thk, slidingco ans usurf from data
     """
@@ -298,12 +305,16 @@ def initialize_optimize(params, state):
             # evalutae th ice flow emulator
             Y = state.iceflow_model(X)
 
-            state.U, state.V = Y_to_UV(params, Y)
-
-            state.U = state.U[0]
-            state.V = state.V[0]
+            U, V = Y_to_UV(params, Y)
 
-            update_2d_iceflow_variables(params, state)
+            U = U[0]
+            V = V[0]
+ 
+            state.ubar = tf.reduce_sum(U * state.vert_weight, axis=0)
+            state.vbar = tf.reduce_sum(V * state.vert_weight, axis=0)
+            
+            state.uvelsurf = U[-1]
+            state.vvelsurf = V[-1]
 
             state.velsurf = tf.stack(
                 [state.uvelsurf, state.vvelsurf], axis=-1
@@ -568,14 +579,14 @@ def initialize_optimize(params, state):
     plt.close("all")
 
     np.savetxt(
-        os.path.join(params.working_dir, "costs.dat"),
+        "costs.dat",
         np.stack(state.costs),
         fmt="%.10f",
         header="        COST_U        COST_H      COST_D       COST_S       REGU_H       REGU_S          HPO           COSTGLEN ",
     )
 
     np.savetxt(
-        os.path.join(params.working_dir, "rms_std.dat"),
+        "rms_std.dat",
         np.stack(
             [
                 state.rmsthk,
@@ -594,19 +605,21 @@ def initialize_optimize(params, state):
     )
 
     os.system(
-        "echo rm " + os.path.join(params.working_dir, "rms_std.dat") + " >> clean.sh"
+        "echo rm " + "rms_std.dat" + " >> clean.sh"
     )
     os.system(
-        "echo rm " + os.path.join(params.working_dir, "costs.dat") + " >> clean.sh"
+        "echo rm " + "costs.dat" + " >> clean.sh"
     )
 
+    # Flag so we can check if initialize was already called
+    state.optimize_initializer_called = True
 
-def update_optimize(params, state):
-    pass
 
+def update(params, state):
+    pass
 
-def finalize_optimize(params, state):
 
+def finalize(params, state):
     if params.iflo_save_model:
         save_iceflow_model(params, state) 
 
@@ -699,7 +712,7 @@ def _update_ncdf_optimize(params, state, it):
 
     if it == 0:
         nc = Dataset(
-            os.path.join(params.working_dir, "optimize.nc"),
+            "optimize.nc",
             "w",
             format="NETCDF4",
         )
@@ -733,18 +746,10 @@ def _update_ncdf_optimize(params, state, it):
 
         nc.close()
 
-        os.system(
-            "echo rm "
-            + os.path.join(params.working_dir, "optimize.nc")
-            + " >> clean.sh"
-        )
+        os.system( "echo rm " + "optimize.nc" + " >> clean.sh" )
 
     else:
-        nc = Dataset(
-            os.path.join(params.working_dir, "optimize.nc"),
-            "a",
-            format="NETCDF4",
-        )
+        nc = Dataset("optimize.nc", "a", format="NETCDF4", )
 
         d = nc.variables["iterations"][:].shape[0]
 
@@ -762,7 +767,7 @@ def _output_ncdf_optimize_final(params, state):
     """
 
     nc = Dataset(
-        os.path.join(params.working_dir, params.opti_save_result_in_ncdf),
+        params.opti_save_result_in_ncdf,
         "w",
         format="NETCDF4",
     )
@@ -791,7 +796,7 @@ def _output_ncdf_optimize_final(params, state):
 
     os.system(
         "echo rm "
-        + os.path.join(params.working_dir, params.opti_save_result_in_ncdf)
+        + params.opti_save_result_in_ncdf
         + " >> clean.sh"
     )
 
@@ -813,12 +818,12 @@ def _plot_cost_functions(params, state, costs):
     plt.ylim(0, 1)
     plt.legend()
 
-    plt.savefig(os.path.join(params.working_dir, "convergence.png"), pad_inches=0)
+    plt.savefig("convergence.png", pad_inches=0)
     plt.close("all")
 
     os.system(
         "echo rm "
-        + os.path.join(params.working_dir, "convergence.png")
+        + "convergence.png"
         + " >> clean.sh"
     )
 
@@ -994,15 +999,10 @@ def _update_plot_inversion(params, state, i):
             clear_output(wait=True)
             display(state.fig)
     else:
-        plt.savefig(
-            os.path.join(params.working_dir, "resu-opti-" + str(i).zfill(4) + ".png"),
-            pad_inches=0,
-        )
+        plt.savefig("resu-opti-" + str(i).zfill(4) + ".png", pad_inches=0 )
         plt.close("all")
 
-        os.system(
-            "echo rm " + os.path.join(params.working_dir, "*.png") + " >> clean.sh"
-        )
+        os.system( "echo rm " + "*.png" + " >> clean.sh" )
 
 
 def _update_plot_inversion_simple(params, state, i):
@@ -1099,13 +1099,13 @@ def _update_plot_inversion_simple(params, state, i):
             display(state.fig)
     else:
         plt.savefig(
-            os.path.join(params.working_dir, "resu-opti-" + str(i).zfill(4) + ".png"),
+            "resu-opti-" + str(i).zfill(4) + ".png",
             pad_inches=0,
         )
         plt.close("all")
 
         os.system(
-            "echo rm " + os.path.join(params.working_dir, "*.png") + " >> clean.sh"
+            "echo rm " + "*.png" + " >> clean.sh"
         )
 
 
diff --git a/igm/modules/preproc/optimize_v1/__init__.py b/igm/modules/preproc/optimize_v1/__init__.py
new file mode 100644
index 0000000..d75bb4c
--- /dev/null
+++ b/igm/modules/preproc/optimize_v1/__init__.py
@@ -0,0 +1,6 @@
+from .optimize_v1 import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/optimize_v1.md b/igm/modules/preproc/optimize_v1/optimize_v1.md
similarity index 100%
rename from igm/modules/preproc/optimize_v1.md
rename to igm/modules/preproc/optimize_v1/optimize_v1.md
diff --git a/igm/modules/preproc/optimize_v1.py b/igm/modules/preproc/optimize_v1/optimize_v1.py
similarity index 99%
rename from igm/modules/preproc/optimize_v1.py
rename to igm/modules/preproc/optimize_v1/optimize_v1.py
index 71a03fc..80d3c7b 100755
--- a/igm/modules/preproc/optimize_v1.py
+++ b/igm/modules/preproc/optimize_v1/optimize_v1.py
@@ -14,10 +14,16 @@
 from netCDF4 import Dataset
 
 from igm.modules.utils import *
-from igm.modules.process.iceflow_v1 import *
+from igm.modules.process.iceflow_v1.iceflow_v1 import initialize as initialize_iceflow_v1
 
-
-def params_optimize_v1(parser):
+def params(parser):
+    parser.add_argument(
+        "--working_dir",
+        type=str,
+        default="",
+        help="Working directory (default empty string)",
+    )
+    
     parser.add_argument(
         "--opti_vars_to_save",
         type=list,
@@ -169,7 +175,7 @@ def params_optimize_v1(parser):
     )
 
 
-def initialize_optimize_v1(params, state):
+def initialize(params, state):
     initialize_iceflow_v1(params, state)
 
     ###### PERFORM CHECKS PRIOR OPTIMIZATIONS
@@ -697,11 +703,11 @@ def initialize_optimize_v1(params, state):
     )
 
 
-def update_optimize_v1(params, state):
+def update(params, state):
     pass
 
 
-def finalize_optimize_v1(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/preproc/pretraining/__init__.py b/igm/modules/preproc/pretraining/__init__.py
new file mode 100644
index 0000000..a9bc514
--- /dev/null
+++ b/igm/modules/preproc/pretraining/__init__.py
@@ -0,0 +1,6 @@
+from .pretraining import (
+	params,
+	initialize,
+	finalize,
+	update
+)
\ No newline at end of file
diff --git a/igm/modules/preproc/pretraining.md b/igm/modules/preproc/pretraining/pretraining.md
similarity index 100%
rename from igm/modules/preproc/pretraining.md
rename to igm/modules/preproc/pretraining/pretraining.md
diff --git a/igm/modules/preproc/pretraining.py b/igm/modules/preproc/pretraining/pretraining.py
similarity index 95%
rename from igm/modules/preproc/pretraining.py
rename to igm/modules/preproc/pretraining/pretraining.py
index 5d52e53..6d10120 100755
--- a/igm/modules/preproc/pretraining.py
+++ b/igm/modules/preproc/pretraining/pretraining.py
@@ -10,15 +10,16 @@
 import os
 import xarray
 
-from igm.modules.process.iceflow import *
+from igm.modules.process.iceflow import params as params_iceflow
+from igm.modules.process.iceflow.iceflow import *
 from igm.modules.utils import *
 
+ 
+def params(parser):
+    
+    # dependency on iceflow parameters...
+    params_iceflow(parser)
 
-def dependency_pretraining():
-    return ["iceflow"]
-
-
-def params_pretraining(parser):
     parser.add_argument(
         "--data_dir",
         type=str,
@@ -96,7 +97,7 @@ def params_pretraining(parser):
     )
 
 
-def initialize_pretraining(params, state):
+def initialize(params, state):
     state.direct_name = (
         "pinnbp"
         + "_"
@@ -117,11 +118,11 @@ def initialize_pretraining(params, state):
         str(params.iflo_dim_arrhenius) + "_" + str(int(params.iflo_new_friction_param))
     )
 
-    os.makedirs(os.path.join(params.working_dir, state.direct_name), exist_ok=True)
+    os.makedirs( state.direct_name, exist_ok=True)
 
     os.system(
         "echo rm -r "
-        + os.path.join(params.working_dir, state.direct_name)
+        + state.direct_name
         + " >> clean.sh"
     )
 
@@ -157,11 +158,11 @@ def initialize_pretraining(params, state):
     train_iceflow_emulator(params, state, subdatasetpath_train)
 
 
-def update_pretraining(params, state):
+def update(params, state):
     pass
 
 
-def finalize_pretraining(params, state):
+def finalize(params, state):
     pass
 
 
@@ -232,7 +233,7 @@ def compute_solutions(params, state):
         )
 
         # define path
-        path = os.path.join(params.working_dir, state.direct_name, code)
+        path = os.path.join(state.direct_name, code)
         os.makedirs(path, exist_ok=True)
 
         fig = plt.figure(figsize=(10, 10))
@@ -273,9 +274,10 @@ def train_iceflow_emulator(params, state, trainingset, augmentation=True):
     if os.path.exists("model0.h5"):
         state.iceflow_model = tf.keras.models.load_model("model0.h5", compile=False)
     else:
-        state.iceflow_model = getattr(igm, params.iflo_network)(
-            params, nb_inputs, nb_outputs
-        )
+        if params.iflo_network=='cnn':
+            state.iceflow_model = cnn(params, nb_inputs, nb_outputs)
+        elif params.iflo_network=='unet':
+            state.iceflow_model = unet(params, nb_inputs, nb_outputs)
 
     state.iceflow_model.summary(line_length=130)
 
@@ -394,9 +396,7 @@ def train_iceflow_emulator(params, state, trainingset, augmentation=True):
             optimizer.lr = params.iflo_retrain_emulator_lr * (0.9 ** (epoch / 100))
 
         if epoch % (params.epochs // 5) == 0:
-            pp = os.path.join(
-                params.working_dir, state.direct_name, "model-" + str(epoch) + ".h5"
-            )
+            pp = os.path.join( state.direct_name, "model-" + str(epoch) + ".h5" )
             state.iceflow_model.save(pp)
 
         if epoch % params.freq_test == 0:
@@ -424,7 +424,7 @@ def train_iceflow_emulator(params, state, trainingset, augmentation=True):
                     + str(int(par[4]))
                 )
 
-                path = os.path.join(params.working_dir, state.direct_name, code)
+                path = os.path.join(state.direct_name, code)
 
                 YP = state.iceflow_model(X)
 
@@ -491,7 +491,7 @@ def train_iceflow_emulator(params, state, trainingset, augmentation=True):
     plt.legend()
     plt.tight_layout()
     plt.savefig(
-        os.path.join(params.working_dir, state.direct_name, "MISFIT.png"), pad_inches=0
+        os.path.join(state.direct_name, "MISFIT.png"), pad_inches=0
     )
     plt.close("all")
 
@@ -518,13 +518,13 @@ def train_iceflow_emulator(params, state, trainingset, augmentation=True):
     plt.legend()
     plt.tight_layout()
     plt.savefig(
-        os.path.join(params.working_dir, state.direct_name, "MISFIT_CO.png"),
+        os.path.join(state.direct_name, "MISFIT_CO.png"),
         pad_inches=0,
     )
     plt.close("all")
 
     state.iceflow_model.save(
-        os.path.join(params.working_dir, state.direct_name, "model.h5")
+        os.path.join(state.direct_name, "model.h5")
     )
 
 
diff --git a/igm/modules/process/__init__.py b/igm/modules/process/__init__.py
index cb01f05..a18a39e 100755
--- a/igm/modules/process/__init__.py
+++ b/igm/modules/process/__init__.py
@@ -2,3 +2,5 @@
 
 # Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
 # Published under the GNU GPL (Version 3), check at the LICENSE file
+
+from . import iceflow
\ No newline at end of file
diff --git a/igm/modules/process/avalanche/__init__.py b/igm/modules/process/avalanche/__init__.py
new file mode 100644
index 0000000..c5f32f1
--- /dev/null
+++ b/igm/modules/process/avalanche/__init__.py
@@ -0,0 +1,6 @@
+from .avalanche import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/avalanche.md b/igm/modules/process/avalanche/avalanche.md
similarity index 100%
rename from igm/modules/process/avalanche.md
rename to igm/modules/process/avalanche/avalanche.md
diff --git a/igm/modules/process/avalanche.py b/igm/modules/process/avalanche/avalanche.py
similarity index 95%
rename from igm/modules/process/avalanche.py
rename to igm/modules/process/avalanche/avalanche.py
index 43baead..1cfc2b5 100755
--- a/igm/modules/process/avalanche.py
+++ b/igm/modules/process/avalanche/avalanche.py
@@ -8,7 +8,7 @@
 import time
 
 
-def params_avalanche(parser):
+def params(parser):
     parser.add_argument(
         "--avalanche_update_freq",
         type=float,
@@ -24,12 +24,12 @@ def params_avalanche(parser):
     )
 
 
-def initialize_avalanche(params, state):
+def initialize(params, state):
     state.tcomp_avalanche = []
     state.tlast_avalanche = tf.Variable(params.time_start, dtype=tf.float32)
 
 
-def update_avalanche(params, state):
+def update(params, state):
     if (state.t - state.tlast_avalanche) >= params.avalanche_update_freq:
         if hasattr(state, "logger"):
             state.logger.info("Update AVALANCHE at time : " + str(state.t.numpy()))
@@ -120,5 +120,5 @@ def update_avalanche(params, state):
         state.tcomp_avalanche[-1] *= -1
 
 
-def finalize_avalanche(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/clim_oggm/__init__.py b/igm/modules/process/clim_oggm/__init__.py
new file mode 100644
index 0000000..db1ea82
--- /dev/null
+++ b/igm/modules/process/clim_oggm/__init__.py
@@ -0,0 +1,6 @@
+from .clim_oggm import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/clim_oggm.md b/igm/modules/process/clim_oggm/clim_oggm.md
similarity index 100%
rename from igm/modules/process/clim_oggm.md
rename to igm/modules/process/clim_oggm/clim_oggm.md
diff --git a/igm/modules/process/clim_oggm.py b/igm/modules/process/clim_oggm/clim_oggm.py
similarity index 93%
rename from igm/modules/process/clim_oggm.py
rename to igm/modules/process/clim_oggm/clim_oggm.py
index 8542d81..c2d5d7e 100755
--- a/igm/modules/process/clim_oggm.py
+++ b/igm/modules/process/clim_oggm/clim_oggm.py
@@ -13,7 +13,7 @@
 from igm.modules.utils import interp1d_tf
 
 
-def params_clim_oggm(parser):
+def params(parser):
     # CLIMATE PARAMETERS
     parser.add_argument(
         "--clim_oggm_update_freq",
@@ -51,10 +51,10 @@ def params_clim_oggm(parser):
     )
 
 
-def initialize_clim_oggm(params, state):
+def initialize(params, state):
     # load the given parameters from the json file
     
-    with open(os.path.join(params.working_dir, params.oggm_RGI_ID, "mb_calib.json"), "r") as json_file:
+    with open(os.path.join(params.oggm_RGI_ID, "mb_calib.json"), "r") as json_file:
         jsonString = json_file.read()
 
     oggm_mb_calib = json.loads(jsonString)
@@ -67,7 +67,7 @@ def initialize_clim_oggm(params, state):
 
     # load climate data from netcdf file climate_historical.nc
     nc = Dataset(
-        os.path.join(params.working_dir, params.oggm_RGI_ID, "climate_historical.nc")
+        os.path.join(params.oggm_RGI_ID, "climate_historical.nc")
     )
 
     time = np.squeeze(nc.variables["time"]).astype("float32")  # unit : year
@@ -116,7 +116,7 @@ def initialize_clim_oggm(params, state):
 
     if params.clim_oggm_clim_trend_array == []:
         state.climpar = np.loadtxt(
-            os.path.join(params.working_dir, params.clim_oggm_file),
+            params.clim_oggm_file,
             skiprows=1,
             dtype=np.float32,
         )
@@ -128,7 +128,7 @@ def initialize_clim_oggm(params, state):
     np.random.seed(params.clim_oggm_seed_par)  # fix the seed
 
 
-def update_clim_oggm(params, state):
+def update(params, state):
     if (state.t - state.tlast_clim_oggm) >= params.clim_oggm_update_freq:
         if hasattr(state, "logger"):
             state.logger.info("update climate at time : " + str(state.t.numpy()))
@@ -181,5 +181,5 @@ def update_clim_oggm(params, state):
         state.tcomp_clim_oggm[-1] *= -1
 
 
-def finalize_clim_oggm(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/enthalpy/__init__.py b/igm/modules/process/enthalpy/__init__.py
new file mode 100644
index 0000000..db0493f
--- /dev/null
+++ b/igm/modules/process/enthalpy/__init__.py
@@ -0,0 +1,6 @@
+from .enthalpy import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/enthalpy.md b/igm/modules/process/enthalpy/enthalpy.md
similarity index 100%
rename from igm/modules/process/enthalpy.md
rename to igm/modules/process/enthalpy/enthalpy.md
diff --git a/igm/modules/process/enthalpy.py b/igm/modules/process/enthalpy/enthalpy.py
similarity index 99%
rename from igm/modules/process/enthalpy.py
rename to igm/modules/process/enthalpy/enthalpy.py
index 585dadd..7855580 100755
--- a/igm/modules/process/enthalpy.py
+++ b/igm/modules/process/enthalpy/enthalpy.py
@@ -11,7 +11,7 @@
 from igm.modules.utils import *
 
 
-def params_enthalpy(parser):
+def params(parser):
     parser.add_argument(
         "--enth_water_density",
         type=float,
@@ -137,7 +137,7 @@ def params_enthalpy(parser):
     )
 
 
-def initialize_enthalpy(params, state):
+def initialize(params, state):
     Ny, Nx = state.thk.shape
 
     state.basalMeltRate = tf.Variable(tf.zeros_like(state.thk))
@@ -175,7 +175,7 @@ def initialize_enthalpy(params, state):
     assert params.iflo_dim_arrhenius == 3
 
 
-def update_enthalpy(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Update ENTHALPY at time : " + str(state.t.numpy()))
 
@@ -310,7 +310,7 @@ def update_enthalpy(params, state):
     state.tcomp_enthalpy[-1] *= -1
 
 
-def finalize_enthalpy(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/process/flow_dt_thk.py b/igm/modules/process/flow_dt_thk.py
deleted file mode 100755
index 66d42e6..0000000
--- a/igm/modules/process/flow_dt_thk.py
+++ /dev/null
@@ -1,32 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file
-
-from igm.modules.process.iceflow import *
-from igm.modules.process.time import *
-from igm.modules.process.thk import *
-
-
-def params_flow_dt_thk(parser):
-    params_iceflow(parser)
-    params_time(parser)
-    params_thk(parser)
-
-
-def initialize_flow_dt_thk(params, state):
-    initialize_iceflow(params, state)
-    initialize_time(params, state)
-    initialize_thk(params, state)
-
-
-def update_flow_dt_thk(params, state):
-    update_iceflow(params, state)
-    update_time(params, state)
-    update_thk(params, state)
-
-
-def finalize_flow_dt_thk(params, state):
-    finalize_iceflow(params, state)
-    finalize_time(params, state)
-    finalize_thk(params, state)
diff --git a/igm/modules/process/flow_dt_thk/__init__.py b/igm/modules/process/flow_dt_thk/__init__.py
new file mode 100644
index 0000000..ce19d9a
--- /dev/null
+++ b/igm/modules/process/flow_dt_thk/__init__.py
@@ -0,0 +1,6 @@
+from .flow_dt_thk import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/flow_dt_thk.md b/igm/modules/process/flow_dt_thk/flow_dt_thk.md
similarity index 100%
rename from igm/modules/process/flow_dt_thk.md
rename to igm/modules/process/flow_dt_thk/flow_dt_thk.md
diff --git a/igm/modules/process/flow_dt_thk/flow_dt_thk.py b/igm/modules/process/flow_dt_thk/flow_dt_thk.py
new file mode 100755
index 0000000..7c87a27
--- /dev/null
+++ b/igm/modules/process/flow_dt_thk/flow_dt_thk.py
@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+
+# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
+# Published under the GNU GPL (Version 3), check at the LICENSE file
+
+from igm.modules.process.iceflow import iceflow
+from igm.modules.process.time import time
+from igm.modules.process.thk import thk
+
+
+def params(parser):
+    iceflow.params(parser)
+    time.params(parser)
+    thk.params(parser)
+
+
+def initialize(params, state):
+    iceflow.initialize(params, state)
+    time.initialize(params, state)
+    thk.initialize(params, state)
+
+
+def update(params, state):
+    iceflow.update(params, state)
+    time.update(params, state)
+    thk.update(params, state)
+
+
+def finalize(params, state):
+    iceflow.finalize(params, state)
+    time.finalize(params, state)
+    thk.finalize(params, state)
diff --git a/igm/modules/process/gflex/__init__.py b/igm/modules/process/gflex/__init__.py
new file mode 100644
index 0000000..17377c8
--- /dev/null
+++ b/igm/modules/process/gflex/__init__.py
@@ -0,0 +1,6 @@
+from .gflex import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/gflex.md b/igm/modules/process/gflex/gflex.md
similarity index 100%
rename from igm/modules/process/gflex.md
rename to igm/modules/process/gflex/gflex.md
diff --git a/igm/modules/process/gflex.py b/igm/modules/process/gflex/gflex.py
similarity index 97%
rename from igm/modules/process/gflex.py
rename to igm/modules/process/gflex/gflex.py
index cad85dd..dda6f02 100755
--- a/igm/modules/process/gflex.py
+++ b/igm/modules/process/gflex/gflex.py
@@ -8,7 +8,7 @@
 import time
 
 
-def params_gflex(parser):
+def params(parser):
     parser.add_argument(
         "--gflex_update_freq",
         type=float,
@@ -29,7 +29,7 @@ def params_gflex(parser):
     )
 
 
-def initialize_gflex(params, state):
+def initialize(params, state):
     import gflex
     from gflex.f2d import F2D
 
@@ -63,8 +63,8 @@ def initialize_gflex(params, state):
     if not hasattr(state,"tcomp_gflex"):
         state.flex.Te0 = state.flex.Te    
 
-def update_gflex(params, state):
-    import gflex
+def update(params, state):    
+    import igm.modules.process.gflex.gflex as gflex
     from scipy.interpolate import griddata
     
     initialize_gflex(params, state)
@@ -161,5 +161,5 @@ def pad_arrays(params,state):
         state.tcomp_gflex[-1] *= -1
 
 
-def finalize_gflex(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/glerosion/__init__.py b/igm/modules/process/glerosion/__init__.py
new file mode 100644
index 0000000..c413cf9
--- /dev/null
+++ b/igm/modules/process/glerosion/__init__.py
@@ -0,0 +1,6 @@
+from .glerosion import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/glerosion.md b/igm/modules/process/glerosion/glerosion.md
similarity index 100%
rename from igm/modules/process/glerosion.md
rename to igm/modules/process/glerosion/glerosion.md
diff --git a/igm/modules/process/glerosion.py b/igm/modules/process/glerosion/glerosion.py
similarity index 92%
rename from igm/modules/process/glerosion.py
rename to igm/modules/process/glerosion/glerosion.py
index 175badf..7c61852 100755
--- a/igm/modules/process/glerosion.py
+++ b/igm/modules/process/glerosion/glerosion.py
@@ -11,7 +11,7 @@
 from igm.modules.utils import *
 
 
-def params_glerosion(parser):
+def params(parser):
     parser.add_argument(
         "--glerosion_cst",
         type=float,
@@ -34,12 +34,12 @@ def params_glerosion(parser):
     )
 
 
-def initialize_glerosion(params, state):
+def initialize(params, state):
     state.tcomp_glerosion = []
     state.tlast_erosion = tf.Variable(params.time_start)
 
 
-def update_glerosion(params, state):
+def update(params, state):
     if (state.t - state.tlast_erosion) >= params.glerosion_update_freq:
         if hasattr(state, "logger"):
             state.logger.info(
@@ -64,5 +64,5 @@ def update_glerosion(params, state):
         state.tcomp_glerosion[-1] *= -1
 
 
-def finalize_glerosion(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/iceflow/__init__.py b/igm/modules/process/iceflow/__init__.py
new file mode 100644
index 0000000..6fa6162
--- /dev/null
+++ b/igm/modules/process/iceflow/__init__.py
@@ -0,0 +1,6 @@
+from .iceflow import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/iceflow.md b/igm/modules/process/iceflow/iceflow.md
similarity index 100%
rename from igm/modules/process/iceflow.md
rename to igm/modules/process/iceflow/iceflow.md
diff --git a/igm/modules/process/iceflow.py b/igm/modules/process/iceflow/iceflow.py
similarity index 96%
rename from igm/modules/process/iceflow.py
rename to igm/modules/process/iceflow/iceflow.py
index 8b205b6..61b8ae4 100755
--- a/igm/modules/process/iceflow.py
+++ b/igm/modules/process/iceflow/iceflow.py
@@ -51,7 +51,7 @@
 ############################################
 
 
-def params_iceflow(parser):
+def params(parser):
     # type of ice flow computations
     parser.add_argument(
         "--iflo_type",
@@ -131,7 +131,7 @@ def params_iceflow(parser):
         default=False,
         help="save the iceflow emaultor at the end of the simulation",
     )
-     
+
     # vertical discretization
     parser.add_argument(
         "--iflo_Nz",
@@ -294,7 +294,12 @@ def params_iceflow(parser):
     )
 
 
-def initialize_iceflow(params, state):
+def initialize(params, state):
+    # This makes it so that if the user included the optimize module, this intializer will not be called again.
+    # This is due to the fact that the optimize module calls the initialize (and params) function of the iceflow module.
+    if hasattr(state, "optimize_initializer_called"):
+        return
+
     state.tcomp_iceflow = []
 
     # here we initialize variable parmaetrizing ice flow
@@ -354,13 +359,16 @@ def initialize_iceflow(params, state):
             + "_"
             + str(int(params.iflo_new_friction_param))
         )
- 
+
         if params.iflo_pretrained_emulator:
-            
-            if params.iflo_emulator=="":
-                if os.path.exists(importlib_resources.files(emulators).joinpath(direct_name)):
+            if params.iflo_emulator == "":
+                if os.path.exists(
+                    importlib_resources.files(emulators).joinpath(direct_name)
+                ):
                     dirpath = importlib_resources.files(emulators).joinpath(direct_name)
-                    print("Found pretrained emulator in the igm package: " + direct_name)
+                    print(
+                        "Found pretrained emulator in the igm package: " + direct_name
+                    )
                 else:
                     print("No pretrained emulator found in the igm package")
             else:
@@ -369,7 +377,7 @@ def initialize_iceflow(params, state):
                     print("Found pretrained emulator: " + params.iflo_emulator)
                 else:
                     print("No pretrained emulator found ")
-  
+
             fieldin = []
             fid = open(os.path.join(dirpath, "fieldin.dat"), "r")
             for fileline in fid:
@@ -386,9 +394,13 @@ def initialize_iceflow(params, state):
                 params.iflo_Nz - 1
             )
             nb_outputs = 2 * params.iflo_Nz
-            state.iceflow_model = getattr(igm, params.iflo_network)(
-                params, nb_inputs, nb_outputs
-            )
+            # state.iceflow_model = getattr(igm, params.iflo_network)(
+            #     params, nb_inputs, nb_outputs
+            # )
+            if params.iflo_network=='cnn':
+                state.iceflow_model = cnn(params, nb_inputs, nb_outputs)
+            elif params.iflo_network=='unet':
+                state.iceflow_model = unet(params, nb_inputs, nb_outputs)
 
     if not params.iflo_type == "emulated":
         if int(tf.__version__.split(".")[1]) <= 10:
@@ -431,7 +443,7 @@ def initialize_iceflow(params, state):
         _update_iceflow_emulated(params, state)
 
 
-def update_iceflow(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Update ICEFLOW at time : " + str(state.t.numpy()))
 
@@ -453,10 +465,9 @@ def update_iceflow(params, state):
     state.tcomp_iceflow[-1] *= -1
 
 
-def finalize_iceflow(params, state):
-
+def finalize(params, state):
     if params.iflo_save_model:
-        save_iceflow_model(params, state) 
+        save_iceflow_model(params, state)
 
 
 ########################################################################
@@ -1080,7 +1091,7 @@ def _update_iceflow_diagnostic(params, state):
 
         print(ERR)
 
-        with open(os.path.join(params.working_dir, "errors.txt"), "ab") as f:
+        with open("errors.txt", "ab") as f:
             np.savetxt(f, np.expand_dims(ERR, axis=0), delimiter=",", fmt="%5.5f")
 
 
@@ -1180,7 +1191,7 @@ def unet(params, nb_inputs, nb_outputs):
 
 
 def save_iceflow_model(params, state):
-    directory = os.path.join(params.working_dir, "iceflow-model")
+    directory = "iceflow-model"
 
     if os.path.exists(directory):
         shutil.rmtree(directory)
diff --git a/igm/modules/process/iceflow_v1/__init__.py b/igm/modules/process/iceflow_v1/__init__.py
new file mode 100644
index 0000000..31d823f
--- /dev/null
+++ b/igm/modules/process/iceflow_v1/__init__.py
@@ -0,0 +1,6 @@
+from .iceflow_v1 import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/iceflow_v1.md b/igm/modules/process/iceflow_v1/iceflow_v1.md
similarity index 100%
rename from igm/modules/process/iceflow_v1.md
rename to igm/modules/process/iceflow_v1/iceflow_v1.md
diff --git a/igm/modules/process/iceflow_v1.py b/igm/modules/process/iceflow_v1/iceflow_v1.py
similarity index 97%
rename from igm/modules/process/iceflow_v1.py
rename to igm/modules/process/iceflow_v1/iceflow_v1.py
index c585fb6..fab31b8 100755
--- a/igm/modules/process/iceflow_v1.py
+++ b/igm/modules/process/iceflow_v1/iceflow_v1.py
@@ -15,7 +15,7 @@
 from igm import emulators
 
 
-def params_iceflow_v1(parser):
+def params(parser):
     parser.add_argument(
         "--init_strflowctrl",
         type=float,
@@ -56,7 +56,7 @@ def params_iceflow_v1(parser):
     )
 
 
-def initialize_iceflow_v1(params, state):
+def initialize(params, state):
     state.tcomp_iceflow = []
 
     # here we initialize variable parmaetrizing ice flow
@@ -112,7 +112,7 @@ def initialize_iceflow_v1(params, state):
         state.PAD = [[0, 0], [0, 0]]
 
 
-def update_iceflow_v1(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Update ICEFLOW at time : " + str(state.t.numpy()))
 
@@ -166,7 +166,7 @@ def update_iceflow_v1(params, state):
     state.tcomp_iceflow[-1] *= -1
 
 
-def finalize_iceflow_v1(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/process/particles/__init__.py b/igm/modules/process/particles/__init__.py
new file mode 100644
index 0000000..be05c7a
--- /dev/null
+++ b/igm/modules/process/particles/__init__.py
@@ -0,0 +1,6 @@
+from .particles import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/particles.md b/igm/modules/process/particles/particles.md
similarity index 100%
rename from igm/modules/process/particles.md
rename to igm/modules/process/particles/particles.md
diff --git a/igm/modules/process/particles.py b/igm/modules/process/particles/particles.py
similarity index 97%
rename from igm/modules/process/particles.py
rename to igm/modules/process/particles/particles.py
index 17c53b1..fd97d14 100755
--- a/igm/modules/process/particles.py
+++ b/igm/modules/process/particles/particles.py
@@ -12,8 +12,9 @@
 
 from igm.modules.utils import *
 
+from igm.modules.process.particles_v1.particles_v1 import seeding_particles
 
-def params_particles(parser):
+def params(parser):
     parser.add_argument(
         "--part_tracking_method",
         type=str,
@@ -34,7 +35,7 @@ def params_particles(parser):
     )
 
 
-def initialize_particles(params, state):
+def initialize(params, state):
     state.tlast_seeding = -1.0e5000
     state.tcomp_particles = []
 
@@ -53,12 +54,12 @@ def initialize_particles(params, state):
     state.gridseed[::rr, ::rr] = True
 
 
-def update_particles(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Update particle tracking at time : " + str(state.t.numpy()))
 
     if (state.t.numpy() - state.tlast_seeding) >= params.part_frequency_seeding:
-        igm.seeding_particles(params, state)
+        seeding_particles(params, state)
 
         # merge the new seeding points with the former ones
         state.xpos = tf.Variable(tf.concat([state.xpos, state.nxpos], axis=-1))
@@ -192,7 +193,7 @@ def update_particles(params, state):
         state.tcomp_particles[-1] *= -1
 
 
-def finalize_particles(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/process/particles_v1/__init__.py b/igm/modules/process/particles_v1/__init__.py
new file mode 100644
index 0000000..b8a6a58
--- /dev/null
+++ b/igm/modules/process/particles_v1/__init__.py
@@ -0,0 +1,6 @@
+from .particles_v1 import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/particles_v1.md b/igm/modules/process/particles_v1/particles_v1.md
similarity index 100%
rename from igm/modules/process/particles_v1.md
rename to igm/modules/process/particles_v1/particles_v1.md
diff --git a/igm/modules/process/particles_v1.py b/igm/modules/process/particles_v1/particles_v1.py
similarity index 98%
rename from igm/modules/process/particles_v1.py
rename to igm/modules/process/particles_v1/particles_v1.py
index 3639335..ccf942c 100755
--- a/igm/modules/process/particles_v1.py
+++ b/igm/modules/process/particles_v1/particles_v1.py
@@ -12,7 +12,7 @@
 from igm.modules.utils import *
 
 
-def params_particles_v1(parser):
+def params(parser):
     parser.add_argument(
         "--part_tracking_method",
         type=str,
@@ -33,7 +33,7 @@ def params_particles_v1(parser):
     )
 
 
-def initialize_particles_v1(params, state):
+def initialize(params, state):
     state.tlast_seeding = -1.0e5000
     state.tcomp_particles = []
 
@@ -52,7 +52,7 @@ def initialize_particles_v1(params, state):
     state.gridseed[::rr, ::rr] = True
 
 
-def update_particles_v1(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info("Update particle tracking at time : " + str(state.t.numpy()))
 
@@ -226,7 +226,7 @@ def update_particles_v1(params, state):
     state.tcomp_particles[-1] *= -1
 
 
-def finalize_particles_v1(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/process/read_output.py b/igm/modules/process/read_output.py
index b40e6ce..e286ab2 100755
--- a/igm/modules/process/read_output.py
+++ b/igm/modules/process/read_output.py
@@ -50,7 +50,7 @@ def params_read_output(parser):
  
 def initialize_read_output(params, state):
  
-    nc = Dataset(os.path.join(params.working_dir, params.rncd_input_file), "r")
+    nc = Dataset(params.rncd_input_file, "r")
 
     state.x    = np.squeeze(nc.variables["x"]).astype("float32")
     state.y    = np.squeeze(nc.variables["y"]).astype("float32")
diff --git a/igm/modules/process/rockflow/__init__.py b/igm/modules/process/rockflow/__init__.py
new file mode 100644
index 0000000..6a7ca53
--- /dev/null
+++ b/igm/modules/process/rockflow/__init__.py
@@ -0,0 +1,6 @@
+from .rockflow import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/rockflow.md b/igm/modules/process/rockflow/rockflow.md
similarity index 100%
rename from igm/modules/process/rockflow.md
rename to igm/modules/process/rockflow/rockflow.md
diff --git a/igm/modules/process/rockflow.py b/igm/modules/process/rockflow/rockflow.py
similarity index 87%
rename from igm/modules/process/rockflow.py
rename to igm/modules/process/rockflow/rockflow.py
index b4a8899..fe1b071 100755
--- a/igm/modules/process/rockflow.py
+++ b/igm/modules/process/rockflow/rockflow.py
@@ -13,7 +13,7 @@
 from igm.modules.utils import *
 
 
-def params_rockflow(parser):
+def params(parser):
     parser.add_argument(
         "--rock_flow_speed",
         type=float,
@@ -22,11 +22,11 @@ def params_rockflow(parser):
     )
 
 
-def initialize_rockflow(params, state):
+def initialize(params, state):
     pass
 
 
-def update_rockflow(params, state):
+def update(params, state):
     slopsurfx, slopsurfy = compute_gradient_tf(state.usurf, state.dx, state.dx)
 
     slop = getmag(slopsurfx, slopsurfy)
@@ -44,5 +44,5 @@ def update_rockflow(params, state):
     state.V.assign(tf.where(thkexp > 0, state.V, diry))
 
 
-def finalize_rockflow(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/smb_oggm/__init__.py b/igm/modules/process/smb_oggm/__init__.py
new file mode 100644
index 0000000..ff24323
--- /dev/null
+++ b/igm/modules/process/smb_oggm/__init__.py
@@ -0,0 +1,6 @@
+from .smb_oggm import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/smb_oggm.md b/igm/modules/process/smb_oggm/smb_oggm.md
similarity index 100%
rename from igm/modules/process/smb_oggm.md
rename to igm/modules/process/smb_oggm/smb_oggm.md
diff --git a/igm/modules/process/smb_oggm.py b/igm/modules/process/smb_oggm/smb_oggm.py
similarity index 93%
rename from igm/modules/process/smb_oggm.py
rename to igm/modules/process/smb_oggm/smb_oggm.py
index 456cc65..539942b 100755
--- a/igm/modules/process/smb_oggm.py
+++ b/igm/modules/process/smb_oggm/smb_oggm.py
@@ -11,7 +11,7 @@
 import json
 
 
-def params_smb_oggm(parser):
+def params(parser):
     parser.add_argument(
         "--smb_oggm_update_freq",
         type=float,
@@ -32,12 +32,12 @@ def params_smb_oggm(parser):
     )
 
 
-def initialize_smb_oggm(params, state):
+def initialize(params, state):
     state.tcomp_smb_oggm = []
     state.tlast_mb = tf.Variable(-1.0e5000)
 
     # load the given parameters from the json file
-    with open(os.path.join(params.working_dir, params.oggm_RGI_ID, "mb_calib.json"), "r") as json_file:
+    with open(os.path.join(params.oggm_RGI_ID, "mb_calib.json"), "r") as json_file:
         jsonString = json_file.read()
 
     oggm_mb_calib = json.loads(jsonString)
@@ -48,7 +48,7 @@ def initialize_smb_oggm(params, state):
     params.melt_f = oggm_mb_calib["melt_f"]  # unit: mm water / (celcius day)
 
 
-def update_smb_oggm(params, state):
+def update(params, state):
     #    mass balance forced by climate with accumulation and temperature-index melt model
     #    Input:  state.precipitation [Unit: kg * m^(-2) * y^(-1)]
     #            state.air_temp      [Unit: °C           ]
@@ -109,5 +109,5 @@ def update_smb_oggm(params, state):
         state.tcomp_smb_oggm[-1] *= -1
 
 
-def finalize_smb_oggm(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/smb_simple/__init__.py b/igm/modules/process/smb_simple/__init__.py
new file mode 100644
index 0000000..5a08db3
--- /dev/null
+++ b/igm/modules/process/smb_simple/__init__.py
@@ -0,0 +1,6 @@
+from .smb_simple import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/smb_simple.md b/igm/modules/process/smb_simple/smb_simple.md
similarity index 100%
rename from igm/modules/process/smb_simple.md
rename to igm/modules/process/smb_simple/smb_simple.md
diff --git a/igm/modules/process/smb_simple.py b/igm/modules/process/smb_simple/smb_simple.py
similarity index 91%
rename from igm/modules/process/smb_simple.py
rename to igm/modules/process/smb_simple/smb_simple.py
index 56024e7..4357808 100755
--- a/igm/modules/process/smb_simple.py
+++ b/igm/modules/process/smb_simple/smb_simple.py
@@ -10,7 +10,7 @@
 from igm.modules.utils import interp1d_tf
 
 
-def params_smb_simple(parser):
+def params(parser):
     parser.add_argument(
         "--smb_simple_update_freq",
         type=float,
@@ -31,10 +31,10 @@ def params_smb_simple(parser):
     )
 
 
-def initialize_smb_simple(params, state):
+def initialize(params, state):
     if params.smb_simple_array == []:
         state.smbpar = np.loadtxt(
-            os.path.join(params.working_dir, params.smb_simple_file),
+            params.smb_simple_file,
             skiprows=1,
             dtype=np.float32,
         )
@@ -45,7 +45,7 @@ def initialize_smb_simple(params, state):
     state.tlast_mb = tf.Variable(-1.0e5000)
 
 
-def update_smb_simple(params, state):
+def update(params, state):
     # update smb each X years
     if (state.t - state.tlast_mb) >= params.smb_simple_update_freq:
         if hasattr(state, "logger"):
@@ -78,5 +78,5 @@ def update_smb_simple(params, state):
         state.tcomp_smb_simple[-1] *= -1
 
 
-def finalize_smb_simple(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/thk/__init__.py b/igm/modules/process/thk/__init__.py
new file mode 100644
index 0000000..caa3e9d
--- /dev/null
+++ b/igm/modules/process/thk/__init__.py
@@ -0,0 +1,6 @@
+from .thk import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/thk.md b/igm/modules/process/thk/thk.md
similarity index 100%
rename from igm/modules/process/thk.md
rename to igm/modules/process/thk/thk.md
diff --git a/igm/modules/process/thk.py b/igm/modules/process/thk/thk.py
similarity index 94%
rename from igm/modules/process/thk.py
rename to igm/modules/process/thk/thk.py
index 93e3e7a..a08a97d 100755
--- a/igm/modules/process/thk.py
+++ b/igm/modules/process/thk/thk.py
@@ -8,7 +8,7 @@
 
 from igm.modules.utils import compute_divflux_slope_limiter
 
-def params_thk(parser):
+def params(parser):
     parser.add_argument(
         "--thk_slope_type",
         type=str,
@@ -22,7 +22,7 @@ def params_thk(parser):
         help="density of ice divided by density of water",
     )
 
-def initialize_thk(params, state):
+def initialize(params, state):
 
     # define the lower ice surface
     if hasattr(state, "sealevel"):
@@ -35,7 +35,7 @@ def initialize_thk(params, state):
 
     state.tcomp_thk = []
 
-def update_thk(params, state):
+def update(params, state):
     if state.it >= 0:
         if hasattr(state, "logger"):
             state.logger.info(
@@ -69,5 +69,5 @@ def update_thk(params, state):
         state.tcomp_thk[-1] *= -1
 
 
-def finalize_thk(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/time/__init__.py b/igm/modules/process/time/__init__.py
new file mode 100644
index 0000000..b51d13a
--- /dev/null
+++ b/igm/modules/process/time/__init__.py
@@ -0,0 +1,6 @@
+from .time import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/time.md b/igm/modules/process/time/time.md
similarity index 100%
rename from igm/modules/process/time.md
rename to igm/modules/process/time/time.md
diff --git a/igm/modules/process/time.py b/igm/modules/process/time/time.py
similarity index 95%
rename from igm/modules/process/time.py
rename to igm/modules/process/time/time.py
index 6d7748c..6bd6199 100755
--- a/igm/modules/process/time.py
+++ b/igm/modules/process/time/time.py
@@ -9,7 +9,7 @@
 import tensorflow as tf
 
 
-def params_time(parser):
+def params(parser):
     parser.add_argument(
         "--time_start",
         type=float,
@@ -39,7 +39,7 @@ def params_time(parser):
     )
 
 
-def initialize_time(params, state):
+def initialize(params, state):
     state.tcomp_time = []
 
     # Initialize the time with starting time
@@ -62,7 +62,7 @@ def initialize_time(params, state):
     state.saveresult = True
 
 
-def update_time(params, state):
+def update(params, state):
     if hasattr(state, "logger"):
         state.logger.info(
             "Update DT from the CFL condition at time : " + str(state.t.numpy())
@@ -104,5 +104,5 @@ def update_time(params, state):
     state.tcomp_time[-1] *= -1
 
 
-def finalize_time(params, state):
+def finalize(params, state):
     pass
diff --git a/igm/modules/process/vert_flow/__init__.py b/igm/modules/process/vert_flow/__init__.py
new file mode 100644
index 0000000..4ee7783
--- /dev/null
+++ b/igm/modules/process/vert_flow/__init__.py
@@ -0,0 +1,6 @@
+from .vert_flow import (
+    params,
+    initialize,
+    finalize,
+    update
+)
diff --git a/igm/modules/process/vert_flow.md b/igm/modules/process/vert_flow/vert_flow.md
similarity index 100%
rename from igm/modules/process/vert_flow.md
rename to igm/modules/process/vert_flow/vert_flow.md
diff --git a/igm/modules/process/vert_flow.py b/igm/modules/process/vert_flow/vert_flow.py
similarity index 92%
rename from igm/modules/process/vert_flow.py
rename to igm/modules/process/vert_flow/vert_flow.py
index 64929ee..69ea66d 100755
--- a/igm/modules/process/vert_flow.py
+++ b/igm/modules/process/vert_flow/vert_flow.py
@@ -10,18 +10,18 @@
 from igm.modules.utils import *
 
 
-def params_vert_flow(parser):
+def params(parser):
     pass
 
 
-def initialize_vert_flow(params, state):
+def initialize(params, state):
     state.tcomp_vert_flow = []
 
 
 #    state.W = tf.zeros_like(state.U[0])
 
 
-def update_vert_flow(params, state):
+def update(params, state):
     """ """
 
     state.tcomp_vert_flow.append(time.time())
@@ -37,7 +37,7 @@ def update_vert_flow(params, state):
     state.tcomp_vert_flow[-1] *= -1
 
 
-def finalize_vert_flow(params, state):
+def finalize(params, state):
     pass
 
 
diff --git a/igm/modules/utils.py b/igm/modules/utils.py
index 5436e32..f6404c2 100755
--- a/igm/modules/utils.py
+++ b/igm/modules/utils.py
@@ -7,11 +7,19 @@
 This util files provides a large number of usefull function for other modules
 """
 
-import numpy as np
-import os, sys, shutil
-import matplotlib.pyplot as plt
 import tensorflow as tf
 
+# constrains wildcard imports - not necessary but for security reasons - updtae if a utility is added...
+__all__ = [
+    "str2bool",
+    "getmag",
+    "getmag3d",
+    "compute_gradient_tf",
+    "compute_divflux",
+    "interp1d_tf",
+    "complete_data",
+    "interpolate_bilinear_tf"
+]
 
 def str2bool(v):
     return v.lower() in ("true", "1")


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..59aa9ac
--- /dev/null
+++ b/tests/__init__.py
@@ -0,0 +1,3 @@
+import sys
+import os
+sys.path.append(os.getcwd()) 
\ No newline at end of file
diff --git a/tests/enthalpy/test-Hewitt/README.md b/tests/enthalpy/test-Hewitt/README.md
deleted file mode 100755
index a88d385..0000000
--- a/tests/enthalpy/test-Hewitt/README.md
+++ /dev/null
@@ -1,28 +0,0 @@
-[![License badge](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-### <h1 align="center" id="title"> Test of the Enthalpy model </h1>
-
-This setup permits to lead the "ice cap experiment" from Hewitt and Schoof (2017) running ''igm_run''.
-[PoLIM](https://github.com/WangYuzhe/PoLIM-Polythermal-Land-Ice-Model) and Wang and al. (2020) greatly helped to lead this experiment.
-
-# Reference
-
-@article{hewitt2017models,
-  title={Models for polythermal ice sheets and glaciers},
-  author={Hewitt, Ian J and Schoof, Christian},
-  journal={The Cryosphere},
-  volume={11},
-  number={1},
-  pages={541--551},
-  year={2017},
-  publisher={Copernicus GmbH}
-}
-
-@article{wang2020two,
-  title={A two-dimensional, higher-order, enthalpy-based thermomechanical ice flow model for mountain glaciers and its benchmark experiments},
-  author={Wang, Yuzhe and Zhang, Tong and Xiao, Cunde and Ren, Jiawen and Wang, Yanfen},
-  journal={Computers \& geosciences},
-  volume={141},
-  pages={104526},
-  year={2020},
-  publisher={Elsevier}
-}
diff --git a/tests/enthalpy/test-Hewitt/cap_enth_exp.py b/tests/enthalpy/test-Hewitt/cap_enth_exp.py
deleted file mode 100755
index c350be3..0000000
--- a/tests/enthalpy/test-Hewitt/cap_enth_exp.py
+++ /dev/null
@@ -1,133 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file
-
-import os
-import numpy as np
-import matplotlib.pyplot as plt
-import tensorflow as tf
-
-from igm.modules.utils import complete_data, compute_gradient_tf
-from igm.modules.process.enthalpy import vertically_discretize_tf
-import igm
-
-
-def params_cap_enth_exp(parser):
-    igm.params_time(parser)
-    igm.params_iceflow(parser)
-
-
-def initialize_cap_enth_exp(params, state):
-    # Initialize the time with starting time
-    state.t = tf.Variable(float(params.time_start))
-
-    state.it = 0
-
-    state.dt = tf.Variable(float(params.time_step_max))
-
-    state.dt_target = tf.Variable(float(params.time_step_max))
-
-    state.time_save = np.ndarray.tolist(
-        np.arange(params.time_start, params.time_end, params.time_save)
-    ) + [params.time_end]
-
-    state.time_save = tf.constant(state.time_save, dtype="float32")
-
-    state.itsave = 0
-
-    state.saveresult = True
-
-    x = np.arange(-100, 101) * 1000  # make x-axis, lenght 100 km,
-    y = np.arange(-100, 101) * 1000  # make y-axis, lenght 100 km,
-
-    X, Y = np.meshgrid(x, y)
-    R = np.sqrt(X**2 + Y**2)
-
-    topg = np.zeros_like(X)
-    usurf = np.maximum(1500 * (1 - (R / 100000) ** 2), topg)
-    thk = usurf
-
-    state.x = tf.constant(x.astype("float32"))
-    state.y = tf.constant(y.astype("float32"))
-
-    state.topg = tf.Variable(topg.astype("float32"))
-    state.thk = tf.Variable(thk.astype("float32"))
-
-    complete_data(state)
-
-    state.depth, state.dz = vertically_discretize_tf(
-        state.thk, params.iflo_Nz, params.iflo_vert_spacing
-    )
-
-    rho = 910.0  # ice density (g/m^3)
-    g = 9.81  # earth's gravity (m/s^2)
-    A = 78 * 10 ** (-18)  # rate factor (Pa^-3 y^-1)
-    n = 3
-    k = A * (rho * g) ** n  # Generic constant (y^-1 m^-3)
-
-    state.U = tf.Variable(
-        tf.zeros((params.iflo_Nz, state.thk.shape[0], state.thk.shape[1]))
-    )
-    state.V = tf.Variable(
-        tf.zeros((params.iflo_Nz, state.thk.shape[0], state.thk.shape[1]))
-    )
-
-    thk = tf.tile(tf.expand_dims(state.thk, 0), (params.iflo_Nz, 1, 1))
-
-    slopsurfx, slopsurfy = compute_gradient_tf(state.usurf, state.dx, state.dx)
-
-    dsx = tf.tile(tf.expand_dims(slopsurfx, 0), (params.iflo_Nz, 1, 1))
-    dsy = tf.tile(tf.expand_dims(slopsurfy, 0), (params.iflo_Nz, 1, 1))
-
-    ds = (dsx**2 + dsy**2) ** 0.5
-
-    state.U.assign(
-        -2
-        * A
-        * ((rho * g) ** n)
-        * ds ** (n - 1)
-        * dsx
-        * (thk ** (n + 1) - state.depth ** (n + 1))
-    )
-    state.V.assign(
-        -2
-        * A
-        * ((rho * g) ** n)
-        * ds ** (n - 1)
-        * dsy
-        * (thk ** (n + 1) - state.depth ** (n + 1))
-    )
-
-    state.air_temp = tf.expand_dims(tf.ones_like(state.usurf), axis=0) * (-10)
-
-    state.arrhenius = tf.Variable(
-        tf.ones((params.iflo_Nz, state.thk.shape[0], state.thk.shape[1]))
-        * params.iflo_init_arrhenius
-    )
-
-    state.slidingco = tf.Variable(tf.ones_like(state.thk) * params.iflo_init_slidingco)
-
-    igm.define_vertical_weight(params, state)
-
-    igm.update_2d_iceflow_variables(params, state)
-
-
-def update_cap_enth_exp(params, state):
-    state.dt = 1.0
-
-    # modify dt such that times of requested savings are reached exactly
-    if state.time_save[state.itsave + 1] <= state.t + state.dt:
-        state.dt = state.time_save[state.itsave + 1] - state.t
-        state.saveresult = True
-        state.itsave += 1
-    else:
-        state.saveresult = False
-
-    state.t.assign(state.t + state.dt)
-
-    state.it += 1
-
-
-def finalize_cap_enth_exp(params, state):
-    pass
diff --git a/tests/enthalpy/test-Hewitt/params.json b/tests/enthalpy/test-Hewitt/params.json
deleted file mode 100755
index 9f4d398..0000000
--- a/tests/enthalpy/test-Hewitt/params.json
+++ /dev/null
@@ -1,14 +0,0 @@
-{
-  "modules_preproc": [],
-  "modules_process": ["cap_enth_exp","vert_flow","enthalpy"],
-  "modules_postproc": ["plotsection","write_ncdf","print_info"],
-  "iflo_dim_arrhenius": 3,
-  "time_start": 0.0,
-  "time_end": 3000.0,
-  "wncd_vars_to_save":  ["topg", "usurf", "thk", "smb", "frictheat", "T","basalMeltRate","ubar","vbar",
-                          "tillwat", "strainheat","arrhenius","tauc","W"],
-  "time_step_max": 2.0,
-  "time_save": 100.0,
-  "plt2d_live": true,
-  "enth_drain_ice_column": true
-}
diff --git a/tests/enthalpy/test-Hewitt/plotsection.py b/tests/enthalpy/test-Hewitt/plotsection.py
deleted file mode 100755
index f2015de..0000000
--- a/tests/enthalpy/test-Hewitt/plotsection.py
+++ /dev/null
@@ -1,97 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file
-
-import numpy as np
-import os, sys, shutil
-import matplotlib.pyplot as plt
-import matplotlib
-import datetime, time
-import tensorflow as tf
-
-from igm.modules.utils import *
-
-
-def params_plotsection(parser):
-    parser.add_argument(
-        "--plt2d_editor",
-        type=str,
-        default="vs",
-        help="Optimized for VS code (vs) or spyder (sp) for live plot",
-    )
-    parser.add_argument(
-        "--plt2d_live",
-        type=str2bool,
-        default=False,
-        help="Display plots live the results during computation instead of making png",
-    )
-
-
-def initialize_plotsection(params, state):
-    state.extent = [np.min(state.x), np.max(state.x), np.min(state.y), np.max(state.y)]
-
-    if params.plt2d_editor == "vs":
-        plt.ion()  # enable interactive mode
-
-    state.fig = plt.figure(dpi=200)
-    state.ax = state.fig.add_subplot(1, 1, 1)
-    state.ax.axis("off")
-    state.ax.set_aspect("equal")
-
-    os.system(
-        "echo rm "
-        + os.path.join(params.working_dir, "temperature*.png")
-        + " >> clean.sh"
-    )
-
-
-def update_plotsection(params, state):
-    if state.saveresult:
-        abs = tf.tile(tf.expand_dims(state.X, 0), (params.iflo_Nz, 1, 1)).numpy()
-        ele = np.expand_dims(state.thk, 0) - state.depth.numpy()
-        tem = (state.T - 273.15).numpy()
-        ome = state.omega.numpy()
-
-        X = np.ndarray.flatten(abs[:, 100, 100:]) / 25
-        Z = np.ndarray.flatten(ele[:, 100, 100:])
-        T = np.ndarray.flatten(tem[:, 100, 100:])
-        O = np.ndarray.flatten(ome[:, 100, 100:])
-
-        im = state.ax.scatter(
-            X, Z, c=T, s=0.2, vmin=-10, vmax=0, cmap=matplotlib.cm.jet
-        )
-        #        im = state.ax.scatter( X, Z, c=O, s=0.2, vmin=0, vmax=0.03, cmap=matplotlib.cm.viridis )
-
-        #        XX = state.X[100,100:].numpy()
-        #        YY = state.basalMeltRate[100,100:].numpy()
-        #        im = state.ax.plot( XX, YY ,'-k') ; state.already_set_cbar = True
-
-        state.ax.set_title("YEAR : " + str(state.t.numpy()), size=15)
-
-        if not hasattr(state, "already_set_cbar"):
-            state.cbar = plt.colorbar(im, label="temperature")
-            state.already_set_cbar = True
-
-        if params.plt2d_live:
-            if params.plt2d_editor == "vs":
-                state.fig.canvas.draw()  # re-drawing the figure
-                state.fig.canvas.flush_events()  # to flush the GUI events
-            else:
-                from IPython.display import display, clear_output
-
-                clear_output(wait=True)
-                display(state.fig)
-        else:
-            plt.savefig(
-                os.path.join(
-                    params.working_dir,
-                    "temperature-" + str(state.t.numpy()).zfill(4) + ".png",
-                ),
-                bbox_inches="tight",
-                pad_inches=0.2,
-            )
-
-
-def finalize_plotsection(params, state):
-    pass
diff --git a/tests/enthalpy/test-Kleiner/README.md b/tests/enthalpy/test-Kleiner/README.md
deleted file mode 100755
index cd00526..0000000
--- a/tests/enthalpy/test-Kleiner/README.md
+++ /dev/null
@@ -1,37 +0,0 @@
-[![License badge](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0)
-### <h1 align="center" id="title"> Test of the Enthalpy model </h1>
-
-This setup permits to lead experiments A and B from Kleiner and al. (2015), and to verify the results with the analalytical solution. For that purpose, you need to type
-
-```bash
-python expA.py
-python expB.py
-```
-
-We are comparing the analytical solution to the data folder was taken from [PoLIM](https://github.com/WangYuzhe/PoLIM-Polythermal-Land-Ice-Model) described by Wang and al. (2020), which greatly helped to lead thes experiments.
-
-# Reference
-
-@article{kleiner2015enthalpy,
-  title={Enthalpy benchmark experiments for numerical ice sheet models},
-  author={Kleiner, Thomas and R{\"u}ckamp, Martin and Bondzio, Johannes H and Humbert, Angelika},
-  journal={The Cryosphere},
-  volume={9},
-  number={1},
-  pages={217--228},
-  year={2015},
-  publisher={Copernicus GmbH}
-}
-
-@article{wang2020two,
-  title={A two-dimensional, higher-order, enthalpy-based thermomechanical ice flow model for mountain glaciers and its benchmark experiments},
-  author={Wang, Yuzhe and Zhang, Tong and Xiao, Cunde and Ren, Jiawen and Wang, Yanfen},
-  journal={Computers \& geosciences},
-  volume={141},
-  pages={104526},
-  year={2020},
-  publisher={Elsevier}
-}
-
-
-
diff --git a/tests/enthalpy/test-Kleiner/expA.py b/tests/enthalpy/test-Kleiner/expA.py
deleted file mode 100755
index 3e60910..0000000
--- a/tests/enthalpy/test-Kleiner/expA.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file
-
-import scipy.io
-import tensorflow as tf
-import igm
-
-from igm.modules.process.enthalpy import *
-
-# this file is avilable at https://github.com/WangYuzhe/PoLIM-Polythermal-Land-Ice-Model
-# verif = scipy.io.loadmat("sol_analytic/enthA_analy_result.mat")
-
-ttf = 300000.0  # 300000
-dt = 50.0
-
-tim = np.arange(0, ttf, dt) + dt  # to put back to 300000
-
-parser = igm.params_core()
-igm.params_iceflow(parser)
-igm.params_enthalpy = params_enthalpy
-igm.params_enthalpy(parser)
-params = parser.parse_args(args=[])
-params.iflo_Nz = 50
-params.iflo_vert_spacing = 1
-
-params.enth_KtdivKc = 10 ** (-5)  # check this value if ok ?
-params.enth_till_wat_max = 200
-params.enth_drain_rate = 0
-
-thk = tf.Variable(1000 * tf.ones((1, 1)))
-
-depth, dz = vertically_discretize_tf(thk, params.iflo_Nz, params.iflo_vert_spacing)
-
-strainheat = tf.Variable(tf.zeros((params.iflo_Nz, 1, 1)))
-frictheat = tf.Variable(0.0 * tf.ones((1, 1)))
-geoheatflux = tf.Variable(0.042 * tf.ones((1, 1)))  # in W m-2
-tillwat = tf.Variable(0.0 * tf.ones((1, 1)))
-
-# Initial enthalpy field
-T = tf.Variable((-30.0 + 273.15) * tf.ones((params.iflo_Nz, 1, 1)))
-E = tf.Variable(params.enth_ci * (T - 223.15))
-omega = tf.Variable(tf.zeros_like(T))
-w = tf.Variable(tf.zeros_like(T))
-
-
-TB = []
-HW = []
-MR = []
-
-surftemp = tf.Variable((-30 + 273.15) * tf.ones((1, 1)))
-
-for it, t in enumerate(tim):
-    if (t >= 100000.0) & (t < 150000.0):
-        surftemp.assign((-5 + 273.15) * tf.ones((1, 1)))
-    else:
-        surftemp.assign((-30 + 273.15) * tf.ones((1, 1)))
-
-    Tpmp, Epmp = TpmpEpmp_from_depth_tf(
-        depth,
-        params.iflo_gravity_cst,
-        params.iflo_ice_density,
-        params.enth_claus_clape,
-        params.enth_melt_temp,
-        params.enth_ci,
-        params.enth_ref_temp,
-    )
-
-    surfenth = surf_enthalpy_from_temperature_tf(
-        surftemp, params.enth_melt_temp, params.enth_ci, params.enth_ref_temp
-    )
-
-    E, basalMeltRate = compute_enthalpy_basalmeltrate(
-        E,
-        Epmp,
-        dt * params.enth_spy,
-        dz,
-        w,
-        surfenth,
-        geoheatflux,
-        strainheat,
-        frictheat,
-        tillwat,
-        params.iflo_thr_ice_thk,
-        params.enth_ki,
-        params.iflo_ice_density,
-        params.enth_water_density,
-        params.enth_ci,
-        params.enth_ref_temp,
-        params.enth_Lh,
-        params.enth_spy,
-        params.enth_KtdivKc,
-        params.enth_drain_ice_column,
-    )
-
-    T, omega = temperature_from_enthalpy_tf(
-        E, Tpmp, Epmp, params.enth_ci, params.enth_ref_temp, params.enth_Lh
-    )
-
-    tillwat = tillwat + dt * (basalMeltRate - params.enth_drain_rate)
-    tillwat = tf.clip_by_value(tillwat, 0, params.enth_till_wat_max)
-
-    TB.append(T[0] - 273.15)
-    HW.append(tillwat)
-    MR.append(basalMeltRate)
-
-    if it % 100 == 0:
-        print(
-            "time :",
-            t,
-            T[0, 0, 0].numpy(),
-            tillwat[0, 0].numpy(),
-            basalMeltRate[0, 0].numpy(),
-        )
-
-fig = plt.figure(figsize=(8, 8))
-
-plt.subplot(311)
-plt.plot(tim, np.stack(TB)[:, 0, 0])
-plt.ylabel("Temperature")
-
-plt.subplot(312)
-plt.plot(tim, np.stack(MR)[:, 0, 0])
-# plt.plot(verif["basalMelt"][1], verif["basalMelt"][0] / 1000, "-r")
-plt.ylabel("Melt rate")
-
-plt.subplot(313)
-plt.plot(tim, np.stack(HW)[:, 0, 0])
-plt.ylabel("Height Water")
-
-plt.savefig("KleinerExpA.png", pad_inches=0)
-
-plt.close()
diff --git a/tests/enthalpy/test-Kleiner/expB.py b/tests/enthalpy/test-Kleiner/expB.py
deleted file mode 100755
index ccbb7c5..0000000
--- a/tests/enthalpy/test-Kleiner/expB.py
+++ /dev/null
@@ -1,161 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file
-
-import scipy.io
-import tensorflow as tf
-import igm
-
-from igm.modules.process.enthalpy import *
-
-# this file is avilable at https://github.com/WangYuzhe/PoLIM-Polythermal-Land-Ice-Model
-# verif = scipy.io.loadmat("sol_analytic/enthB_analy_result.mat")
-
-ttf = 2000.0  # 300000
-dt = 1.0
-
-tim = np.arange(0, ttf, dt) + dt  # to put back to 300000
-
-parser = igm.params_core()
-igm.params_iceflow(parser)
-igm.params_enthalpy = params_enthalpy
-igm.params_enthalpy(parser)
-params = parser.parse_args(args=[])
-params.iflo_Nz = 100
-params.iflo_vert_spacing = 1
-
-PAGlen = 5.3e-24
-# in Pa^{-3} s^{-1}
-alpha = 4 * np.pi / 180  # inclination angle
-ws = -0.2  # # % [m a-1]
-Prho = params.iflo_ice_density
-Pg = params.iflo_gravity_cst
-
-params.enth_KtdivKc = 10 ** (-5)  # check this value if ok ?
-params.enth_till_wat_max = 200
-params.enth_drain_rate = 0
-
-thk = tf.Variable(200 * tf.ones((1, 1)))
-
-depth, dz = vertically_discretize_tf(thk, params.iflo_Nz, params.iflo_vert_spacing)
-
-strainheat = tf.Variable(
-    2.0
-    * PAGlen
-    * ((Prho * Pg * np.sin(alpha)) ** 4.0)
-    * (depth**4.0)
-    * params.enth_spy
-)
-strainheat = strainheat / params.enth_spy  # strainheatsec in Pa s^{-1} = W m^{-3}
-
-frictheat = tf.Variable(0.0 * tf.ones((1, 1)))
-geoheatflux = tf.Variable(0.0 * tf.ones((1, 1)))  # in W m-2
-tillwat = tf.Variable(0.0 * tf.ones((1, 1)))
-
-# Initial enthalpy field
-T = tf.Variable((-1.5 + 273.15) * tf.ones((params.iflo_Nz, 1, 1)))
-E = tf.Variable(params.enth_ci * (T - 223.15))
-omega = tf.Variable(tf.zeros_like(T))
-w = tf.Variable(ws * tf.ones_like(T)) / params.enth_spy
-
-surftemp = tf.Variable((-3.0 + 273.15) * tf.ones((1, 1)))
-
-for it, t in enumerate(tim):
-    Tpmp, Epmp = TpmpEpmp_from_depth_tf(
-        depth,
-        params.iflo_gravity_cst,
-        params.iflo_ice_density,
-        params.enth_claus_clape,
-        params.enth_melt_temp,
-        params.enth_ci,
-        params.enth_ref_temp,
-    )
-
-    surfenth = surf_enthalpy_from_temperature_tf(
-        surftemp, params.enth_melt_temp, params.enth_ci, params.enth_ref_temp
-    )
-
-    E, basalMeltRate = compute_enthalpy_basalmeltrate(
-        E,
-        Epmp,
-        dt * params.enth_spy,
-        dz,
-        w,
-        surfenth,
-        geoheatflux,
-        strainheat,
-        frictheat,
-        tillwat,
-        params.iflo_thr_ice_thk,
-        params.enth_ki,
-        params.iflo_ice_density,
-        params.enth_water_density,
-        params.enth_ci,
-        params.enth_ref_temp,
-        params.enth_Lh,
-        params.enth_spy,
-        params.enth_KtdivKc,
-        params.enth_drain_ice_column,
-    )
-
-    T, omega = temperature_from_enthalpy_tf(
-        E, Tpmp, Epmp, params.enth_ci, params.enth_ref_temp, params.enth_Lh
-    )
-
-    tillwat = tillwat + dt * (basalMeltRate - params.enth_drain_rate)
-    tillwat = tf.clip_by_value(tillwat, 0, params.enth_till_wat_max)
-
-    if it % 100 == 0:
-        print(
-            "time :",
-            t,
-            T[0, 0, 0].numpy(),
-            tillwat[0, 0].numpy(),
-            basalMeltRate[0, 0].numpy(),
-        )
-
-
-hz = np.arange(0, params.iflo_Nz) / (params.iflo_Nz - 1)
-
-fig = plt.figure(figsize=(6, 8))
-
-plt.subplot(131)
-plt.plot(E[:, 0, 0] / 1000, hz * 200)
-#plt.plot(
-#    verif["enthB_analy_E"] / 1000,
-#    verif["enthB_analy_z"] * 200,
-#    "--r",
-#    marker="o",
-#    markersize=5,
-#    markevery=20,
-#)
-plt.ylabel("Enthalpy (* 1000)")
-
-plt.subplot(132)
-plt.plot(T[:, 0, 0] - 273.15, hz * 200)
-#plt.plot(
-#    verif["enthB_analy_T"] - 273.15,
-#    verif["enthB_analy_z"] * 200,
-#    "--r",
-#    marker="o",
-#    markersize=5,
-#    markevery=20,
-#)
-plt.ylabel("Temperature (Deg Celcius)")
-
-plt.subplot(133)
-plt.plot(omega[:, 0, 0] * 100, hz * 200)
-#plt.plot(
-#    verif["enthB_analy_omega"] * 100,
-#    verif["enthB_analy_z"] * 200,
-#    "--r",
-#    marker="o",
-#    markersize=5,
-#    markevery=20,
-#)
-plt.ylabel("Water content (omega)  (/100)")
-
-plt.savefig("KleinerExpB.png", pad_inches=0)
-
-plt.close()
diff --git a/tests/enthalpy/test-Kleiner/test-TMDA.py b/tests/enthalpy/test-Kleiner/test-TMDA.py
deleted file mode 100755
index d9b58a1..0000000
--- a/tests/enthalpy/test-Kleiner/test-TMDA.py
+++ /dev/null
@@ -1,27 +0,0 @@
-#!/usr/bin/env python3
-
-import tensorflow as tf
-
-from igm.modules.process.enthalpy import *
-
-M = tf.ones((10))
-U = -2 * tf.ones((9))
-L = -2 * tf.ones((9))
-R = tf.Variable([4, 6, 7, 8, 9, 4, 6, 7, 8, 9.03])
-
-RES = solve_TDMA(L, M, U, R)
-
-A = tf.Variable(tf.zeros((10, 10)))
-
-for i in range(10):
-    A[i, i].assign(M[i])
-
-for i in range(9):
-    A[i, i + 1].assign(U[i])
-    A[i + 1, i].assign(L[i])
-
-tf.multiply(A, RES)
-
-RHS = tf.linalg.matvec(A, RES)
-
-print(RHS)
diff --git a/tests/instructed_oggm/run_instructed_oggm.py b/tests/instructed_oggm/run_instructed_oggm.py
deleted file mode 100644
index c0f635f..0000000
--- a/tests/instructed_oggm/run_instructed_oggm.py
+++ /dev/null
@@ -1,59 +0,0 @@
-
-"""
-    This script shows a simple example of using the "instruction" to run the IGM model
-    within the OGGM framework using the IGM_Model2D class. The IGM model is called at
-    each time step to compute the ice velocity and the ice thickness evolution. The
-    ice flow model is called through the python wrapper of the IGM package.
-    
-    Code written by: Julien Jehl, Fabien Maussion, and Guillaume Jouvet
-"""
-
-### general modules 
-import matplotlib.pyplot as plt
-import numpy as np
-import xarray as xr 
-
-### modules related to oggm & igm
-from oggm import cfg, utils, workflow, tasks
-from oggm.cfg import G, SEC_IN_YEAR, SEC_IN_DAY
-from oggm.core.massbalance import LinearMassBalance
-from igm.instructed_oggm import IGM_Model2D
-
-# Choice of a glacier
-rgi_ids = ["RGI60-11.01450"]  # Aletsch
-# rgi_ids = ['RGI60-11.00897']  # Hitereisferner
-
-cfg.initialize(logging_level="WARNING")
-cfg.PATHS["working_dir"] = utils.gettempdir(dirname="OGGM-Distrib", reset=True)
-base_url = "https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/exps/igm_v1/"
-
-gdirs = workflow.init_glacier_directories(
-    rgi_ids, prepro_base_url=base_url, from_prepro_level=2, prepro_border=30
-)
-gdir = gdirs[0]
-
-with xr.open_dataset(gdir.get_filepath("gridded_data")) as ds:
-    ds = ds.load()
-
-thick = ds.consensus_ice_thickness.where(~ds.consensus_ice_thickness.isnull(), 0)
-bed  = ds.topo - thick
-mask = ds.glacier_mask.data == 1
-topo = ds.topo
-
-# Define SMB model
-mb = LinearMassBalance(ela_h=2900.0)
- 
-# Define IGM model
-sdmodel = IGM_Model2D( bed.data, init_ice_thick=thick.data, dx=gdir.grid.dx,
-                       mb_model=mb, y0=0, mb_filter=mask, x=ds.x, y=ds.y)
-
-# Run the model
-ods = sdmodel.run_until_and_store(500, grid=gdir.grid, print_stdout="My run")
-
-# Plot the results
-for i in range(0,ods.ice_thickness.shape[0],10):
-    plt.imshow(ods.ice_thickness[i, :])
-    plt.colorbar()
-    plt.savefig("snapshot" + str(i) + ".png")
-    plt.close()
-    print(i)
diff --git a/tests/params.json b/tests/params.json
new file mode 100755
index 0000000..9362b94
--- /dev/null
+++ b/tests/params.json
@@ -0,0 +1,32 @@
+{
+  "modules_preproc": [
+    "load_ncdf",
+    "track_usurf_obs"
+  ],
+  "modules_process": [
+    "clim_aletsch",
+    "smb_accmelt",
+    "iceflow",
+    "time",
+    "thk",
+    "rockflow",
+    "vert_flow",
+    "particles"
+  ],
+  "modules_postproc": [
+    "write_ncdf",
+    "plot2d",
+    "print_info",
+    "print_comp"
+  ],
+  "time_start": 1880.0,
+  "time_end": 2020.0,
+  "time_save": 5.0,
+  "plt2d_live": true,
+  "iflo_init_slidingco": 0.0595,
+  "weight_accumulation": 1.0,
+  "weight_ablation": 1.25,
+  "part_frequency_seeding": 500,
+  "part_tracking_method": "simple",
+  "part_density_seeding": 1
+}
\ No newline at end of file
diff --git a/tests/params_comments.json b/tests/params_comments.json
new file mode 100755
index 0000000..a2554ca
--- /dev/null
+++ b/tests/params_comments.json
@@ -0,0 +1,41 @@
+#
+
+{
+  "modules_preproc": [
+    "load_ncdf",
+    "track_usurf_obs"
+  ],
+  "modules_process": [
+    "clim_aletsch",
+    "smb_accmelt",
+    "iceflow",
+    "time",
+    "thk",
+    ##
+    "rockflow",
+    "vert_flow",
+    "particles"
+  ],
+  "modules_postproc": [
+    "write_ncdf",
+    "plot2d",
+    "print_info",
+    "print_comp"
+  ],
+  "time_start": 1880.0,
+  "time_end": 2020.0,
+  // this is a comment
+  "time_save": 5.0,
+  "plt2d_live": true,
+  "iflo_init_slidingco": 0.0595,
+  "weight_accumulation": 1.0,
+  "weight_ablation": 1.25,
+  "part_frequency_seeding": 500,
+  "part_tracking_method": "simple",
+  "part_density_seeding": 1
+}
+
+
+
+//
+//
diff --git a/tests/real_glacier/check.py b/tests/real_glacier/check.py
deleted file mode 100755
index 26623a0..0000000
--- a/tests/real_glacier/check.py
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file 
-
-import tensorflow as tf
-import numpy as np
-
-def params_check(parser):
-    pass 
-
-def initialize_check(params,state):
-    pass    
-
-def update_check(params,state):
-    pass
-
-def finalize_check(params, state):
-    
-    modules = [A for A in state.__dict__.keys() if 'tcomp_' in A]
-
-    state.tcomp_all = [ np.sum([np.sum(getattr(state,m)) for m in modules]) ]
-   
-    vol = np.sum(state.thk) * (state.dx**2) / 10**9
-    if (vol<2.8)&(vol>.75):
-       print(":-) test Passed in "+ str(int(state.tcomp_all[0]))+" sec")
-    else:
-       print(":-( test NOT PASSED")
-       
-    
\ No newline at end of file
diff --git a/tests/real_glacier/params.json b/tests/real_glacier/params.json
deleted file mode 100755
index abcaeab..0000000
--- a/tests/real_glacier/params.json
+++ /dev/null
@@ -1,22 +0,0 @@
-{
-  "modules_preproc": ["oggm_shop"],
-  "modules_process": ["smb_simple",
-                      "iceflow",
-                      "time",
-                      "thk",
-                      "vert_flow",
-                      "particles"],
-  "modules_postproc": ["check"],
-  "oggm_RGI_ID": "RGI60-11.01238",
-  "time_start": 2000.0,
-  "time_end": 2100.0,
-  "time_save": 5.0,
-  "smb_simple_array": [ 
-                        ["time", "gradabl", "gradacc", "ela", "accmax"],
-                        [ 2000,      0.009,     0.005,  2700,      2.0],
-                        [ 2100,      0.009,     0.005,  2600,      2.0]
-                      ],
-  "iflo_new_friction_param": true,
-  "iflo_dim_arrhenius": 2,
-  "iflo_init_slidingco": 0.1
-}
diff --git a/tests/run_tests.sh b/tests/run_tests.sh
new file mode 100644
index 0000000..3aaf1e9
--- /dev/null
+++ b/tests/run_tests.sh
@@ -0,0 +1,5 @@
+# will run all tests in the tests folder and ignore depreciated warnings (numpy).
+# This is across all tests - ideally, use a decorator to ignore warnings on specific functions
+# -p is to not create pycached from pytest (important for pypi as deployment will mess up the package)
+# ! Supposed to be newer -B but confirm!
+python -m pytest -W ignore::DeprecationWarning -p no:cacheprovider
\ No newline at end of file
diff --git a/tests/synthetic/check.py b/tests/synthetic/check.py
deleted file mode 100755
index 1b2325b..0000000
--- a/tests/synthetic/check.py
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file 
-
-import tensorflow as tf
-import numpy as np
-
-def params_check(parser):
-    pass 
-
-def initialize_check(params,state):
-    pass    
-
-def update_check(params,state):
-    pass
-
-def finalize_check(params, state):
-    
-    modules = [A for A in state.__dict__.keys() if 'tcomp_' in A]
-
-    state.tcomp_all = [ np.sum([np.sum(getattr(state,m)) for m in modules]) ]
-   
-    vol = np.sum(state.thk) * (state.dx**2) / 10**9
-    if (vol<11.5)&(vol>11.0):
-       print(":-) test Passed in "+ str(int(state.tcomp_all[0]))+" sec (ref 11)")
-    else:
-       print(":-( test NOT PASSED")
-       
-    
diff --git a/tests/synthetic/make_synthetic.py b/tests/synthetic/make_synthetic.py
deleted file mode 100755
index d128802..0000000
--- a/tests/synthetic/make_synthetic.py
+++ /dev/null
@@ -1,42 +0,0 @@
-#!/usr/bin/env python3
-
-# Copyright (C) 2021-2023 Guillaume Jouvet <guillaume.jouvet@unil.ch>
-# Published under the GNU GPL (Version 3), check at the LICENSE file 
-
-import os
-import numpy as np
-import matplotlib.pyplot as plt
-import tensorflow as tf
-
-from igm.modules.utils import complete_data
-
-def params_make_synthetic(parser):
-    pass
-
-
-def initialize_make_synthetic(params, state):
-
-    x = np.arange(0, 100) * 100  # make x-axis, lenght 10 km, resolution 100 m
-    y = np.arange(0, 200) * 100  # make y-axis, lenght 20 km, resolution 100 m
-
-    X, Y = np.meshgrid(x, y)
-
-    topg = 1000 + 0.15 * Y + ((X - 5000) ** 2) / 50000  # define the bedrock topography
-    thk = np.zeros_like(topg)
-
-    state.x = tf.constant(x.astype("float32"))
-    state.y = tf.constant(y.astype("float32"))
-
-    state.topg = tf.Variable(topg.astype("float32"))
-    state.thk = tf.Variable(thk.astype("float32"))
-
-    complete_data(state)
-
-
-def update_make_synthetic(params, state):
-    pass
-
-
-def finalize_make_synthetic(params, state):
-    pass
-    
diff --git a/tests/synthetic/params.json b/tests/synthetic/params.json
deleted file mode 100755
index de0322b..0000000
--- a/tests/synthetic/params.json
+++ /dev/null
@@ -1,17 +0,0 @@
-{
-  "modules_preproc": ["make_synthetic"],
-  "modules_process": ["smb_simple",
-                      "iceflow",
-                      "time",
-                      "thk"
-                    ],
-  "modules_postproc": ["check"],
-  "smb_simple_array": [
-                        ["time", "gradabl", "gradacc", "ela", "accmax"],
-                        [ 2000,      0.009,     0.005,  2900,      2.0],
-                        [ 2100,      0.009,     0.005,  3300,      2.0]
-                      ],
-  "time_start": 2000.0,
-  "time_end": 2100.0,
-  "time_save": 10.0
-}
diff --git a/tests/test_all.sh b/tests/test_all.sh
deleted file mode 100644
index 5ea4cf5..0000000
--- a/tests/test_all.sh
+++ /dev/null
@@ -1,15 +0,0 @@
-# This permits to perform a minimaum aamount of testing
-
-cd real_glacier 
-igm_run
-cd ..
-
-cd synthetic 
-igm_run
-cd ..
-
-cd instructed_oggm 
-python run_instructed_oggm.py
-cd ..
-
-
diff --git a/tests/test_loading_modules.py b/tests/test_loading_modules.py
new file mode 100644
index 0000000..40f12af
--- /dev/null
+++ b/tests/test_loading_modules.py
@@ -0,0 +1,14 @@
+import igm
+import pytest
+from igm.common import validate_module, load_modules_from_directory
+from json import JSONDecodeError
+
+
+# def test_load_valid_module():
+#     pass
+
+# def test_load_custom_module_missing_function():
+#     pass
+
+# def test_load_custom_module_package_missing_function():
+#     pass
diff --git a/tests/test_params.py b/tests/test_params.py
new file mode 100644
index 0000000..678616e
--- /dev/null
+++ b/tests/test_params.py
@@ -0,0 +1,183 @@
+import igm
+import pytest
+
+# from igm.common import validate_module, load_modules_from_directory
+from json import JSONDecodeError
+
+
+def test_load_json_params_with_comments():
+    """Tests that adding comments to the json will not fail with a decode error."""
+
+    igm.get_modules_list("params_comments.json")
+
+
+def test_load_igm_modules():
+    """Tests that the core igm modules LIST loaded from the params.json file are loaded."""
+    modules_dict = igm.get_modules_list("params.json")
+
+    preproc_modules = modules_dict["modules_preproc"]
+    process_modules = modules_dict["modules_process"]
+    postproc_modules = modules_dict["modules_postproc"]
+
+    igm_core_preproc_modules = ["load_ncdf"]
+    igm_core_process_modules = [
+        "iceflow",
+        "time",
+        "thk",
+        "rockflow",
+        "vert_flow",
+        "particles",
+    ]
+    igm_core_postproc_modules = ["write_ncdf", "plot2d", "print_info", "print_comp"]
+
+    assert set(igm_core_preproc_modules).issubset(preproc_modules)
+    assert set(igm_core_process_modules).issubset(process_modules)
+    assert set(igm_core_postproc_modules).issubset(postproc_modules)
+
+
+def test_load_modules():
+    """Tests that the core igm modules LIST AND the custom modules from the params.json file are loaded."""
+    modules_dict = igm.get_modules_list("params.json")
+    assert modules_dict == {
+        "modules_preproc": ["load_ncdf", "track_usurf_obs"],
+        "modules_process": [
+            "clim_aletsch",
+            "smb_accmelt",
+            "iceflow",
+            "time",
+            "thk",
+            "rockflow",
+            "vert_flow",
+            "particles",
+        ],
+        "modules_postproc": ["write_ncdf", "plot2d", "print_info", "print_comp"],
+    }
+
+
+def test_params_core():
+    parser = igm.params_core()
+    params, __ = parser.parse_known_args()
+
+    assert vars(params) == {
+        "param_file": "params.json",
+        "modules_preproc": ["oggm_shop"],
+        "modules_process": ["iceflow", "time", "thk"],
+        "modules_postproc": ["write_ncdf", "plot2d", "print_info"],
+        "logging": False,
+        "logging_level": 30,
+        "logging_file": "",
+        "print_params": True,
+        "gpu_info": False,
+        "gpu_id": 0,
+        "saved_params_filename": "params_saved",
+    }
+
+# def test_params_core_overwrite():
+#     state = igm.State()  # class acting as a dictionary
+#     parser = igm.params_core()
+#     # args = ['--working_dir', 'test', '--gpu_info']  # replace with your arguments
+#     parser = get_parser(parser)
+#     __, params, __ = igm.setup_igm(state=state, parser=parser)
+
+#     assert vars(params) == {
+#         "working_dir": "test",
+#         "param_file": "params.json",
+#         "modules_preproc": ["oggm_shop"],
+#         "modules_process": ["iceflow", "time", "thk"],
+#         "modules_postproc": ["write_ncdf", "plot2d", "print_info"],
+#         "logging": False,
+#         "logging_level": 30,
+#         "logging_file": "",
+#         "print_params": True,
+#         "gpu_info": False,
+#         "gpu_id": 0,
+#         "saved_params_filename": "params_saved",
+#     }
+
+# def get_parser(parser):
+#     params, _ = parser.parse_known_args()
+#     params_dict = vars(params)
+#     # print(params_dict)
+#     params_dict.update({"working_dir": "test"})
+#     # print(params_dict)
+
+#     parser.set_defaults(**params_dict)
+    
+#     # print(parser)
+#     # params = parser.parse_args()
+#     return parser
+
+# def test_params_core_overwrite():
+#     state = igm.State()  # class acting as a dictionary
+#     parser = igm.params_core()
+
+#     params, _ = parser.parse_known_args()
+#     params_dict = vars(params)
+#     print(params_dict)
+#     params_dict.update({"working_dir": "test"})
+#     print(params_dict)
+
+#     parser.set_defaults(**params_dict)
+    
+#     print(parser)
+#     params = parser.parse_args()
+#     # args = ['--arg1', 'value1', '--arg2', 'value2']  # replace with your arguments
+#     # parsed_args = parser(args)
+#     __, params, __ = igm.setup_igm(state=state, parser=parser)
+
+#     assert vars(params) == {
+#         "working_dir": "test",
+#         "param_file": "params.json",
+#         "modules_preproc": ["oggm_shop"],
+#         "modules_process": ["iceflow", "time", "thk"],
+#         "modules_postproc": ["write_ncdf", "plot2d", "print_info"],
+#         "logging": False,
+#         "logging_level": 30,
+#         "logging_file": "",
+#         "print_params": True,
+#         "gpu_info": False,
+#         "gpu_id": 0,
+#         "saved_params_filename": "params_saved",
+#     }
+
+
+# def test_params_json():
+#     parser = igm.params_core()
+#     params, _ = parser.parse_known_args()
+
+#     modules_dict = igm.get_modules_list(params.param_file)
+#     imported_modules = igm.load_modules(modules_dict)
+#     imported_modules = igm.load_dependent_modules(imported_modules)
+
+#     for module in imported_modules:
+#         module.params(parser)
+
+#     core_and_module_params = parser.parse_args()
+#     params = igm.load_user_defined_params(
+#         param_file=core_and_module_params.param_file,
+#         params_dict=vars(core_and_module_params),
+#     )
+
+#     parser.set_defaults(**params)
+#     params = parser.parse_args()
+
+
+# def test_params_cli():
+#     parser = igm.params_core()
+#     params, _ = parser.parse_known_args()
+
+#     modules_dict = igm.get_modules_list(params.param_file)
+#     imported_modules = igm.load_modules(modules_dict)
+#     imported_modules = igm.load_dependent_modules(imported_modules)
+
+#     for module in imported_modules:
+#         module.params(parser)
+
+#     core_and_module_params = parser.parse_args()
+#     params = igm.load_user_defined_params(
+#         param_file=core_and_module_params.param_file,
+#         params_dict=vars(core_and_module_params),
+#     )
+
+#     parser.set_defaults(**params)
+#     params = parser.parse_args()