Index: branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3597)
+++ branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 3600)
@@ -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_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3597)
+++ branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 3600)
@@ -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_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3597)
+++ branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 3600)
@@ -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, 'gdept_0' , gdept_0 )    !    ! stretched system
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 )
+         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 )
       ENDIF
       
Index: branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3597)
+++ branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 3600)
@@ -15,4 +15,5 @@
    !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
    !!            3.3  ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
+   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn adn Furner stretching functio
    !!----------------------------------------------------------------------
 
@@ -28,4 +29,5 @@
    !!       zgr_sco      : s-coordinate
    !!       fssig        : sigma coordinate non-dimensional function
+   !!       fgamma       : Siddorn and Furner stretching function
    !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
    !!---------------------------------------------------------------------
@@ -47,13 +49,24 @@
 
    !                                       !!* Namelist namzgr_sco *
+   LOGICAL  ::   ln_s_sh94   = .false.      ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
+   LOGICAL  ::   ln_s_sf12   = .true.       ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
+   LOGICAL  ::   ln_sigcrit  = .false.      ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 
+   !
    REAL(wp) ::   rn_sbot_min =  300._wp     ! minimum depth of s-bottom surface (>0) (m)
    REAL(wp) ::   rn_sbot_max = 5250._wp     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1)
+   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for transition from sigma to stretched coordinates
+   ! Song and Haidvogel 1994 stretching parameters
    REAL(wp) ::   rn_theta    =    6.00_wp   ! surface control parameter (0<=rn_theta<=20)
    REAL(wp) ::   rn_thetb    =    0.75_wp   ! bottom control parameter  (0<=rn_thetb<= 1)
-   REAL(wp) ::   rn_rmax     =    0.15_wp   ! maximum cut-off r-value allowed (0<rn_rmax<1)
-   LOGICAL  ::   ln_s_sigma  = .false.      ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
-   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter for song and haidvogel stretching
+   REAL(wp) ::   rn_bb       =    0.80_wp   ! stretching parameter 
    !                                        ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
-   REAL(wp) ::   rn_hc       =  150._wp     ! Critical depth for s-sigma coordinates
+   ! Siddorn and Furner stretching parameters
+   REAL(wp) ::   rn_alpha    =    4.4_wp    ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
+   REAL(wp) ::   rn_efold    =    0.0_wp    !  efold length scale for transition to stretched coord
+   REAL(wp) ::   rn_zs       =    1.0_wp    !  depth of surface grid box
+                           !  bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
+   REAL(wp) ::   rn_zb_a     =    0.024_wp  !  bathymetry scaling factor for calculating Zb
+   REAL(wp) ::   rn_zb_b     =   -0.2_wp    !  offset for calculating Zb
 
   !! * Substitutions
@@ -1034,56 +1047,4 @@
    END SUBROUTINE zgr_zps
 
-
-   FUNCTION fssig( pk ) RESULT( pf )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_init  ***
-      !!       
-      !! ** Purpose :   provide the analytical function in s-coordinate
-      !!          
-      !! ** Method  :   the function provide the non-dimensional position of
-      !!                T and W (i.e. between 0 and 1)
-      !!                T-points at integer values (between 1 and jpk)
-      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
-      REAL(wp)             ::   pf   ! sigma value
-      !!----------------------------------------------------------------------
-      !
-      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
-         &     - TANH( rn_thetb * rn_theta                                )  )   &
-         & * (   COSH( rn_theta                           )                      &
-         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
-         & / ( 2._wp * SINH( rn_theta ) )
-      !
-   END FUNCTION fssig
-
-
-   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE eos_init  ***
-      !!
-      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
-      !!
-      !! ** Method  :   the function provides the non-dimensional position of
-      !!                T and W (i.e. between 0 and 1)
-      !!                T-points at integer values (between 1 and jpk)
-      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
-      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
-      REAL(wp)             ::   pf1   ! sigma value
-      !!----------------------------------------------------------------------
-      !
-      IF ( rn_theta == 0 ) then      ! uniform sigma
-         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
-      ELSE                        ! stretched sigma
-         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
-            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
-            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
-      ENDIF
-      !
-   END FUNCTION fssig1
-
-
    SUBROUTINE zgr_sco
       !!----------------------------------------------------------------------
@@ -1110,29 +1071,38 @@
       !!            esigt(k) = fsdsig(k    )
       !!            esigw(k) = fsdsig(k-0.5)
-      !!      This routine is given as an example, it must be modified
-      !!      following the user s desiderata. nevertheless, the output as
+      !!      Three options for stretching are give, and they can be modified
+      !!      following the users requirements. Nevertheless, the output as
       !!      well as the way to compute the model levels and scale factors
-      !!      must be respected in order to insure second order a!!uracy
+      !!      must be respected in order to insure second order accuracy
       !!      schemes.
       !!
-      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
+      !!      The three methods for stretching available are:
+      !! 
+      !!           s_sh94 (Song and Haidvogel 1994)
+      !!                a sinh/tanh function that allows sigma and stretched sigma
+      !!
+      !!           s_sf12 (Siddorn and Furner 2012?)
+      !!                allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !! 
+      !!          s_tanh  (Madec et al 1996)
+      !!                a cosh/tanh function that gives stretched coordinates        
+      !!
       !!----------------------------------------------------------------------
       !
       INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
       INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
-      REAL(wp) ::   zcoeft, zcoefw, zrmax, ztaper   ! temporary scalars
+      REAL(wp) ::   zrmax, ztaper   ! temporary scalars
       !
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
-
-      NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
-      !!----------------------------------------------------------------------
+
+      NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
+                           rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
+     !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('zgr_sco')
       !
       CALL wrk_alloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
-      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
-      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
       !
       REWIND( numnam )                       ! Read Namelist namzgr_sco : sigma-stretching parameters
@@ -1144,19 +1114,23 @@
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namzgr_sco'
-         WRITE(numout,*) '      sigma-stretching coeffs '
-         WRITE(numout,*) '      maximum depth of s-bottom surface (>0)       rn_sbot_max   = ' ,rn_sbot_max
-         WRITE(numout,*) '      minimum depth of s-bottom surface (>0)       rn_sbot_min   = ' ,rn_sbot_min
-         WRITE(numout,*) '      surface control parameter (0<=rn_theta<=20)  rn_theta      = ', rn_theta
-         WRITE(numout,*) '      bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ', rn_thetb
-         WRITE(numout,*) '      maximum cut-off r-value allowed              rn_rmax       = ', rn_rmax
-         WRITE(numout,*) '      Hybrid s-sigma-coordinate                    ln_s_sigma    = ', ln_s_sigma
-         WRITE(numout,*) '      stretching parameter (song and haidvogel)    rn_bb         = ', rn_bb
-         WRITE(numout,*) '      Critical depth                               rn_hc         = ', rn_hc
-      ENDIF
-
-      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
-      esigt3  = 0._wp   ;   esigw3  = 0._wp 
-      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
-      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
+         WRITE(numout,*) '     stretching coeffs '
+         WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
+         WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
+         WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
+         WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
+         WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
+         WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
+         WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
+         WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
+         WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
+         WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
+         WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
+         WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
+         WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
+         WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
+         WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
+         WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
+         WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
+      ENDIF
 
       hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
@@ -1352,139 +1326,28 @@
       ! non-dimensional "sigma" for model level depth at w- and t-levels
 
-      IF( ln_s_sigma ) THEN        ! Song and Haidvogel style stretched sigma for depths
-         !                         ! below rn_hc, with uniform sigma in shallower waters
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-
-               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
-                  DO jk = 1, jpk
-                     gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
-                     gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
-                  END DO
-               ELSE ! shallow water, uniform sigma
-                  DO jk = 1, jpk
-                     gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
-                     gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
-                  END DO
-               ENDIF
-               IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw3 1 jpk    ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
-               !
-               DO jk = 1, jpkm1
-                  esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
-                  esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
-               END DO
-               esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
-               esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
-               !
-               ! Coefficients for vertical depth as the sum of e3w scale factors
-               gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
-               DO jk = 2, jpk
-                  gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
-               END DO
-               !
-               DO jk = 1, jpk
-                  zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
-                  zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-                  gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
-                  gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
-                  gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
-               END DO
-               !
-            END DO   ! for all jj's
-         END DO    ! for all ji's
-
-         DO ji = 1, jpim1
-            DO jj = 1, jpjm1
-               DO jk = 1, jpk
-                  esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
-                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                  esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
-                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                  esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
-                     &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
-                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                  esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
-                     &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                  esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
-                     &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                  !
-                  e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  !
-                  e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-                  e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) )
-               END DO
-            END DO
-         END DO
-
-         CALL lbc_lnk( e3t , 'T', 1._wp )
-         CALL lbc_lnk( e3u , 'U', 1._wp )
-         CALL lbc_lnk( e3v , 'V', 1._wp )
-         CALL lbc_lnk( e3f , 'F', 1._wp )
-         CALL lbc_lnk( e3w , 'W', 1._wp )
-         CALL lbc_lnk( e3uw, 'U', 1._wp )
-         CALL lbc_lnk( e3vw, 'V', 1._wp )
-
-         !
-      ELSE   ! not ln_s_sigma
-         !
-         DO jk = 1, jpk
-           gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
-           gsigt(jk) = -fssig( REAL(jk,wp)        )
-         END DO
-         IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
-         !
-         ! Coefficients for vertical scale factors at w-, t- levels
-!!gm bug :  define it from analytical function, not like juste bellow....
-!!gm        or betteroffer the 2 possibilities....
-         DO jk = 1, jpkm1
-            esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
-            esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
-         END DO
-         esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
-         esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
-
-!!gm  original form
-!!org DO jk = 1, jpk
-!!org    esigt(jk)=fsdsig( FLOAT(jk)     )
-!!org    esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
-!!org END DO
-!!gm
-         !
-         ! Coefficients for vertical depth as the sum of e3w scale factors
-         gsi3w(1) = 0.5_wp * esigw(1)
-         DO jk = 2, jpk
-            gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
-         END DO
-!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
-         DO jk = 1, jpk
-            zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
-            zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
-            gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
-            gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
-            gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
-         END DO
-!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               DO jk = 1, jpk
-                 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-                 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-                 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
-                 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
-                 !
-                 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
-                 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
-                 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
-               END DO
-            END DO
-         END DO
-         !
-      ENDIF ! ln_s_sigma
-
-
+
+!========================================================================
+! Song and Haidvogel  1994 (ln_s_sh94=T)
+! Siddorn and Furner 2012 (ln_sf12=T)
+! or  tanh function       (both false)                    
+!========================================================================
+      IF      ( ln_s_sh94 ) THEN 
+                           CALL s_sh94()
+      ELSE IF ( ln_s_sf12 ) THEN
+                           CALL s_sf12()
+      ELSE                 
+                           CALL s_tanh()
+      ENDIF 
+
+      CALL lbc_lnk( e3t , 'T', 1._wp )
+      CALL lbc_lnk( e3u , 'U', 1._wp )
+      CALL lbc_lnk( e3v , 'V', 1._wp )
+      CALL lbc_lnk( e3f , 'F', 1._wp )
+      CALL lbc_lnk( e3w , 'W', 1._wp )
+      CALL lbc_lnk( e3uw, 'U', 1._wp )
+      CALL lbc_lnk( e3vw, 'V', 1._wp )
+
+      fsdepw(:,:,:) = gdepw (:,:,:)
+      fsde3w(:,:,:) = gdep3w(:,:,:)
       !
       where (e3t   (:,:,:).eq.0.0)  e3t(:,:,:) = 1.0
@@ -1544,5 +1407,5 @@
             &                          ' w ', MAXVAL( fse3w (:,:,:) )
       ENDIF
-      !
+      !  END DO
       IF(lwp) THEN                                  ! selected vertical profiles
          WRITE(numout,*)
@@ -1574,28 +1437,447 @@
       ENDIF
 
-!!gm bug?  no more necessary?  if ! defined key_helsinki
-      DO jk = 1, jpk
+!================================================================================
+! check the coordinate makes sense
+!================================================================================
+      DO ji = 1, jpi
          DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
-                  WRITE(ctmp1,*) 'zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
-                  CALL ctl_stop( ctmp1 )
-               ENDIF
-               IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
-                  WRITE(ctmp1,*) 'zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
-                  CALL ctl_stop( ctmp1 )
-               ENDIF
-            END DO
-         END DO
-      END DO
-!!gm bug    #endif
+
+            IF( hbatt(ji,jj) > 0._wp) THEN
+               DO jk = 1, mbathy(ji,jj)
+                 ! check coordinate is monotonically increasing
+                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
+                    WRITE(ctmp1,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
+                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
+                    CALL ctl_stop( ctmp1 )
+                 ENDIF
+                 ! and check it has never gone negative
+                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
+                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
+                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    CALL ctl_stop( ctmp1 )
+                 ENDIF
+                 ! and check it never exceeds the total depth
+                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
+                    CALL ctl_stop( ctmp1 )
+                 ENDIF
+               END DO
+
+               DO jk = 1, mbathy(ji,jj)-1
+                 ! and check it never exceeds the total depth
+                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                    WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    CALL ctl_stop( ctmp1 )
+                 ENDIF
+               END DO
+
+            ENDIF
+
+         END DO
+      END DO
       !
       CALL wrk_dealloc( jpi, jpj,      zenv, ztmp, zmsk, zri, zrj, zhbat                           )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
+      !
+   END SUBROUTINE zgr_sco
+
+!!======================================================================
+   SUBROUTINE s_sh94()
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_sh94  ***
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
+      !!                mixed S/sigma coordinate
+      !!
+      !! Reference : Song and Haidvogel 1994. 
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
+
+      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
+      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
+
+      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
+      esigt3  = 0._wp   ;   esigw3  = 0._wp 
+      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
+      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
+
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+
+            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
+               DO jk = 1, jpk
+                  gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
+                  gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
+               END DO
+            ELSE ! shallow water, uniform sigma
+               DO jk = 1, jpk
+                  gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
+                  gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
+                  END DO
+            ENDIF
+            !
+            DO jk = 1, jpkm1
+               esigt3(ji,jj,jk  ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
+               esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
+            END DO
+            esigw3(ji,jj,1  ) = 2._wp * ( gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ) )
+            esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) )
+            !
+            ! Coefficients for vertical depth as the sum of e3w scale factors
+            gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1)
+            DO jk = 2, jpk
+               gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
+            END DO
+            !
+            DO jk = 1, jpk
+               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
+               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
+               gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
+            END DO
+           !
+         END DO   ! for all jj's
+      END DO    ! for all ji's
+
+      DO ji = 1, jpim1
+         DO jj = 1, jpjm1
+            DO jk = 1, jpk
+               esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk)     &
+                  &                + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !
+               e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               !
+               e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+               e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
+            END DO
+        END DO
+      END DO
+
       CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
       CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_sco')
-      !
-   END SUBROUTINE zgr_sco
+
+   END SUBROUTINE s_sh94
+
+   SUBROUTINE s_sf12
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_sf12 *** 
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
+      !!                mixed S/sigma/Z coordinate
+      !!
+      !!                This method allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !!
+      !!
+      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      REAL(wp) ::   fsmth          ! smoothing around critical depth
+      REAL(wp) ::   zss, zbb       ! Surface and bottom cell thickness in sigma space
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3           
+
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
+      CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
+
+      gsigw3  = 0._wp   ;   gsigt3  = 0._wp   ;   gsi3w3  = 0._wp
+      esigt3  = 0._wp   ;   esigw3  = 0._wp 
+      esigtu3 = 0._wp   ;   esigtv3 = 0._wp   ;   esigtf3 = 0._wp
+      esigwu3 = 0._wp   ;   esigwv3 = 0._wp
+
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+
+          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
+              
+              zbb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
+                                                     ! could be changed by users but care must be taken to do so carefully
+              zbb = 1.0_wp-(zbb/hbatt(ji,jj))
+            
+              zss = rn_zs / hbatt(ji,jj) 
+              
+              IF (rn_efold /= 0.0_wp) THEN
+                fsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
+              ELSE
+                fsmth = 1.0_wp 
+              ENDIF
+               
+              DO jk = 1, jpk
+                gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
+                gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
+              ENDDO
+              gsigw3(ji,jj,:) = fgamma( gsigw3(ji,jj,:), zbb, zss, fsmth  )
+              gsigt3(ji,jj,:) = fgamma( gsigt3(ji,jj,:), zbb, zss, fsmth  )
+ 
+          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
+
+            DO jk = 1, jpk
+              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
+            END DO
+
+          ELSE  ! shallow water, z coordinates
+
+            DO jk = 1, jpk
+              gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+              gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+            END DO
+
+          ENDIF
+
+          DO jk = 1, jpkm1
+             esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
+             esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
+          END DO
+          esigw3(ji,jj,1  ) = 2.0_wp * (gsigt3(ji,jj,1  ) - gsigw3(ji,jj,1  ))
+          esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
+
+          ! Coefficients for vertical depth as the sum of e3w scale factors
+          gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
+          DO jk = 2, jpk
+             gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
+          END DO
+
+          DO jk = 1, jpk
+             gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigt3(ji,jj,jk)
+             gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsigw3(ji,jj,jk)
+             gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*gsi3w3(ji,jj,jk)
+          END DO
+
+        ENDDO   ! for all jj's
+      ENDDO    ! for all ji's
+
+      DO ji=1,jpi
+        DO jj=1,jpj
+
+          DO jk = 1, jpk
+                esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+                esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) +  &
+                                      hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+                esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+
+             e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*esigt3(ji,jj,jk)
+             e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*esigtu3(ji,jj,jk)
+             e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*esigtv3(ji,jj,jk)
+             e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*esigtf3(ji,jj,jk)
+             !
+             e3w(ji,jj,jk)=hbatt(ji,jj)*esigw3(ji,jj,jk)
+             e3uw(ji,jj,jk)=hbatu(ji,jj)*esigwu3(ji,jj,jk)
+             e3vw(ji,jj,jk)=hbatv(ji,jj)*esigwv3(ji,jj,jk)
+          END DO
+
+        ENDDO
+      ENDDO
+
+      CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3                                      )
+      CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 )
+
+   END SUBROUTINE s_sf12
+
+   SUBROUTINE s_tanh()
+
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE s_tanh*** 
+      !!                     
+      !! ** Purpose :   stretch the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate stretch 
+      !!
+      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
+      !!----------------------------------------------------------------------
+
+      INTEGER  ::   ji, jj, jk           ! dummy loop argument
+      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+
+      DO jk = 1, jpk
+        gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
+        gsigt(jk) = -fssig( REAL(jk,wp)        )
+      END DO
+      IF( nprint == 1 .AND. lwp )   WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
+      !
+      ! Coefficients for vertical scale factors at w-, t- levels
+!!gm bug :  define it from analytical function, not like juste bellow....
+!!gm        or betteroffer the 2 possibilities....
+      DO jk = 1, jpkm1
+         esigt(jk  ) = gsigw(jk+1) - gsigw(jk)
+         esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
+      END DO
+      esigw( 1 ) = 2._wp * ( gsigt(1  ) - gsigw(1  ) ) 
+      esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) )
+      !
+      ! Coefficients for vertical depth as the sum of e3w scale factors
+      gsi3w(1) = 0.5_wp * esigw(1)
+      DO jk = 2, jpk
+         gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
+      END DO
+!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
+      DO jk = 1, jpk
+         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
+         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
+         gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw )
+         gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft )
+      END DO
+!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 1, jpk
+              e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+              e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
+              !
+              e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
+              e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
+              e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
+            END DO
+         END DO
+      END DO
+   END SUBROUTINE s_tanh
+
+   FUNCTION fssig( pk ) RESULT( pf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fssig ***
+      !!       
+      !! ** Purpose :   provide the analytical function in s-coordinate
+      !!          
+      !! ** Method  :   the function provide the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
+      REAL(wp)             ::   pf   ! sigma value
+      !!----------------------------------------------------------------------
+      !
+      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb )  )   &
+         &     - TANH( rn_thetb * rn_theta                                )  )   &
+         & * (   COSH( rn_theta                           )                      &
+         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
+         & / ( 2._wp * SINH( rn_theta ) )
+      !
+   END FUNCTION fssig
+
+
+   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fssig1 ***
+      !!
+      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
+      !!
+      !! ** Method  :   the function provides the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
+      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
+      REAL(wp)             ::   pf1   ! sigma value
+      !!----------------------------------------------------------------------
+      !
+      IF ( rn_theta == 0 ) then      ! uniform sigma
+         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
+      ELSE                        ! stretched sigma
+         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta )              &
+            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
+            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
+      ENDIF
+      !
+   END FUNCTION fssig1
+
+
+   FUNCTION fgamma( pk1, Zbb, Zss, F ) RESULT( gam )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE fgamma  ***
+      !!
+      !! ** Purpose :   provide analytical function for the s-coordinate
+      !!
+      !! ** Method  :   the function provides the non-dimensional position of
+      !!                T and W (i.e. between 0 and 1)
+      !!                T-points at integer values (between 1 and jpk)
+      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
+      !!
+      !!                This method allows the maintenance of fixed surface and or
+      !!                bottom cell resolutions (cf. geopotential coordinates) 
+      !!                within an analytically derived stretched S-coordinate framework.
+      !!
+      !! Reference  :   Siddorn and Furner, in prep
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pk1(jpk)   ! continuous "k" coordinate
+      REAL(wp)                ::   gam(jpk)   ! stretched coordinate
+      REAL(wp), INTENT(in   ) ::   Zbb      ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   Zss      ! surface box depth
+      REAL(wp), INTENT(in   ) ::   F        ! Smoothing parameter
+      REAL(wp)                ::   a1,a2,a3 ! local variables
+      REAL(wp)                ::   n1,n2    ! local variables
+      REAL(wp)                ::   A,B,X    ! local variables
+      integer                 ::   jk
+      !!----------------------------------------------------------------------
+      !
+
+      n1  =  1./(jpk-1.)
+      n2  =  1. -  n1
+
+      a1 = (rn_alpha+2.0_wp)*n1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n1**(rn_alpha+2.0_wp) 
+      a2 = (rn_alpha+2.0_wp)*n2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*n2**(rn_alpha+2.0_wp)
+      a3 = ( n2**3.0_wp - a2)/( n1**3.0_wp - a1)
+     
+      A = Zbb - a3*(Zss-a1)-a2
+      A = A/( n2-0.5_wp*(a2+n2**2.0_wp) - a3*(n1-0.5_wp*(a1+n1**2.0_wp) ) )
+      B = (Zss - a1 - A*( n1-0.5_wp*(a1+n1**2.0_wp ) ) ) / (n1**3.0_wp - a1)
+      X = 1.0_wp-A/2.0_wp-B
+ 
+      DO jk = 1, jpk
+        gam(jk) = A*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+B*pk1(jk)**3.0_wp + X*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
+        gam(jk) = gam(jk)*F+pk1(jk)*(1.0_wp-F)
+      ENDDO 
+
+      !
+   END FUNCTION fgamma
 
    !!======================================================================
Index: branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90
===================================================================
--- branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90	(revision 3597)
+++ branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/par_AMM_12km.h90	(revision 3600)
@@ -19,5 +19,5 @@
       jpidta  = 198,        &  !: first horizontal dimension > or = to jpi
       jpjdta  = 224,        &  !: second                     > or = to jpj
-      jpkdta  = 33,         &  !: number of levels           > or = to jpk
+      jpkdta  = 51,         &  !: number of levels           > or = to jpk
       ! total domain matrix size
       jpiglo  = jpidta,      &  !: first  dimension of global domain --> i
