!parameter inputs via namelist convention.
!(1) Use ' ' (single quotes) for chars;
!(2) integer values are fine for real vars/arrays;
!(3) if multiple entries for a parameter are found, the last one wins - please avoid this
!(4) array inputs follow column major (like FORTRAN) and can spill to multiple lines
!(5) space allowed before/after '='

&CORE
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Core (mandatory) parameters; no defaults
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Pre-processing option. Useful for checking grid violations.
  ipre = 0 !Pre-processor flag (1: on; 0: off)

! Baroclinic/barotropic option. If ibc=0 (baroclinic model), ibtp is not used.
  ibc = 0 !Baroclinic option
  ibtp = 0 

  rnday = 5 !total run time in days
  dt = 150. !Time step in sec

! Grid for WWM (USE_WWM)
  msc2 = 24     !same as msc in .nml ... for consitency check between SCHISM and WWM
  mdc2 = 30     !same as mdc in .nml

! Define # of tracers in tracer modules (if enabled)
  ntracer_gen = 2 !user defined module (USE_GEN)
  ntracer_age = 4 !age calculation (USE_AGE). Must be =2*N where N is # of age tracers
  sed_class = 5 !SED3D (USE_SED)
  eco_class = 27 !EcoSim (USE_ECO): must be between [25,60]

! Global output controls
  nspool = 24 !output step spool
  ihfskip = 288 ! zy; Fei:576 !stack spool; every ihfskip steps will be put into 1_*, 2_*, etc...
/

&OPT
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Optional parameters. The values shown below are default unless otherwise noted
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!-----------------------------------------------------------------------
! 2nd pre-proc flag. If ipre2/=0, code will output some diagnotic outputs and stop.
! These outputs are: drag coefficients (Cdp)
!-----------------------------------------------------------------------
  ipre2 = 0

! Starting time
  start_year = 2023 !int
  start_month = 6 !int
  start_day = 1 !int
  start_hour = 12 !double
  utc_start = 0 !double

!-----------------------------------------------------------------------
! Coordinate option: 1: Cartesian; 2: lon/lat (hgrid.gr3=hgrid.ll in this case,
! and orientation of element is outward of earth)
!-----------------------------------------------------------------------
  ics = 2 !Coordinate option

!-----------------------------------------------------------------------
! Hotstart option. 0: cold start; 1: hotstart with time reset to 0; 2: 
! continue from the step in hotstart.nc
!-----------------------------------------------------------------------
  ihot = 1  ! Fei: 2

!-----------------------------------------------------------------------
! Equation of State type used
! ieos_type=0: UNESCO 1980 (nonlinear); =1: linear function of T ONLY, i.e. 
! \rho=eos_b+eos_a*T, where eos_a<=0 in kg/m^3/C
!-----------------------------------------------------------------------
  ieos_type = 0
  ieos_pres = 0 !used only if ieos_type=0. 0: without pressure effects 
  eos_a = -0.1 !needed if ieos_type=1; should be <=0 
  eos_b = 1001. !needed if ieos_type=1

!-----------------------------------------------------------------------
! Main ramp option
!-----------------------------------------------------------------------
! zy: Fei commented out this line: nramp = 0 !ramp-up option (1: on; 0: off)
! nramp = 0  ! tested, MUST have!; zy rm'ed leading "!"; ramp-up option (1: on; 0: off)
  dramp = 0. ! zy- Fei original: 1. !needed if nramp=1; ramp-up period in days

! zy: Fei commented out this line: nrampbc = 0 !ramp-up flag for baroclinic force
! nrampbc = 0  ! added back "!"; zy rm'ed leading "!"; ramp-up flag for baroclinic force
  drampbc = 0. ! zy; Fei's original: 4. !not used if nrampbc=0

!-----------------------------------------------------------------------
! Method for momentum advection. 0: ELM; 1: upwind (not quite working yet)
!-----------------------------------------------------------------------
  iupwind_mom = 0

!-----------------------------------------------------------------------
! Methods for computing velocity at nodes. 
! If indvel=0, conformal linear shape function is used; if indvel=1, averaging method is used.
! For indvel=0, a stabilization method is needed (see below). 
!-----------------------------------------------------------------------
  indvel = 0 !default: 1
 
!-----------------------------------------------------------------------
! 2 stabilization methods, mostly for indvel=0.
! (1) Horizontal viscosity option. ihorcon=0: no viscosity is used; =1: Lapacian;
! =2: bi-harmonic. If ihorcon=1, horizontal viscosity _coefficient_ (<=1/8, related
! to diffusion number) is given in hvis_coef0, and the diffusion # 
! is problem dependent; [0.001-1/8] seems to work well.
! If ihorcon=2, diffusion number is given by hvis_coef0 (<=0.025).
! If indvel=1, no horizontal viscosity is needed. 
! (2) Shapiro filter (see below)
!
! For non-eddying regime applications (nearshore, estuary, river), two easiest options are: 
! (1) indvel=1, ishapiro=ihorcon=0, and any choices of inter_mom; 
! (2) indvel=0, ishapiro=1 (shapiro0=0.5), ihorcon=inter_mom=0.
! For applications that include eddying regime, refer to the manual.
!-----------------------------------------------------------------------
  ihorcon = 2
  hvis_coef0 = 0.025 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1
!  cdh = 0.01 !needed only if ihorcon/=0; land friction coefficient - not active yet

!-----------------------------------------------------------------------
! 2nd stabilization method via Shapiro filter. This should normally be used 
! if indvel=ihorcon=0. To transition between eddying/non-eddying regimes, use
! indvel=0, ihorcon/=0, and ishapiro=-1 (shapiro.gr3).
!-----------------------------------------------------------------------
  ishapiro = 2 !on/off flag. Default: 0
  shapiro0 = 1.e2 !Shapiro filter strength, needed only if ishapiro=1; max is 0.5 
  niter_shap = 1 !needed of ishapiro/=0 - # of iterations with Shapiro filter. Suggested: 1

!-----------------------------------------------------------------------
! Implicitness factor (0.5<thetai<=1).
!-----------------------------------------------------------------------
  thetai = 1.0

!-----------------------------------------------------------------------
! If WWM is used, set coupling/decoupling flag. Not used if USE_WWM is distabled in Makefile
!       0: decoupled so 2 models will run independently; 
!       1: full coupled (elevation, vel, and wind are all passed to WWM); 
!       2: elevation and currents in wwm, no wave force in SCHISM;  
!       3: no elevation and no currents in wwm, wave force in SCHISM;
!       4: elevation but no currents in wwm, wave force in SCHISM;
!       5: elevation but no currents in wwm, no wave force in SCHISM;
!       6: no elevation but currents in wwm, wave force in SCHISM;
!       7: no elevation but currents in wwm, no wave force in SCHISM;
! Note that all these parameters must be present in this file (even though not used).
!-----------------------------------------------------------------------
  icou_elfe_wwm = 0 
  nstep_wwm = 1  !call WWM every this many time steps. If /=1, consider using quasi-steady mode in WWM
  iwbl = 0 !wave boundary layer formulation (used only if USE_WMM and 
           !icou_elfe_wwm/=0 and nchi=1. If icou_elfe_wwm=0, set iwbl=0): 
           !1-modified Grant-Madsen formulation; 2-Soulsby (1997)
  hmin_radstress = 1. !min. total water depth used only in radiation stress calculation [m]
!  nrampwafo = 0       !ramp-up option for the wave forces (1: on; 0: off)
  drampwafo = 1.      !needed if nrampwafo=1; ramp-up period in days
  turbinj = 0.15      !% of depth-induced wave breaking energy injected in turbulence (default: 0.15 (15%), as proposed by Feddersen, 2012)
  fwvor_advxy_stokes = 1 !  --> Stokes drift advection (xy), Coriolis
  fwvor_advz_stokes = 1  !  --> Stokes drift advection (z) , Coriolis
  fwvor_gradpress = 1    !  --> Pressure term
  fwvor_breaking = 1     !  --> Wave breaking
  cur_wwm = 0            ! Coupling current in WWM
                         ! 0: surface layer current
                         ! 1: depth-averaged current
                         ! 2: computed according to Kirby and Chen (1989)
  wafo_obcramp = 0       ! Ramp on wave forces at open boundary (1: on / 0: off)
                         !  --> this option requires the input file wafo_ramp.gr3
                         !      which defines the ramp value (between 0 and 1)
                         !      at nodes over the whole domain

!-----------------------------------------------------------------------
! Bed deformation option (0: off; 1: vertical deformation only; 2: 3D bed deformation). 
! If imm=1, bdef.gr3 is needed; if imm=2, user needs to update depth info etc
! in the code (not working for ics=2 yet).
!-----------------------------------------------------------------------
  imm = 0
  ibdef = 10 !needed if imm=1; # of steps used in deformation

!-----------------------------------------------------------------------
! Reference latitude for beta-plane approximation when ncor=1 (not used if ics=2)
!-----------------------------------------------------------------------
  slam0 = -77.07  !lon - not really used
  sfea0 = 24 !lat

!-----------------------------------------------------------------------
! Option to deal with under resolution near steep slopes in deeper depths
! 0: use h[12,_bcc below; /=0: use hw_* below
!-----------------------------------------------------------------------
  iunder_deep = 0 

!-----------------------------------------------------------------------
! Baroclinicity calculation in off/nearshore with iunder_deep=ibc=0. 
! The 'below-bottom' gradient
! is zeroed out if h>=h2_bcc (i.e. like Z) or uses const extrap
! (i.e. like terrain-following) if h<=h1_bcc(<h2_bcc) (and linear
! transition in between based on local depth)
!-----------------------------------------------------------------------
  h1_bcc = 50. ![m]
  h2_bcc = 100. ![m]; >h1_bcc

!-----------------------------------------------------------------------
! Hannah-Wright-like ratio & depth used to account for under-resolution
! in a b-clinic model. Used only if ibc=0 and iunder_deep/=0.
! The b-clinic force at prism centers is calculated with a reconstruction
! method in horizontal that has a stencil of an element and its adjacent elements.
! If the depths change is too much between the elem and its adjacent elem
! the under-resolution occurs (with steep bottom slope) and b-clinic force
! needs to be zeroed out below the (higher) bottom, specifically, if
! max(2 depths)>=hw_depth and abs(diff(2 depths))>=hw_ratio*max(2 depths).
!-----------------------------------------------------------------------
  hw_depth = 1.e6 !threshold depth in [m]
  hw_ratio = 0.5 !ratio

!-----------------------------------------------------------------------
! Hydraulic model option. If ihydraulics/=0, hydraulics.in 
! is required. This option cannot be used with non-hydrostatic model.
!-----------------------------------------------------------------------
  ihydraulics = 0

!-----------------------------------------------------------------------
! Point sources/sinks option (0: no; 1: on). If =1, needs source_sink.in (list of elements),
! vsource,th, vsink.th, and msource.th. Source/sinks can be specified at an elem more
! than once, and the code will accumulate the volumes, but for mass conc, one of
! the specified values will be taken (no summation).
!-----------------------------------------------------------------------
  if_source = 1
! nramp_ss = 1 !needed if if_source=1; ramp-up flag for source/sinks
  dramp_ss = 0  ! zy; Fei original: 3 !needed if if_source=1; ramp-up period in days
  meth_sink = 1 !options to treat sinks @ dry elem

!-----------------------------------------------------------------------
! Horizontal diffusivity option. if ihdif=1, horizontal diffusivity is given in hdif.gr3
!-----------------------------------------------------------------------
  ihdif = 0 

!-----------------------------------------------------------------------
! Bottom friction. 
!           nchi=0: drag coefficients specified in drag.gr3; nchi=-1: Manning's 
!           formulation (even for 3D prisms) with n specified in manning.gr3. 
!           nchi=1: bottom roughness (in meters) specified in rough.gr3 (and in this case, negative
!           or 0 depths in rough.gr3 indicate time-independent Cd, not roughness!).
!           Cd is calculated using the log law, when dzb>=dzb_min; when dzb<dzb_min,
!           Cd=Cdmax*exp[dzb_decay*(1-dzb/dzb_min)], where Cdmax=Cd(dzb=dzb_min),
!           and dzb_decay (<=0) is a decay const specified below. We recommend dzb_decay=0
!           and may remove this constant in the future.
!           If iwbl/=0, nchi must =1.
!-----------------------------------------------------------------------
  nchi = 0
  dzb_min = 0.5 !needed if nchi=1; min. bottom boundary layer thickness [m].
  !dzb_decay = 0. !needed if nchi=1; a decay const. [-]. should =0
  hmin_man = 1. !needed if nchi=-1: min. depth in Manning's formulation [m]

!-----------------------------------------------------------------------
! Coriolis. If ncor=-1, specify "rlatitude" (in degrees); if ncor=0,
! specify Coriolis parameter in "coricoef"; if ncor=1, model uses
! lat/lon in hgrid.ll for beta-plane approximation if ics=1, and in this case,
! the latitude specified in CPP projection ('sfea0') is used. If ncor=1 and ics=2,
! Coriolis is calculated from local latitude, and 'sfea0' is not used.
!-----------------------------------------------------------------------
  ncor = 1 !should usually be 1 if ics=2
  rlatitude = 46 !if ncor=-1
  coricoef = 0 !1.e-4 !if ncor=0

!-----------------------------------------------------------------------
! Elevation initial condition flag for cold start only. If ic_elev=1, elev.ic (in *.gr3 format) is needed
! to specify the initial elevations; otherwise elevation is initialized to 0 everywhere 
!-----------------------------------------------------------------------
  ic_elev = 0

!-----------------------------------------------------------------------
! Elevation boundary condition ramp-up flag. =0: ramp up from 0; =1: ramp up from
! elev. values read in from elev.ic or hotstart.nc - if neither is present, from 0.
! This flag is mainly used to start the simulation from non-zero elev.
! The ramp-up period is same as 'dramp' below.
!-----------------------------------------------------------------------
  nramp_elev = 0

!-----------------------------------------------------------------------
! Optional inverse barometric effects on the elev. b.c.
! If inv_atm_bnd=1, the elev.'s at boundary are corrected by the difference
! between the actual atmos. pressure and a reference pressure (prmsl_ref below)
!-----------------------------------------------------------------------
  inv_atm_bnd = 0 !0: off; 1: on
  prmsl_ref = 101325. !reference atmos. pressure on bnd [Pa]

!-----------------------------------------------------------------------
! Initial condition for T,S. This value only matters for ihot=0 (cold start).
! If flag_ic(1:2)=1, the initial T,S field is read in from temp.ic and salt.ic (horizontally varying).
! If flag_ic(1:2)=2, the initial T,S field is read in from ts.ic (vertical varying).
! If ihot=0 && flag_ic(1:2)=2 || ibcc_mean=1, ts.ic is used for removing mean density profile.
! flag_ic(1) must =flag_ic(2)
!-----------------------------------------------------------------------
  flag_ic(1) = 0 !T
  flag_ic(2) = 0 !S

! initial conditions for other tracers.
! 1: needs inputs [MOD]_hvar_[1,2,...].ic ('1...' is tracer id); format of each file is similar to salt.ic;
!    i.e. horizontally varying i.c. is used for each tracer.
! 2: needs [MOD]_vvar_[1,2,...].ic. Format of each file (for each tracer in tis MOD) is similar to ts.ic
!    (i.e. level #, z-coord., tracer value). Verically varying i.c. is used for each tracer.
! 0: model sets own i.c. (EcoSim; TIMOR)
  flag_ic(3) = 1 !GEN (user defined module)
!  flag_ic(4) = 1 !Age
  flag_ic(5) = 1 !SED3D
  flag_ic(6) = 1 !EcoSim
  flag_ic(7) = 1 !ICM
  flag_ic(8) = 1 !CoSINE
  flag_ic(9) = 1 !FIB
  flag_ic(10) = 1 !TIMOR
  flag_ic(11) = 1 !FABM

!  ic_GEN = 1 !user defined module
!  ic_AGE = 1 !Age
!  ic_SED = 1 !SED3D
!  ic_ECO = 1 !EcoSim
!  ic_ICM = 1 !ICM
!  ic_COS = 1 !CoSINE
!  ic_FIB = 1 !FIB

!-----------------------------------------------------------------------
! Settling vel [m/s] for GEN module (positive downward)
!-----------------------------------------------------------------------
  gen_wsett = 1.e-4

!-----------------------------------------------------------------------
! Mean T,S profile option. If ibcc_mean=1 (or ihot=0 and flag_ic(1)=2), mean profile
! is read in from ts.ic, and will be removed when calculating baroclinic force.
! No ts.ic is needed if ibcc_mean=0.
!-----------------------------------------------------------------------
  ibcc_mean = 0 

!-----------------------------------------------------------------------
! Max. horizontal velocity magnitude, used mainly to prevent problem in 
! bulk aerodynamic module
!-----------------------------------------------------------------------
  rmaxvel = 4.0

!-----------------------------------------------------------------------
!  Following parameters control backtracking
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!  min. vel for invoking btrack and for abnormal exit in quicksearch
!-----------------------------------------------------------------------
  velmin_btrack = 1.e-4
!-----------------------------------------------------------------------
! Nudging factors for starting side/node - add noise to avoid underflow
! The starting location is nudged to: old*(1-btrack_nudge)+btrack_nudge*centroid
! Suggested value: btrack_nudge=9.013e-3
!-----------------------------------------------------------------------
  btrack_nudge= 9.013e-3 

!-----------------------------------------------------------------------
! Behavior when trajectory hits open bnd. If ibtrack_openbnd=0, slide with
! tangential vel; otherwise, stop and exit btrack (recommended?)
!-----------------------------------------------------------------------
  !ibtrack_openbnd = 0 !Default: 0

!-----------------------------------------------------------------------
! Wetting and drying. 
! - if ihhat=1, \hat{H} is made non-negative to enhance robustness near 
! wetting and drying; if ihhat=0, no retriction is imposed for this quantity. 
! - inunfl=0 is used for normal cases and inunfl=1 is used for more accurate wetting
! and drying if grid resolution is sufficiently fine.
! - if shorewafo=1, we impose radiation stress R_s = g*grad(eta) (balance between radiation stress
! gradients and the barotropic gradients) at the numerical shoreline (boundary between
! dry and wet elements). This option ensures that the shallow depth in dry elements does not
! create unphysical and very high wave forces at the shoreline (advised for morphodynamics runs).
!-----------------------------------------------------------------------
  ihhat = 1 
  inunfl = 0
  h0 = 1.e-4     !min. water depth for wetting/drying [m]
  shorewafo = 0 !Matters only if USE_WWM

!-----------------------------------------------------------------------
! Solver options
! USE_PETSC controls the solver type. If it's diabled, the default JCG 
! solver is used. If it's enabled, use PetSc lib. Some of the parameters
! have different meanings under these 2 options. Also with PetSc one can
! use cmd line options to choose solver etc.
!-----------------------------------------------------------------------
  moitn0 = 50 !output spool for solver info; used only with JCG
  mxitn0 = 1500 !max. iteration allowed
  rtol0 = 1.e-12 !error tolerance

!-----------------------------------------------------------------------
! Advection (ELM) option. If nadv=1, backtracking is done using Euler method; 
! nadv=2, using 2nd order Runge-Kutta; if nadv=0, advection in momentum 
! is turned off/on in adv.gr3 (the depths=0,1, or 2 also control methods 
! in backtracking as above). dtb_max/min are the max/min steps allowed -
! actual step is calculated adaptively based on local gradient.
!-----------------------------------------------------------------------
  nadv = 1
  dtb_max = 30. !in sec
  dtb_min = 10.

!-----------------------------------------------------------------------
! If inter_mom=0, linear interpolation is used for velocity at foot of char. line.
! If inter_mom=1 or -1, Kriging is used, and the choice of covariance function is
! specified in 'kr_co'. If inter_mom=1, Kriging is applied to whole domain;
! if inter_mom=-1, the regions where Kriging is used is specified in krvel.gr3 
! (depth=0: no kriging; depth=1: with kriging). 
!-----------------------------------------------------------------------
  inter_mom = 0 
  kr_co = 2 !not used if inter_mom=0

!-----------------------------------------------------------------------
! Transport method. 
! If itr_met=1, upwind method is used. If itr_met>=2, TVD or WENO method is used
! on an element/prism if the total depth (at all nodes of the elem.)>=h_tvd and the flag in
! tvd.prop = 1 for the elem. (tvd.prop is required in this case); 
! otherwise upwind is used for efficiency. 
! itr_met=3 (horizontal TVD) or 4 (horizontal WENO): implicit TVD in the vertical dimension. 
! Also if itr_met==3 and h_tvd>=1.e5, some parts of the code are bypassed for efficiency
! Controls for WENO are not yet in place
!-----------------------------------------------------------------------
  itr_met = 4 !Default: 1
  h_tvd = 5. !used only if itr_met>=2; cut-off depth (m) 
  !If itr_met=3 or 4, need the following 2 tolerances of convergence. The convergence
  !is achieved when sqrt[\sum_i(T_i^s+1-T_i^s)^2]<=eps1_tvd_imp*sqrt[\sum_i(T_i^s)^2]+eps2_tvd_imp
  eps1_tvd_imp = 1.e-4 !suggested value is 1.e-4, but for large suspended load, need to use a smaller value (e.g. 1.e-9)
  eps2_tvd_imp = 1.e-14  

  !Optional hybridized ELM transport for efficiency
  ielm_transport = 1 !1: turn on. Used only with itr_met>=3
  max_subcyc = 10 !used only if ielm_transport/=0. Max # of subcycling per time step in transport allowed


  !if itr_met = 4, the following parameters are needed
  !if itr_met=4 and ipre=1, diagnostic outputs are generated for weno accuracy and stencil quality, 
  !  see subroutine weno_diag in src/Hydro/misc_subs.F90 for details
  ip_weno = 2   !order of accuracy: 0- upwind; 1- linear polynomial, 2nd order; 2- quadratic polynomial, 3rd order
  courant_weno=0.5 !Courant number for weno transport
  nquad = 2  !number of quad points on each side, nquad= 1 or 2
  ntd_weno = 1 !order of temporal discretization: (1) Euler (default); (3): 3rd-order Runge-Kutta (only for benchmarking)
  epsilon1 = 1.e-15   !coefficient for 2nd order weno smoother
  epsilon2 = 1.e-6  !1st coefficient for 3rd order weno smoother (larger values are more prone to numerical dispersion
                     !, 1.e-10 should be fairly safe, recommended values: 1.e-8 ~ 1.e-6 (from the real applications so far)
  i_prtnftl_weno = 0 !option for writing nonfatal errors on invalid temp. or salinity for density: (0) off; (1) on.

  !Inactive at the moment:
  epsilon3 = 1.e-25  !2nd coefficient for 3rd order weno smoother (inactive at the moment)
  !Elad filter has not been implemented yet; preliminary tests showed it might not be necessary
  ielad_weno = 0      !ielad, if ielad=1, use ELAD method to suppress dispersion (inactive at the moment)
  small_elad = 1.e-4  !small (inactive at the moment)

!-----------------------------------------------------------------------
! Atmos. option. nws=3 is reserved for coupling with atmospheric model.
! If nws=0, no atmos. forcing is applied. If nws=1, atmos.
! variables are read in from wind.th. If nws=2, atmos. variables are
! read in from sflux_ files.
! If nws=4, ascii format is used for wind and atmos. pressure at each node (see source code).
! If nws>0, 'iwindoff' can be used to scale wind speed (with windfactor.gr3).
!
! Stress calculation:
! If nws=2, ihconsv=1 and iwind_form=0, the stress is calculated from heat exchange
! routine; 
! Otherwise if iwind_form=-1, the stress is calculated from Pond & Pichard formulation;
! if iwind_form=1, Hwang (2018) formulation (Cd tapers off at high wind).
! If WWM is enabled and icou_elfe_wwm>0 and iwind_form=-2, stress is overwritten by WWM.
! In all cases, if USE_ICE the stress in ice-covered portion is calculated by ICE routine.
!-----------------------------------------------------------------------
  nws = 2
  wtiminc = 3600. !time step for atmos. forcing. Default: same as dt

! (Fei's original with leading "!")  nrampwind = 1 !ramp-up option for atmos. forcing
! nrampwind = 0  ! zy added back !; zy uncommented above line (Fei's original)
  drampwind = 0.  ! zy; Fei original: 2. !needed if nrampwind/=0; ramp-up period in days


  !iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3
  iwind_form = -1 !needed if nws/=0
  !If impose_net_flux/=0 and nws=2, read in net _surface_ heat flux as var 'dlwrf' 
  !(Downward Long Wave) in sflux_rad (solar radiation is still used separately),
  !and if PREC_EVAP is on, also read in net P-E as 'prate' (Surface Precipitation Rate) in sflux_prc.
!  impose_net_flux = 0

!-----------------------------------------------------------------------
! Heat and salt exchange. isconsv=1 needs ihconsv=1; ihconsv=1 needs nws=2.
! If isconsv=1, need to compile with precip/evap module turned on.
!-----------------------------------------------------------------------
  ihconsv = 1 !heat exchange option
  isconsv = 1 !evaporation/precipitation model

! Min total water depth for applying heat/salt exchange (i.e. off if the local water depth <hmin_airsea_ex)
  hmin_airsea_ex = 0.5 ![m]

!-----------------------------------------------------------------------
! Turbulence closure.
!-----------------------------------------------------------------------
  itur = 3 !Default: 0
  dfv0 = 1.e-2 !needed if itur=0
  dfh0 = 1.e-4 !needed if itur=0
  mid = 'KL' !needed if itur=3,5. Use KE if itur=5 
  stab = 'KC' !needed if itur=3 or 5. Use 'GA' if turb_met='MY'; otherwise use 'KC'. 
  xlsc0 = 0.1 !needed if itur=3 or 5. Scale for surface & bottom mixing length (>0)

!-----------------------------------------------------------------------
! Sponge layer for elevation and vel.
! If inu_elev=0, no relaxation is applied to elev.
! If inu_elev=1, relax. constants are specified in elev_nudge.gr3
! and applied to eta=0 (thus a depth=0 means no relaxation).
! Similarly for inu_uv (with input uv_nudge.gr3)
!-----------------------------------------------------------------------
  inu_elev = 0
  inu_uv = 0

!-----------------------------------------------------------------------
! Nudging options for tracers. If inu_[MOD]=0, no nudging is used. If inu_tr=1,
! nudge to initial condition according to relaxation constants specified.
! If inu_tr=2, nudge to values in [MOD]_nu.nc (with step 'step_nu_tr').
! The relaxation constants = [horizontal relax (specified in [MOD]_nudge.gr3) +
! vertical relax] times dt, where vertical relax is a linear function of 
! vnh[1,2] and vnf[1,2], and [MOD] are tracer model names.
! Code will ignore junk values (<=-99) inside [MOD]_nu.nc, so 1 way to avoid 
! nudging for a tracer is to set its nudged values to -9999
!-----------------------------------------------------------------------
  inu_tr(1) = 2 !T
  inu_tr(2) = 2 !S
  inu_tr(3) = 0 !GEN
  inu_tr(4) = 0 !Age
  inu_tr(5) = 0 !SED3D
  inu_tr(6) = 0 !EcoSim 
  inu_tr(7) = 0 !ICM 
  inu_tr(8) = 0 !CoSINE 
  inu_tr(9) = 0 !FIB
  inu_tr(10) = 0 !TIMOR 
  inu_tr(11) = 0 !FABM 

  vnh1 = 400 !vertical nudging depth 1
  vnf1 = 0. !vertical relax \in [0,1]
  vnh2 = 500 !vertical nudging depth 2 (must >vnh1)
  vnf2 = 0. !vertical relax

  step_nu_tr = 21600.  ! zy; Fei's original: 86400. !time step [sec] in all [MOD]_nu.nc (for inu_[MOD]=2)

!-----------------------------------------------------------------------
! Cut-off depth for cubic spline interpolation near bottom when computing horizontal gradients
! e.g. using hgrad_nodes() (radiation stress, and gradients of qnon and qhat in non-hydro model). 
! If depth > h_bcc1 ('deep'),
! a min. (e.g. max bottom z-cor for the element) is imposed in the spline and so a more
! conservative method is used without extrapolation beyond bottom; 
! otherwise constant extrapolation below bottom is used.
!-----------------------------------------------------------------------
  h_bcc1 = 100. !h_bcc1

!-----------------------------------------------------------------------
! Dimensioning parameters for inter-subdomain btrack. 
! If error occurs like 'bktrk_subs: overflow' or 'MAIN: nbtrk > mxnbt'
! gradually increasing these will solve the problem
!-----------------------------------------------------------------------
  s1_mxnbt = 0.5
  s2_mxnbt = 3.5

!-----------------------------------------------------------------------
! Flag for harmonic analysis for elevation. If used , need to turn on cpp flags
! in Makefile first. Otherwise set it to 0.
!-----------------------------------------------------------------------
  iharind = 0

!-----------------------------------------------------------------------
! Conservation check option. If iflux=1, some fluxes are computed
! in regions specified in fluxflag.prop (regional number from -1 to an arbitrary integer).
! in output flux.dat, positive means flux from region n to region n-1 (n>=1)
! output file flux.dat need to be saved before continuing hotstart
!-----------------------------------------------------------------------
!  iflux = 1
  iflux = 0

!-----------------------------------------------------------------------
! Test flags for debugging. These flags should be turned off normally.
!-----------------------------------------------------------------------
! Williamson test #5 (zonal flow over an isolated mount); if
! on, ics must =2
!-----------------------------------------------------------------------
  izonal5 = 0 !"0" - no test; otherwise on

!-----------------------------------------------------------------------
! Rotating Gausshill test with stratified T,S (1: on; 0: off)
! Surface T,S read in from *.ic; code generates stratification
!-----------------------------------------------------------------------
  ibtrack_test = 0

!-----------------------------------------------------------------------
! Rouse profile test (1: on; 0: off)
! If on, must turn on USE_TIMOR
!-----------------------------------------------------------------------
  irouse_test = 0

!-----------------------------------------------------------------------
! Flag to choose FIB model for bacteria decay (used with USE_FIB)
! flag_fib = 1 - Constant decay rate (/day) in .gr3 format
!                (kkfib_1.gr3 and kkfib_2.gr3)
! flag_fib = 2 - Decay rate computed from Canteras et al., 1995
! flag_fib = 3 - Decay rate computed from Servais et al., 2007
!----------------------------------------------------------------------
  flag_fib = 1

!----------------------------------------------------------------------
! Marsh model parameters (only if USE_MARSH is on)
!----------------------------------------------------------------------
  slr_rate = 120. !sea-level rise rate in mm/yr, times morphological acceleration if used

!----------------------------------------------------------------------
! Vegetation model
! If isav=1, need 3 extra inputs: (1) sav_D.gr3 (depth is stem diameter in meters);
! (2) sav_N.gr3 (depth is # of stems per m^2);
! (3) sav_h.gr3 (height of canopy in meters). 
! If one of these depths=0 at a node, the code will set all to 0. 
! If USE_MARSH is on and isav=1, all .gr3 must have constant depths!
!----------------------------------------------------------------------
  isav = 0 !on/off flag

!----------------------------------------------------------------------
! Coupling step with ICE module.
!----------------------------------------------------------------------
  nstep_ice = 1 !call ice module every nstep_ice steps of SCHISM


!----------------------------------------------------------------------
! Specify level #'s if age module is invoked (USE_AGE), for 1st half of tracers only
!----------------------------------------------------------------------
  level_age = 9, -999 !defauly: -999 (all levels)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Physical constants
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!-----------------------------------------------------------------------
! Earth's radii at pole and equator (to define an ellipsoid)
!-----------------------------------------------------------------------
  rearth_pole = 6378206.4
  rearth_eq = 6378206.4

!-----------------------------------------------------------------------
! Specific heat of water (C_p) in J/kg/K
!-----------------------------------------------------------------------
  shw = 4184.d0

!-----------------------------------------------------------------------
! Reference water density for Boussinesq approximation
!-----------------------------------------------------------------------
  rho0 = 1000.d0 !kg/m^3

!-----------------------------------------------------------------------
! Fraction of vertical flux closure adjustment applied at surface, then subtracted
! from all vertical fluxes. This is currently done for T,S only
! 0.0 <= vclose_surf_frac < 1.0
! 1: fully from surface (i.e. no correction as before); 0: fully from bottom 
!-----------------------------------------------------------------------
  vclose_surf_frac=1.0

!-----------------------------------------------------------------------
! Option to enforce strict mass conservation for each tracer model (only works with itr_met=3,4)
! At moment the scheme has not accounted for bottom 'leaking' (e.g. in SED), 
! so iadjust_mass_consv0(5) must =0
!-----------------------------------------------------------------------
  iadjust_mass_consv0(1)=0 !T
  iadjust_mass_consv0(2)=0 !S
  iadjust_mass_consv0(3)=0 !GEN
  iadjust_mass_consv0(4)=0 !AGE
  iadjust_mass_consv0(5)=0 !SED3D (code won't allow non-0 for this module)
  iadjust_mass_consv0(6)=0 !EcoSim
  iadjust_mass_consv0(7)=0 !ICM
  iadjust_mass_consv0(8)=0 !CoSiNE
  iadjust_mass_consv0(9)=0 !Feco
  iadjust_mass_consv0(10)=0 !TIMOR
  iadjust_mass_consv0(11)=0 !FABM
/

&SCHOUT
!-----------------------------------------------------------------------
! Output section - all optional. Values shown are default unless otherwise stated.
!-----------------------------------------------------------------------
! Option for hotstart outputs
!-----------------------------------------------------------------------
  nhot = 1 !1: output *_hotstart every 'hotout_write' steps
  nhot_write = 288 ! zy; Fei original: 1152 !must be a multiple of ihfskip if nhot=1

!-----------------------------------------------------------------------
! Station output option. If iout_sta/=0, need output skip (nspool_sta) and
! a station.in. If ics=2, the cordinates in station.in must be in lon., lat,
! and z (positive upward; not used for 2D variables). 
!-----------------------------------------------------------------------
  iout_sta = 1
  nspool_sta = 2 !zy; Fei original: 12 !needed if iout_sta/=0; mod(nhot_write,nspool_sta) must=0

!-----------------------------------------------------------------------
! Global output options
! The variable names that appear in nc output are shown in {}
!-----------------------------------------------------------------------
  iof_hydro(1) = 1 !0: off; 1: on - elev. [m]  {elev}
  iof_hydro(2) = 0 !air pressure [Pa]  {air_pressure}
  iof_hydro(3) = 0 !air temperature [C] {air_temperature}
  iof_hydro(4) = 0 !Specific humidity [-] {specific_humidity}
  iof_hydro(5) = 0 !solar (shortwave) radiation [W/m/m] {solar_radiation}
  iof_hydro(6) = 0 !sensible flux (positive upward) [W/m/m]  {sensible_flux}
  iof_hydro(7) = 0 !latent heat flux (positive upward) [W/m/m] {latent_heat}
  iof_hydro(8) = 0 !upward longwave radiation (positive upward) [W/m/m] {upward_longwave}
  iof_hydro(9) = 0 !downward longwave radiation (positive downward) [W/m/m] {downward_longwave}
  iof_hydro(10) = 0 !total flux=-flsu-fllu-(radu-radd) [W/m/m] {total_heat_flux}
  iof_hydro(11) = 1 !evaporation rate [kg/m/m/s] {evaporation}
  iof_hydro(12) = 1 !precipitation rate [kg/m/m/s] {precipitation}
  iof_hydro(13) = 0 !Bottom stress vector [kg/m/s^2(Pa)] {bottom_stress}
  iof_hydro(14) = 1 !wind velocity vector [m/s] {wind_speed}
  iof_hydro(15) = 1 !wind stress vector [m^2/s/s] {wind_stress}
  iof_hydro(16) = 1 !depth-averaged vel vector [m/s] {dahv}
  iof_hydro(17) = 0 !vertical velocity [m/s] {vertical_velocity}
  iof_hydro(18) = 1 !water temperature [C] {temp}
  iof_hydro(19) = 1 !water salinity [PSU] {salt}
  iof_hydro(20) = 0 !water density [kg/m^3] {water_density}
  iof_hydro(21) = 0 !eddy diffusivity [m^2/s] {diffusivity}
  iof_hydro(22) = 0 !eddy viscosity [m^2/s] {viscosity}
  iof_hydro(23) = 0 !turbulent kinetic energy {TKE}
  iof_hydro(24) = 0 !turbulent mixing length [m] {mixing_length}
!  iof_hydro(25) = 1 !z-coord {zCoordinates}  3D - this flag will be reset by code to be always on
  iof_hydro(26) = 1 !horizontal velocity vector [m/s] {horizontalVelX,Y}  3D vector
  iof_hydro(27) = 0 !vertical vel. @elem [m/s] {wvel_elem}
  iof_hydro(28) = 0 !T @prism centers [C] {temp_elem}
  iof_hydro(29) = 0 !S @prism centers [PSU] {salt_elem}
  iof_hydro(30) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers  {pressure_gradient}

!-----------------------------------------------------------------------
! Outputs from optional modules. Only uncomment these if the USE_* is on
!-----------------------------------------------------------------------
! Outputs from WWM (USE_WWM must be on in Makefile)
!-----------------------------------------------------------------------
!  iof_wwm(1)  = 1 !sig. height (m) {WWM_1}
!  iof_wwm(2)  = 1 !Mean average period (sec) - TM01 {WWM_2}
!  iof_wwm(3)  = 1 !Zero down crossing period for comparison with buoy (s) - TM02 {WWM_3}
!  iof_wwm(4)  = 1 !Average period of wave runup/overtopping - TM10 {WWM_4}
!  iof_wwm(5)  = 1 !Mean wave number (1/m) {WWM_5}
!  iof_wwm(6)  = 1 !Mean wave length (m) {WWM_6}
!  iof_wwm(7)  = 1 !Mean average energy transport direction (degr) - MWD in NDBC? {WWM_9}
!  iof_wwm(8) = 1 !Mean directional spreading (degr) {WWM_10}
!  iof_wwm(9) = 1 !Discrete peak period (sec) - Tp {WWM_11}
!  iof_wwm(10) = 1 !Continuous peak period based on higher order moments (sec) {WWM_12}
!  iof_wwm(11) = 1 !Peak phase vel. (m/s) {WWM_13}
!  iof_wwm(12) = 1 !Peak n-factor {WWM_14}
!  iof_wwm(13) = 1 !Peak group vel. (m/s) {WWM_15}
!  iof_wwm(14) = 1 !Peak wave number {WWM_16}
!  iof_wwm(15) = 1 !Peak wave length {WWM_17}
!  iof_wwm(16) = 1 !Peak (dominant) direction (degr) {WWM_18}
!  iof_wwm(17) = 1 !Peak directional spreading {WWM_19}
!  iof_wwm(18) = 1 !Discrete peak direction (radian?) {WWM_20}
!  iof_wwm(19) = 1 !Orbital vel. (m/s) {WWM_21}
!  iof_wwm(20) = 1 !RMS Orbital vel. (m/s) {WWM_22}
!  iof_wwm(21) = 1 !Bottom excursion period (sec?) {WWM_23}
!  iof_wwm(22) = 1 !Bottom wave period (sec) {WWM_24}
!  iof_wwm(23) = 1 !Uresell number based on peak period {WWM_25}
!  iof_wwm(24) = 0 !Friction velocity (m/s?) {WWM_26}
!  iof_wwm(25) = 1 !Charnock coefficient {WWM_27}
!  iof_wwm(26) = 1 !Rougness length {WWM_28}
!  iof_wwm(27) = 1 !WWM_energy vector {WWM_energy_dir}
!  iof_wwm(28) = 1 !Wave force vector (m.s-2) computed by wwm @side centers and whole levels {wave_force} 

!-----------------------------------------------------------------------
! Tracer module outputs. In most cases, actual # of outputs depends on # of tracers used
!-----------------------------------------------------------------------
! Outputs for user-defined tracer module (USE_GEN)
!-----------------------------------------------------------------------
!  iof_gen(1) = 1 !1st tracer {GEN_1}
!  iof_gen(2) = 1 !2nd tracer {GEN_2}

!-----------------------------------------------------------------------
! Outputs for (age). Indices from "1" to "ntracer_age/2"; [days]
!-----------------------------------------------------------------------
!  iof_age(1) = 1 {AGE_1}
!  iof_age(2) = 1 {AGE_2}

!-----------------------------------------------------------------------
! Specific outputs in SED3D (USE_SED must be on in Makefile;
! otherwise these are not needed)
!-----------------------------------------------------------------------
!  iof_sed(1) = 1 !bottom depth _change_ from init. condition (m) {SED_depth_change}
!  iof_sed(2) = 1 ! Bed median grain size in the active layer (mm) {SED_D50}
!  iof_sed(3) = 1 ! Bottom shear stress (Pa) {SED_bed_stress}
!  iof_sed(4) = 1 ! Bottom roughness lenghth (mm) {SED_bed_roughness}
!  iof_sed(5) = 1 !total suspended concentration (g/L) {SED_TSC}
!  iof_sed(6) = 1 ! total bed thickness @elem (m) {bed_thickness}
!  iof_sed(7) = 1 ! total bed age over all layers @elem (sec) {bed_age}
!  iof_sed(8) = 1 ! Sediment transport roughness length @elem (m) (z0st_elem) {z0st}
!  iof_sed(9) = 1 !current-ripples roughness length @elem (m) (z0cr_elem) {z0cr}
!  iof_sed(10) = 1 !sand-waves roughness length (m) @elem (z0sw_elem) {z0sw}
!  iof_sed(11) = 1 !wave-ripples roughness length @elem (m) (z0wr_elem) {z0wr}

!  Example of using 2 classes
!  iof_sed(12) = 1 !conc. of 1st class (one output need by each class) [g/L] {SED3D_1}
!  iof_sed(13) = 1 !Bedload transport rate vector (kg.m-1.s-1) for 1st tracer (one output need by tracer) {SED_bdld_1}
!  iof_sed(14) = 1 !Bed fraction 1st tracer (one output need by each class) [-] {SED_bedfrac_1}
!  iof_sed(15) = 1 !conc. of 2nd class {SED3D_2}
!  iof_sed(16) = 1 !Bedload transport of 2nd class {SED_bdld_2}
!  iof_sed(17) = 1 !Bed fraction of 2nd class {SED_bedfrac_3}

!-----------------------------------------------------------------------
! EcoSim outputs 
!-----------------------------------------------------------------------
!  iof_eco(1) = 1 {ECO_1}

!-----------------------------------------------------------------------
! ICM outputs 
!-----------------------------------------------------------------------
!  iof_icm(1) = 1 !Chlorophyll {ICM_Chl}
!  iof_icm(2) = 1 !PH values (ICM_PH on) {ICM_pH}
!  iof_icm(3) = 1 !ICM primary production @elem [gC/m^3/day] {ICM_PrmPrdt}
!  iof_icm(4) = 1 !ICM totoal inorganic nitrogen (DIN) @elem [gN/m^3] {ICM_DIN}
!  iof_icm(5) = 1 !ICM paticulate organic nitrogen (PON) @elem [gN/m^3] {ICM_PON}
!!  sediment flux @elem
!  iof_icm(6) = 1 !ICM bed sediment flux arrays: SED_BENDOC (output name:ICM_SED_BENDOC) @elem [gC/(m^2 day)] {ICM_SED_BENDOC}
!  iof_icm(7) = 1 !ICM bed sediment flux arrays: SED_BENNH4 (output name:ICM_SED_BENNH4) @elem [gC/(m^2 day)] {ICM_SED_BENNH4}
!  iof_icm(8) = 1 !ICM bed sediment flux arrays: SED_BENNO3 (output name:ICM_SED_BENNO3)@elem [gC/(m^2 day)] {ICM_SED_BENNO3}
!  iof_icm(9) = 1 !ICM bed sediment flux arrays: SED_BENPO4 (output name:ICM_SED_BENPO4) @elem [gC/(m^2 day)] {ICM_SED_BENPO4}
!  iof_icm(10) = 1 !ICM bed sediment flux arrays: SED_BENCOD (output name:ICM_SED_BENCOD) @elem [gC/(m^2 day)] {ICM_SED_BENCOD}
!  iof_icm(11) = 1 !ICM bed sediment flux arrays: SED_BENDO (output name:ICM_SED_BENDO) @elem [gC/(m^2 day)] {ICM_SED_BENDO}
!  iof_icm(12) = 1 !ICM bed sediment flux arrays: SED_BENSA (output name:ICM_SED_BENSA) @elem [gC/(m^2 day)] {ICM_SED_BENSA}
!!  sav @elem
!  iof_icm(13) = 1 !ICM SAV leaf biomass @elem [gC/m^3] (k=1 is surface) {ICM_lfsav}
!  iof_icm(14) = 1 !ICM SAV stem biomass @elem [gC/m^3] {ICM_stsav}
!  iof_icm(15) = 1 !ICM SAV root biomass @elem [gC/m^3] {ICM_rtsav}
!  iof_icm(16) = 1 !ICM SAV total leaf biomass @elem [gC/m^2] {ICM_tlfsav}
!  iof_icm(17) = 1 !ICM SAV total stem biomass @elem [gC/m^2] {ICM_tstsav}
!  iof_icm(18) = 1 !ICM SAV total root biomass @elem [gC/m^2] {ICM_trtsav}
!  iof_icm(19) = 1 !ICM SAV canopy height @elem [m] {ICM_hcansav}
!!  sediment conc @elem
!  iof_icm(20) = 1 !bottom NH4 conc {ICM_CNH4}
!  iof_icm(21) = 1 !bottom NO3 conc {ICM_CNH3}
!  iof_icm(22) = 1 !bottom P conc {ICM_CPIP}
!  iof_icm(23) = 1 !bottom Si conc {ICM_CPOS}
!  iof_icm(24) = 1 !bottom CH4 conc {ICM_CCH4}
!  iof_icm(25) = 1 !bottom SO4 conc {ICM_CSO4}
!  iof_icm(26) = 1 !bottom H2S conc {ICM_CH2S}
!  iof_icm(27) = 1 !bottom PON g1 conc {ICM_SEDPON1}
!  iof_icm(28) = 1 !bottom PON g2 conc {ICM_SEDPON2}
!  iof_icm(29) = 1 !bottom PON g3 conc {ICM_SEDPON3}
!  iof_icm(30) = 1 !bottom POP g1 conc {ICM_SEDPOP1}
!  iof_icm(31) = 1 !bottom POP g2 conc {ICM_SEDPOP2}
!  iof_icm(32) = 1 !bottom POP g3 conc {ICM_SEDPOP3}
!  iof_icm(33) = 1 !bottom POC g1 conc {ICM_SEDPOC1}
!  iof_icm(34) = 1 !bottom POC g2 conc {ICM_SEDPOC2}
!  iof_icm(35) = 1 !bottom POC g3 conc {ICM_SEDPOC3}
!  iof_icm(36) = 1 !erosion flux H2S {ICM_EROH2S}
!  iof_icm(37) = 1 !ersoion flux LPOC {ICM_EROLPOC}
!  iof_icm(38) = 1 !ersoion flux RPOC {ICM_ERORPOC}
!!  water column rates @elem
!  iof_icm(39) = 1 !DO consumption {ICM_DO_consumption}
!  iof_icm(40) = 1 !PB growth #1 {ICM_GP1}
!  iof_icm(41) = 1 !PB growth #2 {ICM_GP2}
!  iof_icm(42) = 1 !PB growth #3 {ICM_GP3}
!  water column foundamental tracers @nodes
!  iof_icm(43) = 0 !Zoo. #1 {ICM_1}
!  iof_icm(44) = 0 !Zoo. #2 {ICM_2}
!  iof_icm(45) = 1 !phyto #1 {ICM_3}
!  iof_icm(46) = 1 !phyto #2 {ICM_4}
!  iof_icm(47) = 1 !phyto #3 {ICM_5}
!  iof_icm(48) = 1 !RPOC {ICM_6}
!  iof_icm(49) = 1 !LPOC {ICM_7}
!  iof_icm(50) = 1 !DOC {ICM_8}
!  iof_icm(51) = 1 !RPON {ICM_9}
!  iof_icm(52) = 1 !LPON {ICM_10}
!  iof_icm(53) = 1 !DON {ICM_11}
!  iof_icm(54) = 1 !NH4 {ICM_12}
!  iof_icm(55) = 1 !NO3 {ICM_13}
!  iof_icm(56) = 1 !RPOP {ICM_14}
!  iof_icm(57) = 1 !LPOP {ICM_15}
!  iof_icm(58) = 1 !DOP {ICM_16}
!  iof_icm(59) = 1 !PO4t {ICM_17}
!  iof_icm(60) = 0 !Si- biogenic {ICM_18}
!  iof_icm(61) = 0 !available Si {ICM_19}
!  iof_icm(62) = 1 !COD: Chemical oxygen demand {ICM_20}
!  iof_icm(63) = 1 !DO {ICM_21}
!! following requires ICM_PH on
!  iof_icm(64) = 1 !TIC {ICM_22}
!  iof_icm(65) = 1 !ALK {ICM_23}
!  iof_icm(66) = 1 !CA {ICM_24}
!  iof_icm(67) = 1 !CACO3 {ICM_25}

!-----------------------------------------------------------------------
! CoSINE outputs
!-----------------------------------------------------------------------
!  iof_cos(1) = 1 ! {COS_1}
!  iof_cos(2) = 1 ! {COS_2}
!  iof_cos(3) = 1 ! {COS_3}
!  iof_cos(4) = 1 ! {COS_4}
!  iof_cos(5) = 1 ! {COS_5}
!  iof_cos(6) = 1 ! {COS_6}
!  iof_cos(7) = 1 ! {COS_7}
!  iof_cos(8) = 1 ! {COS_8}
!  iof_cos(9) = 1 ! {COS_9}
!  iof_cos(10) = 1 ! {COS_10}
!  iof_cos(11) = 1 ! {COS_11}
!  iof_cos(12) = 1 ! {COS_12}
!  iof_cos(13) = 1 ! {COS_13}

!-----------------------------------------------------------------------
! Fecal indicating bacteria module
!-----------------------------------------------------------------------
!  iof_fib(1) = 1 ! {FIB_1}

!-----------------------------------------------------------------------
! Specific outputs in SED2D (USE_SED2D must be on in Makefile;
! otherwise these are not needed)
!-----------------------------------------------------------------------
!  iof_sed2d(1)  = 1 !bottom depth _change_ from init. condition (m) {SED2D_depth_change}
!  iof_sed2d(2)  = 1 !drag coefficient used in transport formulae SED2D_Cd{}
!  iof_sed2d(3) = 1 !Courant number (b.qtot.dt / h.dx) {SED2D_cflsed}
!  iof_sed2d(4)    = 1 !Top layer d50 (m) {SED2D_d50}
!  iof_sed2d(5)   = 1 !total transport rate vector (kg/m/s) {SED2D_total_transport}
!  iof_sed2d(6)   = 1 !suspended tranport rate vector (kg/m/s) {SED2D_susp_load}
!  iof_sed2d(7)   = 1 !bedload transport rate vector (kg/m/s) {SED2D_bed_load}
!  iof_sed2d(8)    = 1 !time averaged total transport rate vector (kg/m/s) {SED2D_average_transport}
!  iof_sed2d(9)  = 1 !bottom slope vector (m/m); negative uphill {SED2D_bottom_slope}
!  iof_sed2d(10) = 1 !Total roughness length @elem (m) (z0eq) {z0eq}
!  iof_sed2d(11) = 1 !current-ripples roughness length @elem (m) (z0cr) {z0cr}
!  iof_sed2d(12) = 1 !sand-waves roughness length @elem (m) (z0sw) {z0sw}
!  iof_sed2d(13) = 1 !wave-ripples roughness length @elem (m) (z0wr) {z0wr}
 
!-----------------------------------------------------------------------
!  marsh flags (USE_MARSH on)
!-----------------------------------------------------------------------
!  iof_marsh(1) = 1 ! {marsh_flag}

!-----------------------------------------------------------------------
! Ice module outputs (if USE_ICE is on)
!-----------------------------------------------------------------------
!  iof_ice(1)= 1 !ice advective velcoity vector [m/s] {ICE_velocity}
!  iof_ice(2)= 1 !strain rate @ elem [1/sec] {ICE_strain_rate}
!  iof_ice(3)= 1 !net heat flux to ocean (>0 warm up SST) [W/m/m] {ICE_net_heat_flux}
!  iof_ice(4)= 1 !net fresh water flux to ocean (>0 freshens up SSS) [kg/s/m/m] {ICE_fresh_water_flux}
!  iof_ice(5)= 1 !ice temperature [C] at air-ice interface {ICE_top_T}

!  iof_ice(6)= 1 !ice volume [m] {ICE_tracer_1}
!  iof_ice(7)= 1 !ice concentration [-] {ICE_tracer_2}
!  iof_ice(8)= 1 !snow volume [m] {ICE_tracer_3}

!-----------------------------------------------------------------------
! Analysis module outputs (USE_ANALYSIS)
!-----------------------------------------------------------------------
  iof_ana(1) = 0 !x-component of \nabla air_pres /\rho_0 [m/s/s] {ANA_air_pres_grad_x}
  iof_ana(2) = 0 !y-component of \nabla air_pres /\rho_0 [m/s/s] {ANA_air_pres_grad_y}
  iof_ana(3) = 0  !\alpha*g*\nabla \Psi [m/s/s] (gradient of tidal potential) {ANA_tide_pot_grad_x}
  iof_ana(4) = 0  !\alpha*g*\nabla \Psi [m/s/s] {ANA_tide_pot_grad_y}
  iof_ana(5) = 0 !\nabla \cdot (\mu \nabla u) [m/s/s] (horizontal viscosity) {ANA_hor_viscosity_x}
  iof_ana(6) = 0 !\nabla \cdot (\mu \nabla u) [m/s/s] {ANA_hor_viscosity_y}
  iof_ana(7) = 0 !-g/rho0* \int_z^\eta dr_dx dz  [m/s/s] (b-clinic gradient) {ANA_bclinic_force_x}
  iof_ana(8) = 0 !-g/rho0* \int_z^\eta dr_dy dz  [m/s/s] {ANA_bclinic_force_y}
  iof_ana(9) = 0 !d (\nu du/dz)/dz [m/s/s] - no vegetation effects (vertical viscosity) {ANA_vert_viscosity_x}
  iof_ana(10) = 0 !d (\nu dv/dz)/dz [m/s/s] - no vegetation effects {ANA_vert_viscosity_y}
  iof_ana(11) = 0 !(u \cdot \nabla) u [m/s/s] (momentum advection) {ANA_mom_advection_x}
  iof_ana(12) = 0 !(u \cdot \nabla) u [m/s/s] {ANA_mom_advection_y}
  iof_ana(13) = 0 !gradient Richardson number [-] {ANA_Richardson}
  iof_ana(14) = 1 !min time step at each element over all 
                   !subcycles in horizontal transport solver [s]   {ANA_transport_min_dt_elem}
/