Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Biblio/Biblio.bib	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Biblio/Biblio.bib	(revision 3507)
@@ -2524,4 +2524,11 @@
 }
 
+@ARTICLE{Siddorn_Furner_OM12,
+  author = {J. Siddorn and R. Furner},
+  title = {An analytical stretching function that combines the best attributes of geopotential and terrain-following vertical coordinates},
+  journal = OM,
+  year = {2012},
+  pages = {submitted},
+}
 @ARTICLE{Simmons_al_OM04,
   author = {H. L. Simmons and S. R. Jayne and L. C. {St. Laurent} and A. J. Weaver},
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Chapters/Chap_DOM.tex
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Chapters/Chap_DOM.tex	(revision 3507)
@@ -502,5 +502,5 @@
 time-variation of the free surface so that the transformation is time dependent: 
 $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step 
-bathymetry or $s$-coordinate (hybride and partial step coordinates have not 
+bathymetry or $s$-coordinate (hybrid and partial step coordinates have not 
 yet been tested in NEMO v2.3). 
 
@@ -743,4 +743,5 @@
 levels are defined from the product of a depth field and either a stretching 
 function or its derivative, respectively:
+
 \begin{equation} \label{DOM_sco_ana}
 \begin{split}
@@ -749,16 +750,18 @@
 \end{split}
 \end{equation}
+
 where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point 
 location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea 
 surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean 
 depth, since a mixed step-like and bottom-following representation of the 
-topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided 
-(\rou{zgr\_sco} routine, see \mdl{domzgr}) $h$ is a smooth envelope bathymetry 
-and steps are used to represent sharp bathymetric gradients.
-
-A new flexible stretching function, modified from \citet{Song_Haidvogel_JCP94} is provided as an example:
+topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f).
+The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
+
+Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application.
+A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is provided as an example:
+
 \begin{equation} \label{DOM_sco_function}
 \begin{split}
-z	&= h_c +( h-h_c)\;c s)	\\
+z	&= h_c +( h-h_c)\;c s	\\
 c(s)	&=  \frac{ \left[	  \tanh{ \left( \theta \, (s+b) \right)} 
 	  	 			- \tanh{ \left(  \theta \, b      \right)}  \right]}
@@ -766,17 +769,59 @@
 \end{split}
 \end{equation}
-where $h_c$ is the thermocline depth and $\theta$ and $b$ are the surface and 
+
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]    \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
+\caption{  \label{Fig_sco_function}   
+Examples of the stretching function applied to a seamount; from left to right: 
+surface, surface and bottom, and bottom intensified resolutions}
+\end{center}   \end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+where $h_c$ is the critical depth (\np{rn\_hc}) total depth at which the coordinate transitions from pure $\sigma$ to the stretched coordinate,  and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and 
 bottom control parameters such that $0\leqslant \theta \leqslant 20$, and 
 $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom 
 increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
 
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!tb]    \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
-\caption{  \label{Fig_sco_function}   
-Examples of the stretching function applied to a sea mont; from left to right: 
-surface, surface and bottom, and bottom intensified resolutions}
-\end{center}   \end{figure}
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that:
+
+\begin{equation}
+z = \gamma\left(h+\zeta\right) \quad \text{ with } \quad 0 \leq \gamma \leq 1
+\label{eq:z}
+\end{equation}
+
+The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
+
+\begin{equation} \label{DOM_gamma_deriv}
+\gamma= A\left(\sigma-\frac{1}{2}\left(\sigma^{2}+f\left(\sigma\right)\right)\right)+B\left(\sigma^{3}-f\left(\sigma\right)\right)+f\left(\sigma\right)
+\end{equation}
+
+Where:
+\begin{equation} \label{DOM_gamma}
+f\left(\sigma\right)=\left(\alpha+2\right)\sigma^{\alpha+1}-\left(\alpha+1\right)\sigma^{\alpha+2}
+\end{equation}
+
+This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth:
+
+\begin{equation} \label{DOM_zb}
+Z_b= h a + b
+\end{equation}
+
+where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively.
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!ht]
+   \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/FIG_DOM_compare_coordinates_surface.pdf}
+        \caption{A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.}
+    \label{fig_compare_coordinates_surface}
+\end{figure}
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell depths in real space. This is not to be confused with the hybrid schemes that superimpose geopotential coordinates on terrain following coordinates thus creating a non-analytical vertical coordinate that therefore may suffer from large gradients in the vertical resolutions. This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.
+
+As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than the critical depth $h_c$. In this example two options are available in depths shallower than $h_c$, with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$.
+
+Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as large slopes lead to hydrostatic consistency. A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, and is output as part of the model mesh file at the start of the run.
 
 % -------------------------------------------------------------------------------------------------------------
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Namelist/namzgr_sco
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Namelist/namzgr_sco	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/DOC/TexFiles/Namelist/namzgr_sco	(revision 3507)
@@ -2,11 +2,24 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20)
-   rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1)
-   rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1)
-   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3507)
@@ -27,8 +27,8 @@
    cn_exp      =  "AMM12"  !  experience name 
    nn_it000    =       1   !  first time step
-   nn_itend    =     576   !  last  time step (std 1 day = 576)
+   nn_itend    =    2880   !  last  time step (std 1 day = 288)
    nn_date0    =  20070101 !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
    nn_leapy    =       1   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.  !  start from rest (F) or from a restart file (T)
+   ln_rstart   =  .true.  !  start from rest (F) or from a restart file (T)
    nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T
                            !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
@@ -37,7 +37,7 @@
    cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input)
    cn_ocerst_out = "restart"   !  suffix of ocean restart name (output)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
-   nn_stock    =     576   !  frequency of creation of a restart file (modulo referenced to 1)
-   nn_write    =      12   !  frequency of write in the output file   (modulo referenced to nit000)
+   nn_istate   =       1   !  output the initial state (1) or not (0)
+   nn_stock    =     2880  !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_write    =     144   !  frequency of write in the output file   (modulo referenced to nit000) 
    ln_dimgnnn  = .false.   !  DIMG file format: 1 file for all processors (F) or by processor (T)
    ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%)
@@ -65,17 +65,18 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   ln_s_sh94   = .false.   !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
-   ln_s_sf12   = .true.    !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
-   ln_sigcrit  = .true.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
    rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
    rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_hc       =   50.0    !  critical depth for transition to stretched coordinates
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
    rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
-                           !  SH94 stretching coefficients 
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
    rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
-   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1)
-   rn_bb       =    0.8    !  stretching with SH94 s-sigma
-                           !  SF12 stretching coefficient 
-   rn_alpha    =    4.4    !  stretching with SH94 s-sigma
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
    rn_efold    =    0.0    !  efold length scale for transition to stretched coord
    rn_zs       =    1.0    !  depth of surface grid box
@@ -83,4 +84,6 @@
    rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
    rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
 !-----------------------------------------------------------------------
@@ -94,5 +97,5 @@
    rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
                            !
-   rn_rdt      =  150.     !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_rdt      =  300.     !  time step for the dynamics (and tracer if nn_acc=0)
    nn_baro     =   30      !  number of barotropic time step            ("key_dynspg_ts")
    rn_atfp     =    0.1    !  asselin time filter parameter
@@ -404,5 +407,5 @@
     cn_mask_file = ''                     !  name of mask file (if ln_mask_file=.TRUE.)
     nn_dyn2d      =  2                    !  boundary conditions for barotropic fields
-    nn_dyn2d_dta  =  2                    !  = 0, bdy data are equal to the initial state
+    nn_dyn2d_dta  =  3                    !  = 0, bdy data are equal to the initial state
                                           !  = 1, bdy data are read in 'bdydata   .nc' files
                                           !  = 2, use tidal harmonic forcing data from files
@@ -412,5 +415,5 @@
                            !  = 1, bdy data are read in 'bdydata   .nc' files
     nn_tra        =  1                    !  boundary conditions for T and S
-    nn_tra_dta    =  0                    !  = 0, bdy data are equal to the initial state
+    nn_tra_dta    =  1                    !  = 0, bdy data are equal to the initial state
                            !  = 1, bdy data are read in 'bdydata   .nc' files
     nn_rimwidth  = 10                      !  width of the relaxation zone
@@ -429,5 +432,5 @@
    bn_v3d  =    'amm12_bdyV_u3d' ,         24        , 'vomecrty' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
    bn_tem  =    'amm12_bdyT_tra' ,         24        , 'votemper' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
-   bn_tem  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
+   bn_sal  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
    cn_dir  =    'bdydta/'
    ln_full_vel = .false.
@@ -606,7 +609,7 @@
    ln_hpg_zco  = .false.   !  z-coordinate - full steps                   
    ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation)
-   ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation)
+   ln_hpg_sco  = .false.    !  s-coordinate (standard jacobian formulation)
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
-   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
+   ln_hpg_prj  = .true.   !  s-coordinate (Pressure Jacobian scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
                                  !           centered      time scheme  (F)
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/AMM12_PISCES/EXP00/namelist	(revision 3507)
@@ -27,8 +27,8 @@
    cn_exp      =  "AMM12"  !  experience name 
    nn_it000    =       1   !  first time step
-   nn_itend    =     576   !  last  time step (std 1 day = 576)
+   nn_itend    =    2880   !  last  time step (std 1 day = 288)
    nn_date0    =  20070101 !  date at nit_0000 (format yyyymmdd) used if ln_rstart=F or (ln_rstart=T and nn_rstctl=0 or 1)
    nn_leapy    =       1   !  Leap year calendar (1) or not (0)
-   ln_rstart   =  .false.  !  start from rest (F) or from a restart file (T)
+   ln_rstart   =  .true.  !  start from rest (F) or from a restart file (T)
    nn_rstctl   =       0   !  restart control => activated only if ln_rstart = T
                            !    = 0 nn_date0 read in namelist ; nn_it000 : read in namelist
@@ -37,7 +37,7 @@
    cn_ocerst_in  = "restart"   !  suffix of ocean restart name (input)
    cn_ocerst_out = "restart"   !  suffix of ocean restart name (output)
-   nn_istate   =       0   !  output the initial state (1) or not (0)
-   nn_stock    =     576   !  frequency of creation of a restart file (modulo referenced to 1)
-   nn_write    =      12   !  frequency of write in the output file   (modulo referenced to nit000)
+   nn_istate   =       1   !  output the initial state (1) or not (0)
+   nn_stock    =     2880  !  frequency of creation of a restart file (modulo referenced to 1)
+   nn_write    =     144   !  frequency of write in the output file   (modulo referenced to nit000) 
    ln_dimgnnn  = .false.   !  DIMG file format: 1 file for all processors (F) or by processor (T)
    ln_mskland  = .false.   !  mask land points in NetCDF outputs (costly: + ~15%)
@@ -65,13 +65,27 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =   10.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 7000.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
    rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
-   rn_thetb    =    1.00   !  bottom control parameter  (0<=thetb<= 1)
-   rn_rmax     =    0.30   !  maximum cut-off r-value allowed (0<r_max<1)
-   ln_s_sigma  = .true.    !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
-/
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
+/
+
 !-----------------------------------------------------------------------
 &namdom        !   space and time domain (bathymetry, mesh, timestep)
@@ -79,10 +93,10 @@
    nn_bathy    =    1      !  compute (=0) or read (=1) the bathymetry file
    nn_closea    =   0      !  remove (=0) or keep (=1) closed seas and lakes (ORCA)
-   nn_msh      =    0      !  create (=1) a mesh file or not (=0)
+   nn_msh      =    1      !  create (=1) a mesh file or not (=0)
    rn_hmin     =   -3.     !  min depth of the ocean (>0) or min number of ocean level (<0)
    rn_e3zps_min=   20.     !  partial step thickness is set larger than the minimum of
    rn_e3zps_rat=    0.1    !  rn_e3zps_min and rn_e3zps_rat*e3t, with 0<rn_e3zps_rat<1
                            !
-   rn_rdt      =  150.     !  time step for the dynamics (and tracer if nn_acc=0)
+   rn_rdt      =  300.     !  time step for the dynamics (and tracer if nn_acc=0)
    nn_baro     =   30      !  number of barotropic time step            ("key_dynspg_ts")
    rn_atfp     =    0.1    !  asselin time filter parameter
@@ -102,6 +116,6 @@
    !
    cn_dir        = './'     !  root directory for the location of the runoff files
-   ln_tsd_init   = .false.  !  Initialisation of ocean T & S with T &S input data (T) or not (F)
-   ln_tsd_tradmp = .false.  !  damping of ocean T & S toward T &S input data (T) or not (F)
+   ln_tsd_init   = .false.   !  Initialisation of ocean T & S with T &S input data (T) or not (F)
+   ln_tsd_tradmp = .false.   !  damping of ocean T & S toward T &S input data (T) or not (F)
 /
 !!======================================================================
@@ -302,5 +316,5 @@
                            !  or to SSS only (=1) or no damping term (=0)
    rn_dqdt     =   -40.    !  magnitude of the retroaction on temperature   [W/m2/K]
-   rn_deds     =  -166.67  !  magnitude of the damping on salinity   [mm/day]
+   rn_deds     =  -27.7    !  magnitude of the damping on salinity   [mm/day]
    ln_sssr_bnd =   .true.  !  flag to bound erp term (associated with nn_sssr=2)
    rn_sssr_bnd =   4.e0    !  ABS(Max/Min) value of the damping erp term [mm/day]
@@ -596,7 +610,7 @@
    ln_hpg_zco  = .false.   !  z-coordinate - full steps                   
    ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation)
-   ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation)
+   ln_hpg_sco  = .false.    !  s-coordinate (standard jacobian formulation)
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
-   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
+   ln_hpg_prj  = .true.   !  s-coordinate (Pressure Jacobian scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
                                  !           centered      time scheme  (F)
@@ -875,5 +889,4 @@
                            !  0 < n : debug section number n
 /
-
 !!======================================================================
 !!            ***  Observation & Assimilation namelists ***
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3507)
@@ -65,12 +65,25 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20)
-   rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1)
-   rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1)
-   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
 !-----------------------------------------------------------------------
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 3507)
@@ -66,12 +66,25 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20)
-   rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1)
-   rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1)
-   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
 !-----------------------------------------------------------------------
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3507)
@@ -65,12 +65,25 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20)
-   rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1)
-   rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1)
-   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
 !-----------------------------------------------------------------------
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3507)
@@ -65,12 +65,25 @@
 &namzgr_sco    !   s-coordinate or hybrid z-s-coordinate
 !-----------------------------------------------------------------------
-   rn_sbot_min =  300.     !  minimum depth of s-bottom surface (>0) (m)
-   rn_sbot_max = 5250.     !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-   rn_theta    =    6.0    !  surface control parameter (0<=rn_theta<=20)
-   rn_thetb    =    0.75   !  bottom control parameter  (0<=rn_thetb<= 1)
-   rn_rmax     =    0.15   !  maximum cut-off r-value allowed (0<rn_max<1)
-   ln_s_sigma  = .false.   !  hybrid s-sigma coordinates
-   rn_bb       =    0.8    !  stretching with s-sigma
-   rn_hc       =  150.0    !  critical depth with s-sigma 
+   ln_s_sh94   = .true.    !  Song & Haidvogel 1994 hybrid S-sigma   (T)|
+   ln_s_sf12   = .false.   !  Siddorn & Furner 2012 hybrid S-z-sigma (T)| if both are false the NEMO tanh stretching is applied
+   ln_sigcrit  = .false.   !  use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
+                           !  stretching coefficients for all functions
+   rn_sbot_min =   10.0    !  minimum depth of s-bottom surface (>0) (m)
+   rn_sbot_max = 7000.0    !  maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   rn_hc       =  150.0    !  critical depth for transition to stretched coordinates
+                        !!!!!!!  Envelop bathymetry
+   rn_rmax     =    0.3    !  maximum cut-off r-value allowed (0<r_max<1)
+                        !!!!!!!  SH94 stretching coefficients  (ln_s_sh94 = .true.)
+   rn_theta    =    6.0    !  surface control parameter (0<=theta<=20)
+   rn_bb       =    0.8    !  stretching with SH94 s-sigma   
+                        !!!!!!!  SF12 stretching coefficient  (ln_s_sf12 = .true.)
+   rn_alpha    =    4.4    !  stretching with SF12 s-sigma
+   rn_efold    =    0.0    !  efold length scale for transition to stretched coord
+   rn_zs       =    1.0    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   rn_zb_a     =    0.024  !  bathymetry scaling factor for calculating Zb
+   rn_zb_b     =   -0.2    !  offset for calculating Zb
+                        !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above]
+   rn_thetb    =    1.0    !  bottom control parameter  (0<=thetb<= 1) 
 /
 !-----------------------------------------------------------------------
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3507)
@@ -174,4 +174,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff     !: interface depth between stretching at  V--F
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1              !: Maximum grid stiffness ratio
 
    !!----------------------------------------------------------------------
@@ -294,5 +295,5 @@
          &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
          &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
-         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) )
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) )
 
       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     &
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3507)
@@ -36,4 +36,5 @@
    USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
    USE timing          ! Timing
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
 
    IMPLICIT NONE
@@ -84,4 +85,5 @@
                              CALL dom_zgr      ! Vertical mesh and bathymetry
                              CALL dom_msk      ! Masks
+      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
       IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh
       !
@@ -322,4 +324,67 @@
    END SUBROUTINE dom_ctl
 
+   SUBROUTINE dom_stiff
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_stiff  ***
+      !!                     
+      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
+      !!
+      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
+      !!                Save the maximum in the vertical direction
+      !!                (this number is only relevant in s-coordinates)
+      !!
+      !!                Haney, R. L., 1991: On the pressure gradient force
+      !!                over steep topography in sigma coordinate ocean models. 
+      !!                J. Phys. Oceanogr., 21, 610???619.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk 
+      REAL(wp) ::   zrxmax
+      REAL(wp), DIMENSION(4) :: zr1
+      !!----------------------------------------------------------------------
+      rx1(:,:) = 0.e0
+      zrxmax   = 0.e0
+      zr1(:)   = 0.e0
+      
+      DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO jk = 1, jpkm1
+               zr1(1) = umask(ji-1,jj  ,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji-1,jj  ,jk  )  & 
+                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1)) &
+                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji-1,jj  ,jk  )  &
+                    &                         -gdepw(ji  ,jj  ,jk+1)-gdepw(ji-1,jj  ,jk+1) + rsmall) )
+               zr1(2) = umask(ji  ,jj  ,jk) *abs( (gdepw(ji+1,jj  ,jk  )-gdepw(ji  ,jj  ,jk  )  &
+                    &                         +gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw(ji+1,jj  ,jk  )+gdepw(ji  ,jj  ,jk  )  &
+                    &                         -gdepw(ji+1,jj  ,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(3) = vmask(ji  ,jj  ,jk) *abs( (gdepw(ji  ,jj+1,jk  )-gdepw(ji  ,jj  ,jk  )  &
+                    &                         +gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1)) &
+                    &                        /(gdepw(ji  ,jj+1,jk  )+gdepw(ji  ,jj  ,jk  )  &
+                    &                         -gdepw(ji  ,jj+1,jk+1)-gdepw(ji  ,jj  ,jk+1) + rsmall) )
+               zr1(4) = vmask(ji  ,jj-1,jk) *abs( (gdepw(ji  ,jj  ,jk  )-gdepw(ji  ,jj-1,jk  )  &
+                    &                         +gdepw(ji  ,jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1)) &
+                    &                        /(gdepw(ji  ,jj  ,jk  )+gdepw(ji  ,jj-1,jk  )  &
+                    &                         -gdepw(ji,  jj  ,jk+1)-gdepw(ji  ,jj-1,jk+1) + rsmall) )
+               zrxmax = MAXVAL(zr1(1:4))
+               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax)
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk( rx1, 'T', 1. )
+
+      zrxmax = MAXVAL(rx1)
+
+      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
+         WRITE(numout,*) '~~~~~~~~~'
+      ENDIF
+
+   END SUBROUTINE dom_stiff
+
+
+
    !!======================================================================
 END MODULE domain
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3507)
@@ -172,6 +172,6 @@
             
       IF( ln_sco ) THEN                                         ! s-coordinate
-         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )         !    ! depth
-         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 
+         CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )
+         CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )
          CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )
          CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )
@@ -187,7 +187,8 @@
          CALL iom_rstput( 0, 0, inum4, 'e3v', e3v )
          CALL iom_rstput( 0, 0, inum4, 'e3w', e3w )
+         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
          !
          CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept )    !    ! stretched system
-         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 
+         CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw )
       ENDIF
       
Index: /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3506)
+++ /branches/2012/dev_r3435_UKMO7_SCOORDS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3507)
@@ -1727,7 +1727,5 @@
       !! ** Purpose :   stretch the s-coordinate system
       !!
-      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
-      !!                mixed S/sigma/Z coordinate
-      !!
+      !! ** Method  :   s-coordinate stretch 
       !!
       !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
