MODULE domzgr
   !!==============================================================================
   !!                       ***  MODULE domzgr   ***
   !! Ocean initialization : domain initialization
   !!==============================================================================

   !!----------------------------------------------------------------------
   !!   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
   !!       zgr_bat_ctl  : check the bathymetry files
   !!       zgr_z        : reference z-coordinate 
   !!       zgr_zco      : z-coordinate 
   !!       zgr_zps      : z-coordinate with partial steps
   !!       zgr_sco      : s-coordinate
   !!---------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE in_out_manager  ! I/O manager
   USE lib_mpp         ! distributed memory computing library
   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
   USE closea
   USE solisl
   USE ini1d           ! initialization of the 1D configuration

   IMPLICIT NONE
   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)


 
   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
   !!----------------------------------------------------------------------

CONTAINS       

   SUBROUTINE dom_zgr
      !!----------------------------------------------------------------------
      !!                ***  ROUTINE dom_zgr  ***
      !!                   
      !! ** Purpose :  set the depth of model levels and the resulting 
      !!      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

      NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
      !!----------------------------------------------------------------------

      ! Read Namelist nam_zgr : vertical coordinate'
      ! ---------------------
      REWIND ( numnam )
      READ   ( numnam, nam_zgr )

      ! Parameter control and print
      ! ---------------------------
      ! Control print
      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'dom_zgr : vertical coordinate'
         WRITE(numout,*) '~~~~~~~'
         WRITE(numout,*) '          Namelist nam_zgr : set vertical coordinate'
         WRITE(numout,*) '             z-coordinate - full steps      ln_zco = ', ln_zco
         WRITE(numout,*) '             z-coordinate - partial steps   ln_zps = ', ln_zps
         WRITE(numout,*) '             s- or hybrid z-s-coordinate    ln_sco = ', ln_sco
      ENDIF

      ! Check Vertical coordinate options
      ioptio = 0
      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( 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' )
      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


   SUBROUTINE zgr_z
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE zgr_z  ***
      !!                   
      !! ** Purpose :   set the depth of model levels and the resulting 
      !!      vertical scale factors.
      !!
      !! ** Method  :   z-coordinate system (use in all type of coordinate)
      !!        The depth of model levels is defined from an analytical
      !!      function the derivative of which gives the scale factors.
      !!        both depth and scale factors only depend on k (1d arrays).
      !!              w-level: gdepw_0  = fsdep(k)
      !!                       e3w_0(k) = dk(fsdep)(k)     = fse3(k)
      !!              t-level: gdept_0  = fsdep(k+0.5)
      !!                       e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
      !!
      !! ** 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
      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
      !!----------------------------------------------------------------------

      ! Set variables from parameters
      ! ------------------------------
       zkth = ppkth       ;   zacr = ppacr
       zdzmin = ppdzmin   ;   zhmax = pphmax

      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
      !
       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 = ppa1 ;       za0 = ppa0 ;          zsur = ppsur
       ENDIF


      ! Parameter print
      ! ---------------
      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates'
         WRITE(numout,*) '    ~~~~~~~'
         IF(  ppkth == 0. ) THEN              
              WRITE(numout,*) '            Uniform grid with ',jpk-1,' layers'
              WRITE(numout,*) '            Total depth    :', zhmax
              WRITE(numout,*) '            Layer thickness:', zhmax/(jpk-1)
         ELSE
            IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
               WRITE(numout,*) '         zsur, za0, za1 computed from '
               WRITE(numout,*) '                 zdzmin = ', zdzmin
               WRITE(numout,*) '                 zhmax  = ', zhmax
            ENDIF
            WRITE(numout,*) '           Value of coefficients for vertical mesh:'
            WRITE(numout,*) '                 zsur = ', zsur
            WRITE(numout,*) '                 za0  = ', za0
            WRITE(numout,*) '                 za1  = ', za1
            WRITE(numout,*) '                 zkth = ', zkth
            WRITE(numout,*) '                 zacr = ', zacr
         ENDIF
      ENDIF


      ! Reference z-coordinate (depth - scale factor at T- and W-points)
      ! ======================
      IF( ppkth == 0.e0 ) THEN            !  uniform vertical grid       

         za1 = zhmax / FLOAT(jpk-1) 
         DO jk = 1, jpk
            zw = FLOAT( jk )
            zt = FLOAT( jk ) + 0.5
            gdepw_0(jk) = ( zw - 1 ) * za1
            gdept_0(jk) = ( zt - 1 ) * za1
            e3w_0  (jk) =  za1
            e3t_0  (jk) =  za1
         END DO

      ELSE

         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
         WRITE(numout,*)
         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
         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

   END SUBROUTINE zgr_z


   SUBROUTINE zgr_bat
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE zgr_bat  ***
      !! 
      !! ** Purpose :   set bathymetry both in levels and meters
      !!
      !! ** Method  :   read or define mbathy and bathy arrays
      !!       * level bathymetry:
      !!      The ocean basin geometry is given by a two-dimensional array,
      !!      mbathy, which is defined as follow :
      !!            mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
      !!                              at t-point (ji,jj).
      !!                            = 0  over the continental t-point.
      !!                            = -n over the nth island t-point.
      !!      The array mbathy is checked to verified its consistency with
      !!      model option. in particular:
      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
      !!                  along closed boundary.
      !!            mbathy must be cyclic IF jperio=1.
      !!            mbathy must be lower or equal to jpk-1.
      !!            isolated ocean grid points are suppressed from mbathy
      !!                  since they are only connected to remaining
      !!                  ocean through vertical diffusion.
      !!      ntopo=-1 :   rectangular channel or bassin with a bump 
      !!      ntopo= 0 :   flat rectangular channel or basin 
      !!      ntopo= 1 :   mbathy is read in 'bathy_level.nc' NetCDF file
      !!                   bathy  is read in 'bathy_meter.nc' NetCDF file
      !!      C A U T I O N : mbathy will be modified during the initializa-
      !!      tion phase to become the number of non-zero w-levels of a water
      !!      column, with a minimum value of 1.
      !!
      !! ** 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
      !!----------------------------------------------------------------------

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) '    zgr_bat : defines level and meter bathymetry'
      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 ) THEN                        ! flat basin

            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,*) '         bathymetry field: flat basin'

            idta(:,:) = jpkm1                            ! flat basin 
            zdta(:,:) = gdepw_0(jpk)
            h_oce     = gdepw_0(jpk)

         ELSE                                         ! bump
            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)
            IF(lwp) WRITE(numout,*) '            bump characteristics: '
            IF(lwp) WRITE(numout,*) '               bump center (i,j)   = ', ii_bump, ii_bump
            IF(lwp) WRITE(numout,*) '               bump height         = ', h_bump , ' meters'
            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 ji = 1, jpidta
                  zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
                  zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
                  zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) )
               END DO
            END DO

            ! idta :
            IF( ln_sco ) THEN      ! s-coordinate (zsc       ): idta()=jpk
               idta(:,:) = jpkm1
            ELSE                   ! z-coordinate (zco or zps): step-like topography
               idta(:,:) = jpkm1
               DO jk = 1, jpkm1
                  DO jj = 1, jpjdta
                     DO ji = 1, jpidta
                        IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) )   idta(ji,jj) = jk
                     END DO
                  END DO
               END DO
            ENDIF
         ENDIF

         ! set 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
            idta( :    ,jpjdta) =  0                ;      zdta( :    ,jpjdta) =  0.e0
         ELSEIF( jperio == 2 ) THEN
            idta( :    , 1    ) = idta( : ,  3  )   ;      zdta( :    , 1    ) = zdta( : ,  3  )
            idta( :    ,jpjdta) = 0                 ;      zdta( :    ,jpjdta) =  0.e0
            idta( 1    , :    ) = 0                 ;      zdta( 1    , :    ) =  0.e0
            idta(jpidta, :    ) = 0                 ;      zdta(jpidta, :    ) =  0.e0
         ELSE
            ih = 0                                  ;      zh = 0.e0
            IF( ln_sco )   ih = jpkm1               ;      IF( ln_sco )   zh = h_oce
            idta( :    , 1    ) = ih                ;      zdta( :    , 1    ) =  zh
            idta( :    ,jpjdta) = ih                ;      zdta( :    ,jpjdta) =  zh
            idta( 1    , :    ) = ih                ;      zdta( 1    , :    ) =  zh
            idta(jpidta, :    ) = ih                ;      zdta(jpidta, :    ) =  zh
         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
            DO ji = 1, nlci
               mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
               bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
            END DO
         END DO

         !                                            ! =============== !
      ELSEIF( ntopo == 1 ) THEN                       !   read in file  !
         !                                            ! =============== !

         IF( ln_zco )   THEN
            CALL iom_open ( 'bathy_level.nc', inum )   ! Level bathymetry
            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
            CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
            CALL iom_close (inum)
         ENDIF
         !                                            ! =============== !
      ELSE                                            !      error      !
         !                                            ! =============== !
         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
            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
               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?"'
#endif
      ! ===========
      ! Zoom domain 
      ! ===========

      IF( lzoom )   CALL zgr_bat_zoom

      ! ================
      ! Bathymetry check
      ! ================

      CALL zgr_bat_ctl

   END SUBROUTINE zgr_bat


   SUBROUTINE zgr_bat_zoom
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE zgr_bat_zoom  ***
      !!
      !! ** Purpose : - Close zoom domain boundary if necessary
      !!              - Suppress Med Sea from ORCA R2 and R05 arctic zoom
      !!
      !! ** Method  : 
      !!
      !! ** 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

      ! Configuration specific domain modifications
      ! (here, ORCA arctic configuration: suppress Med Sea)
      IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
         SELECT CASE ( jp_cfg )
         !                                        ! =======================
         CASE ( 2 )                               !  ORCA_R2 configuration
            !                                     ! =======================
            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
            !                                     ! =======================
         CASE ( 05 )                              !  ORCA_R05 configuration
            !                                     ! =======================
            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
         !
         mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0   ! zero over the Med Sea boxe
         !
      ENDIF

   END SUBROUTINE zgr_bat_zoom


   SUBROUTINE zgr_bat_ctl
      !!----------------------------------------------------------------------
      !!                    ***  ROUTINE zgr_bat_ctl  ***
      !!
      !! ** Purpose :   check the bathymetry in levels
      !!
      !! ** Method  :   The array mbathy is checked to verified its consistency
      !!      with the model options. in particular:
      !!            mbathy must have at least 1 land grid-points (mbathy<=0)
      !!                  along closed boundary.
      !!            mbathy must be cyclic IF jperio=1.
      !!            mbathy must be lower or equal to jpk-1.
      !!            isolated ocean grid points are suppressed from mbathy
      !!                  since they are only connected to remaining
      !!                  ocean through vertical diffusion.
      !!      C A U T I O N : mbathy will be modified during the initializa-
      !!      tion phase to become the number of non-zero w-levels of a water
      !!      column, with a minimum value of 1.
      !!
      !! ** 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
      !!----------------------------------------------------------------------

      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) '    zgr_bat_ctl : check the bathymetry'
      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'

      ! ================
      ! Bathymetry check
      ! ================

      IF( .NOT. lk_c1d )   THEN

         ! 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,:)
               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

         ! 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
         DO ji = 1, jpi
            ikmax = MAX( ikmax, mbathy(ji,jj) )
         END DO
      END DO
      !!! test a faire:   ikmax = MAX( mbathy(:,:) )   ???

      IF( ikmax > jpkm1 ) THEN
         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
         IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
      ELSE IF( ikmax < jpkm1 ) THEN
         IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
         IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
      ENDIF

      IF( lwp .AND. nprint == 1 ) THEN
         WRITE(numout,*)
         WRITE(numout,*) ' bathymetric field '
         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

   END SUBROUTINE zgr_bat_ctl


   SUBROUTINE zgr_zco
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zgr_zco  ***
      !!
      !! ** Purpose :   define the z-coordinate system
      !!
      !! ** 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
            fsdept(:,:,jk) = gdept_0(jk)
            fsdepw(:,:,jk) = gdepw_0(jk)
            fsde3w(:,:,jk) = gdepw_0(jk)
            fse3t (:,:,jk) = e3t_0(jk)
            fse3u (:,:,jk) = e3t_0(jk)
            fse3v (:,:,jk) = e3t_0(jk)
            fse3f (:,:,jk) = e3t_0(jk)
            fse3w (:,:,jk) = e3w_0(jk)
            fse3uw(:,:,jk) = e3w_0(jk)
            fse3vw(:,:,jk) = e3w_0(jk)
         END DO
      ENDIF

   END SUBROUTINE zgr_zco



#  include "domzgr_zps.h90"


#if ! defined key_zco
   !!----------------------------------------------------------------------
   !!   NOT 'key_zco' :                                    3D gdep. and e3.
   !!----------------------------------------------------------------------

   SUBROUTINE zgr_sco
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zgr_sco  ***
      !!                     
      !! ** Purpose :   define the s-coordinate system
      !!
      !! ** Method  :   s-coordinate
      !!         The depth of model levels is defined as the product of an
      !!      analytical function by the local bathymetry, while the vertical
      !!      scale factors are defined as the product of the first derivative
      !!      of the analytical function by the bathymetry.
      !!      (this solution save memory as depth and scale factors are not
      !!      3d fields)
      !!          - Read bathymetry (in meters) at t-point and compute the
      !!         bathymetry at u-, v-, and f-points.
      !!            hbatu = mi( hbatt )
      !!            hbatv = mj( hbatt )
      !!            hbatf = mi( mj( hbatt ) )
      !!          - Compute gsigt, gsigw, esigt, esigw from an analytical
      !!         function and its derivative given as statement function.
      !!            gsigt(k) = fssig (k+0.5)
      !!            gsigw(k) = fssig (k    )
      !!            esigt(k) = fsdsig(k+0.5)
      !!            esigw(k) = fsdsig(k    )
      !!      This routine is given as an example, it must be modified
      !!      following the user s desiderata. nevertheless, the output as
      !!      well as the way to compute the model levels and scale factors
      !!      must be respected in order to insure second order a!!uracy
      !!      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
      REAL(wp) ::   fssig, fsdsig, pfkk

      INTEGER  ::   iip1, ijp1, iim1, ijm1
      REAL(wp) ::   &
         fssigt, fssigw, zcoeft, zcoefw,   &
         zrmax, ztaper

      REAL(wp), DIMENSION(jpi,jpj) ::   &
         zenv, ztmp, zmsk, zri, zrj, zhbat

      NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
      !!----------------------------------------------------------------------

      ! a. Sigma stretching coefficient
       fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))
       fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &
                    * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))

      ! b. Vertical derivative of sigma stretching coefficient
 
!!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ...
      !!----------------------------------------------------------------------

      ! Read Namelist nam_zgr_sco : sigma-stretching parameters
      ! -------------------------
      REWIND( numnam )
      READ  ( numnam, nam_zgr_sco )

      IF(lwp) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
         WRITE(numout,*) '~~~~~~~~~~~'
         WRITE(numout,*) '         Namelist nam_zgr_sco'
         WRITE(numout,*) '            sigma-stretching coeffs '
         WRITE(numout,*) '            maximum depth of s-bottom surface (>0)    sbot_max = ' ,sbot_max
         WRITE(numout,*) '            minimum depth of s-bottom surface (>0)    sbot_min = ' ,sbot_min
         WRITE(numout,*) '            surface control parameter (0<=theta<=20)  theta    = ', theta
         WRITE(numout,*) '            bottom control parameter  (0<=thetb<= 1)  thetb    = ', thetb
         WRITE(numout,*) '            maximum cut-off r-value allowed           r_max    = ' , r_max
      ENDIF

      ! ???
      ! ------------------
      hift(:,:) = sbot_min
      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
      ! =============================
      ! 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
            zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
         END DO
      END DO

      jl = 0
      zrmax = 1.e0
      !                                                     ! ================ !
      DO WHILE ( jl <= 10000 .AND. zrmax > r_max )          !  Iterative loop  !
         !                                                  ! ================ !
         jl = jl + 1

         zrmax = 0.e0
         zmsk(:,:) = 0.e0
 
         DO jj = 1, nlcj
            DO ji = 1, nlci
               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
               zri(ji,jj) = ABS( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
               zrj(ji,jj) = ABS( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
               zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
               IF( zri(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
               IF( zri(ji,jj) > r_max )   zmsk(iip1,jj  ) = 1.0
               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,jj  ) = 1.0
               IF( zrj(ji,jj) > r_max )   zmsk(ji  ,ijp1) = 1.0
            END DO
         END DO
         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. )
         DO jj = 1, nlcj
            DO ji = 1, nlci
                zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )    
            END DO
         END DO
 
         WRITE(numout,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )

         DO jj = 1, nlcj
            DO ji = 1, nlci
               iip1 = MIN( ji+1, nlci )     ! last  line (ji=nlci)
               ijp1 = MIN( jj+1, nlcj )     ! last  raw  (jj=nlcj)
               iim1 = MAX( ji-1,  1  )      ! first line (ji=nlci)
               ijm1 = MAX( jj-1,  1  )      ! first raw  (jj=nlcj)
               IF( zmsk(ji,jj) == 1.0 ) THEN
                  ztmp(ji,jj) =   (                                                                                   &
             &      zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1)   &
             &    + zenv(iim1,jj  )*zmsk(iim1,jj  ) + zenv(ji,jj  )*    2.e0      + zenv(iip1,jj  )*zmsk(iip1,jj  )   &
             &    + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1)   &
             &                    ) / (                                                                               &
             &                      zmsk(iim1,ijp1) +               zmsk(ji,ijp1) +                 zmsk(iip1,ijp1)   &
             &    +                 zmsk(iim1,jj  ) +                   2.e0      +                 zmsk(iip1,jj  )   &
             &    +                 zmsk(iim1,ijm1) +               zmsk(ji,ijm1) +                 zmsk(iip1,ijm1)   &
             &                        )
               ENDIF
            END DO
         END DO
 
         DO jj = 1, nlcj
            DO ji = 1, nlci
               IF( zmsk(ji,jj) == 1.0 )   zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
            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. )
         DO jj = 1, nlcj
            DO ji = 1, nlci
               IF( zenv(ji,jj) == 0.e0 )   zenv(ji,jj) = ztmp(ji,jj)
            END DO
         END DO
         !                                                  ! ================ !
      END DO                                                !     End loop     !
      !                                                     ! ================ !

      ! save envelop bathymetry 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(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
      ! --------------------------------------------

      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) )
        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. )
      DO jj = 1, jpj
         DO ji = 1, jpi
            IF( hbatu(ji,jj) == 0.e0 ) THEN
               IF( zhbat(ji,jj) == 0.e0 )   hbatu(ji,jj) = sbot_min
               IF( zhbat(ji,jj) /= 0.e0 )   hbatu(ji,jj) = zhbat(ji,jj)
            ENDIF
         END DO
      END DO

      zhbat(:,:) = hbatv(:,:)
      CALL lbc_lnk( hbatv, 'V', 1. )
      DO jj = 1, jpj
         DO ji = 1, jpi
            IF( hbatv(ji,jj) == 0.e0 ) THEN
               IF( zhbat(ji,jj) == 0.e0 )   hbatv(ji,jj) = sbot_min
               IF( zhbat(ji,jj) /= 0.e0 )   hbatv(ji,jj) = zhbat(ji,jj)
            ENDIF
         END DO
      END DO

      zhbat(:,:) = hbatf(:,:)
      CALL lbc_lnk( hbatf, 'F', 1. )
      DO jj = 1, jpj
         DO ji = 1, jpi
            IF( hbatf(ji,jj) == 0.e0 ) THEN
               IF( zhbat(ji,jj) == 0.e0 )   hbatf(ji,jj) = sbot_min
               IF( zhbat(ji,jj) /= 0.e0 )   hbatf(ji,jj) = zhbat(ji,jj)
            ENDIF
         END DO
      END DO


!!bug:  key_helsinki a verifer
      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )

      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 hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
      ENDIF
!! helsinki

      ! 2. Computation of gsig and esig fields
      ! --------------------------------------

      ! 2.1 Coefficients for model level depth at w- and t-levels

!!bug gm : change the def of statment function....
      DO jk = 1, jpk
        gsigw(jk) = -fssigw(FLOAT(jk))
        gsigt(jk) = -fssigt(FLOAT(jk))
      END DO

      IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk    ', gsigw(1), gsigw(jpk)

!!org DO jk = 1, jpk
!!org    gsigw(jk) = -fssig ( FLOAT(jk)    )
!!org    gsigt(jk) = -fssig ( FLOAT(jk)+0.5)
!!org END DO

      ! 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)

!!org DO jk = 1, jpk
!!org    esigw(jk)=fsdsig( FLOAT(jk))
!!org    esigt(jk)=fsdsig( FLOAT(jk)+0.5)
!!org END DO

      ! 2.3 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)
      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)
      DO jj = 1, jpj
         DO ji = 1, jpi
            DO jk = 1, jpk
              e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
              e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
              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))
              e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
            END DO
         END DO
      END DO

! HYBRID
      DO jj = 1, jpj
         DO ji = 1, jpi
            DO jk = 1, jpkm1
               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
               IF( scobot(ji,jj) == 0.e0             )   mbathy(ji,jj) = 0
            END DO
         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
         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,*)
         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)
         DO jj = mj0(20), mj1(20)
            DO ji = mi0(20), mi1(20)
               WRITE(numout,*)
               WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
               WRITE(numout,*) ' ~~~~~~  --------------------'
               WRITE(numout,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
         DO jj = mj0(74), mj1(74)
            DO ji = mi0(100), mi1(100)
               WRITE(numout,*)
               WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
               WRITE(numout,*) ' ~~~~~~  --------------------'
               WRITE(numout,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
      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'
               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'
               ENDIF
            END DO
         END DO
      END DO
!!bug gm    #endif

   END SUBROUTINE zgr_sco

#else
   !!----------------------------------------------------------------------
   !!   Default option :                                      Empty routine
   !!----------------------------------------------------------------------
   SUBROUTINE zgr_sco
   END SUBROUTINE zgr_sco
#endif


   !!======================================================================
END MODULE domzgr