Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6579)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6593)
@@ -28,10 +28,10 @@
    USE par_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
+   USE usrdef         ! User defined routine
+!   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
    !
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
-   USE usrdef        ! User defined routine
 
    IMPLICIT NONE
@@ -57,62 +57,37 @@
       !!      the Coriolis factor (in s-1).
       !!
-      !! ** Method  :   The geographical position of the model grid-points is
-      !!      defined from analytical functions, fslam and fsphi, the deriva-
-      !!      tives of which gives the horizontal scale factors e1,e2.
-      !!      Defining two function fslam and fsphi and their derivatives in 
-      !!      the two horizontal directions (fse1 and fse2), the model grid-
-      !!      point position and scale factors are given by:
-      !!         t-point:
-      !!      glamt(i,j) = fslam(i    ,j    )   e1t(i,j) = fse1(i    ,j    )
-      !!      gphit(i,j) = fsphi(i    ,j    )   e2t(i,j) = fse2(i    ,j    )
-      !!         u-point:
-      !!      glamu(i,j) = fslam(i+1/2,j    )   e1u(i,j) = fse1(i+1/2,j    )
-      !!      gphiu(i,j) = fsphi(i+1/2,j    )   e2u(i,j) = fse2(i+1/2,j    )
-      !!         v-point:
-      !!      glamv(i,j) = fslam(i    ,j+1/2)   e1v(i,j) = fse1(i    ,j+1/2)
-      !!      gphiv(i,j) = fsphi(i    ,j+1/2)   e2v(i,j) = fse2(i    ,j+1/2)
-      !!            f-point:
-      !!      glamf(i,j) = fslam(i+1/2,j+1/2)   e1f(i,j) = fse1(i+1/2,j+1/2)
-      !!      gphif(i,j) = fsphi(i+1/2,j+1/2)   e2f(i,j) = fse2(i+1/2,j+1/2)
-      !!      Where fse1 and fse2 are defined by:
-      !!         fse1(i,j) = ra * rad * SQRT( (cos(phi) di(fslam))**2
-      !!                                     +          di(fsphi) **2 )(i,j)
-      !!         fse2(i,j) = ra * rad * SQRT( (cos(phi) dj(fslam))**2
-      !!                                     +          dj(fsphi) **2 )(i,j)
-      !!
-      !!        The coriolis factor is given at z-point by:
-      !!                     ff = 2.*omega*sin(gphif)      (in s-1)
-      !!
-      !!        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 grid-point position and
-      !!      horizontal scale factors must be respected in order to insure
-      !!      second order accuracy schemes.
-      !!
-      !! N.B. If the domain is periodic, verify that scale factors are also
-      !!      periodic, and the coriolis term again.
-      !!
-      !! ** Action  : - define  glamt, glamu, glamv, glamf: longitude of t-, 
-      !!                u-, v- and f-points (in degre)
-      !!              - define  gphit, gphiu, gphiv, gphit: latitude  of t-,
-      !!               u-, v-  and f-points (in degre)
-      !!        define e1t, e2t, e1u, e2u, e1v, e2v, e1f, e2f: horizontal
-      !!      scale factors (in meters) at t-, u-, v-, and f-points.
-      !!        define ff: coriolis factor at f-point
-      !!
-      !! References :   Marti, Madec and Delecluse, 1992, JGR
-      !!                Madec, Imbard, 1996, Clim. Dyn.
+      !! ** Method  :   Controlled by ln_read_cfg logical
+      !!              =T : all needed arrays are read in mesh_mask.nc file 
+      !!              =F : user-defined configuration, all needed arrays 
+      !!                   are computed in usr-def_hgr subroutine 
+      !!
+      !!                If Coriolis factor is neither read nor computed (iff=0)
+      !!              it is computed from gphit assuming that the mesh is
+      !!              defined on the sphere :
+      !!                   ff = 2.*omega*sin(gphif)      (in s-1)
+      !!    
+      !!                If u- & v-surfaces are neither read nor computed (ie1e2u_v=0)
+      !!              (i.e. no use of reduced scale factors in some straits)
+      !!              they are computed from e1u, e2u, e1v and e2v as:
+      !!                   e1e2u = e1u*e2u   and   e1e2v = e1v*e2v  
+      !!    
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
+      !!              - define Coriolis parameter at f-point                   (in 1/s)
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define associated horizontal metrics at t-, u-, v- and f-points
+      !!                (inverse of scale factors 1/e1 & 1/e2, surface e1*e2, ratios e1/e2 & e2/e1)
       !!----------------------------------------------------------------------
       INTEGER  ::   ji, jj                   ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers
-      INTEGER  ::   ijeq                     ! index of equator T point (used in case 4)
-      REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars
-      REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      -
-      REAL(wp) ::   zphi0, zbeta, znorme     !
-      REAL(wp) ::   zarg, zf0, zminff, zmaxff
-      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
-      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
-      INTEGER  ::   isrow                    ! index for ORCA1 starting row
-      INTEGER  ::   ie1e2u_v                 ! fag for u- & v-surface read in coordinate file or not
+      !INTEGER  ::   ii0, ii1, ij0, ij1, iff  ! temporary integers
+      !INTEGER  ::   ijeq                     ! index of equator T point (used in case 4)
+      !REAL(wp) ::   zti, zui, zvi, zfi       ! local scalars
+      !REAL(wp) ::   ztj, zuj, zvj, zfj       !   -      -
+      !REAL(wp) ::   zphi0, zbeta, znorme     !
+      !REAL(wp) ::   zarg, zf0, zminff, zmaxff
+      !REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
+      !REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
+      !INTEGER  ::   isrow                    ! index for ORCA1 starting row
+      INTEGER  ::   ie1e2u_v                 ! flag for u- & v-surfaces 
+      INTEGER  ::   iff                      ! flag for Coriolis parameter
       !!----------------------------------------------------------------------
       !
@@ -122,193 +97,25 @@
          WRITE(numout,*)
          WRITE(numout,*) 'dom_hgr : define the horizontal mesh from ithe following par_oce parameters '
-         WRITE(numout,*) '~~~~~~~      type of horizontal mesh           jphgr_msh = ', jphgr_msh
-         WRITE(numout,*) '             position of the first row and     ppglam0  = ', ppglam0
-         WRITE(numout,*) '             column grid-point (degrees)       ppgphi0  = ', ppgphi0
-         WRITE(numout,*) '             zonal      grid-spacing (degrees) ppe1_deg = ', ppe1_deg
-         WRITE(numout,*) '             meridional grid-spacing (degrees) ppe2_deg = ', ppe2_deg
-         WRITE(numout,*) '             zonal      grid-spacing (meters)  ppe1_m   = ', ppe1_m  
-         WRITE(numout,*) '             meridional grid-spacing (meters)  ppe2_m   = ', ppe2_m  
-      ENDIF
-      !
-      !
-      SELECT CASE( jphgr_msh )   !  type of horizontal mesh  
-      !
-      CASE ( 0 )                     !==  read in coordinate.nc file  ==!
-         !
+         WRITE(numout,*) '~~~~~~~   '
+      ENDIF
+      !
+      !
+      IF( ln_read_cfg ) THEN        !==  read in mesh_mask.nc file  ==!
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          curvilinear coordinate on the sphere read in "coordinate" file'
-         !
-         ie1e2u_v = 0             ! set to unread e1e2u and e1e2v
-         iff = 0                  ! set to unread iff
-         !
-         CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask
-         !
-         IF( ie1e2u_v == 0 ) THEN      ! e1e2u and e1e2v have not been read: compute them
-            !                          ! e2u and e1v does not include a reduction in some strait: apply reduction
-            e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
-            e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
-         ENDIF
-         !
-      CASE ( 1 )                     !==  geographical mesh on the sphere with regular (in degree) grid-spacing  ==!
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere with regular grid-spacing'
-         IF(lwp) WRITE(numout,*) '          given by ppe1_deg and ppe2_deg' 
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - 1 + njmpp - 1 )
-               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - 1 + njmpp - 1 )
-               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - 1 + njmpp - 1 ) + 0.5
-               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - 1 + njmpp - 1 ) + 0.5
-         ! Longitude
-               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
-               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
-               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
-               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
-         ! Latitude
-               gphit(ji,jj) = ppgphi0 + ppe2_deg * ztj
-               gphiu(ji,jj) = ppgphi0 + ppe2_deg * zuj
-               gphiv(ji,jj) = ppgphi0 + ppe2_deg * zvj
-               gphif(ji,jj) = ppgphi0 + ppe2_deg * zfj
-         ! e1
-               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
-               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
-               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
-               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
-         ! e2
-               e2t(ji,jj) = ra * rad * ppe2_deg
-               e2u(ji,jj) = ra * rad * ppe2_deg
-               e2v(ji,jj) = ra * rad * ppe2_deg
-               e2f(ji,jj) = ra * rad * ppe2_deg
-            END DO
-         END DO
-         !
-      CASE ( 2:3 )                   !==  f- or beta-plane with regular grid-spacing  ==!
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          f- or beta-plane with regular grid-spacing'
-         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m' 
-         !
-         ! Position coordinates (in kilometers)
-         !                          ==========
-         glam0 = 0._wp
-         gphi0 = - ppe2_m * 1.e-3
-         !
-#if defined key_agrif 
-         IF ( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN    ! for EEL6 configuration only
-            IF( .NOT. Agrif_Root() ) THEN
-              glam0  = Agrif_Parent(glam0) + (Agrif_ix())*Agrif_Parent(ppe1_m) * 1.e-3
-              gphi0  = Agrif_Parent(gphi0) + (Agrif_iy())*Agrif_Parent(ppe2_m) * 1.e-3
-              ppe1_m = Agrif_Parent(ppe1_m)/Agrif_Rhox()
-              ppe2_m = Agrif_Parent(ppe2_m)/Agrif_Rhoy()          
-            ENDIF
-         ENDIF
-#endif         
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               glamt(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 )       )
-               glamu(ji,jj) = glam0 + ppe1_m * 1.e-3 * ( REAL( ji - 1 + nimpp - 1 ) + 0.5 )
-               glamv(ji,jj) = glamt(ji,jj)
-               glamf(ji,jj) = glamu(ji,jj)
-               !
-               gphit(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 )       )
-               gphiu(ji,jj) = gphit(ji,jj)
-               gphiv(ji,jj) = gphi0 + ppe2_m * 1.e-3 * ( REAL( jj - 1 + njmpp - 1 ) + 0.5 )
-               gphif(ji,jj) = gphiv(ji,jj)
-            END DO
-         END DO
-         !
-         ! Horizontal scale factors (in meters)
-         !                              ======
-         e1t(:,:) = ppe1_m      ;      e2t(:,:) = ppe2_m
-         e1u(:,:) = ppe1_m      ;      e2u(:,:) = ppe2_m
-         e1v(:,:) = ppe1_m      ;      e2v(:,:) = ppe2_m
-         e1f(:,:) = ppe1_m      ;      e2f(:,:) = ppe2_m
-         !
-      CASE ( 4 )                     !==  geographical mesh on the sphere, isotropic MERCATOR type  ==!
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          geographical mesh on the sphere, MERCATOR type'
-         IF(lwp) WRITE(numout,*) '          longitudinal/latitudinal spacing given by ppe1_deg'
-         IF ( ppgphi0 == -90 ) CALL ctl_stop( ' Mercator grid cannot start at south pole !!!! ' )
-         !
-         !  Find index corresponding to the equator, given the grid spacing e1_deg
-         !  and the (approximate) southern latitude ppgphi0.
-         !  This way we ensure that the equator is at a "T / U" point, when in the domain.
-         !  The formula should work even if the equator is outside the domain.
-         zarg = rpi / 4. - rpi / 180. * ppgphi0 / 2.
-         ijeq = ABS( 180./rpi * LOG( COS( zarg ) / SIN( zarg ) ) / ppe1_deg )
-         IF(  ppgphi0 > 0 )  ijeq = -ijeq
-         !
-         IF(lwp) WRITE(numout,*) '          Index of the equator on the MERCATOR grid:', ijeq
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zti = REAL( ji - 1 + nimpp - 1 )         ;   ztj = REAL( jj - ijeq + njmpp - 1 )
-               zui = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zuj = REAL( jj - ijeq + njmpp - 1 )
-               zvi = REAL( ji - 1 + nimpp - 1 )         ;   zvj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
-               zfi = REAL( ji - 1 + nimpp - 1 ) + 0.5   ;   zfj = REAL( jj - ijeq + njmpp - 1 ) + 0.5
-         ! Longitude
-               glamt(ji,jj) = ppglam0 + ppe1_deg * zti
-               glamu(ji,jj) = ppglam0 + ppe1_deg * zui
-               glamv(ji,jj) = ppglam0 + ppe1_deg * zvi
-               glamf(ji,jj) = ppglam0 + ppe1_deg * zfi
-         ! Latitude
-               gphit(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* ztj ) )
-               gphiu(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zuj ) )
-               gphiv(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zvj ) )
-               gphif(ji,jj) = 1./rad * ASIN ( TANH( ppe1_deg *rad* zfj ) )
-         ! e1
-               e1t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
-               e1u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
-               e1v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
-               e1f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
-         ! e2
-               e2t(ji,jj) = ra * rad * COS( rad * gphit(ji,jj) ) * ppe1_deg
-               e2u(ji,jj) = ra * rad * COS( rad * gphiu(ji,jj) ) * ppe1_deg
-               e2v(ji,jj) = ra * rad * COS( rad * gphiv(ji,jj) ) * ppe1_deg
-               e2f(ji,jj) = ra * rad * COS( rad * gphif(ji,jj) ) * ppe1_deg
-            END DO
-         END DO
-         !
-      CASE ( 5 )                   !==  beta-plane with regular grid-spacing and rotated domain ==! (GYRE configuration)
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          beta-plane with regular grid-spacing and rotated domain (GYRE configuration)'
-         IF(lwp) WRITE(numout,*) '          given by ppe1_m and ppe2_m'
-         !
-         IF(lwp) THEN
-            WRITE(numout,*) '               control print of value of ln_read_cfg     = ', ln_read_cfg
-         ENDIF
-         IF (ln_read_cfg) THEN 
-            IF(lwp) WRITE(numout,*) '          ln_read_cfg read input files'
-            ie1e2u_v = 0                  ! set to unread e1e2u and e1e2v
-            iff = 0                       ! set to unread ff
-            !
-            CALL hgr_read( ie1e2u_v, iff )     ! read the coordinate.nc file !SF to be changed in mesh_mask
-            !
-            IF( ie1e2u_v == 0 ) THEN   ! e1e2u and e1e2v have not been read: compute them
-               !                       ! e2u and e1v does not include a reduction in some strait: apply reduction
-               e1e2u (:,:) = e1u(:,:) * e2u(:,:)
-               e1e2v (:,:) = e1v(:,:) * e2v(:,:)
-            ENDIF
-         ELSE
-            IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module'
-            CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    &
-            &                 glamt  , glamu , glamv , glamf ,    &
-            &                 gphit  , gphiu , gphiv , gphif ,    &
-            &                 e1t    , e1u   , e1v   , e1f   ,    &
-            &                 e2t    , e2u   , e2v   , e2f   ,    &
-            &                 ie1e2u_v  )
-         !
-         ENDIF
-         !
-      CASE DEFAULT
-         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
-         CALL ctl_stop( ctmp1 )
-         !
-      END SELECT
-      
+         IF(lwp) WRITE(numout,*) '          read horizontal mesh in "mesh_mask" file'
+         !
+         CALL hgr_read( ie1e2u_v, iff )     
+         !
+      ELSE    
+         IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module'
+         CALL usr_def_hgr( nbench , jp_cfg, iff   , ff    ,    &
+         &                 glamt  , glamu , glamv , glamf ,    &
+         &                 gphit  , gphiu , gphiv , gphif ,    &
+         &                 e1t    , e1u   , e1v   , e1f   ,    &
+         &                 e2t    , e2u   , e2v   , e2f   ,    &
+         &                 ie1e2u_v  )
+         !
+      ENDIF
+      !
       ! associated horizontal metrics
       ! -----------------------------
@@ -321,5 +128,5 @@
       e1e2t (:,:) = e1t(:,:) * e2t(:,:)   ;   r1_e1e2t(:,:) = 1._wp / e1e2t(:,:)
       e1e2f (:,:) = e1f(:,:) * e2f(:,:)   ;   r1_e1e2f(:,:) = 1._wp / e1e2f(:,:)
-      IF( jphgr_msh /= 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
+      IF( ie1e2u_v == 0 ) THEN               ! e1e2u and e1e2v have not been set: compute them
          e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
          e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
@@ -330,61 +137,12 @@
       e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
       e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
-
-
+      
       ! ================= !
       !  Coriolis factor  !
       ! ================= !
-
-      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
-      !
-      CASE ( 0, 1, 4 )               ! mesh on the sphere
+      IF ( iff == 0 ) THEN       ! Coriolis parameter has not been defined: compute it on the sphere 
          !
          ff(:,:) = 2. * omega * SIN( rad * gphif(:,:) ) 
          !
-      CASE ( 2 )                     ! f-plane at ppgphi0 
-         !
-         ff(:,:) = 2. * omega * SIN( rad * ppgphi0 )
-         !
-         IF(lwp) WRITE(numout,*) '          f-plane: Coriolis parameter = constant = ', ff(1,1)
-         !
-      CASE ( 3 )                     ! beta-plane
-         !
-         zbeta   = 2. * omega * COS( rad * ppgphi0 ) / ra                       ! beta at latitude ppgphi0
-         zphi0   = ppgphi0 - REAL( jpjglo/2) * ppe2_m / ( ra * rad )           ! latitude of the first row F-points
-         !
-#if defined key_agrif
-         IF( cp_cfg == 'eel' .AND. jp_cfg == 6 ) THEN       ! for EEL6 configuration only
-            IF( .NOT.Agrif_Root() ) THEN
-              zphi0 = ppgphi0 - REAL( Agrif_Parent(jpjglo)/2)*Agrif_Parent(ppe2_m) / (ra * rad)
-            ENDIF
-         ENDIF
-#endif         
-         zf0     = 2. * omega * SIN( rad * zphi0 )                              ! compute f0 1st point south
-         !
-         ff(:,:) = ( zf0  + zbeta * gphif(:,:) * 1.e+3 )                        ! f = f0 +beta* y ( y=0 at south)
-         !
-         IF(lwp) THEN
-            WRITE(numout,*) 
-            WRITE(numout,*) '          Beta-plane: Beta parameter = constant = ', ff(nldi,nldj)
-            WRITE(numout,*) '          Coriolis parameter varies from ', ff(nldi,nldj),' to ', ff(nldi,nlej)
-         ENDIF
-         IF( lk_mpp ) THEN 
-            zminff=ff(nldi,nldj)
-            zmaxff=ff(nldi,nlej)
-            CALL mpp_min( zminff )   ! min over the global domain
-            CALL mpp_max( zmaxff )   ! max over the global domain
-            IF(lwp) WRITE(numout,*) '          Coriolis parameter varies globally from ', zminff,' to ', zmaxff
-         END IF
-         !
-      END SELECT
-
-
-      ! Control of domain for symetrical condition
-      ! ------------------------------------------
-      ! The equator line must be the latitude coordinate axe
-
-      IF( nperio == 2 ) THEN
-         znorme = SQRT( SUM( gphiu(:,2) * gphiu(:,2) ) ) / REAL( jpi )
-         IF( znorme > 1.e-13 ) CALL ctl_stop( ' ===>>>> : symmetrical condition: rerun with good equator line' )
       ENDIF
       !
@@ -403,6 +161,6 @@
       USE iom
       !!
-      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0)
-      INTEGER, INTENT( inout ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0)
+      INTEGER, INTENT( out ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in mesh_mask file (=1) or not (=0)
+      INTEGER, INTENT( out ) ::   kff        ! fag: kff read in mesh_mask file (=1) or not (=0)
       !
       INTEGER ::   inum   ! temporary logical unit
@@ -415,5 +173,4 @@
       ENDIF
       !
-      !SF CALL iom_open( 'coordinates', inum )
       CALL iom_open( 'mesh_mask', inum )
       !
