Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6972)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6973)
@@ -106,6 +106,4 @@
       CALL dom_msk( ik_top, ik_bot )   ! Masks
       !
-      IF( ln_sco )  CALL dom_stiff     ! Maximum stiffness ratio/hydrostatic consistency
-      !
       DO jj = 1, jpj                   ! depth of the iceshelves
          DO ji = 1, jpj
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6972)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6973)
@@ -171,5 +171,5 @@
       CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
       
-      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor
+      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )       !    ! coriolis factor
       CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
       
@@ -180,19 +180,22 @@
       CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
       zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )   !    ! nb of ocean T-points
             
-      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0 )         !    ! scale factors
-      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0 )
-      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0 )
-      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0 )
-      !
-      CALL dom_stiff( zprt )
-      CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
-      !
-      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d )  !    ! stretched system
-      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d )
-      CALL iom_rstput( 0, 0, inum, 'gdept_0', gdept_0, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'gdepw_0', gdepw_0, ktype = jp_r8 )
-      
+      !															             ! vertical mesh
+      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0, ktype = jp_r8  )    !    ! scale factors
+      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0, ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0, ktype = jp_r8  )
+      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0, ktype = jp_r8  )
+      !
+      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d , ktype = jp_r8 )  ! stretched system
+      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gdept_0'  , gdept_0  , ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gdepw_0'  , gdepw_0  , ktype = jp_r8 )
+      !
+      IF( ln_sco ) THEN                                         ! s-coordinate stiffness
+         CALL dom_stiff( zprt )
+         CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
+   	ENDIF
+      !
       !                                     ! ============================
       CALL iom_close( inum )                !        close the files 
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/dom_oce.f90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/dom_oce.f90	(revision 6972)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/dom_oce.f90	(revision 6973)
@@ -222,5 +222,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff       !: interface depth between stretching at v--f
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio
 !!gm end
 
@@ -362,5 +361,5 @@
          &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
          &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
-         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1(jpi,jpj) , STAT=ierr(8) )
+         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , STAT=ierr(8) )
 
       ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       &
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domain.f90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domain.f90	(revision 6972)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domain.f90	(revision 6973)
@@ -20,5 +20,4 @@
    !!   dom_nam        : read and contral domain namelists
    !!   dom_ctl        : control print for the ocean domain
-   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
    !!----------------------------------------------------------------------
    USE oce             ! ocean variables
@@ -88,5 +87,4 @@
                      CALL dom_zgr               ! Vertical mesh and bathymetry
                      CALL dom_msk               ! Masks
-      IF( ln_sco )   CALL dom_stiff             ! Maximum stiffness ratio/hydrostatic consistency
       !
       ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1)   ! Reference ocean thickness
@@ -370,65 +368,4 @@
 
 
-   SUBROUTINE dom_stiff
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_stiff  ***
-      !!                     
-      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
-      !!
-      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
-      !!                Save the maximum in the vertical direction
-      !!                (this number is only relevant in s-coordinates)
-      !!
-      !!                Haney, R. L., 1991: On the pressure gradient force
-      !!                over steep topography in sigma coordinate ocean models. 
-      !!                J. Phys. Oceanogr., 21, 610???619.
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk 
-      REAL(wp) ::   zrxmax
-      REAL(wp), DIMENSION(4) ::   zr1
-      !!----------------------------------------------------------------------
-      rx1(:,:) = 0._wp
-      zrxmax   = 0._wp
-      zr1(:)   = 0._wp
-      !
-      DO ji = 2, jpim1
-         DO jj = 2, jpjm1
-            DO jk = 1, jpkm1
-               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
-                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
-                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
-                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
-               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
-                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
-                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
-                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
-               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
-                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
-                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
-                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
-               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
-                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
-                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
-                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
-               zrxmax = MAXVAL( zr1(1:4) )
-               rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( rx1, 'T', 1. )
-      !
-      zrxmax = MAXVAL( rx1 )
-      !
-      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
-         WRITE(numout,*) '~~~~~~~~~'
-      ENDIF
-      !
-   END SUBROUTINE dom_stiff
-
-
    SUBROUTINE cfg_write
       !!----------------------------------------------------------------------
@@ -537,11 +474,6 @@
       CALL iom_rstput( 0, 0, inum, 'top    level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF)
       !
-      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio (Not required anyway)
-        !SF  CALL dom_stiff( z2d ) !commented because at compilation error:  The
-                                   !number of actual arguments cannot be greater than the number of dummy
-                                   !arguments.   [DOM_STIFF]
-                                   !CALL dom_stiff( z2d )
-                                   ! --------------^ compilation aborted for
-                                   !/workgpfs/rech/gzi/rgzi011/commit-simplif2-TOOLS/NEMOGCM/CONFIG/TEST/BLD/ppsrc/nemo/domain.f90
+      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio  (Not required anyway)
+         CALL dom_stiff( z2d )
          CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
       ENDIF
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domwri.f90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domwri.f90	(revision 6972)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/TOOLS/DOMAINcfg/src/domwri.f90	(revision 6973)
@@ -13,4 +13,5 @@
    !!   dom_wri        : create and write mesh and mask file(s)
    !!   dom_uniq       : identify unique point of a grid (TUVF)
+   !!   dom_stiff      : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)
    !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
@@ -27,20 +28,10 @@
    PUBLIC   dom_wri              ! routine called by inidom.F90
    PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90
-   !! * Substitutions
-   !!----------------------------------------------------------------------
-   !!                   ***  vectopt_loop_substitute  ***
-   !!----------------------------------------------------------------------
-   !! ** purpose :   substitute the inner loop start/end indices with CPP macro
-   !!                allow unrolling of do-loop (useful with vector processors)
-   !!----------------------------------------------------------------------
+   PUBLIC   dom_stiff            ! routine called by inidom.F90
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $ 
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: domwri.F90 6689 2016-06-13 14:24:43Z flavoni $ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -286,5 +277,4 @@
          CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
          CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
-         CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 )             !    ! Max. grid stiffness ratio
          !
          CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
@@ -292,4 +282,6 @@
          CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )     
          CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
+	      CALL dom_stiff( zprt )
+   	   CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )       !    ! Max. grid stiffness ratio
       ENDIF
       
@@ -416,4 +408,72 @@
    END SUBROUTINE dom_uniq
 
+
+   SUBROUTINE dom_stiff( px1 )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_stiff  ***
+      !!                     
+      !! ** Purpose :   Diagnose maximum grid stiffness/hydrostatic consistency
+      !!
+      !! ** Method  :   Compute Haney (1991) hydrostatic condition ratio
+      !!                Save the maximum in the vertical direction
+      !!                (this number is only relevant in s-coordinates)
+      !!
+      !!                Haney, 1991, J. Phys. Oceanogr., 21, 610-619.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL ::   px1   ! stiffness
+      !
+      INTEGER  ::   ji, jj, jk 
+      REAL(wp) ::   zrxmax
+      REAL(wp), DIMENSION(4) ::   zr1
+      REAL(wp), DIMENSION(jpi,jpj) ::   zx1
+      !!----------------------------------------------------------------------
+      zx1(:,:) = 0._wp
+      zrxmax   = 0._wp
+      zr1(:)   = 0._wp
+      !
+      DO ji = 2, jpim1
+         DO jj = 2, jpjm1
+            DO jk = 1, jpkm1
+!!gm   remark: dk(gdepw) = e3t   ===>>>  possible simplification of the following calculation....
+!!             especially since it is gde3w which is used to compute the pressure gradient
+!!             furthermore, I think gdept_0 should be used below instead of w point in the numerator
+!!             so that the ratio is computed at the same point (i.e. uw and vw) ....
+               zr1(1) = ABS(  ( gdepw_0(ji  ,jj,jk  )-gdepw_0(ji-1,jj,jk  )               & 
+                    &          +gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) )             &
+                    &       / ( gdepw_0(ji  ,jj,jk  )+gdepw_0(ji-1,jj,jk  )               &
+                    &          -gdepw_0(ji  ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall )  ) * umask(ji-1,jj,jk)
+               zr1(2) = ABS(  ( gdepw_0(ji+1,jj,jk  )-gdepw_0(ji  ,jj,jk  )               &
+                    &          +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) )             &
+                    &       / ( gdepw_0(ji+1,jj,jk  )+gdepw_0(ji  ,jj,jk  )               &
+                    &          -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji  ,jj,jk+1) + rsmall )  ) * umask(ji  ,jj,jk)
+               zr1(3) = ABS(  ( gdepw_0(ji,jj+1,jk  )-gdepw_0(ji,jj  ,jk  )               &
+                    &          +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj+1,jk  )+gdepw_0(ji,jj  ,jk  )               &
+                    &          -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj  ,jk+1) + rsmall )  ) * vmask(ji,jj  ,jk)
+               zr1(4) = ABS(  ( gdepw_0(ji,jj  ,jk  )-gdepw_0(ji,jj-1,jk  )               &
+                    &          +gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) )             &
+                    &       / ( gdepw_0(ji,jj  ,jk  )+gdepw_0(ji,jj-1,jk  )               &
+                    &          -gdepw_0(ji,jj  ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall )  ) * vmask(ji,jj-1,jk)
+               zrxmax = MAXVAL( zr1(1:4) )
+               zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zx1, 'T', 1. )
+      !
+      IF( PRESENT( px1 ) )    px1 = zx1
+      !
+      zrxmax = MAXVAL( zx1 )
+      !
+      IF( lk_mpp )   CALL mpp_max( zrxmax ) ! max over the global domain
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax
+         WRITE(numout,*) '~~~~~~~~~'
+      ENDIF
+      !
+   END SUBROUTINE dom_stiff
+
    !!======================================================================
 END MODULE domwri
