diff --git a/runs/iesm-ngc4008/exp.iesm-ngc4008.run b/runs/iesm-ngc4008/exp.iesm-ngc4008.run
new file mode 100755
index 0000000000000000000000000000000000000000..d4b77da4253d3e89338fdd5deac355812e1536ef
--- /dev/null
+++ b/runs/iesm-ngc4008/exp.iesm-ngc4008.run
@@ -0,0 +1,549 @@
+#! /bin/ksh
+#=============================================================================
+#SBATCH --account=p72044
+#SBATCH --partition=skylake_0096
+#SBATCH --qos=skylake_0096
+###SBATCH --qos=skylake_0096_devel
+#SBATCH --job-name=iesm-ngc4008
+#SBATCH --nodes=3 
+###SBATCH --nodes=1 
+#SBATCH --ntasks-per-node=48
+#SBATCH --ntasks-per-core=1
+#SBATCH --output=/gpfs/data/fs72044/avoigt_teach/climlab_s2025/msc-climate-modeling-lab-s2025-code-4-students/runs/iesm-ngc4008/LOG.exp.iesm-ngc4008.run.%j.o
+#SBATCH --error=/gpfs/data/fs72044/avoigt_teach/climlab_s2025/msc-climate-modeling-lab-s2025-code-4-students/runs/iesm-ngc4008/LOG.exp.iesm-ngc4008.run.%j.o
+#SBATCH --exclusive
+#SBATCH --time=02:30:00
+###SBATCH --time=00:05:00
+#SBATCH --mail-user=aiko.voigt@univie.ac.at
+#SBATCH --mail-type=BEGIN,END,FAIL
+
+set -x # debugging command: enables a mode of the shell where all executed commands are printed to the terminal
+ulimit -s unlimited # unsets limits for RAM
+
+# MPI variables
+# -------------
+no_of_nodes=3
+mpi_procs_pernode=48
+((mpi_total_procs=no_of_nodes * mpi_procs_pernode))
+
+# manual fix for mpi pinning with intel mpirun on vsc
+# see https://wiki.vsc.ac.at/doku.php?id=doku:vsc5quickstart#intel_mpi
+export I_MPI_PIN_RESPECT_CPUSET=0
+
+#
+# blocking length
+# ---------------
+nproma=16
+
+#=============================================================================
+# Input variables:
+
+# SIMULATION NAME
+EXP=iesm-ngc4008
+
+ICONFOLDER=/gpfs/data/fs72044/avoigt_teach/climlab_s2025/icon-esm-univie       # DIRECTORY OF ICON MODEL CODE
+RUNSCRIPTDIR=/gpfs/data/fs72044/avoigt_teach/climlab_s2025/msc-climate-modeling-lab-s2025-code-4-students/runs/iesm-ngc4008/
+basedir=$ICONFOLDER # icon base directory
+
+. ${ICONFOLDER}/run/add_run_routines
+
+# experiment directory, with plenty of space, create if new
+EXPDIR=/gpfs/data/fs72044/avoigt_teach/experiments/s2025/${EXP}
+if [ ! -d ${EXPDIR} ] ;  then
+  mkdir -p ${EXPDIR}
+fi
+#
+ls -ld ${EXPDIR}
+if [ ! -d ${EXPDIR} ] ;  then
+    mkdir ${EXPDIR}
+fi
+ls -ld ${EXPDIR}
+
+cd $EXPDIR
+
+#=================================================================================
+
+#-----------------------------------------------------------------------------
+# global timing
+# NOTE: WE NEED TO START IN YEAR 2021 AS SST AND SIC START IN 2020, AND THE MODEL LOOKS 1 YEAR BACKWARD FOR THESE FIELDS
+initial_date="2021-01-01"
+final_date="2050-01-01"
+start_date=$initial_date
+end_date=$final_date
+y0=${start_date%%-*}
+yN=${end_date%%-*}
+
+# restart intervals
+restart_interval="P2Y"
+checkpoint_interval="P1Y"
+
+file_interval="P1Y"
+
+############################################################
+#
+# NO FURTHER CHANGES TO THE DIRECTORIES AND SIMULATION NAME
+# SHOULD BE NEEDED BELOW THIS LINE
+#
+############################################################
+
+#-----------------------------------------------------------------------------
+# Provide input files
+# $Id: format.tmpl 9264 2021-06-21 21:24:57Z m221078 $
+#
+# [files]
+
+# [files.atmosphere]
+data_dir=/gpfs/data/fs72044/avoigt_teach/ICON-inputdata/amip-VSC4
+data_dir_ngc4008=/gpfs/data/fs72044/avoigt_teach/ICON-inputdata/ngc4008_for_iconesm
+
+# [files.atmosphere.mapped]
+grid_dir=$data_dir/grid
+ln -sfv $grid_dir/icon_grid_0013_R02B04_G.nc icon_grid_G.nc
+
+# [files.atmosphere.mapped.initial]
+initial_dir=$data_dir/initial_condition
+ln -sfv $initial_dir/ifs2icon_1979010100_R02B04_G.nc ifs2icon.nc
+
+# [files.atmosphere.mapped.ozone]
+ozone_dir=$data_dir_ngc4008/ozone
+for((yr = y0 + -1; yr <= yN + 1; ++yr)) 
+do
+    ln -sfv $ozone_dir/bc_ozone_ssp370_${yr}.remapcon_grid_0013_R02B04_G.nc bc_ozone_${yr}.nc
+done
+
+# [files.atmosphere.mapped.ocean_surface]
+ocean_surface_dir=$data_dir_ngc4008/sst_and_seaice
+ln -sfv $ocean_surface_dir/siconcbcs.ngc4008.2020-2050.r2b4.nc bc_sic.nc
+ln -sfv $ocean_surface_dir/tosbcs.ngc4008.2020-2050.r2b4.nc bc_sst.nc
+
+# [files.atmosphere.mapped.aerosols]
+aerosols_dir=$data_dir/aerosol
+ln -sfv $aerosols_dir/bc_aeropt_kinne_lw_b16_coa.nc .
+ln -sfv $aerosols_dir/bc_aeropt_kinne_sw_b14_coa.nc .
+ln -sfv $aerosols_dir/bc_aeropt_kinne_sw_b14_fin_1850.nc bc_aeropt_kinne_sw_b14_fin.nc
+
+# [files.atmosphere.model]
+model_dir=/gpfs/data/fs72044/avoigt_teach/climlab_s2025/icon-esm-univie
+
+# [files.atmosphere.model.data]
+ln -sfv $model_dir/data/lsdata.nc .
+ln -sfv $model_dir/data/ECHAM6_CldOptProps.nc .
+cp $data_dir_ngc4008/aerosols/MACv2.0-SP-merged-historical-and-SSP3-70_v1.nc MACv2.0-SP_v1.nc
+
+# [files.atmosphere.model.run]
+run_dir=$model_dir/run
+cp -fv $run_dir/dict.iconam.mpim dict.txt
+
+# [files.atmosphere.independent]
+ln -sfv $data_dir_ngc4008/greenhouse_gases/greenhouse_ssp370.nc bc_greenhouse_gases.nc
+
+# [files.atmosphere.independent.volcano_aerosols]
+for((yr = y0 + -1; yr <= yN + 1; ++yr)) 
+do
+    label=${yr}
+        ((yr >= 2015)) && label=2014
+    volcano_aerosols_dir=$data_dir/aerosol
+    ln -sfv $volcano_aerosols_dir/bc_aeropt_cmip6_volc_lw_b16_sw_b14_${label}.nc \
+        bc_aeropt_cmip6_volc_lw_b16_sw_b14_${yr}.nc
+done # offsets
+
+# [files.atmosphere.independent.solar_irradiance]
+solar_irradiance_dir=$data_dir/solar_radiation
+ln -sfv $solar_irradiance_dir/swflux_14band_cmip6_1850-2299-v3.2.nc \
+    bc_solar_irradiance_sw_b14.nc
+
+# [files.land]
+land_dir=$data_dir/land
+
+# [files.land.mapped]
+ln -sfv $land_dir/ic_land_soil_1976.nc ic_land_soil.nc
+ln -sfv $land_dir/bc_land_frac_11pfts_1976.nc bc_land_frac.nc
+ln -sfv $land_dir/bc_land_phys_1976.nc bc_land_phys.nc
+ln -sfv $land_dir/bc_land_soil_1976.nc bc_land_soil.nc
+ln -sfv $land_dir/bc_land_sso_1976.nc bc_land_sso.nc
+
+# [files.land.hydro]
+hydro_dir=$land_dir
+# preliminary test version
+ln -sfv $hydro_dir/hdpara_r2b4_0013_0035_v3.nc bc_land_hd.nc
+
+# [files.land.model]
+model_dir=$basedir/externals/jsbach/data
+ln -sfv $model_dir/lctlib_nlct21.def .
+
+
+#-----------------------------------------------------------------------------
+# automatic restart setup
+# set some default values and derive some run parameteres
+restart=${restart:=".false."}
+restartSemaphoreFilename='isRestartRun.sem'
+#AUTOMATIC_RESTART_SETUP:
+if [ -f ${restartSemaphoreFilename} ]; then
+  restart=.true.
+  #  do not delete switch-file, to enable restart after unintended abort
+  #[[ -f ${restartSemaphoreFilename} ]] && rm ${restartSemaphoreFilename}
+fi
+#END AUTOMATIC_RESTART_SETUP
+
+# wait 5min to let GPFS finish the write operations
+if [ "x$restart" != 'x.false.' -a "x$submit" != 'x' ]; then
+  if [ x$(df -T ${EXPDIR} | cut -d ' ' -f 2) = gpfs ]; then
+    sleep 10;
+  fi
+fi
+
+#
+# create ICON master, coupling and model namelists
+# ------------------------------------------------
+# For a complete list see Namelist_overview and Namelist_overview.pdf
+#
+
+cat > icon_master.namelist << EOF
+&master_nml
+    lrestart =  ${restart}
+/
+&master_time_control_nml
+    calendar = 'proleptic gregorian'
+    checkpointtimeintval = '$checkpoint_interval'
+    restarttimeintval = '$restart_interval'
+    experimentstartdate = '2021-01-01' ! TODO: hack to reproduce result, NOTE: WE NEED TO START IN YEAR 2021 AS SST AND SIC START IN 2020, AND THE MODEL LOOKS 1 YEAR BACKWARD FOR THESE FIELDS
+    experimentstopdate = '$final_date'
+/
+&master_model_nml ! 'atmo'
+    model_name = 'atmo'
+    model_namelist_filename = 'NAMELIST_atm'
+    model_type = 1
+/
+&jsb_control_nml
+    is_standalone = .false.
+/
+&jsb_model_nml
+    model_id = 1
+    model_name = 'JSBACH'
+    model_shortname = 'jsb'
+    model_description = 'JSBACH land surface model'
+    model_namelist_filename = 'NAMELIST_lnd'
+/
+EOF
+
+#-----------------------------------------------------------------------------
+# II. ATMOSPHERE and LAND
+#-----------------------------------------------------------------------------
+#
+# atmosphere namelist
+# -------------------
+cat > NAMELIST_atm << EOF
+&parallel_nml
+    nproma = $nproma
+    num_io_procs = 0
+    num_prefetch_proc = 0
+    pio_type = 1 !1=default, assync I/O, 2=experimental CDI, 0=nothing?
+/
+&grid_nml
+    dynamics_grid_filename = 'icon_grid_G.nc'
+/
+&run_nml
+    num_lev = 47 ! number of full levels
+    modeltimestep = 'PT15M'
+    ltestcase = .false. ! run testcase
+    ldynamics = .true. ! dynamics
+    ltransport = .true. ! transport
+    iforcing = 2 ! 0: none, 1: HS, 2: ECHAM, 3: NWP
+    output = 'nml'
+    msg_level = 8 ! level of details report during integration
+    restart_filename = 'AMIP-VSC4_restart_atm_<rsttime>.nc'
+    activate_sync_timers = .true.
+/
+&extpar_nml
+    itopo = 1 ! 1: read topography from the grid file
+    itype_lwemiss = 0
+/
+&initicon_nml
+    init_mode = 2 ! 2: initialize from IFS analysis
+    ifs2icon_filename = 'ifs2icon.nc'
+/
+&nonhydrostatic_nml
+    ndyn_substeps = 8 ! dtime/dt_dyn
+    damp_height = 50000. ! [m]
+    rayleigh_coeff = 0.1000 ! set to 0.1001 for rerun with little change
+    vwind_offctr = 0.2
+    divdamp_fac = 0.004
+/
+&interpol_nml
+    rbf_scale_mode_ll = 1
+/
+&sleve_nml
+    min_lay_thckn = 40. ! [m]
+    top_height = 83000. ! [m]
+    stretch_fac = 0.9
+    decay_scale_1 = 4000. ! [m]
+    decay_scale_2 = 2500. ! [m]
+    decay_exp = 1.2
+    flat_height = 16000. ! [m]
+/
+&diffusion_nml
+/
+&transport_nml
+    ihadv_tracer = 52, 2, 2
+    itype_hlimit = 3, 4, 4
+    ivadv_tracer = 3, 3, 3
+    tracer_names = 'hus', 'clw', 'cli'
+/
+&echam_phy_nml
+    ! domain 1
+    ! atmospheric physics ("" = never)
+    echam_phy_config(1)%dt_rad = 'PT90M'
+    echam_phy_config(1)%dt_vdf = 'PT15M'
+    echam_phy_config(1)%dt_cnv = 'PT15M'
+    echam_phy_config(1)%dt_cld = 'PT15M'
+    echam_phy_config(1)%dt_gwd = 'PT15M'
+    echam_phy_config(1)%dt_sso = 'PT15M'
+    ! atmospheric chemistry ("" = never)
+    echam_phy_config(1)%dt_mox = 'PT15M'
+    ! surface (true or false)
+    echam_phy_config(1)%ljsb = .true.
+    echam_phy_config(1)%lamip = .true.
+    echam_phy_config(1)%lice = .true.
+    echam_phy_config(1)%lmlo = .false.
+    echam_phy_config(1)%llake = .true.
+/
+&echam_rad_nml
+    ! domain 1
+    echam_rad_config(1)%isolrad = 1
+    echam_rad_config(1)%irad_h2o = 1
+    echam_rad_config(1)%irad_co2 = 4  
+    echam_rad_config(1)%irad_ch4 = 4
+    echam_rad_config(1)%irad_n2o = 4
+    echam_rad_config(1)%irad_o3 = 8 
+    echam_rad_config(1)%irad_o2 = 2
+    echam_rad_config(1)%irad_cfc11 = 4 
+    echam_rad_config(1)%irad_cfc12 = 4 
+    echam_rad_config(1)%irad_aero = 18
+/
+&echam_gwd_nml
+/
+&echam_sso_nml
+/
+&echam_vdf_nml
+/
+&echam_cnv_nml
+/
+&echam_cld_nml
+/
+&echam_cov_nml
+/
+&ccycle_nml
+/
+&sea_ice_nml
+/
+! Parameters for all output files
+! -------------------------------
+&io_nml
+    output_nml_dict = 'dict.txt'
+    netcdf_dict = 'dict.txt'
+    itype_pres_msl = 4
+    ! restart_file_type = 5
+    ! restart_write_mode = 'joint procs multifile' ! not useful in r2b4 setup
+    ! lnetcdf_flt64_output = .true. ! 64 bit output in all files
+    ! lkeep_in_sync = .true. ! sync after each timestep
+    write_initial_state = .false.
+/
+&dbg_index_nml
+    idbg_mxmn = 0 ! initialize MIN/MAX debug output
+    idbg_val = 0 ! initialize one cell debug output
+    idbg_slev = 1 ! initialize start level for debug output
+    idbg_elev = 2 ! initialize end level for debug output
+    dbg_lat_in = 30.0 ! latitude location of one cell debug output
+    dbg_lon_in = -30.0 ! longitude location of one cell debug output
+    str_mod_tst = 'InterFaceOce' ! define modules to print out in debug mode
+/
+
+! Define output files
+! -------------------
+!
+! 3-dimensional files include 'ps' and 'pfull' to allow the vertical
+! interpolation to pressure levels by cdo ap2pl.
+
+/
+! Standard AMIP output...
+&output_nml ! 'atm_3d'
+    output_filename = '${EXP}_atm_3d'
+    filename_format = '<output_filename>_<levtype_l>_<datetime2>'
+    filetype = 5
+    remap = 0
+    operation = 'mean'
+    output_grid = .TRUE.
+    output_start = '${initial_date}'
+    output_end = '${final_date}'
+    output_interval = 'P1M'
+    file_interval = '$file_interval' !'P10YT1S'
+    include_last = .false.
+    ml_varlist = 'zg', 'ps', 'pfull', 'rho', 'ta', 'ua', 'va', 'wap', 'hus',
+                 'clw', 'cli', 'hur', 'cl'
+/
+! Standard AMIP output...
+&output_nml ! 'atm_2d'
+    output_filename = '${EXP}_atm_2d'
+    filename_format = '<output_filename>_<levtype_l>_<datetime2>'
+    filetype = 5
+    remap = 0
+    operation = 'mean'
+    output_grid = .TRUE.
+    output_start = '${initial_date}'
+    output_end = '${final_date}'
+    output_interval = 'P1M'
+    file_interval = '$file_interval' !'P10YT1S'
+    include_last = .false.
+    ml_varlist = 'orog', 'ps', 'psl', 'cosmu0', 'rsdt', 'rsut', 'rsutcs',
+                 'rlut', 'rlutcs', 'rsds', 'rsdscs', 'rlds', 'rldscs', 'rsus',
+                 'rsuscs', 'rlus', 'ts', 'sic', 'sit', 'albedo', 'clt', 'prlr',
+                 'prls', 'prcr', 'prcs', 'pr', 'prw', 'cllvi', 'clivi', 'hfls',
+                 'hfss', 'evspsbl', 'tauu', 'tauv', 'tauu_sso', 'tauv_sso',
+                 'diss_sso', 'sfcwind', 'uas', 'vas', 'tas', 'dew2', 'ptp'
+/
+! As an example, included output for rsdt on lat-lon grid for every time step for January 1979
+&output_nml ! 'rsdt_latlon_everytimestep'
+    output_filename = '${EXP}_rsdt_latlon_everytimestep'
+    filename_format = '<output_filename>_<levtype_l>_<datetime2>'
+    filetype = 5
+    remap            = 1                         ! 1: latlon,  0: native grid
+    reg_lat_def      = -90.0,0.5,90.0
+    reg_lon_def      = -180.0,0.5,180.0
+    operation = 'false'
+    output_grid = .TRUE.
+    output_start = '1979-01-01'
+    output_end = '1979-02-01'
+    output_interval = 'PT15M'
+    file_interval = 'P1D'
+    include_last = .false.
+    ml_varlist = 'rsdt'
+/
+EOF
+
+
+# jsbach namelist
+# ---------------
+
+cat > NAMELIST_lnd << EOF
+&jsb_model_nml
+    usecase = 'jsbach_pfts'
+    use_lakes = .true.
+/
+&jsb_seb_nml
+    bc_filename = 'bc_land_phys.nc'
+    ic_filename = 'ic_land_soil.nc'
+/
+&jsb_rad_nml
+    use_alb_veg_simple = .false. ! if jsbach_lite
+    bc_filename = 'bc_land_phys.nc'
+    ic_filename = 'ic_land_soil.nc'
+/
+&jsb_turb_nml
+    bc_filename = 'bc_land_phys.nc'
+    ic_filename = 'ic_land_soil.nc'
+/
+&jsb_sse_nml
+    l_heat_cap_map = .false.
+    l_heat_cond_map = .false.
+    l_heat_cap_dyn = .false.
+    l_heat_cond_dyn = .false.
+    l_snow = .true.
+    l_dynsnow = .true.
+    l_freeze = .true.
+    l_supercool = .true.
+    bc_filename = 'bc_land_soil.nc'
+    ic_filename = 'ic_land_soil.nc'
+/
+&jsb_hydro_nml
+    l_organic = .false.
+    bc_filename = 'bc_land_soil.nc'
+    ic_filename = 'ic_land_soil.nc'
+    bc_sso_filename = 'bc_land_sso.nc'
+/
+&jsb_assimi_nml
+    active = .true. ! if jsbach_pfts
+/
+&jsb_pheno_nml
+    scheme = 'logrop' ! 'climatology' if jsbach_lite
+    bc_filename = 'bc_land_phys.nc'
+    ic_filename = 'ic_land_soil.nc'
+/
+&jsb_carbon_nml
+    active = .true.
+    bc_filename = 'bc_land_carbon.nc'
+    ic_filename = 'ic_land_carbon.nc'
+    read_cpools = .false.
+/
+&jsb_fuel_nml
+    active = .true.
+    fuel_algorithm = 1
+/
+&jsb_disturb_nml
+    active = .false.
+    ic_filename = 'ic_land_soil.nc'
+    bc_filename = 'bc_land_phys.nc'
+    fire_algorithm = 1
+    windbreak_algorithm = 1
+    lburn_pasture = .false.
+/
+EOF
+
+
+## setup for status check & restart
+final_status_file=${EXPDIR}/${EXP}.final_status
+
+## Copy icon executable to working directory
+cp -p $ICONFOLDER/bin/icon ./icon.exe
+##
+
+## Start model
+date
+ulimit -s unlimited
+
+ldd icon.exe
+
+START="/gpfs/opt/sw/skylake/spack-0.19.0/opt/spack/linux-almalinux8-skylake_avx512/intel-2021.7.1/intel-oneapi-mpi-2021.7.1-fzg6q4xcj7efjmce3cuqa2b7cum5d3po/mpi/2021.7.1/bin/mpiexec -n $mpi_total_procs"
+MODEL=${EXPDIR}/icon.exe
+
+rm -f finish.status
+
+${START} ${MODEL}
+
+if [ -r finish.status ] ; then
+  check_final_status 0 "${START} ${MODEL}"
+else
+  check_final_status -1 "${START} ${MODEL}"
+fi
+
+#-----------------------------------------------------------------------------
+finish_status=`cat finish.status`
+echo $finish_status
+
+#-----------------------------------------------------------------------------
+namelist_list=""
+#-----------------------------------------------------------------------------
+# check if we have to restart, ie resubmit
+#   Note: this is a different mechanism from checking the restart
+if [ $finish_status = "RESTART" ]; then
+  echo "restart next experiment..."
+  this_script="${RUNSCRIPTDIR}/exp.${EXP}.run"
+  echo 'this_script: ' "$this_script"
+  # note that if ${restartSemaphoreFilename} does not exist yet, then touch will create it
+  touch ${restartSemaphoreFilename}
+  cd ${RUNSCRIPTDIR}
+  sbatch exp.${EXP}.run
+else
+  [[ -f ${restartSemaphoreFilename} ]] && rm ${restartSemaphoreFilename}
+fi
+
+#-----------------------------------------------------------------------------
+
+cd ${RUNSCRIPTDIR}
+
+#-----------------------------------------------------------------------------
+#
+
+echo "============================"
+echo "Script run successfully: ${finish_status}"
+echo "============================"
+#-----------------------------------------------------------------------------