#!/bin/bash

# ICON
#
# ------------------------------------------
# Copyright (C) 2004-2026, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
# Contact information: icon-model.org
# See AUTHORS.TXT for a list of authors
# See LICENSES/ for license information
# SPDX-License-Identifier: BSD-3-Clause
# ------------------------------------------

# ----------------------------------------------------------------------------
#
# Coupled ocean-atmosphere run with arbitrary grids with hydrological discharge (HD) through YAC.
#
# Notes: - turn on debug stdout with namelists dbg_index_nml
#        - PBS options: NE=1-64, CPE=8
#        - sample script for new users to provide a quick start for running
#          a full featured atm/oce/hd coupled model with 80/40km resolution
#          and pyICON output for postprocssing (NOT for buildbot use)
#
# Martin Koehler, Roland Wirth, Trang Van Pham, Stephan Lorenz, Rene Redler, Spring 2024
#
# ----------------------------------------------------------------------------
#PBS -S /bin/bash
#PBS -q sx_norm
#PBS -v NE=32,CPE=8          # NE: number VE nodes, CPE: number cores on each VE (8!)
#PBS --venode=${NE}          # number VEs requested
#PBS --venum-lhost=1         # number VEs per logical host (partitioning, 2 OK)
#PBS -l cpunum_job=3         # number x86 cores per logical host (1ini+#IO)
#PBS -l coresz_prc=0
#PBS -T necmpi_hydra         # when compiling with mpi/2.21.0 (see VE config script)
#PBS --use-hca=2             # number of HCA per logical host (2)
#PBS -l elapstim_req=24:00:00
##PBS -o /hpc/uwork/mkoehler/wq/coupledR2B5-R2B6.%s
#PBS -j o
# ----------------------------------------------------------------------------

module load mpi/3.5.0
set -x
ulimit -s unlimited
ulimit -s 2097152
ulimit -c 0

. /etc/profile


# environment variables for the experiment and the target system
# --------------------------------------------------------------
export EXPNAME="couple_001"
job_name="run_aR2B5_oR2B6_coupled"                          # name of run script
RUNSCRIPTDIR=/hpc/uhome/${USER}/icon/run-icon/coupled       # directory of run script
ICONDIR=$HOME/icon/icon-nwp
WORK_DIR=/hpc/uwork/${USER}

pi_control='no'                                             # yes/no PI-control

# start and end date+time of experiment
# -------------------------------------
start_date=${start_date:="1979-01-01T00:00:00Z"}
    end_date=${end_date:="1979-01-11T00:00:00Z"}
#   end_date=${end_date:="1980-01-01T00:00:00Z"}
#   end_date=${end_date:="1990-01-01T00:00:00Z"}
#start_date=${start_date:="1850-01-01T00:00:00Z"}
#    end_date=${end_date:="2350-01-01T00:00:00Z"}


# restart=checkpoint intervals (= file length)
# --------------------------------------------
restart_interval="P1Y"          #  no restart used for NWP
checkpoint_interval="P1Y"       #  write checkpoint
#restart_interval="P100Y"       #  no restart used for NWP
#checkpoint_interval="P100Y"    #  write checkpoint
#restart_interval="P2D"         #  no restart used for NWP
#checkpoint_interval="P2D"      #  write checkpoint

# file interval - not used, so select longer than restart-interval
# (make sure to delete the partial nc files but not the restart)
atm_file_interval="P100Y"
oce_file_interval="P100Y"

# output intervals
# ----------------
#atm_output_interval="PT450S"   # "PT450S"
#oce_output_interval="PT30M"    # "PT2H"
atm_output_interval="P1D"       # "PT30M"
oce_output_interval="P1D"       # "PT1H"

# levels: L64, L72 OR L128SMT
LEV=L72

#=============================================================================
#
# OpenMP environment variables
# ----------------------------
export OMP_NUM_THREADS=1         # openMPI deactivated - best (needs compiling)
export ICON_THREADS=1
export OMP_SCHEDULE=dynamic,1
export OMP_DYNAMIC="false"
export OMP_STACKSIZE=200M
export NMPI_PROGINF=DETAIL       # get performance details (Jens-Olaf Beismann)

# Run information (reserve last NE for HD)
let NE1=${NE}-1
let PPN=${CPE}*${NE1}/${OMP_NUM_THREADS}         # no. of MPI procs. per NQS job
echo "PPN etc. $NN $NE $CPE $PPN"

#
# MPI variables
# -------------

# NEC has PPN=NE*CPE mpi processes (OMP_NUM_THREADS=1)
# they will be devided in ocean and atmosphere

mpi_total_procs=${PPN}

# at least one node for ocean and atmosphere each is necessary
# default: half of requested nodes for ocean

#((mpi_oce_procs=mpi_total_procs/2))
mpi_oce_procs=192       #112 #448
mpi_hd_procs=1
let NE2=${NE1}-${mpi_hd_procs}
((mpi_atm_procs=PPN-mpi_oce_procs))

echo $mpi_total_procs ${mpi_atm_procs} ${mpi_oce_procs} ${mpi_hd_procs}

nproma_atm=752
nproma_oce=752

# asynchronous diagnostic output processes
# ----------------------------------------

# Note that "mpi_atm_io_procs" must match the number of output files
mpi_atm_io_procs=2      # >0 for atmosphere plus land (not working for monitoring)
mpi_oce_io_procs=0      # >0 for ocean is not working yet

#=============================================================================


# load local setting and profile, if existing
# -------------------------------------------
if [ -a ../setting ]
then
  echo "Load Setting"
  . ../setting
fi
if [[ -a  /etc/profile ]]
then
	. /etc/profile
fi

#=============================================================================

# directories with absolute paths
# -------------------------------
rundir=${ICONDIR}/run
basedir=${ICONDIR}
export basedir

# experiment directory, with plenty of space, create if new
EXPDIR=${WORK_DIR}/run-icon/experiments/${EXPNAME}
if [ ! -d ${EXPDIR} ] ;  then
  mkdir -p ${EXPDIR}
fi
if [ ! -d ${EXPDIR} ] ;  then
    mkdir ${EXPDIR}
fi
ls -ld ${EXPDIR}
check_error $? "${EXPDIR} does not exist?"

cd ${EXPDIR}

rm -rf stdout.* stderr.*      # remove standard output from previous run

# path to model binary, including the executable:
MODEL_VE_orig=${ICONDIR}/build/VE/bin/icon
MODEL_VH_orig=${ICONDIR}/build/VH/bin/icon
#MODEL_HD_orig=${ICONDIR}/build/HD/bin/hd_05.exe    # located in VE/externals/hd/bin/hd_05.exe
MODEL_HD_orig=${ICONDIR}/build/VE/externals/hd/bin/hd_05.exe
MODEL_VE=${EXPDIR}/bin/icon.VE
MODEL_VH=${EXPDIR}/bin/icon.VH
MODEL_HD=${EXPDIR}/bin/icon.HD

# copy binaries to permanent spot during first start to enforce consistent binaries for all restarts
if [ ! -f "$EXPDIR/isRestartRun.sem" ]; then
  mkdir -p bin
  \cp -f ${MODEL_VE_orig} ${MODEL_VE}
  \cp -f ${MODEL_VH_orig} ${MODEL_VH}
  \cp -f ${MODEL_HD_orig} ${MODEL_HD}
fi

# how to submit the next job
# --------------------------
#submit="sbatch"
submit="qsub"

# cdo for post-processing
# -----------------------
cdo="cdo"
cdo_diff="cdo diffn"

# define script functions used in the experiment run script
# ---------------------------------------------------------
. ${ICONDIR}/run/add_run_routines

#=============================================================================
#
# Basic test case for the ICON-NWP coupled to the ICON-O.
#
# This section of the run script contains the specifications of the experiment.
# The specifications are passed by namelist to the program.
# For a complete list see Namelist_overview.pdf
#
# Initial revision: 2021-01-27, S. Lorenz, MPI-M
#
#-----------------------------------------------------------------------------
#
# Basic specifications of the simulation
# --------------------------------------
#
# These variables are set in the header section of the completed run script:
#
# EXPNAME = experiment name
# NPROMA  = array blocking length / inner loop length
#-----------------------------------------------------------------------------

# (0) Basic model configuration
# -----------------------------

atmos_gridID="0030"
atmos_refinement="R02B05"

rad_gridID="0029"
rad_refinement="R02B04"

ocean_gridID="0035"
ocean_refinement="R02B06"
#
#-----------------------------------------------------------------------------

# namelist files
# --------------
atmo_namelist=NAMELIST_${EXPNAME}_atm
ocean_namelist=NAMELIST_${EXPNAME}_ocean

#--------------------------------------------------------------------------------------------------
#
# (1) Define the model time stepping
# ----------------------------------
#
# NWP:                           # 360s or 450s for R2B5
dtime=450                        # NWP atmospheric timestep (s)   (same as in atmTimeStep!!)
oceTimeStep="PT30M"              # ocean time step
atmTimeStep="PT450S"             # atmos time step (for coupler)  (same as dtime!!)
rivTimeStep="P01D"               # for YAC /same as riv_dt)
dtime_riv=86400                  # river timestep

ocecouplingTimeStep="PT30M"      # coupling time step
rivcouplingTimeStep="P01D"       # coupling time step atm->riv->oce
longcouplingTimeStep="P100Y"     # no coupling (long coupling period)

#--------------------------------------------------------------------------------------------------

# (2) unset some setting of create_target_header for mistral

unset OMPI_MCA_coll_fca_enable
unset OMPI_MCA_coll_fca_priority

#--------------------------------------------------------------------------------------------------

# (3) Set variables to configure the experiment:
# ----------------------------------------------

# use alternatively: full restart, init from (ocean) restart, init from climatology or false
#  - set to false if ocean was already initialised
#initialiseOcean="fromRestart"              # first run
initialiseOcean="fromClimatology"           # from Levitus observations (current time, mean last 10 years)
# initialiseOcean="fromRestartwithHamocc"
# initialiseOcean="coldStart"

# use both lines for atm/oce initialisation from other experiment
#initialiseOcean="FALSE"                    # restart run (maybe copy from other run)
#restart=.true.                             # deactivate semaphore mechanism
#restart_exp=slo1325
#restart_yea=4000
restart_exp=mk056_R2B6L72                   # mk099_R2B6L128SMT  mk050_R2B6L64_1000y+1monL72
restart_yea=20200101T00                     # 19790101T00        20000101T00

# set initialiseOcean automatically to false via isRestart.sem
#if [ -d "$basedir/experiments/$EXPNAME" ] ; then
#  if [ -f "$basedir/experiments/$EXPNAME/isRestartRun.sem" ]; then
#    initialiseOcean="FALSE"
#  fi
#fi
if [ -d "$EXPDIR" ] ; then
  if [ -f "$EXPDIR/isRestartRun.sem" ]; then
    initialiseOcean="FALSE"
  fi
fi

# output file selection
# ---------------------
#
# output_<xyz>=yes : yes --> output files for <xyz>, any other value --> no files for <xyz>

# monitoring output switches
output_atm_mon=yes
output_oce_mon=yes
output_oce_moc=yes       #  calc_moc is called every timestep - switch off for high resolutions
output_oce_def=yes       #  annual mean is sufficient
#oce_output_interval_def="P1Y"  # annual mean output for deep ocean
oce_output_interval_def=$oce_output_interval
oce_output_interval_mon=$oce_output_interval

# sea ice, snow, SST and some deeper levels ocean output
output_oce_ice=yes

# standard debug-flux arrays at surface
output_atm_dbg=yes
output_oce_dbg=yes

# output switches for large atm/oce files - set to "yes" if needed for standard quickplots
#  - yes writes much output and needs more computing time at least in high resolution


#--------------------------------------------------------------------------------------------------
# I. coupling section
#--------------------------------------------------------------------------------------------------

if [ $mpi_total_procs -lt 2 ] ; then
  check_error 1 "This setup requires at least 2 mpi processes. Exit"
fi

# I.1 Split the number of total procs and assign to each component
# ----------------------------------------------------------------

oce_min_rank=`expr ${mpi_total_procs} - ${mpi_oce_procs} + ${mpi_atm_io_procs} + 1`
oce_max_rank=`expr ${oce_min_rank} + ${mpi_oce_procs} - 1`
oce_inc_rank=1
atm_min_rank=0
atm_max_rank=`expr ${oce_min_rank} - 1`
atm_inc_rank=1

# I.2 Fill model list
# -------------------

namelist_list[0]="$atmo_namelist"
modelname_list[0]="atmo"
modeltype_list[0]=1
minrank_list[0]=$atm_min_rank
maxrank_list[0]=$atm_max_rank
incrank_list[0]=$atm_inc_rank

namelist_list[1]="$ocean_namelist"
modelname_list[1]="ocean"
modeltype_list[1]=2
minrank_list[1]=$oce_min_rank
maxrank_list[1]=$oce_max_rank
incrank_list[1]=$oce_inc_rank

modelname_list[2]="HD"

# I.3 YAC coupling library configuration
# --------------------------------------

atm_lag=1
oce_lag=1
riv_lag=1

# co2_flux and co2_mixing_ratio are listed as transients
# below but are not configured for coupling in the couples
# section. Their fields are therefore not considered for
# the search not for the data exchange even if yac_fget and
# yac_fput are called for these fields.

cat > coupling_${EXPNAME}.yaml << EOF
definitions:
  components:
    atm: &atm
      component: ${modelname_list[0]}
      grid: icon_atmos_grid
    oce: &oce
      component: ${modelname_list[1]}
      grid: icon_ocean_grid
    riv: &riv
      component: ${modelname_list[2]}
      grid: HD_GRID
  atm2oce: &atm2oce
    source: *atm
    target: *oce
    time_reduction: average
    src_lag: ${atm_lag}
    tgt_lag: ${oce_lag}
  oce2atm: &oce2atm
    source: *oce
    target: *atm
    time_reduction: average
    src_lag: ${oce_lag}
    tgt_lag: ${atm_lag}
  atm2riv: &atm2riv
    source: *atm
    target: *riv
    time_reduction: average
    src_lag: ${atm_lag}
    tgt_lag: ${riv_lag}
  riv2oce: &riv2oce
    source: *riv
    target: *oce
    time_reduction: average
    src_lag: ${riv_lag}
    tgt_lag: ${oce_lag}
  interp_stacks:
    hcsbb_interp_stack: &hcsbb_interp_stack
      interpolation:
        - bernstein_bezier
        - nnn:
            n: 4
            weighted: arithmetic_average
        - fixed:
            user_value: -999.9
    conserv_interp_stack: &conserv_interp_stack
      interpolation:
        - conservative:
            order: 1
            enforced_conservation: false
            partial_coverage: true
            normalisation: fracarea
        - fixed:
            user_value: -999.9
    conserv_interp_dest: &conserv_interp_dest
      interpolation:
        - conservative:
            order: 1
            enforced_conservation: false
            partial_coverage: true
            normalisation: destarea
    spmap_interp_stack: &spmap_interp_stack
      interpolation:
        - source_to_target_map:
            spread_distance: 0.0
            max_search_distance: 0.0
        - fixed:
            user_value: 0.0

timestep_unit: ISO_format
calendar: proleptic-gregorian
coupling:
  - <<: [ *atm2oce, *hcsbb_interp_stack ]
    coupling_period: ${ocecouplingTimeStep}
    field: [surface_downward_eastward_stress,
            surface_downward_northward_stress]
#  - <<: [ *atm2oce, *hcsbb_interp_stack ]
#    coupling_period: ${ocecouplingTimeStep}
#    field: [10m_wind_speed]
  - <<: [ *atm2oce, *conserv_interp_stack ]
    coupling_period: ${ocecouplingTimeStep}
    field: [surface_fresh_water_flux,
            total_heat_flux,
            atmosphere_sea_ice_bundle]
  - <<: [ *oce2atm, *conserv_interp_stack ]
    coupling_period: ${ocecouplingTimeStep}
    field: [sea_surface_temperature,
            ocean_sea_ice_bundle]
  - <<: [ *atm2riv, *conserv_interp_dest ]
    coupling_period: ${rivcouplingTimeStep}
    field: [surface_water_runoff,
            soil_water_runoff]
  - <<: [ *riv2oce, *spmap_interp_stack ]
    coupling_period: ${rivcouplingTimeStep}
    field: river_runoff
EOF


#
# yaml files for yac
# ------------------
#
#add_required_file ${basedir}/run/coupling_${EXPNAME}.yaml     ./coupling.yaml
add_required_file                 coupling_${EXPNAME}.yaml     ./coupling.yaml
#
#-----------------------------------------------------------------------------
# II. ATMOSPHERE
#-----------------------------------------------------------------------------
#
#-----------------------------------------------------------------------------
# NWP model parameters
nlev=90           # nlev = number of full levels
iforcing=3        # 3 for inwp forcing; 0 for no forcing
#
#-----------------------------------------------------------------------------

# link files directly
atmo_grid_target=atmo_${atmos_refinement}_DOM01                      #  needed in name of ifs2icon and extpar data
rad_grid_target=rad_${rad_refinement}_DOM00
atmo_grid_folder=/hpc/rhome/routfox/routfox/icon/grids/public/edzw   # grid directory

add_link_file ${atmo_grid_folder}/icon_grid_${rad_gridID}_${rad_refinement}_R.nc             ${rad_grid_target}.nc
add_link_file ${atmo_grid_folder}/icon_grid_${rad_gridID}_${rad_refinement}_R-grfinfo.nc     ${rad_grid_target}-grfinfo.nc
add_link_file ${atmo_grid_folder}/icon_grid_${atmos_gridID}_${atmos_refinement}_G.nc         ${atmo_grid_target}.nc
add_link_file ${atmo_grid_folder}/icon_grid_${atmos_gridID}_${atmos_refinement}_G-grfinfo.nc ${atmo_grid_target}-grfinfo.nc

# NWP Extpar and ifs2icon specificatons

INDIR=/hpc/uwork/mkoehler/run-icon/coupled/proto2
#add_link_file $INDIR/external_parameter_icon_icon_grid_0013_R02B04_G_tiles.nc           extpar_iconR02B04_DOM01.nc
#add_link_file $INDIR/external_parameter_icon_icon_grid_0013_R02B04_G_tiles_oceanLSM.nc  extpar_iconR02B04_DOM01.nc  # old
#add_link_file $INDIR/external_parameter_icon_icon_grid_0036_R02B04_G_tiles_oceanLSM.nc2 extpar_iconR02B04_DOM01.nc
#add_link_file $INDIR/ExtParSLM_grid_${atmos_gridID}.nc     extpar_atmo_${atmos_refinement}_DOM01.nc  # Stephan
#add_link_file $INDIR/new-extpar-withfract30-35.nc          extpar_atmo_${atmos_refinement}_DOM01.nc  # Stephan
extparfile=icon_extpar_oceLSM_a${atmos_gridID}_${atmos_refinement}_o${ocean_gridID}_${ocean_refinement}_20180625_tiles.nc
add_link_file $INDIR/lsm-extpar/${extparfile}               extpar_atmo_${atmos_refinement}_DOM01.nc  # Martin

#add_link_file $INDIR/nn-0013-ifs2icon_R2B04.nc            ifs2icon_R2B04_DOM01.nc
#add_link_file $INDIR/ifs2icon_0013_r0002_20140601.nc      ifs2icon_R2B04_DOM01.nc
#add_link_file $INDIR/ifs2icon_1979010100_0043_R02B04_G.nc ifs2icon_R2B04_DOM01.nc
add_link_file  $INDIR/ifs2icon_2019010100_0030_R02B05.nc   ifs2icon_R2B05_DOM01.nc

# files needed for radiation

add_link_file ${basedir}/data/ECHAM6_CldOptProps.nc .
add_link_file ${basedir}/data/rrtmg_lw.nc .
ecRad_data_path=${basedir}'/externals/ecrad/data'

# Kinne background aerosols for the year 1850 (irad_aero=12, filename without year)

datadir=/hpc/uwork/icon-sml/edzw/${atmos_gridID}/aerosol_kinne/r0001/
add_link_file ${datadir}/bc_aeropt_kinne_lw_b16_coa.nc      bc_aeropt_kinne_lw_b16_coa.nc
add_link_file ${datadir}/bc_aeropt_kinne_sw_b14_coa.nc      bc_aeropt_kinne_sw_b14_coa.nc

# use exclusively Kinne background aerosols for the year 1850, filename without year (irad_aero=12)
add_link_file ${datadir}/bc_aeropt_kinne_sw_b14_fin_1850.nc bc_aeropt_kinne_sw_b14_fin.nc

# HD input files

HD_DIR=/hpc/uwork/icon-sml
add_link_file $HD_DIR/HD-data/input/05deg/masks_05.nc                   masks.nc
add_link_file $HD_DIR/HD-data/input/05deg/hdpara_vs1_11.nc              hdpara.nc
add_link_file $HD_DIR/HD-data/input/05deg/hdstart_05.nc                 hdstart.nc    # HD restart file !!
add_link_file $HD_DIR/edzw/${atmos_gridID}-${ocean_gridID}/hd/hd_receive_${atmos_gridID}_${atmos_refinement}.nc \
                                                                        hd_receive.nc

#-----------------------------------------------------------------------------

# If the model runs with mpi_procs=1, then this script exits without error:
if [[ ${mpi_total_procs} = 1 ]]; then
  echo "This test must not be started if mpi_procs=1 --> exit test without error"
  check_error 0 "exp.nat_jww_nwp_mpiomp"
  exit 0
fi

# If running in parallel, mpi_procs*no_of_nodes must be bigger than 2:
#n_procs=`expr ${mpi_procs_pernode} \* ${no_of_nodes}`
n_procs=${mpi_total_procs}

if [[ ${n_procs} = 2 ]]; then
  echo "This test needs at least 3 MPI processes --> increase mpi_procs from 2 to 3"
  mpi_procs=3
  START=${START%2}3
fi

#-----------------------------------------------------------------------------
#
# write ICON-NWP atmospheric namelist parameters
# ------------------------
# For a complete list see Namelist_overview and Namelist_overview.pdf

cat > ${atmo_namelist} << EOF
!
&coupling_mode_nml
 coupled_to_ocean        = .TRUE.
 coupled_to_hydrodisc    = .TRUE.
/
!
&parallel_nml
 nproma                  = ${nproma_atm}
 p_test_run              = .false.
 l_test_openmp           = .false.
 l_log_checks            = .true.
 num_io_procs            = ${mpi_atm_io_procs}
 num_restart_procs       = 0
 iorder_sendrecv         = 3
 proc0_shift             = 1
 use_omp_input           = .true.
/
&grid_nml
 dynamics_grid_filename  = '${atmo_grid_target}.nc'
 radiation_grid_filename = '${rad_grid_target}.nc' ! reduced radiattion grid
!radiation_grid_filename =  " "                    ! full radiatation grid
 dynamics_parent_grid_id = 0,1
 lredgrid_phys           = .true.
 lfeedback               = .false.
/
&run_nml
 num_lev                 = ${nlev},${nlev}      ! number of full levels of vertical grid
 dtime                   = ${dtime}             ! [s] timestep in seconds
 ldynamics               = .TRUE.               ! dynamics
 ltransport              = .TRUE.               ! transport
 iforcing                = ${iforcing}          !
 ntracer                 =  5                   ! number of tracers - default 0
 ltestcase               = .false.              ! initialize with real data
 ltimer                  = .true.               !
 timers_level            = 20                   ! 10 useful and cheap
 msg_level               = 10                   ! detailed report during integration (7 or 13, > 10 expensive)
 output                  = 'nml'
 check_uuid_gracefully   = .true.
 Restart_filename        = "${EXPNAME}_restart_atm_<rsttime>.nc"
/
&initicon_nml
 lcouple_ocean_coldstart = .TRUE.               ! true: normal coupled forecasts, false: coupled data assimilation
/
&dynamics_nml
 lcoriolis               = .TRUE.
 divavg_cntrwgt          = 0.50
/
&diffusion_nml
 hdiff_order             = 5
 hdiff_efdt_ratio        = 24.0
 hdiff_smag_fac          = 0.025
 lhdiff_vn               = .TRUE.
 lhdiff_temp             = .TRUE.
 hdiff_multfac           = 1.0
 hdiff_tv_ratio          = 1.0
/
&nwp_phy_nml
 inwp_gscp               = 1
 inwp_convection         = 1
 inwp_radiation          = 4
 inwp_cldcover           = 1
 inwp_turb               = 1
 inwp_satad              = 1
 inwp_sso                = 1
 inwp_gwd                = 1
 inwp_surface            = 1
 latm_above_top          = .false.
 efdt_min_raylfric       = 7200.
 itype_z0                = 2
 icapdcycl               = 3
!icpl_aero_conv          = 1
!icpl_aero_gscp          = 1
 icpl_o3_tp              = 1
!lshallow_liquid_only    = .true.
 ithermo_water           = 1          ! 1: T-dependent latent heat in sat-adjust and microphysics
 dt_rad                  = 1800.
 dt_conv                 = 1800.
 dt_sso                  = 1800.
 dt_gwd                  = 1800.
/
&nwp_tuning_nml
 itune_albedo            = 0                         ! somewhat reduced albedo (w.r.t. MODIS data) over Sahara in order to reduce cold bias
 tune_zceff_min          = 0.025                     ! ** default value to be used for R3B7; use 0.025 for R2B6 in order to get similar temperature biases in upper troposphere **
 tune_gkdrag             = 0.075                     ! R2B6: 0.075
 tune_gkwake             = 1.5                       ! R2B6: 1.5
 tune_gfrcrit            = 0.425                     ! R2B6: 0.425
 tune_dust_abs           = 0.
 tune_zvz0i              = 0.85
 tune_box_liq_asy        = 3.25                      ! oper global: 3.0 , oper D2: 3.25, default: 2.5
 tune_box_liq            = 0.05
 tune_rcucov             = 0.075
 tune_rhebc_land         = 0.825
 tune_gust_factor        = 7.0
 tune_eiscrit            = 7.0                       ! estimated inversion strength to turn off shallow
 icpl_turb_clc           = 1
 lcalib_clcov            = .false.                   ! turn off TCC, HCC, MCC, LCC tuning
/
&turbdiff_nml
 tkhmin                  = 0.6
 tkhmin_strat            = 1.0
 tkmmin                  = 0.75
 tkmmin_strat            = 4
 alpha0                  = 0.0123
 alpha0_max              = 0.0335
 alpha1                  = 0.125
 pat_len                 = 750.
 c_diff                  = 0.2
 rlam_heat               = 10.0
 rat_sea                 = 0.8
 ltkesso                 = .true.                    ! SSO dissipation energy used in TKE equation
 frcsmot                 = 0.0                       ! these 2 switches together apply vertical smoothing of the TKE source terms
 imode_frcsmot           = 2                         ! in the tropics (only), which reduces the moist bias in the tropical lower troposphere
 itype_sher              = 3                         ! use horizontal shear production terms with 1/SQRT(Ri) scaling to prevent unwanted side effects
 ltkeshs                 = .true.
 a_hshr                  = 2.0
 icldm_turb              = 1                         ! 2: Gauss clouds for turbulence    1: grid scale clouds
 icldm_tran              = 2                         ! 2: Gauss clouds for surface layer 1: grid scale clouds
/
&transport_nml
 ivadv_tracer            = 3,3,3,3,3
 itype_hlimit            = 3,4,4,4,4,0
 ihadv_tracer            = 32,2,2,2,2,0
/
&interpol_nml
 nudge_zone_width        = 8
 l_intp_c2l              = .true.
 rbf_scale_mode_ll       = 2
!l_mono_c2l              = .true.
!lsq_high_ord            = 3
!support_baryctr_intp    =.true.,
/
&nonhydrostatic_nml
 iadv_rhotheta           = 2
 ivctype                 = 2
 itime_scheme            = 4
 exner_expol             = 0.333
 vwind_offctr            = 0.2                       ! 0.2 for R2B6 and higher resolution, 0.3 for lower resolution
 divdamp_order           = 24                        ! 2 ass, 24 fc
 divdamp_type            = 32                        ! optional: 2 assimilation cycle, 32 forecast
 divdamp_fac             = 0.004                     ! 0.004 for R2B6; recommendation for R3B7: 0.003
 divdamp_trans_start     = 12500
 divdamp_trans_end       = 17500
 igradp_method           = 3
 l_zdiffu_t              = .true.
 thslp_zdiffu            = 0.02
 thhgtd_zdiffu           = 125.
 htop_moist_proc         = 22500.
 hbot_qvsubstep          = 16000.
/
&sleve_nml
 min_lay_thckn           = 20.                       ! lowest level thickness (between half-levels)
 max_lay_thckn           = 400.                      ! maximum layer thickness below htop_thcknlimit
 htop_thcknlimit         = 14000.
 top_height              = 75000.
 stretch_fac             = 0.9
 decay_scale_1           = 4000.
 decay_scale_2           = 2500.
 decay_exp               = 1.2
 flat_height             = 16000.
/
&io_nml
 lflux_avg               = .FALSE.                   ! false: accumulated fluxes
 itype_pres_msl          = 5                         ! (1) 3: IFS-type extrapolation
 itype_rh                = 1                         ! (1) 2: mixed phase (water and ice)
!dt_checkpoint           = 2592000.                  ! 30 days
 inextra_3d              = 2                         ! 3D extra variables
 inextra_2d              = 10                        ! 2D extra variables
!lnetcdf_flt64_output    = .TRUE.                    ! 64 bit output in all files
!lkeep_in_sync           = .TRUE.                    ! sync after each timestep
/
!--------------------------------------------------------------------------------
! lat-lon output
!--------------------------------------------------------------------------------
&output_nml
 output_filename         = '${EXPNAME}_atm_latlon'   ! file name base
 output_start            = "${start_date}"           ! output_start = output_end
 output_end              = "${end_date}"
 output_interval         = "${atm_output_interval}"  ! the output interval and
!operation               = "mean"                    ! only working with remap=0
 file_interval           = "${atm_file_interval}"    ! the file interval
 filetype                = 4                         ! output format: 2=GRIB2, 4=NETCDFv2
 dom                     = -1
 mode                    = 1                         ! 1: forecast mode (relative t-axis)
 include_last            = .TRUE.                    ! flag whether to include the last time step
 filename_format         = '<output_filename>_<levtype>_<datetime2>'
 output_grid             = .TRUE.                    ! flag whether grid information is added to output.
 remap                   = 1                         ! 1: latlon,  0: native grid
 reg_lon_def             = 0.,1.,360.
 reg_lat_def             = -90.,1.,90.
 ml_varlist              = 'clct', 'tqv', 'tqc_dia', 'tqi_dia', 'sp_10m', 't_2m' ,'t_g', 'qv_2m', 'h_ice', 't_ice',
                           'accthb_s','accthb_t','accsob_s','accsob_t','accshfl_s','acclhfl_s','accumfl_s','accvmfl_s',
                           'pres_sfc', 'tot_prec', 't_seasfc', 'fr_seaice', 'frac_t_9',
                           'fr_land', 'fr_lake','frac_t_1', 'accsod_t'
                         ! 'thb_s','thb_t','sob_s','sob_t','shfl_s','lhfl_s',
                         ! 'lsm_switch', 'lsm_ctr_c',
                         ! 'extra_2d1','extra_2d2','extra_2d3','extra_2d4','extra_2d5',
                         ! 'extra_2d6','extra_2d7','extra_2d8','extra_2d9','extra_2d10'
                         ! 't_so', 'extra_3d1',
 pl_varlist              = 'temp', 'u'
 p_levels                = 20000, 50000, 85000
!m_levels                = "5...(nlev+1)"
/
!--------------------------------------------------------------------------------
! native accumulated output
!--------------------------------------------------------------------------------
&output_nml
 output_filename         = '${EXPNAME}_atm_native'   ! file name base
 output_start            = "${start_date}"           ! output_start = output_end
 output_end              = "${end_date}"
 output_interval         = "${atm_output_interval}"  ! the output interval and
!operation               = "mean"                    ! only working with remap=0, not accumulated!
 file_interval           = "${atm_file_interval}"    ! the file interval
 filetype                = 4                         ! output format: 2=GRIB2, 4=NETCDFv2
 dom                     = -1
 mode                    = 1                         ! 1: forecast mode (relative t-axis)
 include_last            = .TRUE.                    ! flag whether to include the last time step
 filename_format         = '<output_filename>_<levtype>_<datetime2>'
 output_grid             = .TRUE.                    ! flag whether grid information is added to output.
 remap                   = 0                         ! 1: latlon,  0: native grid
 ml_varlist              = 'accthb_s','accthb_t','accsob_s','accsob_t','accshfl_s','acclhfl_s',
                           'accumfl_s','accvmfl_s','tot_prec',
                           'lsm_switch', 'lsm_ctr_c', 'accsod_t',
                           'runoff_g', 'runoff_s'
 pl_varlist              = 'temp', 'u'
 p_levels                = 20000, 50000, 85000
/
!--------------------------------------------------------------------------------
! native time-mean output
!--------------------------------------------------------------------------------
&output_nml
 output_filename         = '${EXPNAME}_atm_nat_mn'   ! file name base
 output_start            = "${start_date}"           ! output_start = output_end
 output_end              = "${end_date}"
 output_interval         = "${atm_output_interval}"  ! the output interval and
 file_interval           = "${atm_file_interval}"    ! the file interval
 operation               = "mean"                    ! only working with remap=0
 filetype                = 4                         ! output format: 2=GRIB2, 4=NETCDFv2
 dom                     = -1
 mode                    = 1                         ! 1: forecast mode (relative t-axis)
 include_last            = .TRUE.                    ! flag whether to include the last time step
 filename_format         = '<output_filename>_<levtype>_<datetime2>'
 output_grid             = .TRUE.                    ! flag whether grid information is added to output.
 remap                   = 0                         ! 1: latlon,  0: native grid
 ml_varlist              = 'clct', 'tqv', 'tqc_dia', 'tqi_dia', 'sp_10m', 't_2m' ,'t_g', 'qv_2m', 'h_ice', 't_ice',
                           'pres_sfc', 't_seasfc', 'fr_seaice', 'frac_t_9'
/
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-NWP Output monitoring: global means
!--------------------------------------------------------------------------------
&output_nml
  output_filename        = '${EXPNAME}_atm_mon'
  filetype               =  4                        ! output format: 2=GRIB2, 4=NETCDFv2
  mode                   =  1                        ! 1: forecast mode, 2: climate mode
  output_start           = "${start_date}"           ! output_start = output_end
  output_end             = "${end_date}"
  output_interval        = "${atm_output_interval}"  ! the output interval and
  file_interval          = "${atm_file_interval}"    ! the file interval
  include_last           = .FALSE.
  operation              = 'mean'
  filename_format        = '<output_filename>_<datetime2>'
  ml_varlist             = 'tas_gmean','rsdt_gmean','rsut_gmean','rlut_gmean','radtop_gmean',
                           'prec_gmean','evap_gmean','pme_gmean'
  output_grid            =  .FALSE.
  remap                  =  0
/
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-NWP Output 2D
!--------------------------------------------------------------------------------
&output_nml
  output_filename        = '${EXPNAME}_atm_2d_ml'
  filetype               =  4                        ! output format: 2=GRIB2, 4=NETCDFv2
  mode                   =  2                        ! 1: forecast mode, 2: climate mode
  dom                    =  1                        ! write all domains
  output_start           = "${start_date}"           ! output_start = output_end
  output_end             = "${end_date}"
  output_interval        = "${atm_output_interval}"  ! the output interval and
  file_interval          = "${atm_file_interval}"    ! the file interval
  include_last           = .FALSE.
  operation              = 'mean'
  filename_format        = '<output_filename>_<datetime2>'
  ml_varlist             = 'pres_msl','sod_t','sou_t','sob_t','thb_t','t_s','fr_land',
                           'clct','shfl_s','lhfl_s','qhfl_s','tot_prec_rate',
                           'tqv_dia','tqc_dia','tqi_dia','umfl_s','vmfl_s','sp_10m','t_2m',
                           't_seasfc', 'fr_seaice'
  output_grid            =  .TRUE.
  remap                  =  0
/
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-NWP Output 3D
!--------------------------------------------------------------------------------
&output_nml
  output_filename        = '${EXPNAME}_atm_3d_ml'
  filetype               =  4      ! output format: 2=GRIB2, 4=NETCDFv2
  mode                   =  2      ! 1: forecast mode, 2: climate mode
  dom                    =  1      ! write all domains
  output_start           = "${start_date}"           ! output_start = output_end
  output_end             = "${end_date}"
  output_interval        = "${atm_output_interval}"  ! the output interval and
  file_interval          = "${atm_file_interval}"    ! the file interval
  include_last           = .FALSE.
  operation              = 'mean'
  filename_format        = '<output_filename>_<datetime2>'
  ml_varlist             = 'pres','geopot','temp','u','v','qv','rh','clc','tot_qc_dia','tot_qi_dia'
  output_grid            =  .TRUE.
  remap                  =  0
/
!&dbg_index_nml
!  idbg_mxmn             = 4                         ! initialize MIN/MAX  debug output
!  idbg_val              = 0                         ! initialize one cell debug output
!  idbg_slev             = 1                         ! initialize start level for debug output
!/
&extpar_nml
 itopo                   = 1
 n_iter_smooth_topo      = 1
 heightdiff_threshold    = 3000.
/
&lnd_nml
 ntiles                  = 3
 nlev_snow               = 3
 lmulti_snow             = .false.
 itype_heatcond          = 3
 idiag_snowfrac          = 20
 lsnowtile               = .true.
 lseaice                 = .true.
 llake                   = .true.
 itype_lndtbl            = 4
 itype_evsl              = 4
 itype_trvg              = 3
 itype_root              = 2
 cwimax_ml               = 5.e-4
 c_soil                  = 1.25
 c_soil_urb              = 0.5
 sstice_mode             = 1                         ! 1 for coupled
 lprog_albsi             = .true.
 hice_min                = 0.05                      ! default of nwp sea-ice model
 hice_max                = 10.0   !L64=4.8, L72=10.0 ! limit sea-ice thickness, consistency with seaice_limit!
 itype_snowevap          = 2
/
&meteogram_output_nml
  lmeteogram_enabled     = .false., .false.
  n0_mtgrm               = 0                          ! meteogram initial time step (0 is first step!)
  ninc_mtgrm             = 2, 4 ! meteogram output interval (in terms of time steps)
  ldistributed           = .false., .false.
  stationlist_tot        =  50.0,       8.6,    'Frankfurt-Flughafen',
                            52.260,    14.313,  'Lindenberg_Obs',
                            52.1873,   13.9894, 'Waldstation_Kehrigk',
                            52.1215,   14.1077, 'Falkenberg',
/
EOF


if [[ "$pi_control" == "yes" ]]; then
   cat >> ${atmo_namelist} << EOF
&radiation_nml
 irad_o3                 = 79
 irad_aero               = 12                   ! 6: default, 12: Kinne aerosols (Trang) 1850
 izenith                 = 4                    ! 4: NWP default, 3: no annual cycle
 isolrad                 = 0                ! 1: Coddington Scaling correction (default NWP)
 albedo_type             = 2                    ! Modis albedo
                         !  values for 1850 CE  ! values for 2012
 vmr_co2                 = 284.3e-06            ! 390.e-06
 vmr_ch4                 = 808.2e-09            ! 1800.e-09
 vmr_n2o                 = 273.0e-09            ! 322.0e-09
 vmr_cfc11               = 0.0	                ! 240.e-12
 vmr_cfc12               = 0.0	                ! 532.e-12
 vmr_o2                  = 0.20946
 direct_albedo           = 4
 direct_albedo_water     = 3
 albedo_whitecap         = 1
 ecRad_llw_cloud_scat    = .true.
 ecRad_data_path         = '${ecRad_data_path}'
 /
EOF

else
   cat >> ${atmo_namelist} << EOF
&radiation_nml
 irad_o3                 = 79
 irad_aero               = 6                ! 6: default, 12: Kinne aerosols (Trang)
 izenith                 = 4                ! 4: NWP default, 3: no annual cycle
 isolrad                 = 1                ! 1: Coddington Scaling correction (default NWP)
 albedo_type             = 2                ! Modis albedo
                         ! values for 2012    values for 1850 CE
 vmr_co2                 = 390.e-06         ! 284.3e-06
 vmr_ch4                 = 1800.e-09        ! 808.2e-09
 vmr_n2o                 = 322.0e-09        ! 273.0e-09
 vmr_o2                  = 0.20946
 vmr_cfc11               = 240.e-12         ! 0.0
 vmr_cfc12               = 532.e-12         ! 0.0
 direct_albedo           = 4
 direct_albedo_water     = 3
 albedo_whitecap         = 1
 ecRad_llw_cloud_scat    = .true.
 ecRad_data_path         = '${ecRad_data_path}'
/
EOF

fi


#-----------------------------------------------------------------------------
# add standard atmo_non-hydrostatic_files
. ${rundir}/add_required_atmo_non-hydrostatic_files


#-----------------------------------------------------------------------------
# III. OCEAN and SEA-ICE
#-----------------------------------------------------------------------------

ocean_grid_source=icon_grid_${ocean_gridID}_${ocean_refinement}_O.nc
ocean_grid_target=ocean_${ocean_refinement}_DOM01.nc
#ocean_grid_folder="/work/mh0727/m300732/input_SBH/${ocean_gridID}"
ocean_grid_folder=/hpc/uwork/mkoehler/run-icon/coupled/proto2
#ocean_data_InputFolder="/pool/data/ICON/oes/input/r0004/${ocean_grid_name}"
ocean_data_InputFolder=/hpc/uwork/mkoehler/run-icon/coupled/proto2

# --------------
# ocean namelist
# --------------

#-----------------------------------------------------------------------------
#set the vertical coordinate  : 0=zlev, 1=zstar

VERT_COR=1
minVerticalLevels=2
if [ ${VERT_COR} = 0 ] ; then
  vert_cor_type=0
  select_lhs=2
  l_lhs_direct=.true.
  STRETCH_C=''
else
  vert_cor_type=1
  select_lhs=1
  l_lhs_direct=.false.
  STRETCH_C=",'stretch_c'"     #  output of variable thickness
fi

#-----------------------------------------------------------------------------
#set the vertical number of levels

if [ ${LEV} == L40 ] ; then

    n_zlev=40
    dzlev_m="12.,10.,10.,10.,10.,10.,13.,15.,20.,25.,30.,35.,40.,45.,50.,55.,60.,70.,80.,90.,
         100.,110.,120.,130.,140.,150.,170.,180.,190.,200.,220.,250.,270.,300.,350.,400.,
         450.,500.,500.,600."

    levidx_100m=9
    levidx_200m=12
    levidx_2000m=30

    minVerticalLevels=2

elif [ ${LEV} == L64 ] ; then

    n_zlev=64
    dzlev_m="12.,10.,10.,10.,10.,10.,10.,10.,10.,10.,11.,12.,13.,14.,15.,16.,17.,18.,20.,22.,
         24.,26.,28.,30.,32.,35.,38.,41.,45.,49.,53.,58.,62.,66.,71.,75.,80.,85.,91.,97.,
         104.,111.,118.,125.,132.,138.,145.,152.,160.,167.,175.,182.,188.,195.,201.,208.,
         213.,219.,224.,230.,235.,241.,250.,260."

    levidx_100m=10
    levidx_200m=17
    levidx_2000m=46

    minVerticalLevels=2

elif [ ${LEV} == L128 ] ; then

    n_zlev=128
    dzlev_m="11.0,   9.0,     8.0,   8.0,     8.0,    8.0,    8.0,    8.0,    8.0,    8.0,\
              8.0,    8.0,    8.0,    8.25,   8.5,    8.75,   9.0,   9.25,    9.5,   9.75,\
             10.0,   10.0,   10.0,   10.0,   10.0,   10.0,   10.0,   10.0,   10.0,   10.0,\
             10.5,   11.0,   11.5,   12.0,   12.5,   13.0,   13.5,   14.0,   14.5,   15.0,\
             15.5,   16.0,   16.5,   17.0,   17.5,   18.0,   18.5,   19.0,   19.5,   20.0,\
             20.5,   21.0,   21.5,   22.0,   22.5,   23.0,   23.5,   24.0,   24.5,   25.0,\
             25.5,   26.0,   26.5,   27.0,   28.5,   29.0,   29.5,   30.0,   30.5,   31.0,\
             31.0,   32.0,   33.0,   34.0,   35.0,   36.0,   37.0,   38.0,   39.0,   40.0,\
             42.0,   44.0,   46.0,   48.0,   50.0,   52.0,   54.0,   56.0,   58.0,   60.0,\
             62.0,   64.0,   66.0,   68.0,   70.0,   72.0,   74.0,   76.0,   78.0,   80.0,\
             82.0,   84.0,   86.0,   88.0,   90.0,   92.0,   94.0,   96.0,   98.0,  100.0,\
            102.0,  104.0,  106.0,  108.0,  110.0,  112.0,  114.0,  116.0,  118.0,  200.0,\
            200.0,  200.0,  200.0,  200.0,  200.0,  200.0,  200.0,  200.0"

    levidx_100m=14
    levidx_200m=24
    levidx_2000m=92

    minVerticalLevels=2

elif [ ${LEV} == L128SMT ] ; then

    n_zlev=128
    dzlev_m="2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 3.0, 3.1, 3.2,\
                 3.4, 3.5, 3.7, 3.9, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.3, 5.5,\
                 5.8, 6.0, 6.3, 6.6, 6.9, 7.2, 7.5, 7.8, 8.2, 8.5, 8.9, 9.3,\
                 9.8, 10.2, 10.7, 11.1, 11.5, 11.9, 12.3, 12.7, 13.1, 13.5,\
                14.0, 14.5, 14.9, 15.4, 15.9, 16.5, 17.0, 17.6, 18.2, 18.8,\
                19.4, 20.0, 20.7, 21.4, 22.1, 22.8, 23.6, 24.4, 25.2, 26.0,\
                26.9, 27.8, 28.7, 29.7, 30.6, 31.7, 32.7, 33.8, 34.9, 36.1,\
                37.3, 38.5, 39.8, 41.1, 42.5, 43.9, 45.3, 46.8, 48.4, 50.0,\
                51.7, 53.4, 55.2, 57.0, 58.9, 60.8, 62.9, 64.9, 67.1, 69.3,\
                71.6, 74.0, 76.5, 79.0, 81.6, 84.3, 87.1, 90.0, 93.0, 96.1,\
                99.3, 102.6, 106.0, 109.5, 113.2, 116.9, 120.8, 124.8, 128.9,\
               133.2, 137.6, 142.2, 146.9, 151.8, 156.9, 162.1, 167.4, 173.0,\
               178.7, 184.7, 190.8, 197.1"

    levidx_100m=27
    levidx_200m=37
    levidx_2000m=96

    minVerticalLevels=12

elif [ ${LEV} == L72 ] ; then

    n_zlev=72
!    dzlev_m="2.0,   2.2,   2.5,   2.8,   3.1,   3.5,   3.9,   4.4,   4.9,   5.4,\
!             5.9,   6.4,   7.1,   7.7,   8.4,   9.2,  10.1,  11.0,  12.0,  13.2,\
!            14.4,  15.7,  17.1,  18.7,  20.4,  22.3,  24.3,  26.5,  28.9,  31.5,\
!            34.3,  37.3,  40.6,  43.1,  45.3,  46.8,  48.4,  50.0,  51.7,  53.4,\
!            55.2,  57.0,  58.9,  60.8,  62.9,  66.6,  72.6,  80.6,  90.6, 100.2,\
!           110.0, 120.3, 128.7, 137.4, 146.4, 155.7, 165.2, 174.8, 184.4, 194.1,\
!           203.6, 212.9, 221.9, 230.5, 238.5, 245.9, 252.4, 258.1, 262.8, 266.4,\
!           268.9, 270.1"

    dzlev_m="2.0, 2.2, 2.5, 2.8, 3.1, 3.5, 3.9, 4.4, 4.9, 5.4,\
            5.9, 6.4, 7.1, 7.7, 8.4, 9.2, 10.1, 11.0, 12.0, 13.2,\
           14.4, 15.7, 17.1, 18.7, 20.4, 22.3, 24.3, 26.5, 28.9, 31.5,\
           34.3, 37.3, 40.6, 43.1, 45.3, 46.8, 48.4, 50.0, 51.7, 53.4,\
           55.2, 57.0, 58.9, 60.8, 62.9, 66.6, 72.6, 80.6, 90.6, 100.2,\
          110.0, 120.3, 128.7, 137.4, 146.4, 155.7, 165.2, 174.8, 184.4, 194.1,\
          203.6, 212.9, 221.9, 230.5, 238.5, 245.9, 252.4, 258.1, 262.8, 266.4,\
          268.9, 270.1"

    levidx_100m=19
    levidx_200m=25
    levidx_2000m=55
    minVerticalLevels=10

fi

#-----------------------------------------------------------------------------
# HAMOCC
#

use_hamocc=no
ocean_vertical_levels=$n_zlev

if [  "x${use_hamocc}"  = "xyes" ]; then
lhamocc=".TRUE."
lbgcadv=".TRUE."
nlev_eu=$ocean_vertical_levels
# set nlev_eu to level belonging to approx 500m
if [ "x$ocean_vertical_levels" = "x40" ];then
nlev_eu=19
fi
if [ "x$ocean_vertical_levels" = "x64" ];then
nlev_eu=26
fi
if [ "x$ocean_vertical_levels" = "x20" ];then
nlev_eu=16
fi
else
lhamocc=".FALSE."
lbgcadv=".FALSE."
fi
l_init_bgc=".false." # Hamocc state variables are in restart
if [[ "$initialiseOcean" == "fromRestart" ]] || [[ "$initialiseOcean" == "fromClimatology" ]]; then
  l_init_bgc=".true." # initialize HAMOCC (HAMOCC state variables are not in restart)
fi
# ----------------------------------------------------------------------------
#
# ocean namelist
# --------------

cat > ${ocean_namelist} << EOF
!
&coupling_mode_nml
 coupled_to_atmo             = .TRUE.
 coupled_to_hydrodisc        = .TRUE.
/
&parallel_nml
 nproma                      = ${nproma_oce}
 num_io_procs                = ${mpi_oce_io_procs}
 num_prefetch_proc           = 0
 pio_type                    = 0
 p_test_run                  = .FALSE.
 l_fast_sum                  = .TRUE.
/
&grid_nml
 dynamics_grid_filename      = "${ocean_grid_target}"
 use_dummy_cell_closure      = .TRUE.
 use_duplicated_connectivity = .FALSE.
/
&run_nml
 modelTimeStep               = "${oceTimeStep}"
 output                      = 'nml'                            ! namelist controlled output scheme
 activate_sync_timers        = .TRUE.
 profiling_output            = 1                                ! aggregated: 1; detailed: 2; in files: 3
 msg_timestamp               = .FALSE.
 timers_level                = 10
 debug_check_level           = 1
 Restart_filename            = "${EXPNAME}_restart_oce_<rsttime>.nc"
/
EOF

# output of fixed geometry
cat >> ${ocean_namelist} << EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O fixed geometry
!--------------------------------------------------------------------------------
&output_nml
  output_filename  = "${EXPNAME}_oce_fx"
  filetype         = 5                                ! output format: 2=GRIB2, 4=NETCDFv2, 5=NETCDFv4
  mode             = 1                                ! 1: forecast mode, 2: climate mode
  filename_format  = "<output_filename>_<datetime2>"
  file_interval    = "${oce_file_interval}"           ! interval in ISO-format
  ml_varlist       = 'lsm_ctr_c','lsm_c','lsm_e','surface_cell_sea_land_mask','surface_edge_sea_land_mask',
                     'surface_vertex_sea_land_mask','vertex_bottomLevel','basin_c','regio_c','bottom_thick_c',
                     'bottom_thick_e','column_thick_c','column_thick_e','wet_c','wet_e','wet_halo_zero_c',
                     'wet_halo_zero_e','prism_thick_c','invConstantPrismThickness','prism_volume','prism_thick_e',
                     'prism_thick_flat_sfc_c','prism_thick_flat_sfc_e','inverse prism_thick_c','prism_center_dist_c',
                     'constantPrismCenters_Zdistance','constantPrismCenters_invZdistance','inv_prism_thick_e',
                     'inv_prism_center_dist_c','inv_prism_center_dist_e','depth_CellMiddle'
  output_grid      = .TRUE.
  include_last     = .false.
  output_bounds    =  0., 0., 3600.                   ! start, end, increment
/
EOF

if [[ "$output_oce_def" == "yes" ]]; then
   cat >> ${ocean_namelist} << EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O 3D Output
!--------------------------------------------------------------------------------
&output_nml
  output_filename            = "${EXPNAME}_oce_def"
  filetype                   = 5                                ! output format: 2=GRIB2, 4=NETCDFv2
  mode                       = 1                                ! 1: forecast mode, 2: climate mode
  filename_format            = "<output_filename>_<datetime2>"
  output_start               = "${start_date}"                  ! start in ISO-format
  output_end                 = "${end_date}"                    ! end in ISO-format
  output_interval            = "${oce_output_interval_def}"     ! interval in ISO-format
  file_interval              = "${oce_file_interval}"
  include_last               = .FALSE.
  output_grid                = .FALSE.
  operation                  = "mean"
  ml_varlist                 = 'group:oce_default', 'group:oce_essentials'
/
EOF
fi

if [[ "$output_oce_ice" == "yes" ]]; then
  cat >> ${ocean_namelist} <<EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O ice and snow output
!--------------------------------------------------------------------------------
&output_nml
  output_filename  = "${EXPNAME}_oce_ice"
  filetype         = 4                                ! output format: 2=GRIB2, 4=NETCDFv2, 5=NETCDFv4
  mode             = 1                                ! 1: forecast mode, 2: climate mode
  output_start     = "${start_date}"                  ! start date in ISO-format
  output_end       = "${end_date}"                    ! end date in ISO-format
  output_interval  = "${oce_output_interval}"         ! interval in ISO-format
  file_interval    = "${oce_file_interval}"           ! interval in ISO-format
  operation        = 'mean'                           ! mean over output interval
  output_grid      = .TRUE.
  filename_format  = "<output_filename>_<datetime2>"
  m_levels         = "1...3,9,23,30,34"               ! surface and some subsurface levels only
  ml_varlist       =  'draftave','hi','hs','conc','to','so'
/
EOF
fi

if [[ "$output_oce_mon" == "yes" ]]; then
   cat >> ${ocean_namelist} << EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O Output monitoring
!--------------------------------------------------------------------------------
&output_nml
  output_filename            = "${EXPNAME}_oce_mon"
  filetype                   = 4                                ! output format: 2=GRIB2, 4=NETCDFv2
  mode                       = 1                                ! 1: forecast mode, 2: climate mode
  filename_format            = "<output_filename>_<datetime2>"
  output_start               = "${start_date}"                  ! start in ISO-format
  output_end                 = "${end_date}"                    ! end in ISO-format
  output_interval            = "${oce_output_interval_mon}"
  file_interval              = "${oce_file_interval}"
  include_last               = .FALSE.
  output_grid                = .FALSE.
  operation                  = "mean"
  ml_varlist                 = 'group:ocean_monitor'
/
EOF
fi

if [[ "$output_oce_moc" == "yes" ]]; then
   cat >> ${ocean_namelist} << EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O AMOC
!--------------------------------------------------------------------------------
&output_nml
  output_filename            = "${EXPNAME}_oce_moc"
  filetype                   = 5                                ! output format: 2=GRIB2, 4=NETCDFv2
  mode                       = 1                                ! 1: forecast mode, 2: climate mode
  filename_format            = "<output_filename>_<datetime2>"
  output_start               = "${start_date}"                  ! start in ISO-format
  output_end                 = "${end_date}"                    ! end in ISO-format
  output_interval            = "${oce_output_interval}"
  file_interval              = "${oce_file_interval}"
  include_last               = .FALSE.
  output_grid                = .FALSE.
  operation                  = "mean"
  ml_varlist                 = 'group:ocean_moc'
/
EOF
fi

if [[ "$output_oce_dbg" == "yes" ]]; then
  cat >> ${ocean_namelist} <<EOF
!--------------------------------------------------------------------------------
! OUTPUT pyICON: ICON-O Output debug
!--------------------------------------------------------------------------------
&output_nml
  output_filename  = "${EXPNAME}_oce_dbg"
  filetype         = 5                                ! output format: 2=GRIB2, 4=NETCDFv2, 5=NETCDFv4
  mode             = 1                                ! 1: forecast mode (relative t-axis); 2: climate mode
  filename_format  = "<output_filename>_<datetime2>"
  output_start     = "${start_date}"                  ! start date in ISO-format
  output_end       = "${end_date}"                    ! end date in ISO-format
  output_interval  = "${oce_output_interval}"         ! interval in ISO-format
  file_interval    = "${oce_file_interval}"           ! interval in ISO-format
  output_grid      = .TRUE.
  operation        = 'mean'                           ! mean over output interval
  include_last     = .false.
  m_levels         = '1,2'                            ! surface and subsurface level only
  ml_varlist       =  'Qtop', 'Qbot','Qbot_slow', 'zHeatOceI','heatOceI','heatOceW','draftave','zUnderIce',
                      'mld','condep','zos','u_vint','hi','hs','conc','to','so','u',
                      'HeatFlux_Total','atmos_fluxes_HeatFlux_ShortWave','atmos_fluxes_HeatFlux_LongWave',
                      'HeatFlux_ShortWave','HeatFlux_LongWave','HeatFlux_Sensible','HeatFlux_Latent',
                      'FrshFlux_Runoff','FrshFlux_Precipitation','FrshFlux_Evaporation','FrshFlux_SnowFall',
                      'FrshFlux_TotalOcean','FrshFlux_VolumeIce','totalsnowfall','Wind_Speed_10m'
/
EOF
fi


cat >> ${ocean_namelist} << EOF
&dbg_index_nml
  idbg_mxmn                  = 4                                ! initialize MIN/MAX  debug output
  idbg_val                   = 0                                ! initialize one cell debug output
  idbg_slev                  = 1                                ! initialize start level for debug output
  idbg_elev                  = 3                                ! 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                ='oceanCouplng'                    ! define modules to print out in debug mode
 !str_mod_tst                = 'all'                            ! define modules to print out in debug mode
/
&ocean_dynamics_nml
! zlev or zstar:
  vert_cor_type = ${vert_cor_type}
  minVerticalLevels = ${minVerticalLevels}

! Number of vertical levels
  n_zlev                =  $n_zlev
  dzlev_m(1:${n_zlev})  =  $dzlev_m

  l_edge_based                               = .FALSE.     ! edge- or cell-based mimetic discretization
! l_partial_cells                            = .FALSE.     ! partial bottom cells=TRUE: local varying bottom depth
  select_solver                              = 4           ! 1=gmres_oce_old; 2=ocean_restart_gmres, 3=mixed precisison restart
                                                           ! 4=CG (default) 5=CGJ 6=BiCG 7=GMRES restart (legacy) 8=MINRES
  use_absolute_solver_tolerance              = .TRUE.
  solver_tolerance                           = 1.0E-9      ! 1e-10 in omip-YVF - enlarged to 1e-9 (LL opt)
  select_lhs                             = ${select_lhs}   ! 1=operator based (default) 2=matrix based - 2=error?
  l_lhs_direct                           = ${l_lhs_direct} ! .true.= use lhs implementation directly  .false.= matrix scanner (default)
  solver_FirstGuess                          = 2           ! 0=start from zeros 1=last timestep smoothed 2=last timestep (default)

! solver_max_iter_per_restart                = 14
! solver_max_restart_iterations              = 100         ! outer (restart solver)

  fast_performance_level                     = 200         ! performance level 12: for cell-based; 5: default
  use_continuity_correction                  = .TRUE.      ! height adjustment according to vertical velocity in dynamics
  cfl_check                                  = .FALSE.
  cfl_write                                  = .FALSE.

  i_bc_veloc_top                             = 1
  i_bc_veloc_bot                             = 1           ! 0: (def) bottom friction off, 1: on
/
&ocean_tracer_transport_nml
  flux_calculation_horz                      = 5           ! 1=upwind, 2=central, 3=Lax-Friedrichs,
                                                           ! 4=Miura, 5=FCT with Zalesak limiter (default)
  flux_calculation_vert                      = 7           ! 6=adpo; 7=upwind biased ppm (default); 8=FCT with zalesak limiter
  ! define low and high order methods to be used in
  ! horizontal flux corrected transport methods
  ! (flux_calculation_horz=4,5)
  fct_low_order_flux                         = 1           ! horizontal low  order method: 1=upwind (def), no other implemented
  fct_high_order_flux                        = 5           ! horizontal high order method: 1=upwind, 2=central, 3=lax_friedrichs, 4=miura_order1
  fct_limiter_horz                           = 100         ! zalesak
/
&ocean_horizontal_diffusion_nml
  laplacian_form                             = 1           ! 1=curlcurl-graddiv
  VelocityDiffusion_order                    = 2           ! 21=biharmonic+laplacian (for the laplacian leith)

  BiharmonicViscosity_scaling                = 1
  BiharmonicViscosity_reference              = 3.5E12      ! [m2/s] biharmonic horizontal viscosity coefficient for velocity
  BiharmonicViscosity_background             = 0.0

  HarmonicViscosity_scaling                  = 1
  HarmonicViscosity_reference                = 0.0         ! [m2/s] constant horizontal viscosity coefficient for velocity
  HarmonicViscosity_background               = 0.0

  TracerHorizontalDiffusion_scaling          = 1
  Temperature_HorizontalDiffusion_Background = 0.0
  Temperature_HorizontalDiffusion_Reference  = 0
  Salinity_HorizontalDiffusion_Background    = 0.0
  Salinity_HorizontalDiffusion_Reference     = 0
/
&ocean_vertical_diffusion_nml
  vert_mix_type                              = 2           ! 1: PP; 2: TKE
  PPscheme_type                              = 0           ! 4: current PPscheme - 0: switched off
! vert_mix_type                              = 1           ! 1: PP; 2: TKE
! PPscheme_type                              = 4           ! 4: current PPscheme - 0: switched off

!  Parameters for PP-type vertical mixing
  !velocity_VerticalDiffusion_background      = 5.0E-5      ! [m2/s]  vertical background viscosity coefficient for velocity
  !temperature_VerticalDiffusion_background   = 5.0E-5      ! [m2/s]  vertical background diffusion coefficient for temperature
  !salinity_VerticalDiffusion_background      = 5.0E-5      ! [m2/s]  vertical background diffusion coefficient for salinity
  !tracer_convection_MixingCoefficient        = 0.1         ! max vertical tracer diffusion for convection used in case of instability
  !convection_InstabilityThreshold            = -1.0E-6     ! used in update_ho_params - default=-5e-8
  !richardsonDiffusion_threshold              = 0.0         ! used in update_ho_params - default=+5e-8
  !tracer_RichardsonCoeff                     = 2.0E-3      ! factor for vertical diffusion coefficient in PP scheme
  !velocity_RichardsonCoeff                   = 2.0E-3      ! factor for vertical viscosity coefficient in PP scheme
  !bottom_drag_coeff                          = 3.0E-3      ! default=2.5E-3; active for i_bc_veloc_bot=1
  !use_wind_mixing                            = .TRUE.      ! TRUE: use wind mixing scheme in MPIOM-type pp-scheme
  !lambda_wind                                = 0.03
  !tracer_TopWindMixing                       = 1.0E-5 ! 0.5E-3
  !! velocity_TopWindMixing                     = 1.0E-5 ! 0.5E-3
  !! tracer_TopWindMixing                       = 2.31480E-6  ! eqals 0.5e-3 m2s-1 (MPIOM)
  !velocity_TopWindMixing                     = 2.31480E-6

!  Parameters for TKE-type vertical mixing
  ! cvmix/tke parameters
  c_k        = 0.3
  c_eps      = 0.7
  alpha_tke  = 30.0
  mxl_min    = 1.d-8
  kappaM_min = 0.0
  kappaM_max = 100.0
  cd         = 3.75
  tke_min    = 1.d-6
  tke_mxl_choice = 2
  tke_surf_min = 1.d-4
  only_tke = .true.
  l_lc         = .false.     ! Use Langmuir parameterisation (Axell, 2002)
  clc          = 0.15        ! Factor in the vertical velocity profile of the Langmuir circulation (default clc=0.1 in Axell (2002))
  use_ubound_dirichlet = .false.
  use_lbound_dirichlet = .false.
/
&ocean_GentMcWilliamsRedi_nml
  GMRedi_configuration                       = 1           ! 0=cartesian diffusion; 1=GM-Redi: bolus advection + isopycnal diffusion
  tapering_scheme                            = 1
  GMRedi_usesRelativeMaxSlopes               = .FALSE.
  S_max                                      = 3.0e-5      ! 1.0
  S_d                                        = 1.0e-4      ! 5e-3 to 5e-4
                                                           !
  k_tracer_GM_kappa_parameter                = 400.0       !
  k_tracer_isoneutral_parameter              = 400.0       ! value for cell-based cartesian diffusion - mpiom: 1000/400km = 400/160km
  k_tracer_dianeutral_parameter              = 0.0         ! 1.0E-5  !
                                                           !
  switch_off_diagonal_vert_expl              = .TRUE.
  gmredi_combined_diagnostic                 = .FALSE.
! switch_on_redi_balance_diagnostic          = .FALSE.     ! not yet available in icon-aes-dyamond++
  revert_vertical_recon_and_transposed       = .TRUE.
  slope_calc_via_temperture_salinity         = .TRUE.
  include_slope_squared_implicit             = .TRUE.      ! think of l_with_vert_tracer_diffusion
  switch_on_tapering_horizontal_diffusion    = .TRUE.
/
&ocean_physics_nml
  i_sea_ice                                  = 1           ! 0 = no sea ice; 1 = sea ice model on; default=1
  lhamocc                                    = ${lhamocc}
  lbgcadv                                    = ${lbgcadv}
  eos_type                                   = 2           ! 1=linear EOS,2=(nonlinear, from MPIOM)
/
&sea_ice_nml
  i_ice_therm                                = 1           ! 1=zero-layer (default), 2=Winton, 0/2: not allowed
  i_ice_dyn                                  = 1           ! 1/0=switch on/off AWI ice dynamics
! i_Qio_type                                 = 3           ! 3 (default): energy of whole grid-area used for melting (MPIOM-type)
! use_constant_tfreez                        = .TRUE.      ! default: TRUE
! use_no_flux_gradients                      = .FALSE.     ! default: TRUE
! hmin                                       = 0.05        ! default: 0.05 - minimum sea-ice value for ice
  leadclose_1                                = 0.25        ! default: 0.5 - value of MPIOM: 0.25
  leadclose_2n                               = 0.666       ! default: 0.0 - value of MPIOM: 2/3
 !albedoW_sim                                = 0.10        ! albedo of the ocean used in sea ice model
 !leadclose_1                                = 0.5         ! default: 0.5 - value of MPIOM: 0.25
 !leadclose_2n                               = 0.666       ! default: 0.0 - value of MPIOM: 2/3
 !pstar                                      = 40000.      ! default: 27500.; MPIOM=20000.
/
EOF

cat >> ${ocean_namelist} << EOF
&ocean_forcing_nml
  iforc_oce                                  = 14          ! ocean forcing: 14 from coupling via YAC
  type_surfRelax_Temp                        = -1          ! -1: use net surface heat flux from atmosphere
                                                           !  0: no relaxation used
                                                           !  1: relaxation switched on for reading (init_oce_relax=1)
                                                           !     or some testcases only
  forcing_enable_freshwater                  = .TRUE.      ! enable/disable freshwater flux
  forcing_windstress_u_type                  = 2           ! 0: zero wind stress, 1: read from file, 2: none
  forcing_windstress_v_type                  = 2           ! 0: zero wind stress, 1: read from file, 2: none
                                                           ! salt-change due to internal fluxes only
  atm_pressure_included_in_ocedyn            = .TRUE.      ! use atm. pressure in the ocean dynamics
  atm_pressure_included_in_icedyn            = .TRUE.      ! use atm. pressure in the ice dynamics
  limit_seaice                               = .TRUE.      ! default: TRUE
  seaice_limit                               = 5.0    !0.4 ! limit sea-ice thickness [m] = seaice_limit * ocean top layer thickness
                                                           ! attention: this value [m] needs to be consistent with hice_max
  lfix_salt_content                          = .FALSE.
  limit_elevation                            = .FALSE.      ! true: adjust daily to zero, default=false: free surface
                                                           ! in absense of a river routing .TRUE. keeps ocen surface constant
! lswr_jerlov                                = .FALSE.     ! SW-penetration (Jerlov-radiation, default=T)
/
EOF

if [[ "$initialiseOcean" == "fromRestart" ]] || [[ "$initialiseOcean" == "fromRestartwithHamocc" ]]; then

  cat >> ${ocean_namelist} << EOF
&ocean_initialConditions_nml
  initial_salinity_type                      = 0                ! 0: none, 1: read S from initial_state.nc
  initial_temperature_type                   = 0                ! 0: none, 1: read T from initial_state.nc
  initialize_fromRestart                     = .TRUE.
/
EOF

else

  if [[ "$initialiseOcean" == "fromClimatology" ]]; then
    cat >> ${ocean_namelist} << EOF
&ocean_initialConditions_nml
  initial_salinity_type                      = 1                ! 0: none, 1: read S from initial_state.nc
  initial_temperature_type                   = 1                ! 0: none, 1: read T from initial_state.nc
  initialize_fromRestart                     = .FALSE.
/
EOF
  else
    if [[ "$initialiseOcean" == "coldStart" ]]; then
      cat >> ${ocean_namelist} << EOF
&ocean_initialConditions_nml
  initial_salinity_type                      = 200              ! 0: none, 1: read S from initial_state.nc
  initial_temperature_type                   = 200              ! 0: none, 1: read T from initial_state.nc
  initial_temperature_top                    = 15.0
  initial_temperature_bottom                 = 15.0
  initial_salinity_top                       = 35.0
  initial_salinity_bottom                    = 35.0
  initialize_fromRestart                     = .false.
/
EOF
    else
      cat >> ${ocean_namelist} << EOF
&ocean_initialConditions_nml
  initial_salinity_type                      = 0                ! 0: none, 1: read S from initial_state.nc
  initial_temperature_type                   = 0                ! 0: none, 1: read T from initial_state.nc
  initialize_fromRestart                     = .FALSE.
/
EOF
    fi
  fi
fi

cat >> ${ocean_namelist} << EOF
&ocean_diagnostics_nml
  diagnostics_level                          = 1
  diagnose_for_horizontalVelocity            = .FALSE.
  diagnose_for_heat_content                  = .TRUE.
/
&io_nml
  restart_file_type                          = 5
  write_last_restart                         = .TRUE.
! restart_write_mode                         = "joint procs multifile"  ! not yet available in ocean model
  lnetcdf_flt64_output                       = .FALSE.          ! T: 64 bit output in all files
  lkeep_in_sync                              = .TRUE.           ! sync after each timestep
/
EOF

add_required_file ${basedir}/run/${ocean_namelist} ./

#-----------------------------------------------------------------------------
#
# Ocean initialisation input
#
if [[ "$initialiseOcean" == "fromRestart" ]]; then
 # spun-up for 900 years
 #add_link_file /work/mh0727/m300732/input/0036/ocean/restart/ler0656_restart_oce_30000101T000000Z.nc    restart_ocean_DOM01.nc
# add_link_file /work/mh0727/m300732/icon-oes-hamocc/experiments/ler0863/ler0863_restart_oce_22000101T000000Z.nc init-restart-oce.nc
# add_link_file ${ocean_data_InputFolder}/ler0863_restart_oce_22000101T000000Z.nc init-restart-oce.nc
  add_link_file ${ocean_data_InputFolder}/${restart_exp}_restart_oce_${restart_yea}.nc restart_ocean_DOM01.nc
fi

if [[ "$initialiseOcean" == "fromRestartwithHamocc" ]]; then
  # Restart file like above + 1 year with Hamocc (initialized from basin averages)
#  add_link_file /work/mh0727/m300732/input/0036/ocean/restart/ler0614_restart_oce_10010101T000000Z.nc init-restart-oce.nc
  add_link_file ${ocean_data_InputFolder}/ler0614_restart_oce_10010101T000000Z.nc init-restart-oce.nc
  add_link_file init-restart-oce.nc restart_ocean_DOM01.nc
fi

if [[ "$initialiseOcean" == "fromClimatology" ]]; then
  # for initializing from Levitus:
  #add_link_file ${ocean_data_InputFolder}/R2B4L40_initial_state.nc                                   initial_state.nc
  #add_link_file ${ocean_data_InputFolder}/ts_phc3.0_annual_icon_icon_grid_0036_R02B04_G_L40.nc       initial_state.nc
  #add_link_file ${ocean_data_InputFolder}/ts_phc3.0_annual_icon_grid_${ocean_gridID}_R02B04_G_L40.nc initial_state.nc
  #add_link_file ${ocean_data_InputFolder}/slo1004-initial_state.nc                                   initial_state.nc
  #add_link_file ${ocean_data_InputFolder}/R2B6${LEV}_initial_state-HH.nc                             initial_state.nc
  add_link_file ${ocean_data_InputFolder}/temp/HelmuthHaak/ts_oras5_icon_OceanOnly_Global_IcosSymmetric_0039km_rotatedZ37d_BlackSea_Greenland_modified_sills_srtm30_1min_L72.nc_2010-01-01   initial_state.nc
fi

# Ocean grid
#
add_link_file ${ocean_grid_folder}/${ocean_grid_source} ./${ocean_grid_target}

#
# HAMOCC forcing files
#
if [  "x${use_hamocc}"  = "xyes" ]; then
  # FIXME: this is still a private path
  datadir=/work/mh0727/m300732/input/0036/ocean/hamocc/
  add_link_file $datadir/MAHOWALDDUST_icon_grid_0036_R02B04_O_remapbil1.nc           dust.nc  # iron deposition
  add_link_file $datadir/ndepo_1-0_gr_185001-185012-clim_icon_grid_0036_R02B04_O.nc  nitrogen.nc  # nitrate deposition
fi
#
#-----------------------------------------------------------------------------

if [ $mpi_total_procs -lt `expr $mpi_oce_procs + 1` ] ; then
   echo "Too few mpi_total_procs for requested mpi_oce_procs."
   echo "-> check mpi_total_procs and mpi_oce_procs. Exiting."
   check_error 0
   exit
fi

#-----------------------------------------------------------------------------
#
# directories definition
#
#RUNSCRIPTDIR=${basedir}/run
if [ x$grids_folder = x ] ; then
   HGRIDDIR=${basedir}/grids
else
   HGRIDDIR=$grids_folder
fi

#-----------------------------------------------------------------------------
final_status_file=${RUNSCRIPTDIR}/${job_name}.final_status
rm -f ${final_status_file}

#-----------------------------------------------------------------------------
# set up the model lists if they do not exist
# this works for single model runs
# for coupled runs the lists should be declared explicilty
if [ x$namelist_list = x ]; then
#  minrank_list=(        0           )
#  maxrank_list=(     65535          )
#  incrank_list=(        1           )
  minrank_list[0]=0
  maxrank_list[0]=65535
  incrank_list[0]=1
  if [ x$atmo_namelist != x ]; then
    # this is the atmo model
    namelist_list[0]="$atmo_namelist"
    modelname_list[0]="atmo"
    modeltype_list[0]=1
    run_atmo="true"
  elif [ x$ocean_namelist != x ]; then
    # this is the ocean model
    namelist_list[0]="$ocean_namelist"
    modelname_list[0]="oce"
    modeltype_list[0]=2
  elif [ x$psrad_namelist != x ]; then
    # this is the psrad model
    namelist_list[0]="$psrad_namelist"
    modelname_list[0]="psrad"
    modeltype_list[0]=3
  elif [ x$hamocc_namelist != x ]; then
    # this is the hamocc model
    namelist_list[0]="$hamocc_namelist"
    modelname_list[0]="hamocc"
    modeltype_list[0]=4
  elif [ x$testbed_namelist != x ]; then
    # this is the testbed model
    namelist_list[0]="$testbed_namelist"
    modelname_list[0]="testbed"
    modeltype_list[0]=99
  else
    check_error 1 "No namelist is defined"
  fi
fi


#-----------------------------------------------------------------------------
# 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}
else
  if [[ "$initialiseOcean" == "FALSE" ]]; then
#   add_link_file ${restart_dir}/${restart_exp}/${restart_exp}_restart_atm_${restart_yea}0101T000000Z.nc init-restart-atm.nc
#   add_link_file ${restart_dir}/${restart_exp}/${restart_exp}_restart_oce_${restart_yea}0101T000000Z.nc init-restart-oce.nc

#   add_link_file ${ocean_data_InputFolder}/slo1414_restart_atm_21000101T000000Z.nc  init-restart-atm.nc
#   add_link_file ${ocean_data_InputFolder}/slo1414_restart_oce_21000101T000000Z.nc  init-restart-oce.nc

#   add_link_file init-restart-atm.nc restart_atm_DOM01.nc
#   add_link_file init-restart-oce.nc restart_ocean_DOM01.nc

    add_link_file ${ocean_data_InputFolder}/slo1414_restart_atm_21000101T000000Z.nc restart_atm_DOM01.nc
    add_link_file ${ocean_data_InputFolder}/slo1414_restart_oce_21000101T000000Z.nc restart_ocean_DOM01.nc
  fi
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
# fill some checks

run_atmo=${run_atmo="false"}
if [ x$atmo_namelist != x ]; then
  run_atmo="true"
fi
run_ocean=${run_ocean="false"}
if [ x$ocean_namelist != x ]; then
  run_ocean="true"
fi
run_psrad=${run_psrad="false"}
if [ x$psrad_namelist != x ]; then
  run_psrad="true"
fi
run_hamocc=${run_hamocc="false"}
if [ x$hamocc_namelist != x ]; then
  run_hamocc="true"
fi


#-----------------------------------------------------------------------------
# add grids to required files
all_grids="${atmo_dyn_grids} ${atmo_rad_grids} ${ocean_grids}"
for gridfile in ${all_grids}; do
  #
  gridfile=${gridfile//\'/} # strip all ' in case ' is used to delimit the grid names
  gridfile=${gridfile//\"/} # strip all " in case " is used to delimit the grid names
  gridfile=${gridfile//\,/} # strip all , in case , is used to separate the grid names
  #
  grfinfofile=${gridfile%.nc}-grfinfo.nc
  #
  ls -l ${HGRIDDIR}/$gridfile
  check_error $? "${HGRIDDIR}/$gridfile does not exist."
  add_link_file ${HGRIDDIR}/${gridfile} ./
  if [ -f ${HGRIDDIR}/${grfinfofile} ]; then
    add_link_file ${HGRIDDIR}/${grfinfofile} ./
  fi
done


#-----------------------------------------------------------------------------
# print_required_files
copy_required_files
link_required_files


#-----------------------------------------------------------------------------
# get restart files

if  [ x$restart_atmo_from != "x" ] ; then
  rm -f restart_atm_DOM01.nc
#  ln -s ${basedir}/experiments/${restart_from_folder}/${restart_atmo_from} ${EXPDIR}/restart_atm_DOM01.nc
  cp ${basedir}/experiments/${restart_from_folder}/${restart_atmo_from} cp_restart_atm.nc
  ln -s cp_restart_atm.nc restart_atm_DOM01.nc
  restart=".true."
fi
if  [ x$restart_ocean_from != "x" ] ; then
  rm -f restart_oce.nc
#  ln -s ${basedir}/experiments/${restart_from_folder}/${restart_ocean_from} ${EXPDIR}/restart_oce.nc
  cp ${basedir}/experiments/${restart_from_folder}/${restart_ocean_from} cp_restart_oce_DOM01.nc
  ln -s cp_restart_oce_DOM01.nc restart_oce_DOM01.nc
  restart=".true."
fi
#-----------------------------------------------------------------------------


read_restart_namelists=${read_restart_namelists:=".true."}

#-----------------------------------------------------------------------------
#
# create ICON master namelist
# ------------------------
# For a complete list see Namelist_overview and Namelist_overview.pdf

#-----------------------------------------------------------------------------
# create master_namelist
master_namelist=icon_master.namelist
if [ x$end_date = x ]; then
cat > $master_namelist << EOF
&master_nml
 lrestart            = $restart
 lrestart_write_last = .true.
/
&master_time_control_nml
 experimentStartDate  = "$start_date"
 restartTimeIntval    = "$restart_interval"
 checkpointTimeIntval = "$checkpoint_interval"
/
&time_nml
 is_relative_time = .false.
/
EOF
else
if [ x$calendar = x ]; then
  calendar='proleptic gregorian'
  calendar_type=1
else
  calendar=$calendar
  calendar_type=$calendar_type
fi
cat > $master_namelist << EOF
&master_nml
 lrestart            = $restart
 lrestart_write_last = .true.
 read_restart_namelists = $read_restart_namelists
/
&master_time_control_nml
 calendar             = "$calendar"
 checkpointTimeIntval = "$checkpoint_interval"
 restartTimeIntval    = "$restart_interval"
 experimentStartDate  = "$start_date"
 experimentStopDate   = "$end_date"
/
&time_nml
 is_relative_time = .false.
/
EOF
fi
#-----------------------------------------------------------------------------


#-----------------------------------------------------------------------------
# add model component to master_namelist
add_component_to_master_namelist()
{

  model_namelist_filename="$1"
  model_name=$2
  model_type=$3
  model_min_rank=$4
  model_max_rank=$5
  model_inc_rank=$6

cat >> $master_namelist << EOF
&master_model_nml
  model_name="$model_name"
  model_namelist_filename="$model_namelist_filename"
  model_type=$model_type
  model_min_rank=$model_min_rank
  model_max_rank=$model_max_rank
  model_inc_rank=$model_inc_rank
/
EOF

}
#-----------------------------------------------------------------------------


no_of_models=${#namelist_list[*]}
echo "no_of_models=$no_of_models"

j=0
while [ $j -lt ${no_of_models} ]
do
  add_component_to_master_namelist "${namelist_list[$j]}" "${modelname_list[$j]}" ${modeltype_list[$j]} ${minrank_list[$j]} ${maxrank_list[$j]} ${incrank_list[$j]}
  j=`expr ${j} + 1`
done


#-----------------------------------------------------------------------------
#
# HD Namelist
#
YYYYMMDD_start=${start_date:0:4}${start_date:5:2}${start_date:8:2}
YYYYMMDD_end=${end_date:0:4}${end_date:5:2}${end_date:8:2}
#
## Log output for gridboxes if nhd_diag=100
xlon1=8.5
xlat1=54.5
xlon2=10.5
xlat2=53.5
## Type of Interpolation from input (atmospheric) grid to HD grid
##          0   Input = Output
##          1   using HDMAP routine by Veronika (default)
##          2   0.5 degree to 5 Min./
IMAP=0
## Coupling type: 0=no, 1=no interpolation, 2=interpolation in HD
ICOUPLE=1
## User specific settings that will be put as attributes in the output file
## via namelist HDUSER_CTL in file namelist.hduser
HD_USER="buildbot"
HD_INST="Deutscher Wetterdienst, Germany"
HD_INSTID=""
EMAIL_ADDRESS=""

if (( $ICOUPLE != 2 )) ; then
   lcoup_out=.FALSE.
else
   lcoup_out=.TRUE.
fi

cat > namelist.hd << end_hd_ctl
&HD_CTL
  date_start    = ${YYYYMMDD_start}
  date_end      = ${YYYYMMDD_end}
  delta_time    = ${dtime_riv}      !! 3600 ! 7200 ! 14400 ! 21600 ! 28800 ! 43200 ! 86400
  runoff_file   = ""
  drainage_file = ""
  forcing_freq  = 0                 ! 0: stepwise, 1: daily
  IOUT          = 6                 ! averaging of HD output, 5: monthly, 6: daily
  UFAKRU        = 0.001             ! ICON-NWP: kg m-2 s-1 = mm/s, HD: m/s
  coupling_type = ${ICOUPLE}        !
  coupling_file = ""
  lcoupling_out = ${lcoup_out}      ! Write discharge on ocean grid dep. on IOUT (couling_type 2 only)
  iform_input   = 1                 ! Format Input files: 0 = SRV (Default), 1 = NetCDCF
  lcoupling_atm = .TRUE.
  lcoupling_oce = .TRUE.
  out_expname   = ${EXPNAME}
  out_datapath  = "${EXPDIR}/"
  icpl_mask_tohd= 2                 ! 2: YAC, use hd_receive_mask from file hd_receive.hd
  icpl_sinks    = 1                 ! 1: redistribute water in sinks to mouth boxes
/
end_hd_ctl

cat > namelist.hdset << end_hydrology_ctl
&HYDROLOGY_CTL
  ldebughd          = .TRUE.  !F
  diag_water_budget = .TRUE.  !F    ! additional prints to diagnose global water budget
  locean            = .FALSE.       ! closure of water budget for ocean coupling
  nhd_diag          = 7             ! 7: Elbe outflow diagnostics at 5min resolution
  lhd_highres       = .TRUE.  !F    ! .TRUE.: CALL hd_highres_write --> /output/hd/hd_YYYY_MM_02_hd_highres.nc
  fllog1            = ${xlon1}
  fblog1            = ${xlat1}
  fllog2            = ${xlon2}
  fblog2            = ${xlat2}
  nremap            = ${IMAP}
  lhd_rout          = .FALSE.
! Factors to allow sensitivity studies on 5 Min.
  fk_rfk            = 1.
  fk_lfk            = 1.
  fk_gfk            = 1.
! Discharge dependent riverflow velocity (irf_vel <> 0)
  irf_vel           = 1
  qrf_ref           = 1000.
/
end_hydrology_ctl

cat > namelist.hduser << end_hduser_ctl
&HDUSER_CTL
  hd_user   = "${HD_USER}"
  hd_cont   = "${EMAIL_ADDRESS}"
  hd_inst   = "${HD_INST}"
  hd_instid = "${HD_INSTID}"
/
end_hduser_ctl

# ----------------------------------------------------------------------
# run the model!
# ----------------------------------------------------------------------

rm -f finish_atmo.status
rm -f finish_ocean.status
date

export VE_ERRCTL_ALLOCATE=MSG
export NMPI_PROGINF=YES
export VE_TRACEBACK=VERBOSE
export NMPI_SEPSELECT=3
export GMON_OUT_PREFIX=scal_prof
export VE_FPE_ENABLE=DIV,FOF,INV
export GFORTRAN_UNBUFFERED_PRECONNECTED=y
export NMPI_EXPORT="GFORTRAN_UNBUFFERED_PRECONNECTED"

# to get better debug output -- VERY slow (Jens-Olaf Beismann)
#export VE_ADVANCEOFF=YES

# NE=4, CPE=8, cpunum_job=3
# geometry: ID of MPI processes (see ranks below, oce_min_rank, etc)
# MPI process ID   node   VH/VE    info
# 0                0      VH       ICON initialization process
# 1-16             0-3    VE       ICON atmosphere
# 17-18            0      VH       IO   atmosphere
# 19-34            0-3    VE       ICON ocean

# normal output
#mpirun -v -vh     -node 0        -np 1                    -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
#          -venode -node 0-${NE1} -np ${mpi_atm_procs}     -env OMP_NUM_THREADS 1 ${MODEL_VE} : \
#          -vh     -node 0        -np ${mpi_atm_io_procs}  -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
#          -venode -node 0-${NE1} -np ${mpi_oce_procs}     -env OMP_NUM_THREADS 1 ${MODEL_VE}

# debug output: each processor with separate output in EXPDIR
# time /opt/nec/ve/mpi/3.5.0/bin/runtime/mpirun \
#   -v -vh  -node 0        -np 1                   -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH} : \
#   -venode -node 0-${NE2} -np ${mpi_atm_procs}    -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} : \
#   -vh     -node 0        -np ${mpi_atm_io_procs} -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH} : \
#   -venode -node 0-${NE2} -np ${mpi_oce_procs}    -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} : \
#   -venode -node ${NE1}   -np ${mpi_hd_procs}     -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_HD}

# normal output
# time /opt/nec/ve/mpi/3.5.0/bin/runtime/mpirun \
#   -v -vh  -node 0        -np 1                   -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
#   -venode -node 0-${NE2} -np ${mpi_atm_procs}    -env OMP_NUM_THREADS 1 ${MODEL_VE} : \
#   -vh     -node 0        -np ${mpi_atm_io_procs} -env OMP_NUM_THREADS 1 ${MODEL_VH} : \
#   -venode -node 0-${NE2} -np ${mpi_oce_procs}    -env OMP_NUM_THREADS 1 ${MODEL_VE} : \
#   -venode -node ${NE1}   -np ${mpi_hd_procs}     -env OMP_NUM_THREADS 1 ${MODEL_HD}


### Roland's processor optimisation:

atm_ves=$(( (mpi_atm_procs + CPE - 1)/CPE ))     # Number of VEs needed (rounded up)
split_ve=$(( atm_ves - 1 ))                      # Last VE might be split between Atmo and ocean
split_num_oce=$(( atm_ves*CPE - mpi_atm_procs )) # Remaining empty cores on last VE
excl_ves=$((split_ve+1))-${NE2}
excl_num_oce=$(( mpi_oce_procs - split_num_oce ))
(( split_num_oce > 0 )) && is_split=1 || is_split=''

# debug output: each processor with separate output in EXPDIR
/opt/nec/ve/mpi/3.5.0/bin/runtime/mpirun \
    -v -vh  -node 0                         -np 1                   -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH} :  \
    -venode -node 0-${split_ve} -nnp ${CPE} -np ${mpi_atm_procs}    -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} :  \
    -vh     -node 0                         -np ${mpi_atm_io_procs} -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VH} :  \
    ${is_split:+-venode -node ${split_ve}   -np ${split_num_oce}    -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} :} \
    -venode -node ${excl_ves}               -np ${excl_num_oce}     -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_VE} :  \
    -venode -node ${NE1}                    -np ${mpi_hd_procs}     -env OMP_NUM_THREADS 1 /opt/nec/ve/bin/mpisep.sh ${MODEL_HD}

date

if [ -r finish_atmo.status ] && [ -r finish_ocean.status ] ; then
  check_final_status 0
else
  check_final_status -1
fi
#
#-----------------------------------------------------------------------------
#
#Read the finish_status files and cut off whitespaces.
"${finish_status_atm}"=`cat finish_atmo.status`
finish_status_atm="${finish_status_atm# }"
"${finish_status_oce}"=`cat finish_ocean.status`
finish_status_oce="${finish_status_oce# }"

echo "${finish_status_atm}"
echo "${finish_status_oce}"
echo "============================"
echo "Script run (atmo) successfully: $finish_status_atm"
echo "Script run (oce) successfully: $finish_status_oce"
echo "============================"

#-----------------------------------------------------------------------------
# rm output_schedule_steps*
#-----------------------------------------------------------------------------
# if [[ "x$use_hamocc" = "xyes" ]]; then
#   # store HAMOCC log file
#   strg="$(ls -rt ${EXPNAME}_hamocc_EU*.nc* | tail -1 )"
#   prefx="${EXPNAME}_hamocc_EU_tendencies"
#   foo=${strg##${prefx}}
#   foo=${foo%%.*}
#   bgcout_file="bgcout_${foo}"
#   mv bgcout $bgcout_file
# fi

#-----------------------------------------------------------------------------
namelist_list=""
#-----------------------------------------------------------------------------
# check if we have to restart, ie resubmit
#   Note: this is a different mechanism from checking the restart
if [ "${finish_status_atm}" = "RESTART" ] || [ "${finish_status_oce}" = "RESTART" ]; then
  echo "restart next experiment..."
  this_script="${RUNSCRIPTDIR}/${job_name}"
  echo 'this_script: ' "$this_script"
  touch ${restartSemaphoreFilename}
  cd ${RUNSCRIPTDIR}
  ${submit} ${this_script}          # $run_param_0
else
  [[ -f ${restartSemaphoreFilename} ]] && rm ${restartSemaphoreFilename}
fi
#
#-----------------------------------------------------------------------------
# automatic call/submission of post processing if available
# if [ "x${autoPostProcessing}" = "xtrue" ]; then
#   # check if there is a postprocessing is available
#   cd ${RUNSCRIPTDIR}
#   targetPostProcessingScript="./post.${EXPNAME}.run"
#   [[ -x $targetPostProcessingScript ]] && ${submit} ${targetPostProcessingScript}
#   cd -
# fi

#-----------------------------------------------------------------------------

# cd $RUNSCRIPTDIR

#-----------------------------------------------------------------------------
