Index: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 7646)
@@ -16,15 +16,18 @@
    !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
    !!                                       add optional read of e1e2u & e1e2v
+   !!             -   ! 2016-04  (S. Flavoni, G. Madec) new configuration interface: read or usrdef.F90
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   dom_hgr       : initialize the horizontal mesh 
-   !!   hgr_read      : read "coordinate" NetCDF file 
+   !!   hgr_read      : read horizontal information in the domain configuration file 
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
+   USE par_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE domwri         ! write 'meshmask.nc' & 'coordinate_e1e2u_v.nc' files
+   USE usrdef_hgr     ! User defined routine
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
@@ -33,10 +36,8 @@
    PRIVATE
 
-   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
-
    PUBLIC   dom_hgr   ! called by domain.F90
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -48,66 +49,32 @@
       !!                  ***  ROUTINE dom_hgr  ***
       !!
-      !! ** Purpose :   Compute the geographical position (in degre) of the 
-      !!      model grid-points,  the horizontal scale factors (in meters) and 
-      !!      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.
-      !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj               ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! 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
+      !! ** Purpose :   Read or compute the geographical position (in degrees)  
+      !!      of the model grid-points, the horizontal scale factors (in meters), 
+      !!      the associated horizontal metrics, and the Coriolis factor (in s-1).
+      !!
+      !! ** 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 ::   ie1e2u_v   ! flag for u- & v-surfaces 
+      INTEGER ::   iff        ! flag for Coriolis parameter
       !!----------------------------------------------------------------------
       !
@@ -117,217 +84,49 @@
          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,*) '~~~~~~~   '
+         WRITE(numout,*) '   namcfg : read (=T) or user defined (=F) configuration    ln_read_cfg  = ', ln_read_cfg
+      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
-         !
-         CALL hgr_read( ie1e2u_v )     ! read the coordinate.nc file
-         !
-         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(:,:) 
+         IF(lwp) WRITE(numout,*) '          read horizontal mesh in ', TRIM( cn_domcfg ), ' file'
+         !
+         CALL hgr_read   ( glamt , glamu , glamv , glamf ,   &    ! geographic position (required)
+            &              gphit , gphiu , gphiv , gphif ,   &    !     -        -
+            &              iff   , ff_f  , ff_t  ,           &    ! Coriolis parameter (if not on the sphere)
+            &              e1t   , e1u   , e1v   , e1f   ,   &    ! scale factors (required)
+            &              e2t   , e2u   , e2v   , e2f   ,   &    !    -     -        -
+            &              ie1e2u_v      , e1e2u , e1e2v     )    ! u- & v-surfaces (if gridsize reduction in some straits)
+         !
+      ELSE                          !==  User defined configuration  ==! 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '          User defined horizontal mesh (usr_def_hgr)'
+         !
+         CALL usr_def_hgr( glamt , glamu , glamv , glamf ,   &    ! geographic position (required)
+            &              gphit , gphiu , gphiv , gphif ,   &    !
+            &              iff   , ff_f  , ff_t  ,           &    ! Coriolis parameter  (if domain not on the sphere)
+            &              e1t   , e1u   , e1v   , e1f   ,   &    ! scale factors       (required)
+            &              e2t   , e2u   , e2v   , e2f   ,   &    !
+            &              ie1e2u_v      , e1e2u , e1e2v     )    ! u- & v-surfaces (if gridsize reduction is used in strait(s))
+         !
+      ENDIF
+      !
+      !                             !==  Coriolis parameter  ==!   (if necessary)
+      !
+      IF( iff == 0 ) THEN                 ! Coriolis parameter has not been defined 
+         IF(lwp) WRITE(numout,*) '          Coriolis parameter calculated on the sphere from gphif & gphit'
+         ff_f(:,:) = 2. * omega * SIN( rad * gphif(:,:) )     ! compute it on the sphere at f-point
+         ff_t(:,:) = 2. * omega * SIN( rad * gphit(:,:) )     !    -        -       -    at t-point
+      ELSE
+         IF( ln_read_cfg ) THEN
+            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been read in ', TRIM( cn_domcfg ), ' file'
+         ELSE
+            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been set in usr_def_hgr routine'
          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'
-         !
-         ! Position coordinates (in kilometers)
-         !                          ==========
-         !
-         ! angle 45deg and ze1=106.e+3 / jp_cfg forced -> zlam1 = -85deg, zphi1 = 29degN
-         zlam1 = -85._wp
-         zphi1 =  29._wp
-         ! resolution in meters
-         ze1 = 106000. / REAL( jp_cfg , wp )            
-         ! benchmark: forced the resolution to be about 100 km
-         IF( nbench /= 0 )   ze1 = 106000._wp     
-         zsin_alpha = - SQRT( 2._wp ) * 0.5_wp
-         zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
-         ze1deg = ze1 / (ra * rad)
-         IF( nbench /= 0 )   ze1deg = ze1deg / REAL( jp_cfg , wp )   ! benchmark: keep the lat/+lon
-         !                                                           ! at the right jp_cfg resolution
-         glam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
-         gphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
-         !
-         IF( nprint==1 .AND. lwp )   THEN
-            WRITE(numout,*) '          ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
-            WRITE(numout,*) '          ze1deg', ze1deg, 'glam0', glam0, 'gphi0', gphi0
-         ENDIF
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5
-               zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5
-               !
-               glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-               gphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
-               !
-               glamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-               gphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
-               !
-               glamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-               gphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
-               !
-               glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-               gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
-            END DO
-         END DO
-         !
-         ! Horizontal scale factors (in meters)
-         !                              ======
-         e1t(:,:) =  ze1     ;      e2t(:,:) = ze1
-         e1u(:,:) =  ze1     ;      e2u(:,:) = ze1
-         e1v(:,:) =  ze1     ;      e2v(:,:) = ze1
-         e1f(:,:) =  ze1     ;      e2f(:,:) = ze1
-         !
-      CASE DEFAULT
-         WRITE(ctmp1,*) '          bad flag value for jphgr_msh = ', jphgr_msh
-         CALL ctl_stop( ctmp1 )
-         !
-      END SELECT
-      
-      ! associated horizontal metrics
-      ! -----------------------------
+      ENDIF
+
+      !
+      !                             !==  associated horizontal metrics  ==!
       !
       r1_e1t(:,:) = 1._wp / e1t(:,:)   ;   r1_e2t (:,:) = 1._wp / e2t(:,:)
@@ -338,113 +137,18 @@
       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
-         e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
+      IF( ie1e2u_v == 0 ) THEN               ! u- & v-surfaces have not been defined
+         IF(lwp) WRITE(numout,*) '          u- & v-surfaces calculated as e1 e2 product'
+         e1e2u (:,:) = e1u(:,:) * e2u(:,:)         ! compute them
          e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
-      ENDIF
-      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in both cases
+      ELSE
+         IF(lwp) WRITE(numout,*) '          u- & v-surfaces have been read in "mesh_mask" file:'
+         IF(lwp) WRITE(numout,*) '                     grid size reduction in strait(s) is used'
+      ENDIF
+      r1_e1e2u(:,:) = 1._wp / e1e2u(:,:)     ! compute their invert in any cases
       r1_e1e2v(:,:) = 1._wp / e1e2v(:,:)
       !   
       e2_e1u(:,:) = e2u(:,:) / e1u(:,:)
       e1_e2v(:,:) = e1v(:,:) / e2v(:,:)
-
-      IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN      ! Control print : Grid informations (if not restart)
-         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), e1f(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, 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), e2f  (1,jj), jj = 1, jpj, 10 )
-      ENDIF
-
-
-      ! ================= !
-      !  Coriolis factor  !
-      ! ================= !
-
-      SELECT CASE( jphgr_msh )   ! type of horizontal mesh
-      !
-      CASE ( 0, 1, 4 )               ! mesh 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
-         !
-      CASE ( 5 )                     ! beta-plane and rotated domain (gyre configuration)
-         !
-         zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
-         zphi0 = 15._wp                                                     ! latitude of the first row F-points
-         zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
-         !
-         ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
-         !
-         IF(lwp) THEN
-            WRITE(numout,*) 
-            WRITE(numout,*) '          Beta-plane and rotated domain : '
-            WRITE(numout,*) '          Coriolis parameter varies in this processor 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
+      !
       !
       IF( nn_timing == 1 )  CALL timing_stop('dom_hgr')
@@ -453,50 +157,70 @@
 
 
-   SUBROUTINE hgr_read( ke1e2u_v )
+   SUBROUTINE hgr_read( plamt , plamu , plamv  , plamf  ,   &    ! gridpoints position (required)
+      &                 pphit , pphiu , pphiv  , pphif  ,   &     
+      &                 kff   , pff_f , pff_t  ,            &    ! Coriolis parameter  (if not on the sphere)
+      &                 pe1t  , pe1u  , pe1v   , pe1f   ,   &    ! scale factors       (required)
+      &                 pe2t  , pe2u  , pe2v   , pe2f   ,   &
+      &                 ke1e2u_v      , pe1e2u , pe1e2v     )    ! u- & v-surfaces (if gridsize reduction in some straits)
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE hgr_read  ***
       !!
-      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
-      !!
-      !!----------------------------------------------------------------------
-      USE iom
-      !!
-      INTEGER, INTENT( inout ) ::   ke1e2u_v   ! fag: e1e2u & e1e2v read in coordinate file (=1) or not (=0)
-      !
-      INTEGER ::   inum   ! temporary logical unit
+      !! ** Purpose :   Read a mesh_mask file in NetCDF format using IOM
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs
+      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter read here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point (if found in file)
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors
+      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces read here, =0 otherwise 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v              ! u- & v-surfaces (if found in file)
+      !
+      INTEGER  ::   inum                  ! logical unit
       !!----------------------------------------------------------------------
       !
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'hgr_read : read the horizontal coordinates'
+         WRITE(numout,*) 'hgr_read : read the horizontal coordinates in mesh_mask'
          WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpk = ', jpk
       ENDIF
       !
-      CALL iom_open( 'coordinates', inum )
-      !
-      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr )
-      !
-      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr )
-      !
-      CALL iom_get( inum, jpdom_data, 'e1t'  , e1t  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e1u'  , e1u  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e1v'  , e1v  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e1f'  , e1f  , lrowattr=ln_use_jattr )
-      !
-      CALL iom_get( inum, jpdom_data, 'e2t'  , e2t  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e2u'  , e2u  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
-      CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
+      CALL iom_open( cn_domcfg, inum )
+      !
+      CALL iom_get( inum, jpdom_data, 'glamt', plamt, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamu', plamu, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamv', plamv, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'glamf', plamf, lrowattr=ln_use_jattr )
+      !
+      CALL iom_get( inum, jpdom_data, 'gphit', pphit, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphiu', pphiu, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphiv', pphiv, lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'gphif', pphif, lrowattr=ln_use_jattr )
+      !
+      CALL iom_get( inum, jpdom_data, 'e1t'  , pe1t  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1u'  , pe1u  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1v'  , pe1v  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e1f'  , pe1f  , lrowattr=ln_use_jattr )
+      !
+      CALL iom_get( inum, jpdom_data, 'e2t'  , pe2t  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2u'  , pe2u  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2v'  , pe2v  , lrowattr=ln_use_jattr )
+      CALL iom_get( inum, jpdom_data, 'e2f'  , pe2f  , lrowattr=ln_use_jattr )
+      !
+      IF(  iom_varid( inum, 'ff_f', ldstop = .FALSE. ) > 0  .AND.  &
+         & iom_varid( inum, 'ff_t', ldstop = .FALSE. ) > 0    ) THEN
+         IF(lwp) WRITE(numout,*) '           Coriolis factor at f- and t-points read in ', TRIM( cn_domcfg ), ' file'
+         CALL iom_get( inum, jpdom_data, 'ff_f'  , pff_f  , lrowattr=ln_use_jattr )
+         CALL iom_get( inum, jpdom_data, 'ff_t'  , pff_t  , lrowattr=ln_use_jattr )
+         kff = 1
+      ELSE
+         kff = 0
+      ENDIF
       !
       IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
-         CALL iom_get( inum, jpdom_data, 'e1e2u'  , e1e2u  , lrowattr=ln_use_jattr )
-         CALL iom_get( inum, jpdom_data, 'e1e2v'  , e1e2v  , lrowattr=ln_use_jattr )
+         IF(lwp) WRITE(numout,*) '           e1e2u & e1e2v read in ', TRIM( cn_domcfg ), ' file'
+         CALL iom_get( inum, jpdom_data, 'e1e2u'  , pe1e2u  , lrowattr=ln_use_jattr )
+         CALL iom_get( inum, jpdom_data, 'e1e2v'  , pe1e2v  , lrowattr=ln_use_jattr )
          ke1e2u_v = 1
       ELSE
@@ -505,6 +229,6 @@
       !
       CALL iom_close( inum )
-      
-    END SUBROUTINE hgr_read
+      !
+   END SUBROUTINE hgr_read
     
    !!======================================================================
