diff --git a/README.md b/README.md
index d696f1a..0ef05a5 100644
--- a/README.md
+++ b/README.md
@@ -54,3 +54,4 @@ module avail access-issm
 For users of ACCESS-ISSM model configurations released by ACCESS-NRI, the exact location of the model executables is not required. Model configurations will be updated with new model components when necessary.
 
 For information on contributing your own fixes to the ACCESS-ISSM `spack.yaml`, see the [Hive Docs article](https://docs.access-hive.org.au/models/run-a-model/create-a-prerelease/) on creating a prerelease.
+
diff --git a/config/versions.json b/config/versions.json
index 2e0bd08..97cd365 100644
--- a/config/versions.json
+++ b/config/versions.json
@@ -1,5 +1,5 @@
 {
     "$schema": "https://raw.githubusercontent.com/ACCESS-NRI/schema/main/au.org.access-nri/model/deployment/config/versions/3-0-0.json",
     "spack": "0.22",
-    "spack-packages": "2025.04.001"
+    "spack-packages": "justin/debug"
 }
diff --git a/spack.yaml b/spack.yaml
index 3075665..c3e5443 100644
--- a/spack.yaml
+++ b/spack.yaml
@@ -1,11 +1,14 @@
 spack:
   specs:
-    - access-issm@git.2025.03.0
+    - access-issm@git.2025.06.0
+
   packages:
     issm:
       require:
-        - '@git.2025.04.11'
+        - '@git.2025.06.19'
         - +wrappers
+        - +debug
+        
     # Mark Python 3.9.2 as external so Spack reuses the system module
     python:
       externals:
@@ -27,3 +30,14 @@ spack:
   # Make the solver unify dependencies
   concretizer:
     unify: true
+
+  # Enable generating module files (TCL by default).
+  modules:
+    default:
+      tcl:
+        include:
+          - access-issm
+          - issm
+          - access-triangle
+        projections:
+          access-issm: '{name}/2025.06.0-{hash:7}'
diff --git a/test/Mismip.py b/test/Mismip.py
new file mode 100644
index 0000000..75e0a42
--- /dev/null
+++ b/test/Mismip.py
@@ -0,0 +1,76 @@
+import numpy as np
+from SetIceShelfBC import SetIceShelfBC
+from mismipbasalforcings import mismipbasalforcings
+
+# Creating thickness
+print('      creating thickness')
+bx = -150 - 728.8*(md.mesh.x/300000)**2 + 343.91*(md.mesh.x/300000)**4 - 50.57*(md.mesh.x/300000)**6
+by = (500./(1 + np.exp(-2./4000.*(md.mesh.y - 80000./2 - 24000))) +
+      500./(1 + np.exp( 2./4000.*(md.mesh.y - 80000./2 + 24000))))
+by0 = (500./(1 + np.exp(-2./4000.*(0     - 80000./2 - 24000))) +
+       500./(1 + np.exp( 2./4000.*(0     - 80000./2 + 24000))))
+md.geometry.bed     = np.maximum(bx + by, -720)
+md.geometry.surface = np.maximum(bx + by0 + 100, 10)
+md.geometry.base    = np.maximum(md.geometry.bed, -90)
+md.geometry.thickness = md.geometry.surface - md.geometry.base
+
+# Creating drag
+print('      creating drag')
+md.friction.coefficient = np.sqrt(3.160e6) * np.ones(md.mesh.numberofvertices)
+md.friction.p           =            3 * np.ones(md.mesh.numberofelements)
+md.friction.q           =            0 * np.ones(md.mesh.numberofelements)
+
+# Creating flow law parameter
+print('      creating flow law parameter')
+md.materials.rheology_B   = 1/((6.338e-25)**(1/3)) * np.ones(md.mesh.numberofvertices)
+md.materials.rheology_n   =            3 * np.ones(md.mesh.numberofelements)
+md.materials.rheology_law = 'None'
+
+# Boundary conditions for diagnostic model
+print('      boundary conditions for diagnostic model')
+md = SetIceShelfBC(md, './Front.exp')
+md.mask.ice_levelset[:]   = -1
+md.mask.ocean_levelset[:] = -1
+pos = np.where((md.mesh.x < 640000.1) & (md.mesh.x > 639999.9))[0]
+md.mask.ice_levelset[pos] = 0
+md.stressbalance.spcvx[:]  = np.nan
+md.stressbalance.spcvy[:]  = np.nan
+pos = np.where(((md.mesh.y <  80000.1) & (md.mesh.y >  79999.9)) |
+               ((md.mesh.y <      0.1) & (md.mesh.y >     -0.1)))[0]
+md.stressbalance.spcvy[pos] = 0
+pos2 = np.where((md.mesh.x <      0.1) & (md.mesh.x >     -0.1))[0]
+md.stressbalance.spcvx[pos2] = 0
+md.stressbalance.spcvy[pos2] = 0
+
+# Forcing conditions
+print('      forcing conditions')
+# --- basal melt and buttressing conditions -----------------------
+md.basalforcings = mismipbasalforcings()      #  ← remove md
+
+
+md.basalforcings.meltrate_factor        =   0
+md.basalforcings.threshold_thickness    =  75
+md.basalforcings.upperdepth_melt        = -100
+md.smb.mass_balance                    = 0.3 * np.ones(md.mesh.numberofvertices)
+md.basalforcings.geothermalflux        = 0.5 * np.ones(md.mesh.numberofvertices)
+md.basalforcings.groundedice_melting_rate = 0. * np.ones(md.mesh.numberofvertices)
+
+# Other model parameters
+md.thermal.spctemperature        = np.nan * np.ones(md.mesh.numberofvertices)
+md.groundingline.migration       = 'SubelementMigration'
+md.materials.rho_ice             = 918
+md.materials.rho_water           = 1028
+md.constants.g                   = 9.8
+md.constants.yts                 = 31556926
+md.transient.isthermal           = 0
+md.transient.isgroundingline     = 1
+md.stressbalance.isnewton        = 0
+
+# Initialization
+md.initialization.vx          = np.ones(md.mesh.numberofvertices)
+md.initialization.vy          = np.ones(md.mesh.numberofvertices)
+md.initialization.vz          = np.ones(md.mesh.numberofvertices)
+md.initialization.vel         = np.sqrt(2) * np.ones(md.mesh.numberofvertices)
+md.initialization.pressure    = md.constants.g * md.materials.rho_ice * md.geometry.thickness
+md.initialization.temperature = 273 * np.ones(md.mesh.numberofvertices)
+
diff --git a/test/mismip_driver.py b/test/mismip_driver.py
new file mode 100644
index 0000000..e370a77
--- /dev/null
+++ b/test/mismip_driver.py
@@ -0,0 +1,780 @@
+#!/usr/bin/env python3
+import os, sys
+
+# Ensure ISSM paths are correctly set
+issm_dir = os.getenv('ISSM_DIR')
+if not issm_dir:
+    print("Error: ISSM_DIR environment variable is not set.")
+    sys.exit(1)
+
+import numpy as np
+from triangle import triangle
+from model import *
+from netCDF4 import Dataset
+from InterpFromGridToMesh import InterpFromGridToMesh
+from bamg import bamg
+from xy2ll import xy2ll
+from plotmodel import plotmodel
+from export_netCDF import export_netCDF
+from loadmodel import loadmodel
+from setmask import setmask
+from parameterize import parameterize
+from setflowequation import setflowequation
+from socket import gethostname
+from solve import solve
+from ll2xy import ll2xy
+from BamgTriangulate import BamgTriangulate
+from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d
+from scipy.interpolate import griddata
+from matplotlib import pyplot as plt
+from gadi_spack import gadi
+from organizer import organizer
+from toolkits import toolkits
+from bcgslbjacobioptions import bcgslbjacobioptions
+from SetMOLHOBC import SetMOLHOBC
+
+# 1) Define the steps and model number
+#steps = [11]
+#modelnum = 1
+
+# 2) Choose model name from model number
+if modelnum == 1:
+    modelname = '1km_viscous'
+elif modelnum == 2:
+    modelname = '2km_viscous'
+elif modelnum == 3:
+    modelname = '1km_coulomb'
+elif modelnum == 4:
+    modelname = '2km_coulomb'
+elif modelnum == 5:
+    modelname = '500m_viscous'
+elif modelnum == 6:
+    modelname = '500m_coulomb'
+elif modelnum == 7:
+    modelname = '200m_viscous'
+elif modelnum == 8:
+    modelname = '200m_coulomb'
+else:
+    raise ValueError('Invalid modelnum!')
+
+# Hard-coded parameters
+clustername = 'gadi'
+# loadonly = 0
+# interactive = 0
+printflag = True
+
+
+
+if clustername.lower() == 'gadi':
+
+    cluster = gadi(
+        'name',         'gadi.nci.org.au',
+        'login',        'jh7060',
+        'srcpath',      '/g/data/au88/jh7060/ISSM',
+        'project',      'au88',
+        'numnodes',     1,
+        'cpuspernode',  32,
+        'time',         48*60,  # minutes
+        'codepath',     '/g/data/au88/jh7060/spack/0.22/release/linux-rocky8-x86_64_v4/gcc-13.2.0/issm-4.24-v52nx3pfx7lpqwfldqlpspflk34wz756/bin',
+        'executionpath','/scratch/au88/jh7060/issm_runs'
+    )
+    loadonly = 0
+    lock     = 1
+else:
+    # --- generic single-node/local fallback ---
+    cluster        = generic()            # <-- no args here
+    cluster.name   = gethostname()        # label for logs
+    cluster.np     = 32              # number of MPI processes
+    cluster.mem    = 64              # GB – change to suit your machine
+    cluster.time   = 60                  # wall-clock minutes
+    cluster.queue  = 'local'              # or the default queue on your scheduler
+    lock     = 1
+    loadonly = 1
+
+# Create an organizer
+# - The “organizer” is the same logic as the MATLAB version.
+# - Note that in Python, you typically use Organizer(...) from issm.
+org = organizer(
+     'repository',   './Models_'     + modelname,
+     'prefix',       'mismip_'       + modelname + '_',
+     'steps',        steps,
+     'trunkprefix',  '34;47;2'
+)
+# ------------------------------------------------------------------------------
+#  STEP1: Mesh_generation
+# ------------------------------------------------------------------------------
+if org.perform('Mesh_generation'):
+    model = model()  # Start a blank model
+    if modelnum == 1 or modelnum == 3:
+        md = bamg(model,
+          'domain',       './Domain.exp',
+          'hmax',         1000,
+          'splitcorners', 1)
+
+        print(md.mesh.x.max() -  md.mesh.x.min())
+    elif modelnum == 2 or modelnum == 4:
+        md = bamg(model, domain='./Domain.exp', hmax=2000, splitcorners=1)
+    elif modelnum == 5 or modelnum == 6:
+        md = bamg(model, domain='./Domain.exp', hmax=500, splitcorners=1)
+    elif modelnum == 7 or modelnum == 8:
+        md = bamg(model, domain='./Domain.exp', hmax=200, splitcorners=1)
+    else:
+        raise RuntimeError("Model not supported yet")
+
+    md.miscellaneous.name = 'MISMIP_' + modelname
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP2: Parameterization
+# ------------------------------------------------------------------------------
+if org.perform('Parameterization'):
+    md = org.loadmodel('Mesh_generation')
+    md = setmask(md, '', '')
+    md = parameterize(md, './Mismip.py')
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP3: Transient_Steadystate
+# ------------------------------------------------------------------------------
+if org.perform('Transient_Steadystate'):
+    md = org.loadmodel('Parameterization')
+    md = setflowequation(md, 'SSA', 'all')
+
+    # Coulomb friction for certain models
+    if (modelnum == 3) or (modelnum == 4) or (modelnum == 6):
+        md.friction = frictioncoulomb()
+        md.friction.coefficient = math.sqrt(3.160e6) * np.ones(md.mesh.numberofvertices)
+        md.friction.coefficientcoulomb = math.sqrt(0.5) * np.ones(md.mesh.numberofvertices)
+        md.friction.p = 3 * np.ones(md.mesh.numberofelements)
+        md.friction.q = np.zeros(md.mesh.numberofelements)
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 200000
+    md.settings.output_frequency = 2000
+    md.settings.checkpoint_frequency = 2000
+    md.stressbalance.maxiter = 30
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('solution', True, 'module', True, 'convergence', True)
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tss1'
+    md.settings.solver_residue_threshold = float('nan')
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype,'loadonly', loadonly, 'lock', lock, 'runtimename', False)
+    export_netCDF(md, "./mismip_1km_viscous_Transient_steadystate.nc")
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP4: Transient_Steadystate_remesh (example of specialized re-meshing)
+# ------------------------------------------------------------------------------
+if org.perform('Transient_Steadystate_remesh'):
+    md = org.loadmodel('Parameterization')
+    # Example: load from another directory
+    # The code below tries to load an existing model from the 2km_viscous folder.
+    # Adapt to your local file paths if needed.
+    md2 = loadmodel('Models_2km_viscous/mismip_2km_viscous_Transient_steadystate2.nc')
+    # Overwrite the final solution’s mesh
+    md2.results.TransientSolution[-1].MeshX = md.mesh.x
+    md2.results.TransientSolution[-1].MeshY = md.mesh.y
+    md2.results.TransientSolution[-1].MeshElements = md.mesh.elements
+
+    md = remesh(md2, './Mismip.py')
+    md = setflowequation(md, 'SSA', 'all')
+    if (modelnum == 3) or (modelnum == 4) or (modelnum == 6):
+        md.friction = frictioncoulomb()
+        md.friction.coefficient = math.sqrt(3.160e6) * np.ones(md.mesh.numberofvertices)
+        md.friction.coefficientcoulomb = math.sqrt(0.5) * np.ones(md.mesh.numberofvertices)
+        md.friction.p = 3 * np.ones(md.mesh.numberofelements)
+        md.friction.q = np.zeros(md.mesh.numberofelements)
+
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 200000
+    md.settings.output_frequency = 5000
+    md.settings.checkpoint_frequency = 5000
+    md.stressbalance.maxiter = 10
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('convergence', False, 'solution', False)
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tssr'
+    md.settings.solver_residue_threshold = float('nan')
+    md.settings.waitonlock = 0
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+ 
+
+# ------------------------------------------------------------------------------
+#  STEP5: Transient_steadystate2
+# ------------------------------------------------------------------------------
+if org.perform('Transient_steadystate2'):
+    md = org.loadmodel('Transient_Steadystate')
+
+    md = setflowequation(md, 'SSA', 'all')
+    # Re-initialize from previous solution
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 10000
+    md.settings.output_frequency = 1000
+    md.settings.checkpoint_frequency = 5000
+    md.stressbalance.maxiter = 10
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('solution', True, 'module', True, 'convergence', True)
+    md.cluster = cluster
+    md.settings.waitonlock = 1
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tss2'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype,'loadonly', loadonly, 'lock', lock)#, 'runtimename', False)
+    #export_netCDF(md, "./mismip_1km_viscous_Transient_steadystate_2.nc")
+
+# ------------------------------------------------------------------------------
+#  STEP6: Transient_steadystate3
+# ------------------------------------------------------------------------------
+if org.perform('Transient_steadystate3'):
+    md = org.loadmodel('Transient_steadystate2')
+    md = setflowequation(md, 'SSA', 'all')
+    # Re-initialize from previous solution
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 200000
+    md.settings.output_frequency = 6000
+    md.stressbalance.maxiter = 10
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('solution', True, 'module', True, 'convergence', True)
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tss3'
+    md.settings.solver_residue_threshold = float('nan')
+    md.settings.waitonlock = 0
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype,'loadonly', loadonly, 'lock', lock)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP7: Transient_steadystate4
+# ------------------------------------------------------------------------------
+if org.perform('Transient_steadystate4'):
+    md = org.loadmodel('Transient_steadystate3')
+    md = setflowequation(md, 'SSA', 'all')
+    # Re-initialize from previous solution
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 200000
+    md.settings.output_frequency = 6000
+    md.stressbalance.maxiter = 10
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('convergence', False, 'solution', True)
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tss4'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP8: Transient_steadystate5
+# ------------------------------------------------------------------------------
+if org.perform('Transient_steadystate5'):
+    md = org.loadmodel('Transient_steadystate4')
+    md = setflowequation(md, 'SSA', 'all')
+    # Re-initialize from previous solution
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md.timestepping.time_step = 1
+    md.timestepping.final_time = 200000
+    md.settings.output_frequency = 6000
+    md.stressbalance.maxiter = 10
+    md.stressbalance.abstol = float('nan')
+    md.stressbalance.restol = 1
+    md.verbose = verbose('convergence', False, 'solution', True)
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_Tss5'
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP9: Transient_extrude
+# ------------------------------------------------------------------------------
+if org.perform('Transient_extrude'):
+    md = org.loadmodel('Transient_steadystate3')
+    # Re-initialize from previous solution
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md = md.extrude( 10, 1.1)
+    md = setflowequation(md, 'HO', 'all')
+    md.transient.isthermal = 0
+    md.transient.issmb = 0
+    md.initialization.temperature[:] = 273.0
+
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP10: GlenSSA
+# ------------------------------------------------------------------------------
+if org.perform('GlenSSA'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate') #200k years relaxation
+    else:
+        md = org.loadmodel('Transient_extrude') #600k years relaxation
+
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md = md.collapse()  # collapse 3D model to 2D if extruded
+    md = setflowequation(md, 'SSA', 'all')
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GSSA'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP11: GlenMOLHO
+# ------------------------------------------------------------------------------
+if org.perform('GlenMOLHO'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'VxShear','VyShear','VxBase','VyBase','VxSurface','VySurface',
+        'VxAverage','VyAverage','StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md = md.collapse()
+    md = setflowequation(md, 'MOLHO', 'all')
+
+    # Example solver settings
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    # If you need special solver settings for small meshes:
+   # md.toolkits = toolkits()
+    md.toolkits = toolkits.addoptions(md.toolkits,'StressbalanceAnalysis',bcgslbjacobioptions())
+
+    md = SetMOLHOBC(md)  # Shear boundary conditions for MOLHO?
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GMOLHO'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP12: GlenHO
+# ------------------------------------------------------------------------------
+if org.perform('GlenHO'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md = setflowequation(md, 'HO', 'all')
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'StrainRatexx','StrainRatexy','StrainRateyy','StrainRatexz','StrainRateyz','StrainRatezz',
+        'StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    #md.toolkits = md.addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GHO'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP13: GlenFS
+# ------------------------------------------------------------------------------
+if org.perform('GlenFS'):
+    # If using a 2D mesh (e.g. 500m case), use the steady state model
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        # or extruded model
+        md = org.loadmodel('Transient_steadystate3')
+
+        # or load from GlenHO to reuse geometry
+        # md = org.loadmodel('GlenHO')
+    md.initialization.vx = md.results.TransientSolution[-1].Vx
+    md.initialization.vy = md.results.TransientSolution[-1].Vy
+    md.initialization.vel = md.results.TransientSolution[-1].Vel
+    md.geometry.thickness = md.results.TransientSolution[-1].Thickness
+    md.geometry.base = md.results.TransientSolution[-1].Base
+    md.geometry.surface = md.results.TransientSolution[-1].Surface
+    md.mask.ocean_levelset = md.results.TransientSolution[-1].MaskOceanLevelset
+
+    md = md.extrude( 10, 1.1)
+    md = setflowequation(md, 'HO', 'all')
+    md.transient.isthermal = 0
+    md.transient.issmb = 0
+    #md.transient.ismasstransport = 0
+    md.initialization.temperature[:] = 273.0
+    
+    md = setflowequation(md, 'FS', 'all')
+    md.stressbalance.shelf_dampening = 1
+    md.masstransport.isfreesurface = 1
+    md.transient.isgroundingline = 1  # or 1, depending on your test
+    md.groundingline.migration='Contact' # or 'SoftMigration'
+    md.timestepping.time_step = 0.00001
+    md.timestepping.final_time = 0.0001
+    md.settings.output_frequency = 1
+   # md = md.SetInput(md, 'Bed', md.geometry.bed)
+    md.toolkits = toolkits.addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
+    md.flowequation.fe_FS = 'TaylorHood'
+    md.stressbalance.maxiter = 20
+    md.stressbalance.restol = 0.5
+    md.stressbalance.reltol = 0.001
+    # If you know surface and thickness:
+    #md.geometry.bed = md.geometry.surface - md.geometry.thickness
+
+    #md.geometry.surface = InterpFromGridToMesh(x, y, usrf, md.mesh.x, md.mesh.y, 0);
+    #md.geometry.bed     = InterpFromGridToMesh(x, y, topg, md.mesh.x, md.mesh.y, 0);
+
+
+    # Or if you already interpolated base (from NetCDF):
+#; % Sometimes 'base' is used internally
+
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GFS'
+    print(md.mesh)
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype,'loadonly', loadonly, 'lock', lock, 'runtimename', False)
+    #org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP14: GlenESSA (enhanced SSA with a factor E=5)
+# ------------------------------------------------------------------------------
+if org.perform('GlenESSA'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md = collapse(md)
+    md = setflowequation(md, 'SSA', 'all')
+    md.materials = matenhancedice(md.materials)
+    md.materials.rheology_E = 5.0 * np.ones(md.mesh.numberofvertices)
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_ESSA'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP15: GlenEMOLHO (enhanced MOLHO)
+# ------------------------------------------------------------------------------
+if org.perform('GlenEMOLHO'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'VxShear','VyShear','VxBase','VyBase','VxSurface','VySurface',
+        'VxAverage','VyAverage','StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md = collapse(md)
+    md = setflowequation(md, 'MOLHO', 'all')
+    md.materials = matenhancedice(md.materials)
+    md.materials.rheology_E[:] = 5.0
+
+    md = SetMOLHOBC(md)
+
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GEMOLHO'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+# ------------------------------------------------------------------------------
+#  STEP16: GlenEHO (enhanced HO)
+# ------------------------------------------------------------------------------
+if org.perform('GlenEHO'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md = setflowequation(md, 'HO', 'all')
+    md.materials = matenhancedice(md.materials)
+    md.materials.rheology_E = 5.0 * np.ones(md.mesh.numberofvertices)
+    md.transient.requested_outputs = [
+        'default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+        'StrainRatexx','StrainRatexy','StrainRateyy','StrainRatexz','StrainRateyz','StrainRatezz',
+        'StrainRateeffective',
+        'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset'
+    ]
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.toolkits = md.addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = True
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GEHO'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP17: GlenEFS (enhanced FS)
+# ------------------------------------------------------------------------------
+if org.perform('GlenEFS'):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+
+    md = setflowequation(md, 'FS', 'all')
+    md.materials = matenhancedice(md.materials)
+    md.materials.rheology_E[:] = 5.0
+    md.stressbalance.shelf_dampening = 1
+    md.masstransport.isfreesurface = 1
+    md.groundingline.migration = 'Contact'
+
+    # Example short run
+    md.timestepping.time_step = 0.001
+    md.timestepping.final_time = 0.002
+    md.settings.output_frequency = 1
+
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = True
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_GEFS'
+
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+if org.perform("ESTARSSA"):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+    md = setflowequation(md, 'SSA', 'all')
+    md.materials = matenhancedice(md.materials)
+    md.materials.rheology_Es = 5.0 * np.ones(md.mesh.numberofvertices)
+    md.materials.rheology_Ec= 3.0/8.0 * md,materials.rheology_Es
+    md.transient.requested_outputs = ['default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+		'StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective','LambdaS','Epsprime',
+		'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset']
+    
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_ESTARSSA'
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+if org.perform("ESTARMOLHO"):
+
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+    
+    md.transient.requested_outputs = ['default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+		'VxShear','VyShear','VxBase','VyBase','VxSurface','VySurface','VxAverage','VyAverage',
+		'StrainRatexx','StrainRatexy','StrainRateyy','StrainRateeffective','LambdaS','Epsprime',
+		'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset']
+    md = collapse(md)
+    md = setflowequation(md, 'MOLHO', 'all')
+    md.materials = matestar(md.materials)
+    md.materials.rheology_Es = 5.0 * np.ones(md.mesh.numberofvertices)
+    md.materials.rheology_Ec= 3.0/8.0 * md,materials.rheology_Es
+
+    md.stressbalance.maxiter = 10
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    if res <= 1000:
+        md.toolkits = toolkits()
+        md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
+
+        md.settings.solver_residue_threshold = np.nan  # 11/05/2019
+        md.stressbalance.maxiter = 50                 # 10/24/2019
+        md.stressbalance.restol = 1e-4                # 12/17/2019 original
+        md.stressbalance.reltol = np.nan              # 11/05/2019
+        md.stressbalance.abstol = np.nan              # 11/05/2019
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_ESTARMOLHO'
+    
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# STEP18------------------------------------------------------------------------------
+if org.perform("ESTARHO"):
+    if modelnum == 5 or modelnum == 6:
+        md = org.loadmodel('Transient_Steadystate')
+    else:
+        md = org.loadmodel('Transient_extrude')
+    md = setflowequation(md, 'HO', 'all')
+    md.materials = matestar(md.materials)
+    md.transient.requested_outputs = ['default','IceVolume','IceVolumeAboveFloatation','GroundedArea',
+		'StrainRatexx','StrainRatexy','StrainRateyy','StrainRatexz','StrainRateyz','StrainRatezz',
+		'StrainRateeffective','LambdaS',
+		'MaskOceanLevelset','IceMaskNodeActivation','MaskIceLevelset']
+    md.materials.rheology_Es = 5.0 * np.ones(md.mesh.numberofvertices)
+    md.materials.rheology_Ec= 3.0/8.0 * md,materials.rheology_Es
+
+    md.timestepping.time_step = 1.0/12.0
+    md.timestepping.final_time = 1000
+    md.settings.output_frequency = 600
+    md.toolkits = addoptions(md.toolkits, 'StressbalanceAnalysis', bcgslbjacobioptions())
+    md.stressbalance.maxiter = 30
+    md.stressbalance.restol = 1
+    md.stressbalance.reltol = 0.001
+    md.stressbalance.abstol = float('nan')
+    md.verbose.convergence = False
+    md.cluster = cluster
+    md.miscellaneous.name = 'MISMIP_' + modelname + '_EFS'
+    
+    solutiontype = 'tr'    # 'tr' = transient, 'sb' = stress‑balance, etc.
+    md = solve(md, solutiontype)
+    org.savemodel(md)
+
+# ------------------------------------------------------------------------------
+#  STEP19: (Example) Analysis / Plotting
+# ------------------------------------------------------------------------------
+if org.perform('analyse'):
+    mdgs = org.loadmodel('GlenSSA')
+    mdgm = org.loadmodel('GlenMOLHO')
+    mdgh = org.loadmodel('GlenHO')
+    mdes = org.loadmodel('ESTARSSA')  # e.g., or GlenESSA
+    mdeh = org.loadmodel('ESTARHO')   # e.g., or GlenEHO
+
+    mdgsV = mdgs.results.TransientSolution[-1].Vel
+    mdgmV = mdgm.results.TransientSolution[-1].Vel
+    mdghV = mdgh.results.TransientSolution[-1].Vel
+    mdesV = mdes.results.TransientSolution[-1].Vel
+    mdehV = mdeh.results.TransientSolution[-1].Vel
+
+    plotmodel(mdgs,'data',mdgsV,'data',mdgmV,'data',mdghV(mdeh.mesh.vertexonsurface==1),
+    'data',mdesV,'data',mdehV(mdeh.mesh.vertexonsurface==1),'nlines',5,'ncols',1,'caxis#all',[1,1000])
+    plt.show()
+    # Plot with Python version of plotmodel (if available),
+    # or use your own visualization approach:
+    # plotmodel(mdgs, data=mdgsV, ... ) # This is an example call
+    pass
+
+print("Done with the MISMIP script in Python.")
+
diff --git a/test/run_mismip_first.sh b/test/run_mismip_first.sh
new file mode 100644
index 0000000..3eef7db
--- /dev/null
+++ b/test/run_mismip_first.sh
@@ -0,0 +1,57 @@
+#!/bin/bash
+# run_mismip_first.sh  – ACCESS‑ISSM (vk83) edition
+#
+# Runs the MISMIP workflow just far enough to generate the
+# Glen–vs‑ESTAR comparison plots, skipping the later “enhanced”
+# or full‑Stokes experiments.
+#
+#      Mesh_generation → Parameterization → Transient_Steadystate
+#      (optional) Transient_extrude
+#      GlenSSA / GlenMOLHO / GlenHO
+#      ESTARSSA / ESTARHO
+#      analyse   ← comparison happens here
+#
+# Gadi job directives ---------------------------------------------------------
+#PBS -P au88
+#PBS -q normal
+#PBS -l ncpus=32
+#PBS -l mem=64GB
+#PBS -l walltime=48:00:00
+#PBS -l jobfs=200GB
+#PBS -l wd
+#PBS -N mismip_first
+#PBS -j oe
+#------------------------------------------------------------------------------
+
+set -eu  # fail hard on the first error
+
+# ── 1.  Modules / software stack ─────────────────────────────────────────────
+module purge
+module use /g/data/vk83/modules            # << NEW: expose ACCESS‑ISSM modules
+module load access-issm            #    (or pick a specific version)
+# The access‑issm module sets:
+#   * ISSM_DIR            → /g/data/vk83/...
+#   * PATH / PYTHONPATH   → the matching Python 3.x executable & libs
+#   * LD_LIBRARY_PATH     → PETSc, MUMPS, HDF5, NetCDF, etc.
+
+# ── 2.  ISSM installation root ───────────────────────────────────────────────
+# The module already exports $ISSM_DIR, so the next line is optional and safe
+# to leave out unless you want to override it.
+# export ISSM_DIR=$ISSM_DIR
+
+# ── 3.  Define the “early” organiser steps we wish to execute  ---------------
+early_steps='[3,5,6,9,10,11,12,13,14,15,16,17,18,19]'
+
+# ── 4.  Patch the Python driver on‑the‑fly & run it  -------------------------
+driver=./mismip_driver.py
+tmp=$(mktemp mismipXXXX.py)
+
+cp "$driver" "$tmp"
+sed -i -E "s/^steps\s*=.*/steps = $early_steps/" "$tmp"
+
+echo "Running first‑stage MISMIP workflow …"
+python "$tmp"
+
+rm "$tmp"
+echo "Done – comparison plot and NetCDF outputs are in the Models_* folders."
+
diff --git a/test/run_pps.sh b/test/run_pps.sh
new file mode 100644
index 0000000..07ed7b8
--- /dev/null
+++ b/test/run_pps.sh
@@ -0,0 +1,47 @@
+#!/bin/bash
+# run_pps.sh – Pre‑processing wrapper for ACCESS‑ISSM MISMIP
+
+#PBS -P au88               # ← change to your project code
+#PBS -q normal
+#PBS -l ncpus=8
+#PBS -l mem=32GB
+#PBS -l walltime=04:00:00
+#PBS -l jobfs=50GB
+#PBS -l wd
+#PBS -N mismip_pps
+#PBS -j oe
+
+set -eu
+
+# ── 1.  Modules / software stack ─────────────────────────────────────────────
+module purge
+module use /g/data/vk83/modules            # << NEW: expose ACCESS‑ISSM modules
+module load access-issm            #    (or pick a specific version)
+# The access‑issm module sets:
+#   * ISSM_DIR            → /g/data/vk83/...
+#   * PATH / PYTHONPATH   → the matching Python 3.x executable & libs
+#   * LD_LIBRARY_PATH     → PETSc, MUMPS, HDF5, NetCDF, etc.
+
+# ── 2.  ISSM installation root ───────────────────────────────────────────────
+# The module already exports $ISSM_DIR, so the next line is optional and safe
+# to leave out unless you want to override it.
+# export ISSM_DIR=$ISSM_DIR
+
+export ISSM_DIR=$(spack location -i issm)
+
+# 2. (Optional) tie into a persistent‑session if you started one
+# source ~/.persistent-sessions/bashrc &>/dev/null || true
+
+# 3. Restrict workflow to early steps: Mesh_generation + Parameterization
+early_steps='[1,2]'
+
+driver=./access-issm/examples/mismip/mismip_driver.py
+
+tmp=$(mktemp mismip_ppsXXXX.py)
+cp "$driver" "$tmp"
+sed -i -E "s/^steps\s*=.*/steps = $early_steps/" "$tmp"
+
+python "$tmp"
+rm "$tmp"
+
+echo "PPS finished – inputs in $PBS_O_WORKDIR/Models_*/"