Index: trunk/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1083)
+++ trunk/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 1099)
@@ -4,7 +4,17 @@
    !! Ocean initialization : domain initialization
    !!==============================================================================
-
+   !! History :  OPA  !  1995-12  (G. Madec)  Original code : s vertical coordinate
+   !!                 !  1997-07  (G. Madec)  lbc_lnk call
+   !!                 !  1997-04  (J.-O. Beismann) 
+   !!            8.5  !  2002-09 (A. Bozec, G. Madec)  F90: Free form and module
+   !!             -   !  2002-09 (A. de Miranda)  rigid-lid + islands
+   !!  NEMO      1.0  !  2003-08  (G. Madec)  F90: Free form and module
+   !!             -   !  2005-10  (A. Beckmann)  modifications for hybrid s-ccordinates & new stretching function
+   !!            2.0  !  2006-04  (R. Benshila, G. Madec)  add zgr_zco
+   !!            3.0  !  2008-06  (G. Madec)  insertion of domzgr_zps.h90 & conding style
    !!----------------------------------------------------------------------
-   !!   dom_zgr     : defined the ocean vertical coordinate system
+
+   !!----------------------------------------------------------------------
+   !!   dom_zgr          : defined the ocean vertical coordinate system
    !!       zgr_bat      : bathymetry fields (levels and meters)
    !!       zgr_bat_zoom : modify the bathymetry field if zoom domain
@@ -15,6 +25,6 @@
    !!       zgr_sco      : s-coordinate
    !!       fssig        : sigma coordinate non-dimensional function
+   !!       dfssig       : derivative of the sigma coordinate function    !!gm  (currently missing!)
    !!---------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -22,6 +32,6 @@
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE closea
-   USE solisl
+   USE closea          ! closed seas
+   USE solisl          ! solver - island in rigid-lid
    USE c1d             ! 1D configuration
 
@@ -29,16 +39,13 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dom_zgr        ! called by dom_init.F90
-
-   !! * Module variables
-      REAL(wp) ::          & !!: Namelist nam_zgr_sco
-      sbot_min =  300.  ,  &  !: minimum depth of s-bottom surface (>0) (m)
-      sbot_max = 5250.  ,  &  !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
-      theta    =    6.0 ,  &  !: surface control parameter (0<=theta<=20)
-      thetb    =    0.75,  &  !: bottom control parameter  (0<=thetb<= 1)
-      r_max    =    0.15      !: maximum cut-off r-value allowed (0<r_max<1)
-
-
+   PUBLIC   dom_zgr    ! called by dom_init.F90
+
+!!gm   DOCTOR name for the namelist parameter...
+   !                                 !!! ** Namelist nam_zgr_sco **
+   REAL(wp) ::   sbot_min =  300.     ! minimum depth of s-bottom surface (>0) (m)
+   REAL(wp) ::   sbot_max = 5250.     ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
+   REAL(wp) ::   theta    =    6.0    ! surface control parameter (0<=theta<=20)
+   REAL(wp) ::   thetb    =    0.75   ! bottom control parameter  (0<=thetb<= 1)
+   REAL(wp) ::   r_max    =    0.15   ! maximum cut-off r-value allowed (0<r_max<1)
  
    !! * Substitutions
@@ -46,5 +53,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! $Id:$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -58,26 +67,23 @@
       !!      vertical scale factors.
       !!
-      !! ** Method  reference vertical coordinate
-      !!
-      !! ** Action : 
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for hybrid s-ccordinates
-      !!----------------------------------------------------------------------
-      INTEGER ::   ioptio = 0      ! temporary integer
-
+      !! ** Method  : - reference 1D vertical coordinate (gdep._0, e3._0)
+      !!              - read/set ocean depth and ocean levels (bathy, mbathy)
+      !!              - vertical coordinate (gdep., e3.) depending on the 
+      !!                coordinate chosen :
+      !!                   ln_zco=T   z-coordinate   (forced if lk_zco)
+      !!                   ln_zps=T   z-coordinate with partial steps
+      !!                   ln_zco=T   s-coordinate 
+      !!
+      !! ** Action  :   define gdep., e3., mbathy and bathy
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio = 0   ! temporary integer
+      !!
       NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
       !!----------------------------------------------------------------------
 
-      ! Read Namelist nam_zgr : vertical coordinate'
-      ! ---------------------
-      REWIND ( numnam )
+      REWIND ( numnam )                ! Read Namelist nam_zgr : vertical coordinate'
       READ   ( numnam, nam_zgr )
 
-      ! Parameter control and print
-      ! ---------------------------
-      ! Control print
-      IF(lwp) THEN
+      IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'dom_zgr : vertical coordinate'
@@ -89,48 +95,22 @@
       ENDIF
 
-      ! Check Vertical coordinate options
-      ioptio = 0
+      ioptio = 0                       ! Check Vertical coordinate options
       IF( ln_zco ) ioptio = ioptio + 1
       IF( ln_zps ) ioptio = ioptio + 1
       IF( ln_sco ) ioptio = ioptio + 1
-      IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
-
+      IF ( ioptio /= 1 )   CALL ctl_stop( ' none or several vertical coordinate options used' )
       IF( lk_zco ) THEN
           IF(lwp) WRITE(numout,*) '          z-coordinate with reduced incore memory requirement'
-          IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
+          IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
       ENDIF
 
       ! Build the vertical coordinate system
       ! ------------------------------------
-
                      CALL zgr_z        ! Reference z-coordinate system (always called)
-
                      CALL zgr_bat      ! Bathymetry fields (levels and meters)
-
       IF( ln_zco )   CALL zgr_zco      ! z-coordinate
-
       IF( ln_zps )   CALL zgr_zps      ! Partial step z-coordinate
-
       IF( ln_sco )   CALL zgr_sco      ! s-coordinate or hybrid z-s coordinate
-
-!!bug gm control print:
-      IF( nprint == 1 .AND. lwp )   THEN
-         WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
-         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
-            &                   ' w ',   MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ),  &
-            &                   ' u ',   MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ),  &
-            &                   ' uw',   MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)),   &
-            &                   ' w ',   MINVAL( fse3w(:,:,:) )
-
-         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
-            &                   ' w ',   MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
-         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ),  &
-            &                   ' u ',   MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ),  &
-            &                   ' uw',   MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)),   &
-            &                   ' w ',   MAXVAL( fse3w(:,:,:) )
-      ENDIF
-!!!bug gm
-
+      !
    END SUBROUTINE dom_zgr
 
@@ -153,18 +133,12 @@
       !!
       !! ** Action  : - gdept_0, gdepw_0 : depth of T- and W-point (m)
-      !!              -  e3t_0, e3w_0    : scale factors at T- and W-levels (m)
-      !!
-      !! Reference :
-      !!      Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!              - e3t_0  , e3w_0   : scale factors at T- and W-levels (m)
+      !!
+      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
+      !!----------------------------------------------------------------------
       INTEGER  ::   jk                     ! dummy loop indices
       REAL(wp) ::   zt, zw                 ! temporary scalars
-      REAL(wp) ::   &
-      zsur , za0, za1, zkth, zacr,      &  ! Values set from parameters in
-      zdzmin, zhmax                        ! par_CONFIG_Rxx.h90
+      REAL(wp) ::   zsur, za0, za1, zkth   ! Values set from parameters in
+      REAL(wp) ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
       !!----------------------------------------------------------------------
 
@@ -177,23 +151,18 @@
       !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
       !
-       IF(  ppa1  == pp_to_be_computed  .AND.  &
+      IF(   ppa1  == pp_to_be_computed  .AND.  &
          &  ppa0  == pp_to_be_computed  .AND.  &
          &  ppsur == pp_to_be_computed           ) THEN
-         za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) )          &
-             / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) &
-             &                         *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
-             &                             - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
-
-         za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
-         zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
-
-       ELSE
+         !
+         za1  = (  ppdzmin - pphmax / FLOAT(jpkm1)  )                                                      &
+            & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
+            &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
+         za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
+         zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
+      ELSE
          za1 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
-       ENDIF
-
-
-      ! Parameter print
-      ! ---------------
-      IF(lwp) THEN
+      ENDIF
+
+      IF(lwp) THEN                         ! Parameter print
          WRITE(numout,*)
          WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
@@ -222,5 +191,4 @@
       ! ======================
       IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid       
-
          za1 = zhmax / FLOAT(jpk-1) 
          DO jk = 1, jpk
@@ -232,23 +200,17 @@
             e3t_0  (jk) =  za1
          END DO
-
-      ELSE
-
+      ELSE                                ! Madec & Imbard 1996 function
          DO jk = 1, jpk
             zw = FLOAT( jk )
             zt = FLOAT( jk ) + 0.5
-            gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) )  )
-            gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) )  )
-            e3w_0  (jk) =          za0      + za1        * TANH(      (zw-zkth)/zacr   )
-            e3t_0  (jk) =          za0      + za1        * TANH(      (zt-zkth)/zacr   )
-         END DO
-         gdepw_0(1) = 0.e0      ! force first w-level to be exactly at zero
-
-      ENDIF
-
-      ! Control and  print
-      ! ==================
-
-      IF(lwp) THEN
+            gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
+            gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
+            e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
+            e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
+         END DO
+         gdepw_0(1) = 0.e0                     ! force first w-level to be exactly at zero
+      ENDIF
+
+      IF(lwp) THEN                        ! control print
          WRITE(numout,*)
          WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
@@ -256,10 +218,9 @@
          WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
       ENDIF
-
-      DO jk = 1, jpk
-         IF( e3w_0(jk)  <= 0. .OR. e3t_0(jk)  <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
-         IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
-      END DO
-
+      DO jk = 1, jpk                      ! control positivity
+         IF( e3w_0  (jk) <= 0.e0 .OR. e3t_0  (jk) <= 0.e0 )   CALL ctl_stop( ' e3w or e3t =< 0 '    )
+         IF( gdepw_0(jk) <  0.e0 .OR. gdept_0(jk) <  0.e0 )   CALL ctl_stop( ' gdepw or gdept < 0 ' )
+      END DO
+      !
    END SUBROUTINE zgr_z
 
@@ -298,24 +259,14 @@
       !! ** Action  : - mbathy: level bathymetry (in level index)
       !!              - bathy : meter bathymetry (in meters)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE iom
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jl, jk       ! dummy loop indices
-      INTEGER ::   inum                 ! temporary logical unit
-      INTEGER  ::   &
-         ii_bump, ij_bump, ih           ! bump center position
-      INTEGER , DIMENSION(jpidta,jpjdta) ::   &
-         idta                           ! global domain integer data
-      REAL(wp) ::   &
-         r_bump, h_bump, h_oce,      &  ! bump characteristics 
-         zi, zj, zh                     ! temporary scalars
-      REAL(wp), DIMENSION(jpidta,jpjdta) ::   &
-         zdta                           ! global domain scalar data
+      !!
+      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices
+      INTEGER  ::   inum                      ! temporary logical unit
+      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
+      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics 
+      REAL(wp) ::   zi     , zj     , zh      ! temporary scalars
+      INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data
+      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data
       !!----------------------------------------------------------------------
 
@@ -324,29 +275,23 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~'
 
-      ! ========================================
-      ! global domain level and meter bathymetry (idta,zdta)
-      ! ========================================
-      !                                               ! =============== ! 
-      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          ! defined by hand !
-         !                                            ! =============== !
-
+      !                                               ! ================== ! 
+      IF( ntopo == 0 .OR. ntopo == -1 ) THEN          !   defined by hand  !
+         !                                            ! ================== !
+         !                                            ! global domain level and meter bathymetry (idta,zdta)
+         !
          IF( ntopo == 0 ) THEN                        ! flat basin
-
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'
-
-            idta(:,:) = jpkm1                            ! flat basin 
-            zdta(:,:) = gdepw_0(jpk)
+            idta(:,:) = jpkm1                            ! before last level
+            zdta(:,:) = gdepw_0(jpk)                     ! last w-point depth
             h_oce     = gdepw_0(jpk)
-
-         ELSE                                         ! bump
+         ELSE                                         ! bump centered in the basin
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin with a bump'
-
-            ii_bump = jpidta / 2           ! i-index of the bump center
-            ij_bump = jpjdta / 2           ! j-index of the bump center
-            r_bump  = 50000.e0             ! bump radius (meters)       
-            h_bump  = 2700.e0              ! bump height (meters)
-            h_oce   = gdepw_0(jpk)         ! background ocean depth (meters)
+            ii_bump = jpidta / 2                           ! i-index of the bump center
+            ij_bump = jpjdta / 2                           ! j-index of the bump center
+            r_bump  = 50000.e0                             ! bump radius (meters)       
+            h_bump  = 2700.e0                              ! bump height (meters)
+            h_oce   = gdepw_0(jpk)                         ! background ocean depth (meters)
             IF(lwp) WRITE(numout,*) '            bump characteristics: '
             IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
@@ -354,6 +299,6 @@
             IF(lwp) WRITE(numout,*) '               bump radius         = ', r_bump , ' index'
             IF(lwp) WRITE(numout,*) '            background ocean depth = ', h_oce  , ' meters'
-            ! zdta :
-            DO jj = 1, jpjdta
+            !                                        
+            DO jj = 1, jpjdta                              ! zdta :
                DO ji = 1, jpidta
                   zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
@@ -362,9 +307,8 @@
                END DO
             END DO
-
-            ! idta :
-            IF( ln_sco ) THEN      ! s-coordinate (zsc       ): idta()=jpk
+            !                                              ! idta :
+            IF( ln_sco ) THEN                                   ! s-coordinate (zsc       ): idta()=jpk
                idta(:,:) = jpkm1
-            ELSE                   ! z-coordinate (zco or zps): step-like topography
+            ELSE                                                ! z-coordinate (zco or zps): step-like topography
                idta(:,:) = jpkm1
                DO jk = 1, jpkm1
@@ -377,6 +321,6 @@
             ENDIF
          ENDIF
-
-         ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
+         !                                            ! set GLOBAL boundary conditions 
+         !                                            ! Caution : idta on the global domain: use of jperio, not nperio
          IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
             idta( :    , 1    ) = -1                ;      zdta( :    , 1    ) = -1.e0
@@ -396,12 +340,8 @@
          ENDIF
 
-         ! =======================================
-         ! local domain level and meter bathymetry (mbathy,bathy)
-         ! =======================================
-         
-         mbathy(:,:) = 0                                 ! set to zero extra halo points
-         bathy (:,:) = 0.e0                              ! (require for mpp case)
-         
-         DO jj = 1, nlcj                                 ! interior values
+         !                                            ! local domain level and meter bathymetries (mbathy,bathy)
+         mbathy(:,:) = 0                                   ! set to zero extra halo points
+         bathy (:,:) = 0.e0                                ! (require for mpp case)
+         DO jj = 1, nlcj                                   ! interior values
             DO ji = 1, nlci
                mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
@@ -409,18 +349,17 @@
             END DO
          END DO
-
-         !                                            ! =============== !
-      ELSEIF( ntopo == 1 ) THEN                       !   read in file  !
-         !                                            ! =============== !
-
-         IF( ln_zco )   THEN
-            CALL iom_open ( 'bathy_level.nc', inum )   ! Level bathymetry
+         !
+         !                                            ! ================ !
+      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
+         !                                            ! ================ !
+         !
+         IF( ln_zco )   THEN                          ! zco : read level bathymetry 
+            CALL iom_open( 'bathy_level.nc', inum )  
             CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
             CALL iom_close (inum)
             mbathy(:,:) = INT( bathy(:,:) )
          ENDIF
-
-         IF( ln_zps .OR. ln_sco )   THEN
-            CALL iom_open ( 'bathy_meter.nc', inum )   ! meter bathymetry
+         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry
+            CALL iom_open( 'bathy_meter.nc', inum ) 
             CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close (inum)
@@ -429,41 +368,32 @@
       ELSE                                            !      error      !
          !                                            ! =============== !
-         WRITE(ctmp1,*) '          parameter , ntopo = ', ntopo
+         WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
          CALL ctl_stop( '    zgr_bat : '//trim(ctmp1) )
       ENDIF
-
-      ! =======================
-      ! NO closed seas or lakes
-      ! =======================
-
-      IF( nclosea == 0 ) THEN
-         DO jl = 1, jpncs
+      !
+      !                                               ! =========================== !
+      IF( nclosea == 0 ) THEN                         !   NO closed seas or lakes   !
+         DO jl = 1, jpncs                             ! =========================== !
             DO jj = ncsj1(jl), ncsj2(jl)
                DO ji = ncsi1(jl), ncsi2(jl)
-                  mbathy(ji,jj) = 0                   ! suppress closed seas
-                  bathy (ji,jj) = 0.e0                ! and lakes
-!!bug gm          bathy (ji,jj) = 300.e0                ! and lakes
+                  mbathy(ji,jj) = 0                   ! suppress closed seas and lakes from bathymetry
+                  bathy (ji,jj) = 0.e0                
                END DO
             END DO
          END DO
       ENDIF
-
 #if defined key_orca_lev10
-      ! 10 time the vertical resolution
-      mbathy(:,:) = 10 * mbathy(:,:)
-      IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"'
+      !                                               ! ================= !
+      mbathy(:,:) = 10 * mbathy(:,:)                  !   key_orca_lev10  !
+      !                                               ! ================= !
+      IF( ln_zps .OR. ln_sco )   CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
 #endif
-      ! ===========
-      ! Zoom domain 
-      ! ===========
-
-      IF( lzoom )   CALL zgr_bat_zoom
-
-      ! ================
-      ! Bathymetry check
-      ! ================
-
-      CALL zgr_bat_ctl
-
+      !                                               ! =============== !
+      IF( lzoom        )   CALL zgr_bat_zoom          !   Zoom domain   !
+      !                                               ! =============== !
+
+      !                                               ! =================== !
+      IF( .NOT. lk_c1d )   CALL zgr_bat_ctl           !   Bathymetry check  !
+      !                                               ! =================== !
    END SUBROUTINE zgr_bat
 
@@ -479,25 +409,21 @@
       !!
       !! ** Action  : - update mbathy: level bathymetry (in level index)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!----------------------------------------------------------------------
-      !! * Local variables
+      !!----------------------------------------------------------------------
       INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integers
       !!----------------------------------------------------------------------
-
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '    zgr_bat_zoom : modify the level bathymetry for zoom domain'
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
-
+      !
       ! Zoom domain
       ! ===========
-
+      !
       ! Forced closed boundary if required
-      IF( lzoom_w )   mbathy( mi0(jpizoom):mi1(jpizoom) , :  ) = 0
-      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
-      IF( lzoom_e )   mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
-      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
-
+      IF( lzoom_s )   mbathy(  : , mj0(jpjzoom):mj1(jpjzoom) )      = 0
+      IF( lzoom_w )   mbathy(      mi0(jpizoom):mi1(jpizoom) , :  ) = 0
+      IF( lzoom_e )   mbathy(      mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , :  ) = 0
+      IF( lzoom_n )   mbathy(  : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) )      = 0
+      !
       ! Configuration specific domain modifications
       ! (here, ORCA arctic configuration: suppress Med Sea)
@@ -508,5 +434,4 @@
             !                                     ! =======================
             IF(lwp) WRITE(numout,*) '                   ORCA R2 arctic zoom: suppress the Med Sea'
-
             ii0 = 141   ;   ii1 = 162      ! Sea box i,j indices
             ij0 =  98   ;   ij1 = 110
@@ -515,8 +440,6 @@
             !                                     ! =======================
             IF(lwp) WRITE(numout,*) '                   ORCA R05 arctic zoom: suppress the Med Sea'
-  
             ii0 = 563   ;   ii1 = 642      ! zero over the Med Sea boxe
             ij0 = 314   ;   ij1 = 370 
-
          END SELECT
          !
@@ -524,5 +447,5 @@
          !
       ENDIF
-
+      !
    END SUBROUTINE zgr_bat_zoom
 
@@ -549,15 +472,8 @@
       !! ** Action  : - update mbathy: level bathymetry (in level index)
       !!              - update bathy : meter bathymetry (in meters)
-      !!
-      !! History :
-      !!   9.0  !  03-08  (G. Madec)  Original code
-      !!   9.0  !  05-10  (A. Beckmann)  modifications for s-ccordinates
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   ji, jj, jl           ! dummy loop indices
-      INTEGER ::   &
-         icompt, ibtest, ikmax          ! temporary integers
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zbathy                         ! temporary workspace
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jl                    ! dummy loop indices
+      INTEGER ::   icompt, ibtest, ikmax         ! temporary integers
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbathy   ! temporary workspace
       !!----------------------------------------------------------------------
 
@@ -566,107 +482,88 @@
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
 
-      ! ================
-      ! Bathymetry check
-      ! ================
-
-      IF( .NOT. lk_c1d )   THEN
-
-         ! Suppress isolated ocean grid points
-
+      !                                          ! Suppress isolated ocean grid points
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
+      IF(lwp) WRITE(numout,*)'                   -----------------------------------'
+      icompt = 0
+      DO jl = 1, 2
+         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
+            mbathy( 1 ,:) = mbathy(jpim1,:)           ! local domain is cyclic east-west
+            mbathy(jpi,:) = mbathy(  2  ,:)
+         ENDIF
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ibtest = MAX(  mbathy(ji-1,jj), mbathy(ji+1,jj),   &
+                  &           mbathy(ji,jj-1), mbathy(ji,jj+1)  )
+               IF( ibtest < mbathy(ji,jj) ) THEN
+                  IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
+                     &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
+                  mbathy(ji,jj) = ibtest
+                  icompt = icompt + 1
+               ENDIF
+            END DO
+         END DO
+      END DO
+      IF( icompt == 0 ) THEN
+         IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
+      ELSE
+         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
+      ENDIF
+      IF( lk_mpp ) THEN
+         zbathy(:,:) = FLOAT( mbathy(:,:) )
+         CALL lbc_lnk( zbathy, 'T', 1. )
+         mbathy(:,:) = INT( zbathy(:,:) )
+      ENDIF
+
+      !                                          ! East-west cyclic boundary conditions
+      IF( nperio == 0 ) THEN
+         IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
+         IF( lk_mpp ) THEN
+            IF( nbondi == -1 .OR. nbondi == 2 ) THEN
+               IF( jperio /= 1 )   mbathy(1,:) = 0
+            ENDIF
+            IF( nbondi == 1 .OR. nbondi == 2 ) THEN
+               IF( jperio /= 1 )   mbathy(nlci,:) = 0
+            ENDIF
+         ELSE
+            IF( ln_zco .OR. ln_zps ) THEN
+               mbathy( 1 ,:) = 0
+               mbathy(jpi,:) = 0
+            ELSE
+               mbathy( 1 ,:) = jpkm1
+               mbathy(jpi,:) = jpkm1
+            ENDIF
+         ENDIF
+      ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
+         IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
+         mbathy( 1 ,:) = mbathy(jpim1,:)
+         mbathy(jpi,:) = mbathy(  2  ,:)
+      ELSEIF( nperio == 2 ) THEN
+         IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: nperio = ', nperio
+      ELSE
+         IF(lwp) WRITE(numout,*) '    e r r o r'
+         IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
+         !         STOP 'dom_mba'
+      ENDIF
+
+      ! Set to zero mbathy over islands if necessary  (lk_isl=F)
+      IF( .NOT. lk_isl ) THEN    ! No island
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'                   suppress isolated ocean grid points'
-         IF(lwp) WRITE(numout,*)'                   -----------------------------------'
-
-         icompt = 0
-
-         DO jl = 1, 2
-
-            IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN
-               mbathy( 1 ,:) = mbathy(jpim1,:)
-               mbathy(jpi,:) = mbathy(  2  ,:)
-            ENDIF
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj),   &
-                     &          mbathy(ji,jj-1),mbathy(ji,jj+1) )
-                  IF( ibtest < mbathy(ji,jj) ) THEN
-                     IF(lwp) WRITE(numout,*) ' the number of ocean level at ',   &
-                        &   'grid-point (i,j) =  ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
-                     mbathy(ji,jj) = ibtest
-                     icompt = icompt + 1
-                  ENDIF
-               END DO
-            END DO
-
-         END DO
-         IF( icompt == 0 ) THEN
-            IF(lwp) WRITE(numout,*)'     no isolated ocean grid points'
-         ELSE
-            IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points suppressed'
-         ENDIF
-         IF( lk_mpp ) THEN
+         IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
+         IF(lwp) WRITE(numout,*) '         ----------------------------'
+         !
+         mbathy(:,:) = MAX( 0, mbathy(:,:) )
+         !
+         !  Boundary condition on mbathy
+         IF( .NOT.lk_mpp ) THEN 
+!!gm        !!bug ???  think about it !
+            !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
             zbathy(:,:) = FLOAT( mbathy(:,:) )
             CALL lbc_lnk( zbathy, 'T', 1. )
             mbathy(:,:) = INT( zbathy(:,:) )
          ENDIF
-
-         ! 3.2 East-west cyclic boundary conditions
-
-         IF( nperio == 0 ) THEN
-            IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west',   &
-               ' boundary: nperio = ', nperio
-            IF( lk_mpp ) THEN
-               IF( nbondi == -1 .OR. nbondi == 2 ) THEN
-                  IF( jperio /= 1 )   mbathy(1,:) = 0
-               ENDIF
-               IF( nbondi == 1 .OR. nbondi == 2 ) THEN
-                  IF( jperio /= 1 )   mbathy(nlci,:) = 0
-               ENDIF
-            ELSE
-               IF( ln_zco .OR. ln_zps ) THEN
-                  mbathy( 1 ,:) = 0
-                  mbathy(jpi,:) = 0
-               ELSE
-                  mbathy( 1 ,:) = jpkm1
-                  mbathy(jpi,:) = jpkm1
-               ENDIF
-            ENDIF
-         ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio ==  6 ) THEN
-            IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions',   &
-               ' on mbathy: nperio = ', nperio
-            mbathy( 1 ,:) = mbathy(jpim1,:)
-            mbathy(jpi,:) = mbathy(  2  ,:)
-         ELSEIF( nperio == 2 ) THEN
-            IF(lwp) WRITE(numout,*) '   equatorial boundary conditions',   &
-               ' on mbathy: nperio = ', nperio
-         ELSE
-            IF(lwp) WRITE(numout,*) '    e r r o r'
-            IF(lwp) WRITE(numout,*) '    parameter , nperio = ', nperio
-            !         STOP 'dom_mba'
-         ENDIF
-
-         ! Set to zero mbathy over islands if necessary  (lk_isl=F)
-         IF( .NOT. lk_isl ) THEN    ! No island
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '         mbathy set to 0 over islands'
-            IF(lwp) WRITE(numout,*) '         ----------------------------'
-
-            mbathy(:,:) = MAX( 0, mbathy(:,:) )
-
-            !  Boundary condition on mbathy
-            IF( .NOT.lk_mpp ) THEN 
-               !!bug ???  y reflechir!
-               !   ... mono- or macro-tasking: T-point, >0, 2D array, no slab
-               zbathy(:,:) = FLOAT( mbathy(:,:) )
-               CALL lbc_lnk( zbathy, 'T', 1. )
-               mbathy(:,:) = INT( zbathy(:,:) )
-            ENDIF
-
-         ENDIF
-
       ENDIF
 
       ! Number of ocean level inferior or equal to jpkm1
-
       ikmax = 0
       DO jj = 1, jpj
@@ -675,6 +572,5 @@
          END DO
       END DO
-      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???
-
+!!gm  !!! test to do:   ikmax = MAX( mbathy(:,:) )   ???
       IF( ikmax > jpkm1 ) THEN
          IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
@@ -685,15 +581,13 @@
       ENDIF
 
+      ! control print
       IF( lwp .AND. nprint == 1 ) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' bathymetric field '
+         WRITE(numout,*) ' bathymetric field :   number of non-zero T-levels '
          WRITE(numout,*) ' ------------------'
-         WRITE(numout,*) ' number of non-zero T-levels '
-         CALL prihin( mbathy, jpi, jpj, 1, jpi,   &
-                      1     , 1  , jpj, 1, 3  ,   &
-                      numout )
-         WRITE(numout,*)
-      ENDIF
-
+         CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
+         WRITE(numout,*)
+      ENDIF
+      !
    END SUBROUTINE zgr_bat_ctl
 
@@ -706,11 +600,8 @@
       !!
       !! ** Method  :   set 3D coord. arrays to reference 1D array if lk_zco=T
-      !!
-      !! History :
-      !!        !  06-04  (R. Benshila, G. Madec)  Original code 
       !!----------------------------------------------------------------------
       INTEGER  ::   jk
       !!----------------------------------------------------------------------
-
+      !
       IF( .NOT.lk_zco ) THEN
          DO jk = 1, jpk
@@ -727,5 +618,5 @@
          END DO
       ENDIF
-
+      !
    END SUBROUTINE zgr_zco
 
@@ -787,46 +678,15 @@
       !!         - - - - - - -   gdept, gdepw and e3. are positives
       !!      
-      !!  Reference :
-      !!     Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
-      !!
-      !!  History :
-      !!    8.5  !  02-09 (A. Bozec, G. Madec)  F90: Free form and module
-      !!    9.0  !  02-09 (A. de Miranda)  rigid-lid + islands
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      INTEGER  ::   &
-         ik, it            ! temporary integers
-      
-      REAL(wp) ::   &  
-         ze3tp, ze3wp,    &  ! Last ocean level thickness at T- and W-points
-         zdepwp,          &  ! Ajusted ocean depth to avoid too small e3t
-         zdepth,          &  !    "         "
-         zmax, zmin,      &  ! Maximum and minimum depth
-         zdiff               ! temporary scalar
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zprt  !    "           "
-
-      LOGICAL ::  ll_print                          ! Allow  control print for debugging
-
+      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk       ! dummy loop indices
+      INTEGER  ::   ik, it           ! temporary integers
+      LOGICAL  ::   ll_print         ! Allow  control print for debugging
+      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
+      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
+      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth
+      REAL(wp) ::   zdiff            ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zprt   ! 3D workspace
       !!---------------------------------------------------------------------
-      !! OPA8.5, LODYC-IPSL (2002)
-      !!---------------------------------------------------------------------
-      
-      ! Local variable for debugging
-      ll_print=.FALSE.
-!!!   ll_print=.TRUE.
-      
-      ! Initialization of constant
-      zmax = gdepw_0(jpk) + e3t_0(jpk)
-      zmin = gdepw_0(4)
-      
-      ! Ocean depth
-      IF(lwp .AND. ll_print) THEN 
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
-         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
-      ENDIF
 
       IF(lwp) WRITE(numout,*)
@@ -835,29 +695,30 @@
       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used'
 
+      ll_print=.FALSE.                     ! Local variable for debugging
+!!    ll_print=.TRUE.
+      
+      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)'
+         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
+      ENDIF
+
 
       ! bathymetry in level (from bathy_meter)
       ! ===================
-
-      ! initialize mbathy to the maximum ocean level available
-      mbathy(:,:) = jpkm1
-
-      ! storage of land and island's number (zera and negative values) in mbathy
+      zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
+      zmin = gdepw_0(4)                    ! minimum depth = 3 levels
+
+      mbathy(:,:) = jpkm1                  ! initialize mbathy to the maximum ocean level available
+
+      !                                    ! storage of land and island's number (zera and negative values) in mbathy
+      WHERE( bathy(:,:) <= 0. )   mbathy(:,:) = INT( bathy(:,:) )
+
+      ! bounded value of bathy
+!!gm  bathy(:,:) = MIN(  zmax,  MAX( bathy(:,:), zmin )  )     !!gm : simpler   as zmin is > 0
       DO jj = 1, jpj
          DO ji= 1, jpi
-            IF( bathy(ji,jj) <= 0. )   mbathy(ji,jj) = INT( bathy(ji,jj) )
-         END DO
-      END DO
-
-      ! bounded value of bathy
-      ! minimum depth == 3 levels
-      ! maximum depth == gdepw_0(jpk)+e3t_0(jpk) 
-      ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk)
-      DO jj = 1, jpj
-         DO ji= 1, jpi
-            IF( bathy(ji,jj) <= 0. ) THEN
-               bathy(ji,jj) = 0.e0
-            ELSE
-               bathy(ji,jj) = MAX( bathy(ji,jj), zmin )
-               bathy(ji,jj) = MIN( bathy(ji,jj), zmax )
+            IF( bathy(ji,jj) <= 0. ) THEN   ;   bathy(ji,jj) = 0.e0
+            ELSE                            ;   bathy(ji,jj) = MIN(  zmax,  MAX( bathy(ji,jj), zmin )  )
             ENDIF
          END DO
@@ -877,91 +738,77 @@
       END DO
 
-     ! Scale factors and depth at T- and W-points
-     
-     ! intitialization to the reference z-coordinate
-     DO jk = 1, jpk
-        gdept(:,:,jk) = gdept_0(jk)
-        gdepw(:,:,jk) = gdepw_0(jk)
-        e3t  (:,:,jk) = e3t_0  (jk)
-        e3w  (:,:,jk) = e3w_0  (jk)
-     END DO
-     hdept(:,:) = gdept(:,:,2 )
-     hdepw(:,:) = gdepw(:,:,3 )     
-     
-     ! 
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           ik = mbathy(ji,jj)
-           ! ocean point only
-           IF( ik > 0 ) THEN
-              ! max ocean level case
-              IF( ik == jpkm1 ) THEN
-                 zdepwp = bathy(ji,jj)
-                 ze3tp  = bathy(ji,jj) - gdepw_0(ik)
-                 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
-                 e3t(ji,jj,ik  ) = ze3tp
-                 e3t(ji,jj,ik+1) = ze3tp
-                 e3w(ji,jj,ik  ) = ze3wp
-                 e3w(ji,jj,ik+1) = ze3tp
-                 gdepw(ji,jj,ik+1) = zdepwp
-                 gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
-                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
-                 ! standard case
-              ELSE
-!!alex
-                 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN
-                    gdepw(ji,jj,ik+1) = bathy(ji,jj)
-                 ELSE
-!!alex ctl          write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1)
-                    gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
-                 ENDIF
-!!Alex
-!!Alex           gdepw(ji,jj,ik+1) = bathy(ji,jj)
-!!bug gm  verifier les gdepw_0
-                 gdept(ji,jj,ik  ) =  gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) )   &
-                    &                             * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
-                    &                             / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
-                 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
-                    &                      / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
-                 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
-                    &                * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
-                 !       ... on ik+1
-                 e3w(ji,jj,ik+1) = e3t(ji,jj,ik)
-                 e3t(ji,jj,ik+1) = e3t(ji,jj,ik)
-                 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
-              ENDIF
-           ENDIF
-        END DO
-     END DO
-
-     it = 0
-     DO jj = 1, jpj
-        DO ji = 1, jpi
-           ik = mbathy(ji,jj)
-           ! ocean point only
-           IF( ik > 0 ) THEN
-              ! bathymetry output
-              hdept(ji,jj) = gdept(ji,jj,ik  )
-              hdepw(ji,jj) = gdepw(ji,jj,ik+1)
-              e3tp (ji,jj) = e3t(ji,jj,ik  )
-              e3wp (ji,jj) = e3w(ji,jj,ik  )
-              ! test
-              zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
-              IF( zdiff <= 0. .AND. lwp ) THEN 
-                 it=it+1
-                 WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
-                 WRITE(numout,*) ' bathy = ', bathy(ji,jj)
-                 WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1),   &
-                                 ' zdiff   = ', zdiff
-                 WRITE(numout,*) ' e3tp    = ', e3t(ji,jj,ik  ), ' e3wp    = ', e3w(ji,jj,ik  )
-              ENDIF
-           ENDIF
-        END DO
+      ! Scale factors and depth at T- and W-points
+      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate
+         gdept(:,:,jk) = gdept_0(jk)
+         gdepw(:,:,jk) = gdepw_0(jk)
+         e3t  (:,:,jk) = e3t_0  (jk)
+         e3w  (:,:,jk) = e3w_0  (jk)
+      END DO
+      hdept(:,:) = gdept(:,:,2 )
+      hdepw(:,:) = gdepw(:,:,3 )     
+      ! 
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ik = mbathy(ji,jj)
+            IF( ik > 0 ) THEN               ! ocean point only
+               ! max ocean level case
+               IF( ik == jpkm1 ) THEN
+                  zdepwp = bathy(ji,jj)
+                  ze3tp  = bathy(ji,jj) - gdepw_0(ik)
+                  ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
+                  e3t(ji,jj,ik  ) = ze3tp
+                  e3t(ji,jj,ik+1) = ze3tp
+                  e3w(ji,jj,ik  ) = ze3wp
+                  e3w(ji,jj,ik+1) = ze3tp
+                  gdepw(ji,jj,ik+1) = zdepwp
+                  gdept(ji,jj,ik  ) = gdept_0(ik-1) + ze3wp
+                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
+                  !
+               ELSE                         ! standard case
+                  IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN   ;   gdepw(ji,jj,ik+1) = bathy(ji,jj)
+                  ELSE                                       ;   gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
+                  ENDIF
+!gm Bug?  check the gdepw_0
+                  !       ... on ik
+                  gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   &
+                     &                          * ((gdept_0(      ik  ) - gdepw_0(ik) )   &
+                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ))
+                  e3t  (ji,jj,ik) = e3t_0  (ik) * ( gdepw  (ji,jj,ik+1) - gdepw_0(ik) )   & 
+                     &                          / ( gdepw_0(      ik+1) - gdepw_0(ik) ) 
+                  e3w  (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) )   &
+                     &                  * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
+                  !       ... on ik+1
+                  e3w  (ji,jj,ik+1) = e3t  (ji,jj,ik)
+                  e3t  (ji,jj,ik+1) = e3t  (ji,jj,ik)
+                  gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
+               ENDIF
+            ENDIF
+         END DO
+      END DO
+      !
+      it = 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ik = mbathy(ji,jj)
+            IF( ik > 0 ) THEN               ! ocean point only
+               hdept(ji,jj) = gdept(ji,jj,ik  )
+               hdepw(ji,jj) = gdepw(ji,jj,ik+1)
+               e3tp (ji,jj) = e3t(ji,jj,ik  )
+               e3wp (ji,jj) = e3w(ji,jj,ik  )
+               ! test
+               zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik  )
+               IF( zdiff <= 0. .AND. lwp ) THEN 
+                  it = it + 1
+                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj
+                  WRITE(numout,*) ' bathy = ', bathy(ji,jj)
+                  WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
+                  WRITE(numout,*) ' e3tp  = ', e3t  (ji,jj,ik), ' e3wp  = ', e3w  (ji,jj,ik  )
+               ENDIF
+            ENDIF
+         END DO
       END DO
 
       ! Scale factors and depth at U-, V-, UW and VW-points
-
-      ! initialisation to z-scale factors
-      DO jk = 1, jpk
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
          e3u (:,:,jk)  = e3t_0(jk)
          e3v (:,:,jk)  = e3t_0(jk)
@@ -969,199 +816,96 @@
          e3vw(:,:,jk)  = e3w_0(jk)
       END DO
-
-     ! Computed as the minimum of neighbooring scale factors
-     DO jk = 1,jpk
-        DO jj = 1, jpjm1
-           DO ji = 1, fs_jpim1   ! vector opt.
-              e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
-              e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
-              e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
-              e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
-           END DO
-        END DO
-     END DO
+      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
+               e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
+               e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
+               e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
+            END DO
+         END DO
+      END DO
+      !                                     ! Boundary conditions
+      CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
+      CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3u (:,:,jk) == 0.e0 )   e3u (:,:,jk) = e3t_0(jk)
+         WHERE( e3v (:,:,jk) == 0.e0 )   e3v (:,:,jk) = e3t_0(jk)
+         WHERE( e3uw(:,:,jk) == 0.e0 )   e3uw(:,:,jk) = e3w_0(jk)
+         WHERE( e3vw(:,:,jk) == 0.e0 )   e3vw(:,:,jk) = e3w_0(jk)
+      END DO
+      
+      ! Scale factor at F-point
+      DO jk = 1, jpk                        ! initialisation to z-scale factors
+         e3f (:,:,jk) = e3t_0(jk)
+      END DO
+      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( e3f, 'F', 1. )          ! Boundary conditions
+      !
+      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
+         WHERE( e3f(:,:,jk) == 0.e0 )   e3f(:,:,jk) = e3t_0(jk)
+      END DO
+!!gm  bug ? :  must be a do loop with mj0,mj1
+      ! 
+      e3t(:,mj0(1),:) = e3t(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2
+      e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
+      e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
+      e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
+      e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
+
+      ! Control of the sign
+      IF( MINVAL( e3t  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t   <= 0' )
+      IF( MINVAL( e3w  (:,:,:) ) <= 0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w   <= 0' )
+      IF( MINVAL( gdept(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
+      IF( MINVAL( gdepw(:,:,:) ) <  0.e0 )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw <  0' )
      
-     ! Boundary conditions
-     CALL lbc_lnk( e3u , 'U', 1. )   ;   CALL lbc_lnk( e3uw, 'U', 1. )
-     CALL lbc_lnk( e3v , 'V', 1. )   ;   CALL lbc_lnk( e3vw, 'V', 1. )
-     
-     ! set to z-scale factor if zero (i.e. along closed boundaries)
-     DO jk = 1, jpk
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk)
-              IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk)
-              IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk)
-              IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk)
-           END DO
-        END DO
-     END DO
-     
-     ! Scale factor at F-point
-     
-     ! initialisation to z-scale factors
-     DO jk = 1, jpk
-        e3f (:,:,jk) = e3t_0(jk)
-     END DO
-     
-     ! Computed as the minimum of neighbooring V-scale factors
-     DO jk = 1, jpk
-        DO jj = 1, jpjm1
-           DO ji = 1, fs_jpim1   ! vector opt.
-              e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
-           END DO
-        END DO
-     END DO
-     ! Boundary conditions
-     CALL lbc_lnk( e3f, 'F', 1. )
-     
-     ! set to z-scale factor if zero (i.e. along closed boundaries)
-     DO jk = 1, jpk
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk)
-           END DO
-        END DO
-     END DO
-!!bug ? gm:  must be a do loop with mj0,mj1
-
-     ! we duplicate factor scales for jj = 1 and jj = 2
-     e3t(:,mj0(1),:) = e3t(:,mj0(2),:) 
-     e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 
-     e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 
-     e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 
-     e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 
-
-
-     
-     ! Compute gdep3w (vertical sum of e3w)
-     
-     gdep3w   (:,:,1) = 0.5 * e3w (:,:,1)
-     DO jk = 2, jpk
-        gdep3w   (:,:,jk) = gdep3w   (:,:,jk-1) + e3w (:,:,jk) 
-     END DO
+      ! Compute gdep3w (vertical sum of e3w)
+      gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
+      DO jk = 2, jpk
+         gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 
+      END DO
           
-     ! Control print
- 9600 FORMAT(9x,' level   gdept    gdepw     e3t      e3w   ')
- 9610 FORMAT(10x,i4,4f9.2)
-      IF(lwp .AND. ll_print) THEN
+      !                                               ! ================= !
+      IF(lwp .AND. ll_print) THEN                     !   Control print   !
+         !                                            ! ================= !
          DO jj = 1,jpj
             DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3t(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3t(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3w(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3w(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3u(ji,jj,ik)
-            END DO
-         END DO
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3u(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = e3v(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3v(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik = MAX(mbathy(ji,jj),1) 
-               zprt(ji,jj) = e3f(ji,jj,ik)
-            END DO
-         END DO
-
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr e3f(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               ik =  MAX(mbathy(ji,jj),1)
-               zprt(ji,jj) = gdep3w(ji,jj,ik)
-            END DO
-         END DO
-         WRITE(numout,*)
-         WRITE(numout,*) 'domzgr gdep3w(mbathy)'
-         CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
-
+               ik = MAX( mbathy(ji,jj), 1 )
+               zprt(ji,jj,1) = e3t   (ji,jj,ik)
+               zprt(ji,jj,2) = e3w   (ji,jj,ik)
+               zprt(ji,jj,3) = e3u   (ji,jj,ik)
+               zprt(ji,jj,4) = e3v   (ji,jj,ik)
+               zprt(ji,jj,5) = e3f   (ji,jj,ik)
+               zprt(ji,jj,6) = gdep3w(ji,jj,ik)
+            END DO
+         END DO
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
       ENDIF  
-
-
-      DO jk = 1,jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi 
-               IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
-                     WRITE(numout,*) ' =========           --------------- '
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr.psteps'
-               ENDIF
-               IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
-                     WRITE(numout,*) ' =========        ------------------ '
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr.psteps'
-               ENDIF
-            END DO
-         END DO
-      END DO  
-
-   IF(lwp) THEN
-      WRITE(numout,*) ' e3t lev 21 '
-      CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3w lev 21  '
-      CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3u lev 21  '
-      CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3v lev 21  '
-      CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3f lev 21  '
-      CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3t lev 22 '
-      CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3w lev 22  '
-      CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3u lev 22  '
-      CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3v lev 22  '
-      CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-      WRITE(numout,*) ' e3f lev 22  '
-      CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
-   ENDIF
-
-      ! ===========
-      ! Zoom domain 
-      ! ===========
-
-      IF( lzoom )   CALL zgr_bat_zoom
-
-      ! ================
-      ! Bathymetry check
-      ! ================
-
-      CALL zgr_bat_ctl
-
+       
+      !                                               ! =============== !
+      IF( lzoom )   CALL zgr_bat_zoom                 !   Zoom domain   !
+      !                                               ! =============== !
+
+      !                                               ! =================== !
+      IF( .NOT. lk_c1d )   CALL zgr_bat_ctl           !   Bathymetry check  !
+      !                                               ! =================== !
    END SUBROUTINE zgr_zps
 
@@ -1223,31 +967,19 @@
       !!      schemes.
       !!
-      !! Reference :
-      !!      Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
-      !!
-      !! History :
-      !!        !  95-12  (G. Madec)  Original code : s vertical coordinate
-      !!        !  97-07  (G. Madec)  lbc_lnk call
-      !!        !  97-04  (J.-O. Beismann) 
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
-      !!   9.0  !  05-10  (A. Beckmann)  new stretching function
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk, jl
-      INTEGER  ::   iip1, ijp1, iim1, ijm1
-      REAL(wp) ::   zrmax, ztaper
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zenv, ztmp, zmsk, zri, zrj, zhbat
-
+      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
+      !!----------------------------------------------------------------------
+      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), DIMENSION(jpi,jpj) ::   zenv, ztmp, zmsk    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zri , zrj , zhbat   !  -     -
+      !!
       NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
       !!----------------------------------------------------------------------
 
-      ! Read Namelist nam_zgr_sco : sigma-stretching parameters
-      ! -------------------------
-      REWIND( numnam )
+      REWIND( numnam )                        ! Read Namelist nam_zgr_sco : sigma-stretching parameters
       READ  ( numnam, nam_zgr_sco )
 
-      IF(lwp) THEN
+      IF(lwp) THEN                            ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
@@ -1262,30 +994,19 @@
       ENDIF
 
-      ! ???
-      ! ------------------
-      hift(:,:) = sbot_min
+      hift(:,:) = sbot_min                     ! set the minimum depth for the s-coordinate
       hifu(:,:) = sbot_min
       hifv(:,:) = sbot_min
       hiff(:,:) = sbot_min
-
-
-      ! set maximum ocean depth 
-      ! -----------------------
-       DO jj = 1, jpj
-          DO ji = 1, jpi
-             bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) )
-          END DO
-       END DO
-
-      ! Define the envelop bathymetry
-      ! =============================
+      !                                        ! set maximum ocean depth 
+      bathy(:,:) = MIN( sbot_max, bathy(:,:) )
+
+      !                                        ! =============================
+      !                                        ! Define the envelop bathymetry   (hbatt)
+      !                                        ! =============================
       ! Smooth the bathymetry (if required) 
-
       scosrf(:,:) = 0.e0              ! ocean surface depth (here zero: no under ice-shelf sea)
       scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
-
- 
+      !
       ! use r-value to create hybrid coordinates
- 
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1293,5 +1014,4 @@
          END DO
       END DO
-
       jl = 0
       zrmax = 1.e0
@@ -1300,8 +1020,6 @@
          !                                                  ! ================ !
          jl = jl + 1
-
          zrmax = 0.e0
          zmsk(:,:) = 0.e0
- 
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1319,6 +1037,5 @@
          IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
          ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
-         ztmp(:,:) = zmsk(:,:)
-         CALL lbc_lnk( zmsk, 'T', 1. )
+         ztmp(:,:) = zmsk(:,:)   ;   CALL lbc_lnk( zmsk, 'T', 1. )
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1326,7 +1043,7 @@
             END DO
          END DO
- 
-         WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
-
+         ! 
+         IF(lwp)WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
+         !
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1348,5 +1065,5 @@
             END DO
          END DO
- 
+         ! 
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1354,8 +1071,7 @@
             END DO
          END DO
-
+         !
          ! Apply lateral boundary condition   CAUTION: kept the value when the lbc field is zero
-         ztmp(:,:) = zenv(:,:)
-         CALL lbc_lnk( zenv, 'T', 1. )
+         ztmp(:,:) = zenv(:,:)   ;   CALL lbc_lnk( zenv, 'T', 1. )
          DO jj = 1, nlcj
             DO ji = 1, nlci
@@ -1366,65 +1082,50 @@
       END DO                                                !     End loop     !
       !                                                     ! ================ !
-
-      ! save envelop bathymetry in hbatt
+      !
+      !                                        ! envelop bathymetry saved in hbatt
       hbatt(:,:) = zenv(:,:) 
-
-!!bug gm  :   CAUTION:  tapering hard coded !!!!  orca2 only
-!!bug gm                NOT valid in mpp ===> taper have been changed
-
-       DO jj = 1, jpj
-          DO ji = 1, jpi
-!!bug mpp    ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
-             ztaper = EXP( -(gphit(ji,jj)/8)**2 )
-             hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
-          END DO
-       END DO
-
-      DO jj = 1, jpj
-         zrmax  = EXP( -(gphit(10,jj)/8)**2 )
-         ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
-         IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax
-      END DO
-
-      IF( nprint == 1 .AND. lwp )   THEN
-         WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
-         WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
-      ENDIF
-
-      ! Control print
+      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) 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)**2 )
+               hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
+            END DO
+         END DO
+      ENDIF
+      !
+      IF(lwp) THEN                             ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
+         WRITE(numout,*)
+         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
+         IF( nprint == 1 )   THEN        
+            WRITE(numout,*) ' bathy  MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
+            WRITE(numout,*) ' hbatt  MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
+         ENDIF
+      ENDIF
+
+      !                                        ! ==============================
+      !                                        !   hbatu, hbatv, hbatf fields
+      !                                        ! ==============================
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
-         WRITE(numout,*)
-         CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
-      ENDIF
-
-      ! 1. computation of hbatu, hbatv, hbatf fields
-      ! --------------------------------------------
-
+         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
+      ENDIF
       hbatu(:,:) = sbot_min
       hbatv(:,:) = sbot_min
       hbatf(:,:) = sbot_min
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
-         WRITE(numout,*)
-      ENDIF
-
       DO jj = 1, jpjm1
         DO ji = 1, jpim1
-           hbatu(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji+1,jj  ) )
-           hbatv(ji,jj)= 0.5 *( hbatt(ji  ,jj)+hbatt(ji  ,jj+1) )
-           hbatf(ji,jj)= 0.25*( hbatt(ji  ,jj)+hbatt(ji  ,jj+1)   &
-                               +hbatt(ji+1,jj)+hbatt(ji+1,jj+1) )
+           hbatu(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
+           hbatv(ji,jj) = 0.5 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
+           hbatf(ji,jj) = 0.25* ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
+              &                 + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
         END DO
       END DO
- 
+      ! 
       ! Apply lateral boundary condition
-      ! CAUTION: retain non zero value in the initial file
-      !          this should be OK for orca cfg, not for EEL
-      zhbat(:,:) = hbatu(:,:)
-      CALL lbc_lnk( hbatu, 'U', 1. )
+!!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
+      zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1. )
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1435,7 +1136,5 @@
          END DO
       END DO
-
-      zhbat(:,:) = hbatv(:,:)
-      CALL lbc_lnk( hbatv, 'V', 1. )
+      zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1. )
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1446,7 +1145,5 @@
          END DO
       END DO
-
-      zhbat(:,:) = hbatf(:,:)
-      CALL lbc_lnk( hbatf, 'F', 1. )
+      zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1. )
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1458,5 +1155,4 @@
       END DO
 
-
 !!bug:  key_helsinki a verifer
       hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
@@ -1466,8 +1162,8 @@
 
       IF( nprint == 1 .AND. lwp )   THEN
-         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 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(:,:) )
@@ -1477,50 +1173,46 @@
 !! helsinki
 
-      ! 2. Computation of gsig and esig fields
-      ! --------------------------------------
-
-      ! 2.1 Coefficients for model level depth at w- and t-levels
-
+      !                                            ! =======================
+      !                                            !   s-ccordinate fields     (gdep., e3.)
+      !                                            ! =======================
+      !
+      ! non-dimensional "sigma" for model level depth at w- and t-levels
       DO jk = 1, jpk
         gsigw(jk) = -fssig( FLOAT(jk)-0.5 )
         gsigt(jk) = -fssig( FLOAT(jk)     )
       END DO
-
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)
-
-      ! 2.2 Coefficients for vertical scale factors at w-, t- levels
-
-!!bug gm :  define it from analytical function, not like juste bellow....
-!!          or betteroffer the 2 possibilities....
-      DO jk=2,jpk
-        esigw(jk)=gsigt(jk)-gsigt(jk-1)
-      END DO
-      DO jk=1,jpkm1
-        esigt(jk)=gsigw(jk+1)-gsigw(jk)
-      END DO
-        esigw(1)=esigw(2)
-        esigt(jpk)=esigt(jpkm1)
-
-      ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
-
+      !
+      ! 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 ) = esigw(  2  )
+      esigt(jpk) = esigt(jpkm1)
+
+!!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 * esigw(1)
       DO jk = 2, jpk
-        gsi3w(jk) = gsi3w(jk-1)+ esigw(jk)
-      END DO
-
-!!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
+         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
-!! change the solver stat....
-!!       zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
-!!       gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)
-         gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*(FLOAT(jk)-0.5)/FLOAT(jpkm1))
-!!       zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 
-!!       gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)
-!!       gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw)
-         gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1))
-         gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1))
-      END DO
-
-!!bug gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
+         zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
+         zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 
+         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(:,:)*zcoefw)
+      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
@@ -1530,5 +1222,5 @@
               e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
               e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
-
+              !
               e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
               e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
@@ -1537,6 +1229,6 @@
          END DO
       END DO
-
-! HYBRID
+      !
+      ! HYBRID : 
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1547,32 +1239,44 @@
          END DO
       END DO
-
-      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
-
-
-      ! ===========
-      ! Zoom domain 
-      ! ===========
-
-      IF( lzoom )   CALL zgr_bat_zoom
-
-      ! 2.4 Control print
-
-      IF(lwp) THEN
+      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
+         &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
+
+
+      !                                            ! ===========
+      IF( lzoom )   CALL zgr_bat_zoom              ! Zoom domain 
+      !                                            ! ===========
+
+      !                                            ! =============
+      IF(lwp) THEN                                 ! Control print
+         !                                         ! =============
          WRITE(numout,*) 
          WRITE(numout,*) ' domzgr: vertical coefficients for model level'
-         WRITE(numout,9400)
-         WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk)
-      ENDIF
- 9400 FORMAT(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')
- 9410 FORMAT(10x,i4,5f11.4)
-
-      IF(lwp) THEN
+         WRITE(numout, "(9x,'  level    gsigt      gsigw      esigt      esigw      gsi3w')" )
+         WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
+      ENDIF
+      IF( nprint == 1  .AND. lwp )   THEN         ! min max values over the local domain
+         WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)   ), ' MAX ', MAXVAL( mbathy(:,:) )
+         WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ),   &
+            &                          ' w ', MINVAL( fsdepw(:,:,:) ), '3w '  , MINVAL( fsde3w(:,:,:) )
+         WRITE(numout,*) ' MIN val e3    t ', MINVAL( fse3t (:,:,:) ), ' f '  , MINVAL( fse3f (:,:,:) ),   &
+            &                          ' u ', MINVAL( fse3u (:,:,:) ), ' u '  , MINVAL( fse3v (:,:,:) ),   &
+            &                          ' uw', MINVAL( fse3uw(:,:,:) ), ' vw'  , MINVAL( fse3vw(:,:,:) ),   &
+            &                          ' w ', MINVAL( fse3w (:,:,:) )
+
+         WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ),   &
+            &                          ' w ', MAXVAL( fsdepw(:,:,:) ), '3w '  , MAXVAL( fsde3w(:,:,:) )
+         WRITE(numout,*) ' MAX val e3    t ', MAXVAL( fse3t (:,:,:) ), ' f '  , MAXVAL( fse3f (:,:,:) ),   &
+            &                          ' u ', MAXVAL( fse3u (:,:,:) ), ' u '  , MAXVAL( fse3v (:,:,:) ),   &
+            &                          ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw'  , MAXVAL( fse3vw(:,:,:) ),   &
+            &                          ' w ', MAXVAL( fse3w (:,:,:) )
+      ENDIF
+      !
+      IF(lwp) THEN                                  ! 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,9420)
-         WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk),     &
-                             fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk)
+         WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
+         WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk),     &
+            &                                 fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
          DO jj = mj0(20), mj1(20)
             DO ji = mi0(20), mi1(20)
@@ -1580,7 +1284,7 @@
                WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
                WRITE(numout,*) ' ~~~~~~  --------------------'
-               WRITE(numout,9420)
-               WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     &
-                    &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
+               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
+               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
+                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
             END DO
          END DO
@@ -1590,49 +1294,31 @@
                WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
                WRITE(numout,*) ' ~~~~~~  --------------------'
-               WRITE(numout,9420)
-               WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk),     &
-                    &                 fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
-            END DO
-         END DO
-      ENDIF
-
- 9420 FORMAT(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')
- 9430 FORMAT(10x,i4,4f9.2)
-
-!!bug gm   no more necessary?  if ! defined key_helsinki
+               WRITE(numout,"(9x,' level   gdept    gdepw    gde3w     e3t      e3w  ')")
+               WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk),     &
+                  &                                 fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
+            END DO
+         END DO
+      ENDIF
+
+!!gm bug?  no more necessary?  if ! defined key_helsinki
       DO jk = 1, jpk
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' e r r o r :         e3w or e3t =< 0 '
-                     WRITE(numout,*) ' =========           --------------- '
-                     WRITE(numout,*)
-                     WRITE(numout,*) '             point (i,j,k)= ',ji,jj,jk
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr'
+                  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. .OR. fsdept(ji,jj,jk) < 0. ) THEN
-                  IF(lwp) THEN
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' e r r o r :      gdepw or gdept < 0 '
-                     WRITE(numout,*) ' =========        ------------------ '
-                     WRITE(numout,*)
-                     WRITE(numout,*) '             point (i,j,k)= ', ji, jj, jk
-                     WRITE(numout,*)
-                  ENDIF
-                  STOP 'domzgr'
+                  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
-!!bug gm    #endif
-
+!!gm bug    #endif
+      !
    END SUBROUTINE zgr_sco
 
 #endif
-
 
    !!======================================================================
