Index: /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90
===================================================================
--- /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 5255)
+++ /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90	(revision 5255)
@@ -0,0 +1,837 @@
+PROGRAM SCOORD_GEN
+      !!----------------------------------------------------------------------
+      !!                  ***  PROGRAM scoord_gen  ***
+      !!                     
+      !! ** Purpose :   define the s-coordinate system
+      !!
+      !! ** Method  :   s-coordinate
+      !!         The depth of model levels is defined as the product of an
+      !!      analytical function by the local bathymetry, while the vertical
+      !!      scale factors are defined as the product of the first derivative
+      !!      of the analytical function by the bathymetry.
+      !!      (this solution save memory as depth and scale factors are not
+      !!      3d fields)
+      !!          - Read bathymetry (in meters) at t-point and compute the
+      !!         bathymetry at u-, v-, and f-points.
+      !!            hbatu = mi( hbatt )
+      !!            hbatv = mj( hbatt )
+      !!            hbatf = mi( mj( hbatt ) )
+      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
+      !!         function and its derivative given as function.
+      !!            z_gsigt(k) = fssig (k    )
+      !!            z_gsigw(k) = fssig (k-0.5)
+      !!            z_esigt(k) = fsdsig(k    )
+      !!            z_esigw(k) = fsdsig(k-0.5)
+      !!      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 accuracy
+      !!      schemes.
+      !!
+      !!      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        
+      !!
+      !!----------------------------------------------------------------------
+      !
+      IMPLICIT NONE
+
+      USE utils
+
+     !!----------------------------------------------------------------------
+      !
+      !
+      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
+      READ( numnam, namzgr_sco )
+      CLOSE( numnam )
+      jpk = rn_jpk
+
+      WRITE(numout,*)
+      WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
+      WRITE(numout,*) '~~~~~~~~~~~'
+      WRITE(numout,*) '   Namelist namzgr_sco'
+      WRITE(numout,*) '     stretching coeffs '
+      WRITE(numout,*) '        Number of levels                             rn_jpk        = ',rn_jpk
+      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,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
+      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'
+
+      
+      ! Read in bathy, jpi and jpj from bathy.nc
+      CALL read_bathy()
+
+      !Allocate all other allocatable variables
+      ios = dom_oce_alloc()
+      IF( ios .NE. 0) THEN
+         WRITE(0,*) 'Unable to allocate all arrays'
+         STOP 1
+      ENDIF
+
+      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
+      hifu(:,:) = rn_sbot_min
+      hifv(:,:) = rn_sbot_min
+      hiff(:,:) = rn_sbot_min
+
+      !                                        ! set maximum ocean depth
+      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
+         END DO
+      END DO
+      !                                        ! =============================
+      !                                        ! Define the envelop bathymetry   (hbatt)
+      !                                        ! =============================
+      ! use r-value to create hybrid coordinates
+      zenv(:,:) = bathy(:,:)
+      !
+      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+           IF( bathy(ji,jj) == 0._wp ) THEN
+             iip1 = MIN( ji+1, jpi )
+             ijp1 = MIN( jj+1, jpj )
+             iim1 = MAX( ji-1, 1 )
+             ijm1 = MAX( jj-1, 1 )
+             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
+        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
+               zenv(ji,jj) = rn_sbot_min
+             ENDIF
+           ENDIF
+         END DO
+      END DO
+      ! 
+      ! smooth the bathymetry (if required)
+      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
+      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
+      !
+      jl = 0
+      zrmax = 1._wp
+      !   
+      !     
+      ! set scaling factor used in reducing vertical gradients
+      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
+      !
+      ! initialise temporary evelope depth arrays
+      ztmpi1(:,:) = zenv(:,:)
+      ztmpi2(:,:) = zenv(:,:)
+      ztmpj1(:,:) = zenv(:,:)
+      ztmpj2(:,:) = zenv(:,:)
+      !
+      ! initialise temporary r-value arrays
+      zri(:,:) = 1._wp
+      zrj(:,:) = 1._wp
+      !                                                            ! ================ !
+      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
+         !                                                         ! ================ !
+         jl = jl + 1
+         zrmax = 0._wp
+         ! we set zrmax from previous r-values (zri and zrj) first
+         ! if set after current r-value calculation (as previously)
+         ! we could exit DO WHILE prematurely before checking r-value
+         ! of current zenv
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
+            END DO
+         END DO
+         zri(:,:) = 0._wp
+         zrj(:,:) = 0._wp
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
+               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
+               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
+                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
+               END IF
+               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
+                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
+               END IF
+               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
+               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
+               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
+               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
+            END DO
+         END DO
+         !
+         WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
+         !
+         DO jj = 1, nlcj
+            DO ji = 1, nlci
+               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
+            END DO
+         END DO
+         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
+         !                                                  ! ================ !
+      END DO                                                !     End loop     !
+      !                                                     ! ================ !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
+         END DO
+      END DO
+      !
+      ! Envelope bathymetry saved in hbatt
+      hbatt(:,:) = zenv(:,:) 
+      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
+         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
+               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
+            END DO
+         END DO
+      ENDIF
+      !
+      !                                        ! ==============================
+      !                                        !   hbatu, hbatv, hbatf fields
+      !                                        ! ==============================
+      WRITE(numout,*)
+      WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+      hbatu(:,:) = rn_sbot_min
+      hbatv(:,:) = rn_sbot_min
+      hbatf(:,:) = rn_sbot_min
+      DO jj = 1, jpjm1
+        DO ji = 1, jpim1   ! NO vector opt.
+           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
+           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
+           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
+              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
+        END DO
+      END DO
+      ! 
+      zhbat(:,:) = hbatu(:,:)   
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatu(ji,jj) == 0._wp ) THEN
+               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      zhbat(:,:) = hbatv(:,:)  
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatv(ji,jj) == 0._wp ) THEN
+               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+      zhbat(:,:) = hbatf(:,:)   
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( hbatf(ji,jj) == 0._wp ) THEN
+               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
+               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
+            ENDIF
+         END DO
+      END DO
+
+!!bug:  key_helsinki a verifer
+      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+
+         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
+            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
+         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
+            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
+         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
+            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
+         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
+            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
+!! helsinki
+
+      !                                            ! =======================
+      !                                            !   s-coordinate fields     (gdep., e3.)
+      !                                            ! =======================
+      !
+      ! non-dimensional "sigma" for model level depth at w- and t-levels
+
+!========================================================================
+! 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 
+
+      !
+      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
+      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
+      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
+      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
+      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
+      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
+      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+
+#if defined key_agrif
+      ! Ensure meaningful vertical scale factors in ghost lines/columns
+      ! TODO: How can we deal with this in offline tool? 
+      IF( .NOT. Agrif_Root() ) THEN
+         !  
+         IF((nbondi == -1).OR.(nbondi == 2)) THEN
+            e3u_0(1,:,:) = e3u_0(2,:,:)
+         ENDIF
+         !
+         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
+            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
+         ENDIF
+         !
+         IF((nbondj == -1).OR.(nbondj == 2)) THEN
+            e3v_0(:,1,:) = e3v_0(:,2,:)
+         ENDIF
+         !
+         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
+            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
+         ENDIF
+         !
+      ENDIF
+#endif
+!!
+      ! HYBRID : 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            DO jk = 1, jpk-1
+               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
+            END DO
+            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
+         END DO
+      END DO
+      WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
+
+      ! min max values over the domain
+      WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
+      WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
+         &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
+      WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
+         &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
+         &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
+         &                          ' w ', MINVAL( e3w_0  (:,:,:) )
+
+      WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
+         &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
+      WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
+         &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
+         &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
+         &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
+      
+         ! selected vertical profiles
+      WRITE(numout,*)
+      WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
+      WRITE(numout,*) ' ~~~~~~  --------------------'
+      WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+      WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
+         &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
+      DO jj = mj0(20), mj1(20)
+         DO ji = mi0(20), mi1(20)
+            WRITE(numout,*)
+            WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
+            WRITE(numout,*) ' ~~~~~~  --------------------'
+            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
+         END DO
+      END DO
+      DO jj = mj0(74), mj1(74)
+         DO ji = mi0(100), mi1(100)
+            WRITE(numout,*)
+            WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
+            WRITE(numout,*) ' ~~~~~~  --------------------'
+            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
+            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
+               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
+         END DO
+      END DO
+
+!================================================================================
+! check the coordinate makes sense
+!================================================================================
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+
+            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(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 STOP 1
+                 ENDIF
+                 ! and check it has never gone negative
+                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
+                    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 STOP 1
+                 ENDIF
+                 ! and check it never exceeds the total depth
+                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
+                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
+                    CALL STOP 1
+                 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(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
+                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
+                    CALL STOP 1
+                 ENDIF
+               END DO
+
+            ENDIF
+
+         END DO
+      END DO
+      !
+      !
+
+!!======================================================================
+   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), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+
+      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
+      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
+      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
+      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
+
+      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
+      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
+      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
+      z_esigwu3 = 0._wp   ;   z_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
+                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
+                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
+               END DO
+            ELSE ! shallow water, uniform sigma
+               DO jk = 1, jpk
+                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
+                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
+                  END DO
+            ENDIF
+            !
+            DO jk = 1, jpk-1
+               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
+               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
+            END DO
+            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
+            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
+            !
+            ! Coefficients for vertical depth as the sum of e3w scale factors
+            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
+            DO jk = 2, jpk
+               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
+            END DO
+            !
+            DO jk = 1, jpk
+               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp)
+               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
+               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
+               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
+               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_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
+               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
+                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
+                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !
+               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               !
+               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
+            END DO
+        END DO
+      END DO
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
+      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
+
+   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) ::   zsmth               ! smoothing around critical depth
+      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
+      !
+      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
+      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
+      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
+      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
+
+      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
+      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
+      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
+      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
+
+      DO ji = 1, jpi
+         DO jj = 1, jpj
+
+          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
+              
+              zzb = 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
+              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
+            
+              zzs = rn_zs / hbatt(ji,jj) 
+              
+              IF (rn_efold /= 0.0_wp) THEN
+                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
+              ELSE
+                zsmth = 1.0_wp 
+              ENDIF
+               
+              DO jk = 1, jpk
+                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
+                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
+              ENDDO
+              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
+              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
+ 
+          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
+
+            DO jk = 1, jpk
+              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
+              z_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
+              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
+              z_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, jpk-1
+             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
+             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
+          END DO
+          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
+          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
+
+          ! Coefficients for vertical depth as the sum of e3w scale factors
+          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
+          DO jk = 2, jpk
+             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
+          END DO
+
+          DO jk = 1, jpk
+             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
+             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
+             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
+          END DO
+
+        ENDDO   ! for all jj's
+      ENDDO    ! for all ji's
+
+      DO ji=1,jpi-1
+        DO jj=1,jpj-1
+
+          DO jk = 1, jpk
+                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
+                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
+                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+
+             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
+             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
+             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
+             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
+             !
+             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
+             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
+             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
+          END DO
+
+        ENDDO
+      ENDDO
+      !
+      !                                               ! =============
+
+      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
+      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_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
+
+      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
+      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
+
+      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
+      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
+
+      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
+      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
+
+      DO jk = 1, jpk
+        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
+        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
+      END DO
+      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_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, jpk-1
+         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
+         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
+      END DO
+      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
+      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
+      !
+      ! Coefficients for vertical depth as the sum of e3w scale factors
+      z_gsi3w(1) = 0.5_wp * z_esigw(1)
+      DO jk = 2, jpk
+         z_gsi3w(jk) = z_gsi3w(jk-1) + z_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(jpk-1,wp)
+         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
+         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
+         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
+         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_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_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
+              !
+              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
+              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
+              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
+            END DO
+         END DO
+      END DO
+
+      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
+      DEALLOCATE( z_esigt, z_esigw )
+
+   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(jpk-1) + 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( jpk-1 )
+      ELSE                        ! stretched sigma
+         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
+            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpk-1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
+            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
+      ENDIF
+      !
+   END FUNCTION fssig1
+
+
+   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
+      !!----------------------------------------------------------------------
+      !!                 ***  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)                ::   p_gamma(jpk)   ! stretched coordinate
+      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
+      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
+      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
+      REAL(wp)                ::   za1,za2,za3    ! local variables
+      REAL(wp)                ::   zn1,zn2        ! local variables
+      REAL(wp)                ::   za,zb,zx       ! local variables
+      integer                 ::   jk
+      !!----------------------------------------------------------------------
+      !
+
+      zn1  =  1./(jpk-1.)
+      zn2  =  1. -  zn1
+
+      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
+      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
+      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
+     
+      za = pzb - za3*(pzs-za1)-za2
+      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
+      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
+      zx = 1.0_wp-za/2.0_wp-zb
+ 
+      DO jk = 1, jpk
+        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
+                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
+                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
+        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
+      ENDDO 
+
+      !
+   END FUNCTION fgamma
+
+END PROGRAM SCOORD_GEN
Index: /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90
===================================================================
--- /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 5255)
+++ /branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90	(revision 5255)
@@ -0,0 +1,137 @@
+MODULE utils
+
+   IMPLICIT NONE
+   USE netcdf
+
+   PUBLIC             ! allows the acces to par_oce when dom_oce is used
+   !                  ! exception to coding rules... to be suppressed ???
+
+   PUBLIC dom_oce_alloc
+   PUBLIC read_bathy
+
+
+   !! All coordinates
+   !! ---------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
+
+   !! s-coordinate and hybrid z-s-coordinate
+   !! =----------------======---------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gsigt, gsigw       !: model level depth coefficient at t-, w-levels (analytic)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   gsi3w              !: model level depth coefficient at w-level (sum of gsigw)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   esigt, esigw       !: vertical scale factor coef. at t-, w-levels
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf      !: ocean depth at the vertical of  v--f
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu      !:                                 t--u points (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot     !: ocean surface and bottom topographies 
+   !                                                                           !  (if deviating from coordinate surfaces in HYBRID)
+   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
+
+   !!----------------------------------------------------------------------
+   !! masks, bathymetry
+   !! ---------------------------------------------------------------------
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters - read from file)
+
+   !! Other variables needed by scoord_gen
+   INTEGER  ::   jpi, jpj, jpk            ! Size of the domain - read from bathy or namelist?
+   INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
+   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
+   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
+   REAL(wp) ::   zrmax, ztaper   ! temporary scalars
+   REAL(wp) ::   zrfact
+   !
+   REAL(wp), DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
+   REAL(wp), DIMENSION(:,:  ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
+
+   NAMELIST/namzgr_sco/rn_jpk, 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
+
+
+   INTEGER FUNCTION dom_oce_alloc()
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(12) :: ierr
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+      !
+      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), &
+         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
+         !
+      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         &
+         &      gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) ,                         &
+         &      gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) )
+         !
+         !
+         !
+      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
+         &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     &
+         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
+         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) )
+
+      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(9) )
+      !
+      dom_oce_alloc = MAXVAL(ierr)
+      !
+   END FUNCTION dom_oce_alloc
+
+
+   SUBROUTINE read_bathy()
+     !! Read bathymetry from input netcdf file
+     INTEGER :: var_id
+
+     CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening mesh_mask file' )
+
+     ! Find the size of the input bathymetry
+     CALL dimlen(ncin, 'x', jpi)    
+     CALL dimlen(ncin, 'y', jpj)    
+     
+     ALLOCATE( bathy(jpi, jpj) )
+     
+     ! Read the bathymetry variable from file
+     CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', tmask_id ), 'Cannot get variable ID for bathymetry')
+     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) )
+
+   END SUBROUTINE read_bathy
+
+   SUBROUTINE dimlen( ncid, dimname, len )
+     ! Determine the length of dimension dimname
+     INTEGER, INTENT(in)          :: ncid
+     CHARACTER(LEN=*), INTENT(in) :: dimname
+     INTEGER, INTENT(out)         :: len
+     ! Local variables
+     INTEGER :: id_var, istatus
+
+     id_var = 1
+     CALL check_nf90( nf90_inq_dimid(ncid, dimname, id_var), 'Dimension not found in file')
+     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))
+
+  END SUBROUTINE dimlen
+
+
+   
+   SUBROUTINE write_coord_file()
+
+   END SUBROUTINE write_coord_file
+
+   SUBROUTINE check_nf90( istat, message )
+      !Check for netcdf errors
+      INTEGER, INTENT(in) :: istat
+      CHARACTER(LEN=*), INTENT(in), OPTIONAL :: message
+
+      IF (istat /= nf90_noerr) THEN
+         WRITE(numerr,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
+         IF ( PRESENT(message) ) THEN ; WRITE(numerr,*) message ; ENDIF
+         STOP
+      ENDIF
+
+   END SUBROUTINE check_nf90
+
+
+END MODULE utils
