Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 7217)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 7806)
@@ -17,4 +17,5 @@
    USE lib_mpp         ! distributed memory computing library
 
+   USE iom
    USE domstp          ! domain: set the time-step
 
@@ -73,10 +74,8 @@
 
       CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+      CALL dom_msk      ! Masks
+      CALL dom_hgr      ! Horizontal grid
       CALL dom_zgr      ! Vertical mesh and bathymetry option
-      CALL dom_grd      ! Create a domain file
-
-     !
-      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
-      !        but could be usefull in many other routines
+      !
       e12t    (:,:) = e1t(:,:) * e2t(:,:)
       e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
@@ -91,16 +90,5 @@
       re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
       !
-      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
-      hv(:,:) = 0._wp
-      DO jk = 1, jpk
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
-
       CALL dom_stp      ! Time step
-      CALL dom_msk      ! Masks
       CALL dom_ctl      ! Domain control
 
@@ -178,6 +166,4 @@
       nstocklist = nn_stocklist
       nwrite = nn_write
-
-
       !                             ! control of output frequency
       IF ( nstock == 0 .OR. nstock > nitend ) THEN
@@ -222,4 +208,5 @@
 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namdom )
+
 
       IF(lwp) THEN
@@ -321,10 +308,259 @@
    END SUBROUTINE dom_nam
 
+   SUBROUTINE dom_msk
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean mask informations and defines the interior domain T-mask.
+      !!
+      !! ** Method  :  Read in a file all the arrays generated in routines
+      !!               dommsk:   'mask.nc' file
+      !!              The interior ocean/land mask is computed from tmask
+      !!              setting to zero the duplicated row and lines due to
+      !!              MPP exchange halos, est-west cyclic and north fold
+      !!              boundary conditions.
+      !!
+      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
+      !!               tpol     : ???
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::  inum   ! local integers
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
+      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
+      !
+      !!---------------------------------------------------------------------
+      
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      CALL wrk_alloc( jpi, jpj, zmbk )
+      zmbk(:,:) = 0._wp
+
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mask', inum )
+
+         !                                                         ! masks (inum2) 
+      CALL iom_get( inum, jpdom_data, 'tmask', tmask )
+      CALL iom_get( inum, jpdom_data, 'umask', umask )
+      CALL iom_get( inum, jpdom_data, 'vmask', vmask )
+      CALL iom_get( inum, jpdom_data, 'fmask', fmask )
+
+      CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
+      CALL lbc_lnk( umask, 'U', 1._wp )      
+      CALL lbc_lnk( vmask, 'V', 1._wp )
+      CALL lbc_lnk( fmask, 'F', 1._wp )
+
+#if defined key_c1d
+      ! set umask and vmask equal tmask in 1D configuration
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
+      IF(lwp) WRITE(numout,*) '**********                                                     ********'
+
+      umask(:,:,:) = tmask(:,:,:)
+      vmask(:,:,:) = tmask(:,:,:)
+#endif
+
+#if defined key_degrad
+      CALL iom_get( inum, jpdom_data, 'facvolt', facvol )
+#endif
+
+      CALL iom_get( inum, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
+      mbathy (:,:) = INT( zmbk(:,:) )
+      misfdep(:,:) = 1                                               ! ice shelf case not yet done
+      
+      CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
+
+      !                                     ! ============================
+      !                                     !        close the files 
+      !                                     ! ============================
+
+      !
+      ! Interior domain mask (used for global sum)
+      ! --------------------
+      ssmask(:,:)  = tmask(:,:,1)
+      tmask_i(:,:) = tmask(:,:,1)
+      iif = jpreci                        ! thickness of exchange halos in i-axis
+      iil = nlci - jpreci + 1
+      ijf = jprecj                        ! thickness of exchange halos in j-axis
+      ijl = nlcj - jprecj + 1
+      !
+      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
+      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
+      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
+      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
+      !
+      !                                   ! north fold mask
+      tpol(1:jpiglo) = 1._wp
+      !                                
+      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
+      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
+      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
+         IF( mjg(ijl-1) == jpjglo-1 ) THEN
+            DO ji = iif+1, iil-1
+               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
+            END DO
+         ENDIF
+      ENDIF 
+      !
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
+      ! least 1 wet u point
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector loop
+            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+         END DO
+         DO ji = 1, jpim1      ! NO vector opt.
+            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
+               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
+         END DO
+      END DO
+      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask_i, 'V', 1._wp )
+      CALL lbc_lnk( fmask_i, 'F', 1._wp )
+
+      ! 3. Ocean/land mask at wu-, wv- and w points 
+      !----------------------------------------------
+      wmask (:,:,1) = tmask(:,:,1) ! ????????
+      wumask(:,:,1) = umask(:,:,1) ! ????????
+      wvmask(:,:,1) = vmask(:,:,1) ! ????????
+      DO jk = 2, jpk
+         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
+      END DO
+      !
+      CALL wrk_dealloc( jpi, jpj, zmbk )
+      !
+      CALL iom_close( inum )
+      !
+   END SUBROUTINE dom_msk
+
+   SUBROUTINE zgr_bot_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_bot_level  ***
+      !!
+      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
+      !!
+      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
+      !!
+      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
+      !!                                     ocean level at t-, u- & v-points
+      !!                                     (min value = 1 over land)
+      !!----------------------------------------------------------------------
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
+      !!----------------------------------------------------------------------
+
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
+      !
+      CALL wrk_alloc( jpi, jpj, zmbk )
+      !
+      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
+      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
+      !                                     ! bottom k-index of W-level = mbkt+1
+      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
+         DO ji = 1, jpim1
+            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
+            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
+         END DO
+      END DO
+      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
+      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
+      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zmbk )
+      !
+   END SUBROUTINE zgr_bot_level
+
+   SUBROUTINE dom_hgr
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_hgr  ***
+      !!                   
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean horizontal mesh informations 
+      !!
+      !! ** Method  :   Read in a file all the arrays generated in routines
+      !!                domhgr:   'mesh_hgr.nc' file
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER ::   ji, jj   ! dummy loop indices
+      INTEGER  ::  inum    ! local integers
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mesh_hgr', inum )
+
+      !                                                         ! horizontal mesh (inum3)
+      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
+      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
+      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
+      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
+
+      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
+      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
+      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
+      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
+
+      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
+      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
+      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
+      
+      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
+      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
+      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
+
+      CALL iom_get( inum, jpdom_data, 'ff', ff )
+
+
+      ! Control printing : Grid informations (if not restart)
+      ! ----------------
+
+      IF(lwp .AND. .NOT.ln_rstart ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '          longitude and e1 scale factors'
+         WRITE(numout,*) '          ------------------------------'
+         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
+            glamv(ji,1), glamf(ji,1),   &
+            e1t(ji,1), e1u(ji,1),   &
+            e1v(ji,1), ji = 1, jpi,10)
+
+         WRITE(numout,*)
+         WRITE(numout,*) '          latitude and e2 scale factors'
+         WRITE(numout,*) '          -----------------------------'
+         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
+            &                     gphiv(1,jj), gphif(1,jj),   &
+            &                     e2t  (1,jj), e2u  (1,jj),   &
+            &                     e2v  (1,jj), jj = 1, jpj, 10 )
+      ENDIF
+
+      !                                     ! ============================
+      !                                     !        close the files 
+      !                                     ! ============================
+      CALL iom_close( inum )
+      !
+9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
+            f19.10, 1x, f19.10, 1x, f19.10 )
+   END SUBROUTINE dom_hgr
+
+
    SUBROUTINE dom_zgr
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_zgr  ***
       !!                   
-      !! ** Purpose :  set the depth of model levels and the resulting 
-      !!      vertical scale factors.
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean horizontal mesh informations and/or set the depth of model levels 
+      !!      and the resulting vertical scale factors.
       !!
       !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
@@ -338,8 +574,11 @@
       !! ** Action  :   define gdep., e3., mbathy and bathy
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio = 0   ! temporary integer
-      INTEGER ::   ios
+      INTEGER  ::  ioptio = 0   ! temporary integer
+      INTEGER  ::  inum, ios
+      INTEGER  ::  ji, jj, jk, ik
+      REAL(wp) ::  zrefdep
       !!
       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw
       !!----------------------------------------------------------------------
 
@@ -372,234 +611,74 @@
       IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
 
-   END SUBROUTINE dom_zgr
-
-   SUBROUTINE dom_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_ctl  ***
-      !!
-      !! ** Purpose :   Domain control.
-      !!
-      !! ** Method  :   compute and print extrema of masked scale factors
-      !!
-      !! History :
-      !!   8.5  !  02-08  (G. Madec)    Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
-      INTEGER, DIMENSION(2) ::   iloc      ! 
-      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
-      !!----------------------------------------------------------------------
-
-      ! Extrema of the scale factors
-
-      IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
-      IF(lwp)WRITE(numout,*) '~~~~~~~'
-
-      IF (lk_mpp) THEN
-         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
-         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
-         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
-         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mesh_zgr', inum )
+
+      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
+      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d )
+      IF( ln_zco .OR. ln_zps ) THEN
+         CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
+         CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
+      ENDIF
+
+!!gm BUG in s-coordinate this does not work!
+      ! deepest/shallowest W level Above/Below ~10m
+      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
+      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
+      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
+!!gm end bug
+
+      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_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
+      ENDIF
+
+      DO jk = 1, jpk
+         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
+         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
+      END DO
+
+      IF( lk_vvl ) THEN
+          CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) )
+          ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at  T-points
+          DO jk = 1, jpk
+             ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
+          END DO
       ELSE
-         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-
-         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi1 = iloc(1) + nimpp - 1
-         ijmi1 = iloc(2) + njmpp - 1
-         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi2 = iloc(1) + nimpp - 1
-         ijmi2 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima1 = iloc(1) + nimpp - 1
-         ijma1 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima2 = iloc(1) + nimpp - 1
-         ijma2 = iloc(2) + njmpp - 1
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
-         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
-         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
-         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
-      ENDIF
-
-   END SUBROUTINE dom_ctl
-
-   SUBROUTINE dom_grd
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_grd  ***
-      !!                   
-      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
-      !!      ocean domain informations (mesh and mask arrays). This (these)
-      !!      file(s) is (are) used for visualisation (SAXO software) and
-      !!      diagnostic computation.
-      !!
-      !! ** Method  :   Read in a file all the arrays generated in routines
-      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
-      !!      the vertical coord. used (z-coord, partial steps, s-coord)
-      !!                    nmsh = 1  :   'mesh_mask.nc' file
-      !!                         = 2  :   'mesh.nc' and mask.nc' files
-      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
-      !!                                  'mask.nc' files
-      !!      For huge size domain, use option 2 or 3 depending on your 
-      !!      vertical coordinate.
-      !!
-      !! ** input file : 
-      !!      meshmask.nc  : domain size, horizontal grid-point position,
-      !!                     masks, depth and vertical scale factors
-      !!----------------------------------------------------------------------
-      USE iom
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
-      REAL(wp) ::   zrefdep         ! local real
-      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
-
-      CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw )
-
-      zmbk(:,:) = 0._wp
-
-      SELECT CASE (nmsh)
-         !                                     ! ============================
-         CASE ( 1 )                            !  create 'mesh_mask.nc' file
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
-            CALL iom_open( 'mesh_mask', inum0 )
-
-            inum2 = inum0                                            ! put all the informations
-            inum3 = inum0                                            ! in unit inum0
-            inum4 = inum0
-
-            !                                  ! ============================
-         CASE ( 2 )                            !  create 'mesh.nc' and 
-            !                                  !         'mask.nc' files
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" '
-            CALL iom_open( 'mesh', inum1 )
-            CALL iom_open( 'mask', inum2 )
-
-            inum3 = inum1                                            ! put mesh informations 
-            inum4 = inum1                                            ! in unit inum1 
-
-            !                                  ! ============================
-         CASE ( 3 )                            !  create 'mesh_hgr.nc'
-            !                                  !         'mesh_zgr.nc' and
-            !                                  !         'mask.nc'     files
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" '
-            CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc'
-            CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc'
-            CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc'
-
-            !                                  ! ===========================
-         CASE DEFAULT                          ! return error 
-            !                                  ! mesh has to be provided
-            !                                  ! ===========================
-            CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   &
-            &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' )
-
-         END SELECT
-
-         !                                                         ! masks (inum2) 
-         CALL iom_get( inum2, jpdom_data, 'tmask', tmask )
-         CALL iom_get( inum2, jpdom_data, 'umask', umask )
-         CALL iom_get( inum2, jpdom_data, 'vmask', vmask )
-         CALL iom_get( inum2, jpdom_data, 'fmask', fmask )
-
-         CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
-         CALL lbc_lnk( umask, 'U', 1._wp )      
-         CALL lbc_lnk( vmask, 'V', 1._wp )
-         CALL lbc_lnk( fmask, 'F', 1._wp )
-
-#if defined key_c1d
-         ! set umask and vmask equal tmask in 1D configuration
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
-         IF(lwp) WRITE(numout,*) '**********                                                     ********'
-
-         umask(:,:,:) = tmask(:,:,:)
-         vmask(:,:,:) = tmask(:,:,:)
-#endif
-
-#if defined key_degrad
-         CALL iom_get( inum2, jpdom_data, 'facvolt', facvol )
-#endif
-
-         !                                                         ! horizontal mesh (inum3)
-         CALL iom_get( inum3, jpdom_data, 'glamt', glamt )
-         CALL iom_get( inum3, jpdom_data, 'glamu', glamu )
-         CALL iom_get( inum3, jpdom_data, 'glamv', glamv )
-         CALL iom_get( inum3, jpdom_data, 'glamf', glamf )
-
-         CALL iom_get( inum3, jpdom_data, 'gphit', gphit )
-         CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu )
-         CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv )
-         CALL iom_get( inum3, jpdom_data, 'gphif', gphif )
-
-         CALL iom_get( inum3, jpdom_data, 'e1t', e1t )
-         CALL iom_get( inum3, jpdom_data, 'e1u', e1u )
-         CALL iom_get( inum3, jpdom_data, 'e1v', e1v )
-         
-         CALL iom_get( inum3, jpdom_data, 'e2t', e2t )
-         CALL iom_get( inum3, jpdom_data, 'e2u', e2u )
-         CALL iom_get( inum3, jpdom_data, 'e2v', e2v )
-
-         CALL iom_get( inum3, jpdom_data, 'ff', ff )
-
-         CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
-         mbathy (:,:) = INT( zmbk(:,:) )
-         misfdep(:,:) = 1                                               ! ice shelf case not yet done
-         
-         CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
-         !
          IF( ln_sco ) THEN                                         ! s-coordinate
-            CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt )
-            CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu )
-            CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv )
-            CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf )
+            CALL iom_get( inum, jpdom_data, 'hbatt', hbatt )
+            CALL iom_get( inum, jpdom_data, 'hbatu', hbatu )
+            CALL iom_get( inum, jpdom_data, 'hbatv', hbatv )
+            CALL iom_get( inum, jpdom_data, 'hbatf', hbatf )
             
-            CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
-            CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw )
-            CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w ) 
-            CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt )
-            CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw )
-
-            CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
-            CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
-            CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
-            CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
-
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
+            CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
+            CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw )
+            CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w ) 
+            CALL iom_get( inum, jpdom_unknown, 'esigt', esigt )
+            CALL iom_get( inum, jpdom_unknown, 'esigw', esigw )
+
+            CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
+            CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
+            CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
+            CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
          ENDIF
 
  
          IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
-            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
-            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
             !
-            IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
-               CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
+            IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
             ELSE                                                        ! 2D bottom scale factors
-               CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp )
-               CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp )
+               CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp )
+               CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp )
                !                                                        ! deduces the 3D scale factors
                DO jk = 1, jpk
@@ -633,10 +712,12 @@
             END IF
 
-            IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
-               CALL iom_get( inum4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
+            IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
+               CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
             ELSE                                                           ! 2D bottom depth
-               CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
-               CALL iom_get( inum4, jpdom_data, 'hdepw', zprw )
+               CALL wrk_alloc( jpi, jpj, zprt, zprw )
+               !
+               CALL iom_get( inum, jpdom_data, 'hdept', zprt )
+               CALL iom_get( inum, jpdom_data, 'hdepw', zprw )
                !
                DO jk = 1, jpk                                              ! deduces the 3D depth
@@ -654,4 +735,5 @@
                   END DO
                END DO
+               CALL wrk_dealloc( jpi, jpj, zprt, zprw )
             ENDIF
             !
@@ -659,8 +741,4 @@
 
          IF( ln_zco ) THEN           ! Vertical coordinates and scales factors
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
-            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )
-            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
             DO jk = 1, jpk
                fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
@@ -672,212 +750,70 @@
             END DO
          ENDIF
-
-!!gm BUG in s-coordinate this does not work!
-      ! deepest/shallowest W level Above/Below ~10m
-      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
-      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
-      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
-!!gm end bug
-
-      ! Control printing : Grid informations (if not restart)
-      ! ----------------
-
-      IF(lwp .AND. .NOT.ln_rstart ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '          longitude and e1 scale factors'
-         WRITE(numout,*) '          ------------------------------'
-         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
-            glamv(ji,1), glamf(ji,1),   &
-            e1t(ji,1), e1u(ji,1),   &
-            e1v(ji,1), ji = 1, jpi,10)
-9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
-            f19.10, 1x, f19.10, 1x, f19.10 )
-
-         WRITE(numout,*)
-         WRITE(numout,*) '          latitude and e2 scale factors'
-         WRITE(numout,*) '          -----------------------------'
-         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
-            &                     gphiv(1,jj), gphif(1,jj),   &
-            &                     e2t  (1,jj), e2u  (1,jj),   &
-            &                     e2v  (1,jj), jj = 1, jpj, 10 )
-      ENDIF
-
-
-      IF( nprint == 1 .AND. lwp ) THEN
-         WRITE(numout,*) '          e1u e2u '
-         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         WRITE(numout,*) '          e1v e2v  '
-         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-      ENDIF
-
-      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_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
-      ENDIF
-
-      DO jk = 1, jpk
-         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
-         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
-      END DO
+         !
+      ENDIF
       !                                     ! ============================
       !                                     !        close the files 
       !                                     ! ============================
-      SELECT CASE ( nmsh )
-         CASE ( 1 )                
-            CALL iom_close( inum0 )
-         CASE ( 2 )
-            CALL iom_close( inum1 )
-            CALL iom_close( inum2 )
-         CASE ( 3 )
-            CALL iom_close( inum2 )
-            CALL iom_close( inum3 )
-            CALL iom_close( inum4 )
-      END SELECT
-      !
-      CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
-      !
-   END SUBROUTINE dom_grd
-
-
-   SUBROUTINE zgr_bot_level
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bot_level  ***
-      !!
-      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
-      !!
-      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
-      !!
-      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
-      !!                                     ocean level at t-, u- & v-points
-      !!                                     (min value = 1 over land)
-      !!----------------------------------------------------------------------
-      !
-      INTEGER ::   ji, jj   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
-      !!----------------------------------------------------------------------
-
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
-      !
-      CALL wrk_alloc( jpi, jpj, zmbk )
-      !
-      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
-      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
-      !                                     ! bottom k-index of W-level = mbkt+1
-      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
-         DO ji = 1, jpim1
-            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
-            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
-         END DO
-      END DO
-      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
-      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
-      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
-      !
-      CALL wrk_dealloc( jpi, jpj, zmbk )
-      !
-   END SUBROUTINE zgr_bot_level
-
-   SUBROUTINE dom_msk
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE dom_msk  ***
-      !!
-      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
-      !!
-      !! ** Method  :   The interior ocean/land mask is computed from tmask
-      !!              setting to zero the duplicated row and lines due to
-      !!              MPP exchange halos, est-west cyclic and north fold
-      !!              boundary conditions.
-      !!
-      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
-      !!               tpol     : ???
-      !!----------------------------------------------------------------------
-      !
-      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
-      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
-      !
-      !!---------------------------------------------------------------------
-      
-      CALL wrk_alloc( jpi, jpj, imsk )
-      !
-      ! Interior domain mask (used for global sum)
-      ! --------------------
-      ssmask(:,:)  = tmask(:,:,1)
-      tmask_i(:,:) = tmask(:,:,1)
-      iif = jpreci                        ! thickness of exchange halos in i-axis
-      iil = nlci - jpreci + 1
-      ijf = jprecj                        ! thickness of exchange halos in j-axis
-      ijl = nlcj - jprecj + 1
-      !
-      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
-      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
-      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
-      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
-      !
-      !                                   ! north fold mask
-      tpol(1:jpiglo) = 1._wp
-      !                                
-      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
-      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
-      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
-         IF( mjg(ijl-1) == jpjglo-1 ) THEN
-            DO ji = iif+1, iil-1
-               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
-            END DO
-         ENDIF
-      ENDIF 
-      !
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
-      ! least 1 wet u point
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector loop
-            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
-         END DO
-         DO ji = 1, jpim1      ! NO vector opt.
-            fmask_i(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
-               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
-         END DO
-      END DO
-      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask_i, 'V', 1._wp )
-      CALL lbc_lnk( fmask_i, 'F', 1._wp )
-
-      ! 3. Ocean/land mask at wu-, wv- and w points 
-      !----------------------------------------------
-      wmask (:,:,1) = tmask(:,:,1) ! ????????
-      wumask(:,:,1) = umask(:,:,1) ! ????????
-      wvmask(:,:,1) = vmask(:,:,1) ! ????????
-      DO jk=2,jpk
-         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
-      END DO
-      !
-      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
-         imsk(:,:) = INT( tmask_i(:,:) )
-         WRITE(numout,*) ' tmask_i : '
-         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         WRITE (numout,*)
-         WRITE (numout,*) ' dommsk: tmask for each level'
-         WRITE (numout,*) ' ----------------------------'
-         DO jk = 1, jpk
-            imsk(:,:) = INT( tmask(:,:,jk) )
-            WRITE(numout,*)
-            WRITE(numout,*) ' level = ',jk
-            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, imsk )
-      !
-   END SUBROUTINE dom_msk
+      CALL iom_close( inum )
+      !
+      !
+   END SUBROUTINE dom_zgr
+
+   SUBROUTINE dom_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_ctl  ***
+      !!
+      !! ** Purpose :   Domain control.
+      !!
+      !! ** Method  :   compute and print extrema of masked scale factors
+      !!
+      !! History :
+      !!   8.5  !  02-08  (G. Madec)    Original code
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
+      INTEGER, DIMENSION(2) ::   iloc      ! 
+      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
+      !!----------------------------------------------------------------------
+
+      ! Extrema of the scale factors
+
+      IF(lwp)WRITE(numout,*)
+      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
+      IF(lwp)WRITE(numout,*) '~~~~~~~'
+
+      IF (lk_mpp) THEN
+         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      ELSE
+         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi1 = iloc(1) + nimpp - 1
+         ijmi1 = iloc(2) + njmpp - 1
+         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi2 = iloc(1) + nimpp - 1
+         ijmi2 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima1 = iloc(1) + nimpp - 1
+         ijma1 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima2 = iloc(1) + nimpp - 1
+         ijma2 = iloc(2) + njmpp - 1
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
+         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
+         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
+         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
+      ENDIF
+
+   END SUBROUTINE dom_ctl
 
    !!======================================================================
