Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 2481)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 2482)
@@ -57,5 +57,5 @@
       INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
       REAL(wp) ::   zrefdep         ! local real
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk, zprt, zprw   ! 2D workspace
       !!----------------------------------------------------------------------
 
@@ -220,6 +220,6 @@
                CALL iom_get( inum4, jpdom_data, 'gdepw', gdepw )
             ELSE                                                           ! 2D bottom depth
-               CALL iom_get( inum4, jpdom_data, 'hdept', hdept )
-               CALL iom_get( inum4, jpdom_data, 'hdepw', hdepw )
+               CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
+               CALL iom_get( inum4, jpdom_data, 'hdepw', zprw )
                !
                DO jk = 1, jpk                                              ! deduces the 3D depth
@@ -231,6 +231,6 @@
                      ik = mbkt(ji,jj)
                      IF( ik > 0 ) THEN
-                        gdepw(ji,jj,ik+1) = hdepw(ji,jj)
-                        gdept(ji,jj,ik  ) = hdept(ji,jj)
+                        gdepw(ji,jj,ik+1) = zprw(ji,jj)
+                        gdept(ji,jj,ik  ) = zprt(ji,jj)
                         gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
                      ENDIF
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2481)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2482)
@@ -18,32 +18,33 @@
    !! ----------------------------
    !                                              !!* Namelist namdom : time & space domain *
-   INTEGER , PUBLIC ::   nn_bathy     =  0         !: = 0/1 ,compute/read the bathymetry file
-   REAL(wp), PUBLIC ::   rn_e3zps_min = 5.0_wp     !: miminum thickness for partial steps (meters)
-   REAL(wp), PUBLIC ::   rn_e3zps_rat = 0.1_wp     !: minimum thickness ration for partial steps
-   INTEGER , PUBLIC ::   nn_msh       = 0          !: = 1 create a mesh-mask file
-   INTEGER , PUBLIC ::   nn_acc       = 0          !: = 0/1 use of the acceleration of convergence technique
-   REAL(wp), PUBLIC ::   rn_atfp      = 0.1_wp     !: asselin time filter parameter
+   INTEGER , PUBLIC ::   nn_bathy     =    0       !: = 0/1 ,compute/read the bathymetry file
+   REAL(wp), PUBLIC ::   rn_hmin      =   -3.0_wp  !: minimum ocean depth (>0) or minimum number of ocean levels (<0)
+   REAL(wp), PUBLIC ::   rn_e3zps_min =    5.0_wp  !: miminum thickness for partial steps (meters)
+   REAL(wp), PUBLIC ::   rn_e3zps_rat =    0.1_wp  !: minimum thickness ration for partial steps
+   INTEGER , PUBLIC ::   nn_msh       =    0       !: = 1 create a mesh-mask file
+   INTEGER , PUBLIC ::   nn_acc       =    0       !: = 0/1 use of the acceleration of convergence technique
+   REAL(wp), PUBLIC ::   rn_atfp      =    0.1_wp  !: asselin time filter parameter
    REAL(wp), PUBLIC ::   rn_rdt       = 3600._wp   !: time step for the dynamics (and tracer if nacc=0)
    REAL(wp), PUBLIC ::   rn_rdtmin    = 3600._wp   !: minimum time step on tracers
    REAL(wp), PUBLIC ::   rn_rdtmax    = 3600._wp   !: maximum time step on tracers
    REAL(wp), PUBLIC ::   rn_rdth      =  800._wp   !: depth variation of tracer step
-   INTEGER , PUBLIC ::   nn_baro      = 64         !: number of barotropic time steps (key_dynspg_ts)
-   INTEGER , PUBLIC ::   nn_closea    =  0         !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
-
-   !                                          ! old non-DOCTOR names still used in the model
-   INTEGER , PUBLIC ::   ntopo                !: = 0/1 ,compute/read the bathymetry file
-   REAL(wp), PUBLIC ::   e3zps_min            !: miminum thickness for partial steps (meters)
-   REAL(wp), PUBLIC ::   e3zps_rat            !: minimum thickness ration for partial steps
-   INTEGER , PUBLIC ::   nmsh                 !: = 1 create a mesh-mask file
-   INTEGER , PUBLIC ::   nacc                 !: = 0/1 use of the acceleration of convergence technique
-   REAL(wp), PUBLIC ::   atfp                 !: asselin time filter parameter
-   REAL(wp), PUBLIC ::   rdt                  !: time step for the dynamics (and tracer if nacc=0)
-   REAL(wp), PUBLIC ::   rdtmin               !: minimum time step on tracers
-   REAL(wp), PUBLIC ::   rdtmax               !: maximum time step on tracers
-   REAL(wp), PUBLIC ::   rdth                 !: depth variation of tracer step
-   INTEGER , PUBLIC ::   nclosea              !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
-
-
-   !                                         !!! associated variables
+   INTEGER , PUBLIC ::   nn_baro      =   64       !: number of barotropic time steps (key_dynspg_ts)
+   INTEGER , PUBLIC ::   nn_closea    =    0       !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
+
+   !                                    !! old non-DOCTOR names still used in the model
+   INTEGER , PUBLIC ::   ntopo           !: = 0/1 ,compute/read the bathymetry file
+   REAL(wp), PUBLIC ::   e3zps_min       !: miminum thickness for partial steps (meters)
+   REAL(wp), PUBLIC ::   e3zps_rat       !: minimum thickness ration for partial steps
+   INTEGER , PUBLIC ::   nmsh            !: = 1 create a mesh-mask file
+   INTEGER , PUBLIC ::   nacc            !: = 0/1 use of the acceleration of convergence technique
+   REAL(wp), PUBLIC ::   atfp            !: asselin time filter parameter
+   REAL(wp), PUBLIC ::   rdt             !: time step for the dynamics (and tracer if nacc=0)
+   REAL(wp), PUBLIC ::   rdtmin          !: minimum time step on tracers
+   REAL(wp), PUBLIC ::   rdtmax          !: maximum time step on tracers
+   REAL(wp), PUBLIC ::   rdth            !: depth variation of tracer step
+   INTEGER , PUBLIC ::   nclosea         !: =0 suppress closed sea/lake from the ORCA domain or not (=1)
+
+
+   !                                                  !!! associated variables
    INTEGER , PUBLIC                 ::   neuler  = 0   !: restart euler forward option (0=Euler)
    REAL(wp), PUBLIC                 ::   atfp1         !: asselin time filter coeff. (atfp1= 1-2*atfp)
@@ -150,5 +151,4 @@
    REAL(wp), PUBLIC, DIMENSION(jpk)     ::   gdept_0, gdepw_0   !: reference depth of t- and w-points (m)
    REAL(wp), PUBLIC, DIMENSION(jpk)     ::   e3t_0  , e3w_0     !: reference vertical scale factors at T- and W-pts (m)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hdept  , hdepw     !: ocean bottom depth at T and W points
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   e3tp   , e3wp      !: ocean bottom level thickness at T and W points
 
@@ -163,6 +163,6 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   scosrf, scobot   !: ocean surface and bottom topographies 
    !                                                          !  (if deviating from coordinate surfaces in HYBRID)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hifv  , hiff     !: interface depth between stretching    at  V--F
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hift  , hifu     !: and quasi-uniform spacing                 T--U  points (m)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hifv  , hiff     !: interface depth between stretching at  V--F
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hift  , hifu     !: and quasi-uniform spacing              T--U  points (m)
 
    !!----------------------------------------------------------------------
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 2481)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 2482)
@@ -56,26 +56,25 @@
       !!      if     nmsh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
       !!      if 3 < nmsh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 
-      !!                        corresponding to the depth of the bottom points hdep[tw]
+      !!                        corresponding to the depth of the bottom t- and w-points
       !!      if 6 < nmsh <= 9: write 2D arrays corresponding to the depth and the
-      !!                        thickness of the bottom points hdep[tw] and e3[tw]_ps
+      !!                        thickness (e3[tw]_ps) of the bottom points 
       !!
       !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
       !!                                   masks, depth and vertical scale factors
       !!----------------------------------------------------------------------
-      INTEGER                          ::   inum0    ! temprary units for 'mesh_mask.nc' file
-      INTEGER                          ::   inum1    ! temprary units for 'mesh.nc'      file
-      INTEGER                          ::   inum2    ! temprary units for 'mask.nc'      file
-      INTEGER                          ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
-      INTEGER                          ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
-      INTEGER                          ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt     ! temporary array for bathymetry 
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu    ! 3D depth of U point
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepv    ! 3D depth of V point
-      CHARACTER(len=21)                ::   clnam0   ! filename (mesh and mask informations)
-      CHARACTER(len=21)                ::   clnam1   ! filename (mesh informations)
-      CHARACTER(len=21)                ::   clnam2   ! filename (mask informations)
-      CHARACTER(len=21)                ::   clnam3   ! filename (horizontal mesh informations)
-      CHARACTER(len=21)                ::   clnam4   ! filename (vertical   mesh informations)
-      !!----------------------------------------------------------------------
+      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
+      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
+      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
+      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
+      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
+      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
+      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
+      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
+      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
+      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
+      INTEGER           ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zprt , zprw    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdepu, zdepv   ! 3D workspace
+     !!----------------------------------------------------------------------
 
       IF(lwp) WRITE(numout,*)
@@ -216,10 +215,10 @@
             DO jj = 1,jpj   
                DO ji = 1,jpi
-                  hdept(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
-                  hdepw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
+                  zprt(ji,jj) = gdept(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1)
+                  zprw(ji,jj) = gdepw(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1)
                END DO
             END DO
-            CALL iom_rstput( 0, 0, inum4, 'hdept' , hdept, ktype = jp_r4 )     
-            CALL iom_rstput( 0, 0, inum4, 'hdepw' , hdepw, ktype = jp_r4 ) 
+            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r4 )     
+            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r4 ) 
          ENDIF
          !
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 2481)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 2482)
@@ -236,6 +236,6 @@
          IF( .NOT. ldbletanh ) THEN
             DO jk = 1, jpk
-               zw = FLOAT( jk )
-               zt = FLOAT( jk ) + 0.5_wp
+               zw = REAL( jk , wp )
+               zt = REAL( jk , wp ) + 0.5_wp
                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 ) )  )
@@ -248,12 +248,12 @@
                zt = FLOAT( jk ) + 0.5_wp
                ! Double tanh function
-               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )    &
-     &                                         + za2 * zacr2* LOG ( COSH( (zw-zkth2)/zacr2 ) )  )
-               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )    &
-     &                                         + za2 * zacr2* LOG ( COSH( (zt-zkth2)/zacr2 ) )  )
-               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )    &
-     &                                         + za2        * TANH(       (zw-zkth2)/zacr2   )
-               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )    &
-     &                                         + za2        * TANH(       (zt-zkth2)/zacr2   )
+               gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
+                  &                            + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
+               gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
+                  &                            + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
+               e3w_0  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )    &
+                  &                            + za2        * TANH(       (zw-zkth2) / zacr2 )
+               e3t_0  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )    &
+                  &                            + za2        * TANH(       (zt-zkth2) / zacr2 )
             END DO
          ENDIF
@@ -315,7 +315,7 @@
       INTEGER  ::   inum                      ! temporary logical unit
       INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position
-      INTEGER  ::   ii0, ii1, ij0, ij1        ! indices
+      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices
       REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics 
-      REAL(wp) ::   zi     , zj     , zh      ! temporary scalars
+      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
       INTEGER , DIMENSION(jpidta,jpjdta) ::   idta   ! global domain integer data
       REAL(wp), DIMENSION(jpidta,jpjdta) ::   zdta   ! global domain scalar data
@@ -497,4 +497,16 @@
       ENDIF
       !
+      !                                               ! =========================== !
+      !                                               !     set a minimum depth     !
+      !                                               ! =========================== !
+      IF( rn_hmin < 0._wp ) THEN    ;   ik = INT( rn_hmin ) + 1                                    ! from a nb of level
+      ELSE                          ;   ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
+      ENDIF
+      zhmin = gdepw_0(ik+1)                                                         ! minimum depth = ik+1 w-levels 
+      WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
+      ELSE WHERE                     ;   bathy(:,:) = MIN(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
+      END WHERE
+      IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
+      !
    END SUBROUTINE zgr_bat
 
@@ -797,5 +809,5 @@
       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) ::   zmax             ! Maximum depth
       REAL(wp) ::   zdiff            ! temporary scalar
       REAL(wp) ::   zrefdep          ! temporary scalar
@@ -820,25 +832,7 @@
       ! bathymetry in level (from bathy_meter)
       ! ===================
-      zmax = gdepw_0(jpk) + e3t_0(jpk)     ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
-
-      zrefdep = 25._wp
-      nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 )   ! shallowest W level Below zrefdep
-      IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
-      zmin = gdepw_0(nlbelow)              ! minimum depth = nlbelow-1 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._wp )   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._wp ) THEN   ;   bathy(ji,jj) = 0._wp
-            ELSE                               ;   bathy(ji,jj) = MIN(  zmax,  MAX( bathy(ji,jj), zmin )  )
-            ENDIF
-         END DO
-      END DO
+      zmax = gdepw_0(jpk) + e3t_0(jpk)          ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
+      mbathy(:,:) = jpkm1                       ! initialize mbathy to the maximum ocean level available
+      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
 
       ! Compute mbathy for ocean points (i.e. the number of ocean levels)
@@ -858,6 +852,4 @@
          e3w  (:,:,jk) = e3w_0  (jk)
       END DO
-      hdept(:,:) = gdept(:,:,2 )
-      hdepw(:,:) = gdepw(:,:,3 )     
       ! 
       DO jj = 1, jpj
@@ -905,6 +897,4 @@
             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  )
@@ -924,8 +914,8 @@
       ! Scale factors and depth at U-, V-, UW and VW-points
       DO jk = 1, jpk                        ! initialisation to z-scale factors
-         e3u (:,:,jk)  = e3t_0(jk)
-         e3v (:,:,jk)  = e3t_0(jk)
-         e3uw(:,:,jk)  = e3w_0(jk)
-         e3vw(:,:,jk)  = e3w_0(jk)
+         e3u (:,:,jk) = e3t_0(jk)
+         e3v (:,:,jk) = e3t_0(jk)
+         e3uw(:,:,jk) = e3w_0(jk)
+         e3vw(:,:,jk) = e3w_0(jk)
       END DO
       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
