#!/bin/ksh
#
# "SMS-like" script to run standalone
#
set -exvu
tdir=${TDIR:-unknown}
npx=1
npy=1
date=20060101
ndays=1
ln_grd=.false.
set -A availdates 20060101
dist=no
pro=.true.
npro=1
sla=.false.
nsla=1
sst=.false.
nsst=1
nsbfil=2
lrnorm=.true.
nfilnf=1
nesnow=0
nubsal=1
nubssh=1
nubvel=0
nbarssh=0
nbarvel=0
nmppsum=2
nbalbkg=0
ninnmax=10
ln_sgl=.false.
ln_ref=.false.
ityp=1
zlam=-160.0
zphi=0.0
zdep=25.0
zval=1.0
ln_tst=.true.
compile=${PLATFORM:-unknown}
debug=0
grid=ORCA2_Z31
global=.true.
prefix="nemovar"
lookup=.false.
lookupres=1.0
nvarex=0
nlevnor=-1

dynmem=no

while getopts DK:t:x:y:i:n:rQa:As:Sp:Pf:OM:E:ubI:1:FL:Z:z:v:UVCWTdg:NFBH:R:4m:G option
do
  case $option in
    D) dynmem=yes;;
    K) compile=$OPTARG;;
    t) tdir=$OPTARG;;
    i) date=$OPTARG;;
    n) ndays=$OPTARG;;
    x) npx=$OPTARG;;
    y) npy=$OPTARG;;
    r) ln_grd=.true.;;
    Q) dist=yes;;
    p) pro=.true.
       npro=$OPTARG;;
    P) pro=.false.;;
    a) sla=.true.
       nsla=$OPTARG;;
    A) sla=.false.;;
    s) sst=.true.
       nsst=$OPTARG;;
    S) sst=.false;;
    f) nsbfil=$OPTARG;;
    O) lrnorm=.false.;;
    M) nfilnf=$OPTARG;;
    E) nesnow=$OPTARG;;
    u) nmppsum=1;;
    b) nbalbkg=1;;
    I) ninnmax=$OPTARG;;
    1) ln_sgl=.true.
       ityp=$OPTARG;;
    F) ln_ref=.true.;;
    L) zlam=$OPTARG;;
    Z) zphi=$OPTARG;;
    z) zdep=$OPTARG;;
    v) zval=$OPTARG;;
    U) nubssh=0
       nubsal=0;;
    V) nubvel=1;;
    C) nbarssh=2;;
    W) nbarvel=1;;
    T) ln_tst=.true.;;
    d) debug=1;;
    g) grid=$OPTARG;;
    B) global=.false.;;
    H) prefix=$OPTARG;;
    G) lookup=.true.;;
    R) lookupres=$OPTARG;;
    4) nvarex=1;;
    m) nlevnor=$OPTARG;;
    \?) exit 1;;
  esac
done
#
# Platform for compilation
#
if [[ $compile = "unknown" ]]; then
   print "Using default for platform $ARCH"
   case $ARCH in
     linux) compile="linux_ia32";;
     rs6000) compile="ecgate";;
     ibm_power4) compile="IBM-ECMWF";;
   esac
fi

if [[ $compile = "unknown" ]]; then
   print "No defaults found for platform $ARCH"
   exit 1
fi
#
# Check that we have testdata for the requested date 
#
datefound=0
ndates=${#availdates[@]}
i=0
while((i<ndates))
do
  print "availdates[$i]=${availdates[$i]}"
  if [[ $date = ${availdates[$i]} ]] ; then
     datefound=1
     index=$i
  fi
  ((i=i+1))
done
if [[ $datefound = 0 ]]; then
  print "Test data not available for date $date"
  exit 1
fi
#
# Check the grid and set grid specific namelists
#
if [[ $grid = ORCA2_Z31 ]]; then
  aht0=2000.0
  aeiv0=2000.0
  ahm0=40000.0
  e3zps_min=5.
  tstep=5760
  nn_fsbc=5 
  dqdt=-200.0
  deds=-27.7
  nmax=800
  nfithts=200
  nfithuv=100
  nfithet=200
  flshts=4.0
  flshuv=2.3
  flshet=4.0
  nlev=31
  normal=background.normalization.nc
  n_cla=1
elif [[ $grid = ORCA1_Z42 ]]; then
  aht0=1000.0
  aeiv0=1000.0
  ahm0=10000.0
  e3zps_min=20.
  tstep=3600.0
  dqdt=-200.0
  deds=-33.333333
  nn_fsbc=6
  nmax=2000
  nfithts=400
  nfithuv=200
  nfithet=400
  flshts=2.0
  flshuv=1.15
  flshet=2.0
  nlev=42
  normal=background.normalization.nc
  n_cla=0
elif [[ $grid = ORCA1_Z64 ]]; then
  aht0=1000.0
  aeiv0=1000.0
  ahm0=10000.0
  e3zps_min=4.
  tstep=3600.0
  dqdt=-200.0
  deds=-33.333333
  nn_fsbc=6
  nmax=2000
  nfithts=400
  nfithuv=200
  nfithet=400
  flshts=2.0
  flshuv=1.15
  flshet=2.0
  nlev=64
  normal=background.normalization.nc
  n_cla=0
else
  print "No setup found for grid $grid"
  exit 1
fi
if [[ $nlevnor != -1 ]]; then
  nlev=$nlevnor
fi

#
# If outfreq is still 1 set to ndays
#
nstepday=$((86400/tstep))
nsteps=$((ndays*nstepday))

if [[ -d $tdir/${prefix}_inner_${grid}_${date}_${compile}_${npx}x${npy} ]]
then
  rm -rf $tdir/${prefix}_inner_${grid}_${date}_${compile}_${npx}x${npy}/*
fi

mkdir $tdir/${prefix}_inner_${grid}_${date}_${compile}_${npx}x${npy} || true

cd $tdir/${prefix}_inner_${grid}_${date}_${compile}_${npx}x${npy}

[[ -d ../nemov3_testdata/${grid}_common ]] || exit 1
ln -s ../nemov3_testdata/${grid}_common/* .

[[ -d ../nemov3_testdata/${grid}_${date}_inner ]] || exit 1
ln -s ../nemov3_testdata/${grid}_${date}_inner/* .

rm -f background.normalization.nc
ln -s ../nemov3_testdata/${grid}_${date}_inner/$normal background.normalization.nc

#
# Set link to non-linear trajectory
#

ln -s ../nemo_outer_${grid}_${date}_${compile}_${npx}x${npy}/assim_trj_*.nc .

if [[ $dist = "yes" ]] ; then
  rm -f dist_coast_uvtf.nc
fi

typeset -Z2 inum
if [[ $pro = ".true." ]]; then
  inum=$npro
  ln -s profiles_${inum}_*.nc profiles_01.nc
fi
if [[ $sla = ".true." ]]; then
  inum=$nsla
  ln -s sla_${inum}_*.nc sla_01.nc
fi
if [[ $sst = ".true." ]]; then
  inum=$nsst
  ln -s sst_${inum}_*.nc sst_01.nc
fi

typeset -Z1 ipro
ipro=$npro
if [[ $ipro = 1 ]]; then
  ln_enact=.true.
else
  ln_enact=.false.
fi

date0=$(echo $date |cut -c1-8)

npxy=$((npx*npy))

cat >namelist.nemovar<<EOF
!-----------------------------------------------------------------------
!       nam_ctl      Control prints & Benchmark
!-----------------------------------------------------------------------
!  ln_ctl     trends control print (expensive!)
!  nprint     level of print (0 no print)
!  nbit_cmp   bit comparison mode parameter (0/1): enables bit comparison between 
!             single and multiple processor runs.
!  nabortx    abort option (stop, abort, mpi_abort etc).           
&namctl
   ln_ctl =  .false.
   nprint =       0
   nbit_cmp =     0
   nabortx =      1
/
&nam_mpp_dyndist 
!-----------------------------------------------------------------------
!       nam_mpp_dyndist Massively Parallel Distribution
!  jpni       number of processors following i
!  jpnj       number of processors following j
!  jpnij      jpnij  number of local domains
!-----------------------------------------------------------------------
   jpni       =  $npx
   jpnj       =  $npy
   jpnij      =  $npxy
/
!-----------------------------------------------------------------------
!       namalg       Incremental algorithm parameters
!-----------------------------------------------------------------------
! nvarex      type of incremental variational assimilation
!                = 0 3D-Var (FGAT)
!                = 1 4D-Var
!                = 2 NMC-type method
! noutmax      maximum number of outer iterations for the experiment
! ninnmax      maximum number of inner iterations for the current 
!              outer iter
! noutit       outer loop iteration counter
! nittrjfrq    frequency of the trajectory used in 4D-VAR
! ln_ref       boolean term for using a reference background state rather
!              than a true background state (typically used for single 
!              observation experiments)
!
&namalg
   nvarex  = $nvarex
   noutmax = 1
   ninnmax = $ninnmax
   noutit  = 1
   nittrjfrq = $nstepday
   ln_ref = $ln_ref
/
!-----------------------------------------------------------------------
!       namtst   assimilation test parameters
!-----------------------------------------------------------------------
!
!  ln_tst          master switch for TL/AD tests
!  ln_tst_obsadj   switch for H adjoint tests
!  ln_tst_bkgadj   switch for B adjoint tests
!  ln_tst_nemotam  switch for M adjoint tests
!  ln_tst_grad     switch for gradient test
!  ln_tst_cpd_tam  switch for M adjoint tests (adjoint of the components + step_tam)
!
 &namtst
   ln_tst = $ln_tst
   ln_tst_obsadj  = .TRUE.
   ln_tst_bkgadj  = .TRUE.
   ln_tst_nemotam = .TRUE.
   ln_tst_grad    = .FALSE.
   ln_tst_cpd_tam = .TRUE.
 /
!-----------------------------------------------------------------------
!       namrun   parameters of the run
!-----------------------------------------------------------------------
!  no         job number
!  cexper     experience name for vairmer format
!  ln_rstart  boolean term for restart (true or false)
!  nrstdt     control of the restart timestep:
!                = 0 restart, do not control nit000 in the restart file.
!                = 1 restart, control nit000 in the restart file. Do not
!                    use the date in the restart file (use ndate0 in namelist)
!                = 2 restart, control nit000 in the restart file, use the date
!                    in the restart file. ndate0 in the namelist is ignored.
!  nit000     number of the first time step
!  nitend     number of the last time step
!  ndate0     initial calendar date aammjj
!  nleapy     Leap year calendar (0/1)
!  ninist     initial state output flag (0/1)
!  nstock     frequency of restart file
!  nwrite     frequency of OUTPUT file
!
!  CAUTION: for usual run scripts, logical value of
!  *******  ln_rstart must be .true. or .false.
!                     and NOT .TRUE. or .FALSE.
&namrun
   no         =       0
   cexper     =  "${grid}"
   ln_rstart  = .false.
   nrstdt     =       0
   nit000     =       1
   nitend     = $nsteps
   ndate0     = $date0
   nleapy     =       1
   ninist     =       0
   nstock     = $nsteps
   nwrite     = $nsteps 
/
!-----------------------------------------------------------------------
!       namdom   space and time domain (bathymetry, mesh, timestep)
!-----------------------------------------------------------------------
!  ntopo      = 0/1 ,compute/read the bathymetry file
!               (mbathy, nb of T-ocean levels)
!  e3zps_min  the thickness of the partial step is set larger than the
!  e3zps_rat     the minimum of e3zps_min and e3zps_rat * e3t
!                (N.B. 0<e3zps_rat<1)
!  nmsh       =1 create a mesh file (coordinates, scale factors, masks)
!  nacc       the acceleration of convergence method
!             = 0, no acceleration, rdt = rdttra
!             = 1, acceleration used, rdt < rdttra(k)
!  atfp       asselin time filter parameter
!  rdt        time step for the dynamics (and tracer if nacc=0)
!  rdtmin     minimum time step on tracers
!  rdtmax     maximum time step on tracers
!  rdth       depth variation of tracer time step
!  rdtbt      barotropic time step (for the time splitting algorithm)
!  nclosea    = 0 no closed sea 
!             = 1 closed sea (Black Sea, Caspian Sea, Great US Lakes...) 
&namdom
   ntopo     =     1
   e3zps_min =     ${e3zps_min}
   e3zps_rat =     0.1
   nmsh      =     0
   nacc      =     0
   atfp      =     0.1
   rdt       =  $tstep
   rdtmin    =  $tstep
   rdtmax    =  $tstep
   rdth      =   800.
   rdtbt     =    90.
   nclosea   =     0
/
!!======================================================================
!!            ***  Surface Boundary Condition namelists  ***
!!======================================================================
!!   namsbc        surface boundary condition
!!   namsbc_ana    analytical         formulation
!!   namsbc_flx    flux               formulation
!!   namsbc_clio   CLIO bulk formulea formulation
!!   namsbc_core   CORE bulk formulea formulation
!!   namsbc_cpl    CouPLed            formulation                       ("key_coupled")
!!   namqsr        penetrative solar radiation
!!   namsbc_rnf    river runoffs
!!   namsbc_ssr    sea surface restoring term (for T and/or S)
!!   namalb        albedo parameters
!!======================================================================
!-----------------------------------------------------------------------
&namsbc        !   Surface Boundary Condition (surface module)
!-----------------------------------------------------------------------
   nn_fsbc     = $nn_fsbc  !  frequency of surface boundary condition computation 
                           !               (= the frequency of sea-ice model call)
   ln_ana      = .false.   !  analytical formulation (T => fill namsbc_ana ) 
   ln_flx      = .true.    !  flux formulation       (T => fill namsbc_flx )
   ln_blk_clio = .false.   !  CLIO bulk formulation  (T => fill namsbc_clio) 
   ln_blk_core = .false.   !  CORE bulk formulation  (T => fill namsbc_core) 
   ln_cpl      = .false.   !  Coupled formulation    (T => fill namsbc_cpl )
   nn_ice      = 0         !  =0 no ice boundary condition   ,
                           !  =1 use observed ice-cover      ,
                           !  =2 ice-model used                             ("key_lim3" or "key_lim2)
   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
   ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
   ln_rnf      = .true.    !  runoffs (T => fill namsbc_rnf)
   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S (T => fill namsbc_ssr)
   nn_fwb      = 1         !  FreshWater Budget: =0 unchecked                              , 
                           !                     =1 global mean of e-p-r set to zero at each nn_fsbc time step   ,
                           !                     =2 annual global mean of e-p-r set to zero
/
!-----------------------------------------------------------------------
&namsbc_ssr    !   surface boundary condition : sea surface restoring
!-----------------------------------------------------------------------
   nn_sstr     =     1     !  add a retroaction term in the surface heat       flux (=1) or not (=0)
   nn_sssr     =     1     !  add a damping     term in the surface freshwater flux (=1) or not (=0)
   dqdt        =   $dqdt   !  magnitude of the retroaction on temperature   [W/m2/K]
   deds        =   $deds   !  magnitude of the damping on salinity   [mm/day/psu]
/
!-----------------------------------------------------------------------
!       nam_zgr       vertical coordinate
!-----------------------------------------------------------------------
!  ln_zco     z-coordinate - full steps      (T/F)
!  ln_zps     z-coordinate - partial steps   (T/F)
!  ln_sco     s- or hybrid z-s-coordinate    (T/F)
&nam_zgr
   ln_zco   =  .false.
   ln_zps   =  .true.
   ln_sco   =  .false.
/
!-----------------------------------------------------------------------
!       namcro   cross land advection
!-----------------------------------------------------------------------
!  n_cla   advection between 2 ocean pts separates by land 
&namcla
   n_cla   = $n_cla
/
!-----------------------------------------------------------------------
!       namlbc   lateral momentum boundary condition
!-----------------------------------------------------------------------
!  shlat   lateral boundary condition on velocity
!                   shlat = 0 , free slip
!               0 < shlat < 2 , partial slip
!                   shlat = 2 , no slip
!               2 < shlat     , strong slip
&namlbc
   shlat  =      2.
/
!-----------------------------------------------------------------------
!       nam_mpp      Massively Parallel Processing
!-----------------------------------------------------------------------
!  c_mpi_send         mpi send/recieve type
!                      = 'S'  : standard blocking send
!                      = 'B'  : buffer blocking send
!                      = 'I'  : immediate non-blocking send
&nam_mpp
   c_mpi_send =  'I'
/
!-----------------------------------------------------------------------
!       namtdp   tracer newtonian damping ('key_tradmp')
!-----------------------------------------------------------------------
!  ndmp    type of damping in temperature and salinity 
!          (='latitude', damping poleward of 'ndmp' degrees and function 
!             of the distance-to-coast. Red and Med Seas as ndmp=-1)
!          (=-1 damping only in Med and Red Seas)
!  ndmpf   =1 create a damping.coeff NetCDF file (the 3D damping array)
!  nmldmp  type of damping in the mixed layer
!          (=0 damping throughout the water column)
!	   (=1 no damping in the mixed layer defined by avt >5cm2/s )
!	   (=2 no damping in the mixed layer defined rho<rho(surf)+.01 )
!  sdmp    surface time scale for internal damping (days)
!  bdmp    bottom time scale for internal damping (days)
!  hdmp    depth of transition between sdmp and bdmp (meters)
&namtdp
   ndmp   =    0
   ndmpf  =    0
   nmldmp =    0
   sdmp   = 1095.
   bdmp   = 1095.
   hdmp   =  800.
/
!-----------------------------------------------------------------------
!       nameos   ocean physical parameters
!-----------------------------------------------------------------------
!  neos    type of equation of state and Brunt-Vaisala frequency
!          = 0, UNESCO (formulation of Jackett and McDougall (1994)
!                                         and of McDougall (1987) )
!          = 1, linear: rho(T)   = rau0 * ( 1.028 - ralpha * T )
!          = 2, linear: rho(T,S) = rau0 * ( rbeta * S - ralpha * T )
!                               with rau0=1020 set in parcst routine
!  ralpha  thermal expension coefficient (linear equation of state)
!  rbeta   saline  expension coefficient (linear equation of state)
&nameos
   neos   =      0
   ralpha =  2.e-4
   rbeta  =  0.001
/
!-----------------------------------------------------------------------
!       namzdf   vertical physics
!-----------------------------------------------------------------------
!  ln_zdfevd  enhanced vertical diffusion         (default T)
!  ln_zdfnpc  Non-Penetrative Convection          (default T)
!  avm0       vertical eddy viscosity for the dynamic (m2/s)
!  avt0       vertical eddy diffusivity for tracers (m2/s)
!  avevd      vertical coefficient for enhanced diffusion scheme (m2/s)
!  nevdm      = 0  apply enhanced mixing on tracer only
!             = 1  apply enhanced mixing on both tracer and momentum
!  ln_zdfexp   vertical physics: (=T)  time splitting (T)     (Default=F)
!                               (=F)  euler backward (F)
!  n_zdfexp   number of sub-timestep for time splitting scheme
&namzdf
   ln_zdfevd = .true.
   ln_zdfnpc = .false.
   avm0      = 1.2e-4
   avt0      = 1.2e-5
   avevd     =   100.
   n_evdm    =     1
   ln_zdfexp =  .false.
   n_zdfexp  =     3
/
!-----------------------------------------------------------------------
!       namtke   turbulent eddy kinetic dependent vertical diffusion
!                ( #ifdef "key_zdftke" )
!-----------------------------------------------------------------------
!  ln_rstke    restart with tke from a run without tke (T) or not (F)
!  nn_itke     number of iterative loops if ln_rstke=T
!  rn_ediff    coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
!  rn_ediss    coef. of the Kolmogoroff dissipation
!  rn_ebb      coef. of the surface input of tke
!  rn_efave    boost of the tke diffusion ( avtke=rn_efave*avm )
!  rn_emin     minimum value of tke [m2/s2]
!  rn_emin0    surface minimum value of tke [m2/s2]
!  nn_mxl      mixing length: = 0 bounded by the distance to surface and bottom
!           = 1 bounded by the local vertical scale factor
!           = 2 first vertical derivative of mixing length bounded by 1
!                             = 3 same criteria as case 2 but applied in a different way
!  nn_pdl      Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
!  nn_avb      profile for constant background used on avt & avm (=1) or not (=0)
!  nn_ave      horizontal averaged on avt (=1) or not (=0)
!  ln_mxl0     mixing length scale surface value as function of wind stress (T) or not (F)
!  rn_lmin     interior buoyancy lenght scale minimum value
!  rn_lmin0    surface  buoyancy lenght scale minimum value
!  nn_etau     exponentially deceasing penetration of tke due to internal & intertial waves
!                      = 0 no penetration ( O(2 km) resolution)
!                      = 1 additional tke source
!                      = 2 additional tke source applied only at the base of the mixed layer
!  nn_htau     type of exponential decrease of tke penetration
!                      = 0  constant 10 m length scale
!                      = 1  ???
!                       = 2  ???
!  rn_efr      fraction of surface tke value which penetrates inside the ocean
!  ln_lc       Langmuir cell effect
!  rn_lc       coef. associated to Langmuir cells
!  nn_havtb    horizontal shape for avtb (=1) or not (=0)
&namtke
   ln_rstke = .false.
   nn_itke     =  50
   rn_ediff    =   0.1
   rn_ediss    =   0.7
   rn_ebb      =   3.75
   rn_efave    =   1.
   rn_emin     =   1.e-6
   rn_emin0    =   1.e-4
   nn_mxl      =   2
   nn_pdl      =   1
   nn_avb      =   0
   nn_ave      =   1
   ln_mxl0     = .false.
   rn_lmin     =   0.4
   rn_lmin0    =   0.4
   nn_etau     =   0
   nn_htau     =   2 
   rn_efr      =   0.05
   ln_lc       = .false.
   rn_lc       =   0.15
   nn_havtb    =   0
/
!-----------------------------------------------------------------------
!       namddm   double diffusive mixing parameterization
!-----------------------------------------------------------------------
!   avts    maximum avs for dd mixing 
!   hsbfr   heat/salt buoyancy flux ratio
&namddm
      avts  = 1.e-4
      hsbfr = 1.6
/
!-----------------------------------------------------------------------
!       nambfr   bottom friction
!-----------------------------------------------------------------------
!  nbotfr  type of bottom friction 
!                  nbotfr = 0 , no slip 
!                  nbotfr = 1 , linear friction
!                  nbotfr = 2 , nonlinear friction
!                  nbotfr = 3 , free slip
!  bfri1   bottom drag coefficient (linear case)
!  bfri2   bottom drag coefficient (non linear case)
!  bfeb2   bottom turbulent kinetic energy  (m^2/s^2)
&nambfr
   nbotfr =       1
   bfri1  =   4.e-4
   bfri2  =   1.e-3
   bfeb2  =  2.5e-3
/
!-----------------------------------------------------------------------
!       namqsr   penetrative solar radiation
!-----------------------------------------------------------------------
!  ln_traqsr : penetrative solar radiation (T) or not (F)     (Default=T)
!  rabs       fraction of qsr associated with xsi1
!  xsi1       first depth of extinction
!  xsi2       second depth of extinction
&namqsr
   ln_traqsr = .true.
   rabs     =   0.58
   xsi1     =   0.35
   xsi2     =   23.0
/
!-----------------------------------------------------------------------
!       namsol   elliptic solver / island / free surface 
!-----------------------------------------------------------------------
!  nsolv     elliptic solver (=1 preconditioned conjugate gradient: pcg)
!                            (=2 successive-over-relaxation: sor)
!                            (=3 FETI: fet, all require "key_feti" defined)
!                            (=4 sor with extra outer halo)
!  nsol_arp  absolute/relative (0/1) precision convergence test
!  nmin      minimum of iterations for the SOR solver
!  nmax      maximum of iterations for the SOR solver
!  nmod      frequency of test for the SOR solver
!  eps       absolute precision of the solver
!  resmax    absolute precision for the SOR solver
!  sor       optimal coefficient for SOR solver
!  epsisl    absolute precision on stream function solver
!  nmisl     maximum pcg iterations for island
!  rnu       strength of the additional force used in free surface b.c.
&namsol
   nsolv     =      2
   nsol_arp  =      1
   nmin      =    300
   nmax      =  $nmax
   nmod      =     10
   eps       =  1.E-6
   resmax    = 1.E-10
   sor       =   1.92
   epsisl    = 1.e-10
   nmisl     =   4000
   rnu       =     1.
/
!-----------------------------------------------------------------------
!       nam_dynadv   option of physics/algorithm (not control by CPP keys)
!-----------------------------------------------------------------------
!  ln_dynadv_vec      vector form flag
!  ln_dynadv_cen2     flux form - 2nd order centered scheme    (default T)
!  ln_dynadv_ubs      flux form - 3rd order UBS scheme         (default F)
&nam_dynadv   
   ln_dynadv_vec   = .TRUE.
   ln_dynadv_cen2  = .FALSE.
   ln_dynadv_ubs   = .FALSE.
/  
!-----------------------------------------------------------------------
!       nam_dynhpg   Hydrostatic pressure gradient option
!-----------------------------------------------------------------------
!  type of pressure gradient scheme (choose one only!)
!     ln_hpg_zco    z-coordinate - full steps                   (default T)
!     ln_hpg_zps    z-coordinate - partial steps (interpolation)
!     ln_hpg_sco    s-coordinate (standard jacobian formulation)
!     ln_hpg_hel    s-coordinate (helsinki modification)
!     ln_hpg_wdj    s-coordinate (weighted density jacobian)
!     ln_hpg_djc    s-coordinate (Density Jacobian with Cubic polynomial)
!     ln_hpg_rot    s-coordinate (ROTated axes scheme)
!  parameters
!    gamm          weighting coefficient (wdj scheme)
&nam_dynhpg
   ln_hpg_zco = .false.
   ln_hpg_zps = .true.
   ln_hpg_sco = .false.
   ln_hpg_hel = .false.
   ln_hpg_wdj = .false.
   ln_hpg_djc = .false.
   ln_hpg_rot = .false.
   gamm       = 0.e0
/
!-----------------------------------------------------------------------
&nam_dynldf    !   lateral diffusion on momentum
!-----------------------------------------------------------------------
!                               !  Type of the operator : 
   ln_dynldf_lap    =  .true.   !     laplacian operator         
   ln_dynldf_bilap  =  .false.  !     bilaplacian operator    
!                               !  Direction of action  : 
   ln_dynldf_level  =  .false.  !     iso-level               
   ln_dynldf_hor    =  .true.   !     horizontal (geopotential)        (require "key_ldfslp" in s-coord.)
   ln_dynldf_iso    =  .false.  !     iso-neutral                      (require "key_ldfslp")
                                !  Coefficient
   ahm0    = $ahm0              !     horizontal eddy viscosity   [m2/s]
   ahmb0   =     0.             !     background eddy viscosity for ldf_iso [m2/s]
/
!-----------------------------------------------------------------------
&nam_dynvor    !   option of physics/algorithm (not control by CPP keys)
!-----------------------------------------------------------------------
   ln_dynvor_ene = .false. !  enstrophy conserving scheme  
   ln_dynvor_ens = .true.  !  energy conserving scheme    
   ln_dynvor_mix = .false. !  mixed scheme               
   ln_dynvor_een = .false. !  energy & enstrophy scheme  
/
!!======================================================================
!!             Tracers & Dynamics vertical physics namelists
!!======================================================================
!!       namzdf        vertical physics
!!       namnpc        non penetrative convection                       
!!       namric        richardson number dependent vertical mixing      ("key_zdfric"      )
!!       namtke        TKE dependent vertical mixing                    ("key_zdftke"      )
!!       namkpp        KPP dependent vertical mixing                    ("key_zdfkpp"      )
!!       namddm        double diffusive mixing parameterization         ("key_zdfddm"      )
!!======================================================================

!-----------------------------------------------------------------------
&namzdf        !   vertical physics
!-----------------------------------------------------------------------
   avm0        =   1.2e-4  !  vertical eddy viscosity   [m2/s]          (background Kz if not "key_zdfcst")
   avt0        =   1.2e-5  !  vertical eddy diffusivity [m2/s]          (background Kz if not "key_zdfcst")
   ln_zdfnpc   = .false.   !  convection: Non-Penetrative algorithm (T) or not (F)
   ln_zdfevd   = .true.    !  convection: enhanced vertical diffusion (T) or not (F)    
   avevd       = 100.      !  vertical coefficient for enhanced diffusion scheme [m2/s]
   n_evdm      =   1       !  enhanced mixing apply on tracer (=0) or on tracer and momentum (=1)
   ln_zdfexp   =  .false.  !  split explicit (T) or implicit (F) time stepping
   n_zdfexp    =   3       !  number of sub-timestep for ln_zdfexp=T
/
!-----------------------------------------------------------------------
&nam_traldf    !   lateral diffusion scheme for tracer 
!-----------------------------------------------------------------------
!                               !  Type of the operator : 
   ln_traldf_lap    =  .true.   !     laplacian operator       
   ln_traldf_bilap  =  .false.  !     bilaplacian operator     
                                !  Direction of action  :
   ln_traldf_level  =  .false.  !     iso-level                
   ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
   ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
!                               !  Coefficient
   aht0        =  $aht0         !     horizontal eddy diffusivity for tracers [m2/s]
   ahtb0       =     0.         !     background eddy diffusivity for ldf_iso [m2/s]
   aeiv0       =  $aeiv0        !     eddy induced velocity coefficient [m2/s]    (require "key_traldf_eiv")
/
!-----------------------------------------------------------------------
!       nam_dynvor   option of physics/algorithm (not control by CPP keys)
!-----------------------------------------------------------------------
!  ln_dynvor_ens   vorticity trends: enstrophy conserving scheme (default T)
!  ln_dynvor_ene      "         "  : energy conserving scheme    (default F)
!  ln_dynvor_mix      "         "  : mixed scheme                (default F)
!  ln_dynvor_een      "         "  : energy & enstrophy scheme   (default F)
&nam_dynvor
   ln_dynvor_ene = .FALSE.
   ln_dynvor_ens = .TRUE.
   ln_dynvor_mix = .FALSE.
   ln_dynvor_een = .FALSE.
//
!                               !  Type of the operator : 
   ln_traldf_lap    =  .true.   !     laplacian operator       
   ln_traldf_bilap  =  .false.  !     bilaplacian operator     
                                !  Direction of action  :
   ln_traldf_level  =  .false.  !     iso-level                
   ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
   ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
!                               !  Coefficient
   aht0        =  $aht0         !     horizontal eddy diffusivity for tracers [m2/s]
   ahtb0       =     0.         !     background eddy diffusivity for ldf_iso [m2/s]
   aeiv0       =  $aeiv0        !     eddy induced velocity coefficient [m2/s]    (require "key_traldf_eiv")
/
!  ln_dynvor_ens   vorticity trends: enstrophy conserving scheme (default T)
!  ln_dynvor_ene      "         "  : energy conserving scheme    (default F)
!  ln_dynvor_mix      "         "  : mixed scheme                (default F)
!  ln_dynvor_een      "         "  : energy & enstrophy scheme   (default F)
&nam_dynvor
   ln_dynvor_ene = .FALSE.
   ln_dynvor_ens = .TRUE.
   ln_dynvor_mix = .FALSE.
   ln_dynvor_een = .FALSE.
/
!-----------------------------------------------------------------------
!       namobs    observation usage switch
!-----------------------------------------------------------------------
! ln_prf     Logical switch for profiles
! ln_t3d     Logical switch for T profile observations
! ln_s3d     Logical switch for S profile observations
! nprosets   Number of profile data sets
! ln_sla     Logical switch for SLA observations
! nslasets   Number of SLA data sets
! ln_sst     Logical switch for SST observations
! nsstsets   Number of SST data sets
! n1dint     Type of vertical interpolation method
!             0 = Linear intepolation.
!             1 = Cubic spline interpolation.
! n2dint     Type of horizontal interpolation method
!             0 = Distance-weighted interpolation
!             1 = Distance-weighted interpolation (small angle)
!             2 = Bilinear interpolation (geographical grid)
!             3 = Bilinear remapping interpolation (general grid)
!             4 = Polynomial interpolation
! ln_nea     Logical switch to reject observations near land
! ln_grd     Redo grid search
! ln_sgl     Single observation experiment
! sstmaxlat  Cut-off latitude for SST observations
&namobs
    ln_prf = $pro
    ln_t3d = $pro
    ln_s3d = $pro
    nprofsets = 1
    ln_sla = $sla
    nslasets = 1
    ln_sst = $sst
    nsstsets = 1
    n1dint = 1
    n2dint = 3
    ln_nea = .FALSE.
    ln_grd = $ln_grd
    ln_sgl = $ln_sgl
    sstmaxlat = 50.0
    ln_grid_global = $global
    ln_grid_search_lookup = $lookup
    grid_search_res = $lookupres
    ln_enact(1) = ${ln_enact}
/
!-----------------------------------------------------------------------
!       namoer    observation error parameters
!-----------------------------------------------------------------------
! nsdt3dp  Flag for type of in situ T observation error
!          standard deviation parameterization
!          = 0 constant
!          = 1 analytical function of depth
!              (fit to Ingleby and Huddleston values)
!          = 2 proportional to the background temperature 
!              error standard deviations
!          = 3 gridded field read in from file and 
!              interpolated to observation points
! nsds3dp  Flag for type of in situ S observation error
!          = 0 constant
!          = 1 analytical function of depth
!              (fit to Ingleby and Huddleston values)
!          = 2 proportional to the background salinity 
!              error standard deviations
!          = 3 gridded field read in from file and 
!              interpolated to observation points
! nsdslap  Flag type of SLA standard error deviation 
!          (0 = constant, 1 = READ 2d field from .nc file)
! nsdsstp  Flag type of SST standard error deviation 
!          (0 = constant, 1 = read 2d field from .nc file)
! rfact3d  Temperature observation standard error deviation 
!          inflation factor at the coast (representativeness error)
! rdist3d  Distance to the coastlines (in metres) at which
!          the temperature obs. error std. deviations are increased
! rfacs3d  Salinity observation standard error deviation 
!          inflation factor at the coast (representativeness error)
! rdiss3d  Distance to the coastlines (in metres) at which
!          the salinity obs. error std. deviations are increased
! rfacsla  Sea level anomoly observation standard error deviation 
!          inflation factor at the coast (representativeness error)
! rdissla  Distance to the coastlines (in metres) at which
!          the Sea level anomoly obs. error std. deviations are increased
! rfacsst  Sea surface temperature observation standard error deviation 
!          inflation factor at the coast (representativeness error)
! rdissst  Distance to the coastlines (in metres) at which
!          the sea surface temperature obs. error std. deviations are increased
! scat3d   Scaling factor used if in situ T standard 
!          error deviation is proportional to the background
!          temperature error standard deviations 
!          (sigo = scat3d x sigb) (if nsdt3dp = 2)
! sdt3dob  Standard error deviation for in situ T DATA 
!          (constant value used if nsdt3dp = 0)
! scas3d   Scaling factor used if in situ S standard 
!          error deviation is proportional to the background
!          salinity error standard deviations 
!          (sigo = scas3d x sigb) (if nsds3dp = 2)
! sds3dob  Standard error deviation for in situ S DATA 
!          (constant value used if nsds3dp = 0)
! sdslaob  Std error dev for along-track altimeter SLA DATA
! sdsstob  Std error dev for SST gridded DATA
! ln_wricsterr  Write constant errors in file for later tuning!
! ln_t3doercoa  Reduce T observation errors near coast lines
! ln_s3doercoa  Reduce S observation errors near coast lines
! ln_slaoercoa  Reduce SLA observation errors near coast lines
! ln_sstoercoa  Reduce SST observation errors near coast lines
!-----------------------------------------------------------------------
&namoer
  nsdt3dp = 1
  nsds3dp = 1
  nsdslap = 0
  nsdsstp = 0
  rfact3d = 2.0
  rdist3d = 300.e3
  rfacs3d = 2.0
  rdiss3d = 300.e3
  rfacsla = 2.0
  rdissla = 300.e3
  rfacsst = 2.0
  rdissst = 300.e3
  scat3d  = 1.0
  sdt3dob = 1.0
  scas3d  = 1.0
  sds3dob = 0.18
  sdslaob = 0.05
  sdsstob = 0.2
  ln_wricsterr = .false.
  ln_t3doercoa = .true.
  ln_s3doercoa = .true.
  ln_slaoercoa = .true.
  ln_sstoercoa = .true.
/
!
!-----------------------------------------------------------------------
!       nambkg    background state parameters
!-----------------------------------------------------------------------
!
!  lrbksd  boolean term for reading background error standard 
!          deviations from file (true or false)
!  nperbk  perturbed background state flag
!           = 1 perturb the background state using the square
!               root of the background error covariance matrix
!           = 0 no perturbation
!  nseedb  seed for the random number generator (set to a large number)
!  nupcov  flag for updating reference in the change of control variable
!           = 1 update the reference state
!           = 0 no update
!  nsdbkgts type of parametrisation scheme for the default background
!           T and S error standard deviations
!            0 = constant
!            1 = normalize vertical gradient by local maximum
!            default = multiply vertical gradient by displacement error
!  nfilsd  flag for filtering standard deviations
!            0 = no filtering
!            1 = apply the horizontal diffusion filter used
!                in the correlation model
!            2 = apply a horizontal 2-grid point Shapiro filter
!  ndtdzm  flag for type of numerical scheme used to compute
!          the derivative dT/dz in the standard deviation
!          parameterisation
!            0 = use finite differences
!            1 = use a cubic spline
!  nbkgran flag for computing randomization estimates of the
!          implied background error standard deviations
!           0 = no computation
!           1 = estimate the diagonal elements of B
!           2 = estimate the diagonal elements of B and H B H^{T}
!           3 = estimate the diagonal elements of B, M B M^{T}, 
!               and H M B M^{T} H^{T}
!  nbkgens ensemble member size for background error standard deviation
!          computation using randomization
!  deltaz  vertical displacement error for the default parametrisation
!          scheme for the T/S standard error deviations
!  wbkgsc  background relative to observations scale factor
!  sdbkgu  background u-velocity  standard error deviation
!  sdbkgv  background v-velocity  standard error deviation
!  sdbkgeeq background SSH standard error deviation       
!           * defines the value at the equator for the parametrised 
!             unbalanced (barotropic) values
!  sdbkgeex background SSH standard error deviation
!           * defines the value in the extratropics for the parametrised 
!             unbalanced (barotropic) values
!  sdelatc  transition latitude where background SSH unbalanced standard 
!           error deviations are set to sdbkgeex 
!  sdbkgt  background temperature standard error deviation
!          * for the parametrised values it defines the maximum allowed
!            value in each profile
!  sdbkgtml background temperature standard error deviation
!          * for the parametrised values it defines the minimum allowed
!            value in the mixed layer
!  sdbkgtdo background temperature standard error deviation
!          * for the parametrised values it defines the minimum allowed
!            value in the deep ocean
!  sdbkgs  background salinity    standard error deviation
!          * for the parametrised values it defines the maximum allowed
!            value in each profile
!  sdbkgsdo background salinity   standard error deviation
!          * for the parametrised values it defines the minimum allowed
!            value in the deep ocean
!
 &nambkg
   lrbksd = .false.
   nperbk = 0,
   nseedb = 934056,
   nupcov = 1,
   nsdbkgts = 2,
   nfilsd = 1,
   ndtdzm = 0,
   nbkgran = 1,
   nbkgens = 100,
   deltaz = 10.0,
   wbkgsc = 1.0,
   sdbkgu = 0.05,
   sdbkgv = 0.05,
   sdbkgeeq = 0.0,
   sdbkgeex = 0.01,
   sdelatc = 20.0,
   sdbkgt = 1.5,
   sdbkgtml = 0.5,
   sdbkgtdo = 0.07,
   sdbkgs = 0.25,
   sdbkgsdo = 0.01,
 /
!
!-----------------------------------------------------------------------
!       nambal    background multivariate balance operator
!-----------------------------------------------------------------------
!
!  nubsal  flag for retaining unbalanced part of salinity (0/1 = off/on)
!  nubssh  flag for retaining unbalanced part of SSH      (0/1 = off/on)
!  nubvel  flag for retaining unbalanced part of velocity (0/1 = off/on)
!  nbasalw flag for reducing salinity balance near coasts (0/1 = off/on)
!  nbasshw flag for reducing SSH      balance near coasts (0/1 = off/on)
!  nbavelw flag for reducing velocity balance near coasts (0/1 = off/on)
!  nbaveld flag for reducing velocity balance at depth    (0/1 = off/on)
!  nbarssh flag for type of baroclinic SSH balance formulation
!           0 = baroclinic SSH approximated by dynamic height
!           1 = baroclinic SSH is computed from the rigid-lid
!               approximation assuming a constant depth
!           2 = baroclinic SSH is computed from the rigid-lid
!               approximation assuming a variable depth, and
!               the resulting elliptic equation is solved iteratively
!  nsolbal solver for the baroclinic SSH elliptic equation
!           1 = preconditioned conjugate gradient PCG
!           2 = successive-over-relaxation SOR
!           3 = FETI
!           4 = successive-over-relaxation SOR with extra halo
!  nsol_arp_fs absolute or relative precision criterion for the
!              SSH balance elliptic equation solver
!              0 = absolute
!              1 = relative
!  nmin_fs     minimum number of iterations for the SSH balance elliptic eq.
!  nmax_fs     maximum number of iterations for the SSH balance elliptic eq.
!  nmod_fs     frequency of solver convergence test
!  eps_fs      solver precision of the SSH balance elliptic eqn.
!  resmax_fs   SOR solver precision of the SSH balance elliptic eqn.
!  sor_fs      optimal coefficient for SOR solver
!  nbarvel flag for type of baroclinic velocity balance formulation
!           0 = baroclinic pressure gradient is integrated from a
!               reference level (level of no motion) to the surface
!               (i.e., velocity is zero below the reference level)
!           1 = baroclinic pressure gradient is integrated over the
!               entire water column
!  nbalbkg flag for computing balanced/unbalanced parts
!          of the background state                        (0/1 = off/on)
!  baldis  distance (in metres) from the coast from which the balance
!          will be gradually reduced to zero (at the coast)
!  baldep  e-folding depth scale for reducing the velocity balance
!          below the surface; reference level for the dynamic height
!          computation (nbarssh = 0); level of no motion when
!          nbarvel = 0
!  nabfil  flag for applying a filter to the local
!          expansion coefficients used in the density balance
!           0 = no filtering
!           1 = apply the horizontal diffusion filter 
!               used in the correlation model
!           2 = apply a horizontal 2-grid point Shapiro filter
!  nmxlcr  mixed-layer depth criterion for T-S balance
!           1   = vertical diffusion coefficient criterion
!           default = density criterion
!  ntsmix  flag for a T-S balance in the mixed layer
!  nsbfil  flag for applying a filter to the local
!          reference coefficient in the salinity balance
!           0 = no filtering
!           1 = apply the horizontal and vertical diffusion filter 
!               used in the correlation model
!           2 = apply a horizontal 2-grid point Shapiro filter
!  ndsdtm  flag for type of numerical scheme used to compute dS/dT
!          in the T-S balance
!           0 = finite differences
!           1 = cubic spline
!  dsdtmax maximum allowed value of |dS/dT| in the T-S balance
!  betals  length scale (in degrees) of the function used to smooth
!          the f-plane and beta-plane geostrophic velocities and
!          of the function used in the pressure correction term
!  nprossh flag for projecting SSH increment to have a zero global mean
!                                                         (0/1 = off/on)
!  nwribal flag for writing in NetCDF format the various coefficients
!          (dS/dT, alpha, beta) used in the balance opts  (0/1 = off/on)
!
 &nambal
   nubsal = $nubsal,
   nubssh = $nubssh,
   nubvel = $nubvel,
   nbasalw = 0,
   nbasshw = 0,
   nbavelw = 1,
   nbaveld = 0,
   nbarssh = $nbarssh,
   nsolbal = 2,
   nsol_arp_fs = 1,
   nmin_fs = 500,
   nmax_fs = 10000,
   nmod_fs = 10,
   eps_fs = 1.e-6,
   resmax_fs = 1.e-14,
   sor_fs = 1.92,
   nbarvel = $nbarvel,
   nbalbkg = $nbalbkg,
   baldis = 200.e3,
   baldep = 1500.0,
   nabfil = 2,
   nmxlcr = 1,
   ntsmix = 0,
   nsbfil = $nsbfil,
   ndsdtm = 1,
   dsdtmax = 1.0,
   betals = 1.55,
   nprossh = 0,
   nwribal = 1,
 /
!
!-----------------------------------------------------------------------
!       namcor    background correlation operator parameters
!-----------------------------------------------------------------------
!  lrnorm   boolean term for reading filter normalization factors
!           from file (true or false)
!  nfilnf   method for computing filter normalization factors:
!             = 0 exact
!             = 1 randomization
!  nlvnow   number of levels on which to compute the filter 
!           normalization factors with the exact method
!  nesnow   size of current ensemble for computing the filter
!           normalization factors with the randomization method
!  ndifsch  Numerical scheme for horizontal (T, S, U, V, SSH)
!           and vertical (T, S, U, V)
!           0 = explicit
!           1 = implicit
!  nfithts  total number of iterations for horizontal 2D filter on T-S
!           (default), or
!           total number of iterations for isopycnal  3D filter on T-S
!           (if 'key_bkgcor_filtsiso') 
!  nfithuv  total number of iterations for horizontal 2D filter on u-v
!           (default), or
!           total number of iterations for isopycnal  3D filter on u-v
!           (if 'key_bkgcor_filuviso') 
!  nfithet  total number of iterations for horizontal 2D filter on SSH
!  nfitzt   total number of iterations for vertical   1D filter on T
!  nfitzs   total number of iterations for vertical   1D filter on S
!  nfitzu   total number of iterations for vertical   1D filter on u
!  nfitzv   total number of iterations for vertical   1D filter on v
!  flshts   correlation horizontal length scale for T and S (specified
!           in degrees then converted to metres)
!  flshuv   correlation horizontal length scale for u and v (specified
!           in degrees then converted to metres)
!  flshet   correlation horizontal length scale for SSH     (specified
!           in degrees then converted to metres)
!  flszt    correlation vertical   length scale for T (in metres)
!  flszs    correlation vertical   length scale for S (in metres)
!  flszu    correlation vertical   length scale for u (in metres)
!  flszv    correlation vertical   length scale for v (in metres)
!  fzscml   stretching factor for the vertical scales within the 
!           mixed layer defined such that 
!            vertical length scale = fzscml x mixed-layer depth
!           (   #if defined "key_bkgcor_filztscoef3d"
!            or #if defined "key_bkgcor_filzuvcoef3d")
!  fzscmr   stretching factor for the vertical scales below the 
!           mixed layer defined such that 
!            vertical length scale = fzscmr x local vertical grid-size
!  shflhts  strength of lateral "shear" along coast for filter on T-S
!  shflhuv  strength of lateral "shear" along coast for filter on u-v 
!                 shfhlt(v) = 0 , "free slip"
!             0 < shfhlt(v) < 2 , "partial slip"
!                 shfhlt(v) = 2 , "no slip"
!             2 < shfhlt(v)     , "strong slip"
!
 &namcor
   lrnorm  = $lrnorm,
   lnorrep = .true.
   nfilnf  = $nfilnf,
   nlvnow  = $nlev,
   nesnow  = $nesnow,
   ndifsch = 0,0,
   nfitht = $nfithts,
   nfiths = $nfithts,
   nfithv = $nfithuv,
   nfithu = $nfithuv,
   nfithet = $nfithet,
   nfitzt  = 10,
   nfitzs  = 10,
   nfitzu  = 10,
   nfitzv  = 10,
   flshs  = $flshts,
   flshv  = $flshuv,
   flsht  = $flshts,
   flshu  = $flshuv,
   flshet  = $flshet,
   flszt   = 10.0,
   flszs   = 10.0,
   flszu   = 10.0,
   flszv   = 10.0,
   fzscml  = 0.75,
   fzscmr  = 1.0,
 /
!
!-----------------------------------------------------------------------
!       namopp    optimization parameters
!-----------------------------------------------------------------------
!
!  repsg   precision of stopping criterion                 
!  toltst  accepted tolerance for the relative error       
!          between the cost/gradient norm estimated by
!          the minimization routine and those computed by
!          the simulator
!  northo  flag for invoking orthogonalization             
!          (0/1 = off/on)
!  npcond  flag for type of preconditioner
!            0 = no preconditioner
!            1 = LBFGS
!            2 = Ritz
!  nppvmx  maximum number of preconditioning vectors
!           - Ritz vectors (Ritz)
!           - vector pairs of gradient and control differences (LBFGS)
!  npairs  number of allowed preconditioning vector pairs 
!          per inner iteration
!  npurge  flag for emptying the preconditioner after each outer iter.
!          (0/1 = off/on)
!  reperr  precision for eigenpair preconditioner          
!  npress  output printing control                         
!            npress = 2 (typical) for CG
!  nptest  flag for testing the divergence projection    
!          (0/1 = off/on)
!  nopdia  flag for computing diagnostics in model space 
!          (0/1 = off/on)
!  nmppsum summation method to use (1=simple, 2=reproducible)
!  
 &namopp
   repsg  = 1.e-10,
   toltst = 1.e-2,
   northo = 1,
   npcond = 0,
   nppvmx = 0,
   npairs = 0,
   npurge = 0,
   reperr = 1.e-5,
   npress = 5,
   nptest = 0,
   nopdia = 0,
   nmppsum= $nmppsum
 /
!-----------------------------------------------------------------------
!       namsgl    single observation usage switch
!-----------------------------------------------------------------------
 &namsgl
   ln_inc = .true.
   ityp   = $ityp
   istp   =  1
   zlam   = $zlam
   zphi   = $zphi
   zval   = $zval
   zdep   = $zdep
 /
EOF

NPROC=$((npx*npy))
typeset -Z4 p
p=0
while ((p<$NPROC))
do
  ln -s geothermal_heating.nc geothermal_heating_${p}.nc
  ln -s coordinates.nc coordinates_${p}.nc
  ((p=p+1))
done
