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 6410)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6422)
@@ -16,4 +16,7 @@
    !!            3.7  ! 2015-09  (G. Madec, S. Flavoni) add cell surface and their inverse
    !!                                       add optional read of e1e2u & e1e2v
+   !!             -   ! 2016-04  (S. Flavoni) change configuration's interface:
+   !!                            read file or CALL usr_def module to compute 
+   !!                            horizontal grid (example given for GYRE)
    !!----------------------------------------------------------------------
 
@@ -23,4 +26,5 @@
    !!----------------------------------------------------------------------
    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
@@ -29,9 +33,10 @@
    USE lib_mpp        ! MPP library
    USE timing         ! Timing
+   USE usrdef        ! User defined routine
 
    IMPLICIT NONE
    PRIVATE
 
-   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
+   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in namelist !SF modify some routines????
 
    PUBLIC   dom_hgr   ! called by domain.F90
@@ -99,15 +104,15 @@
       !!                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 !
+      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  ::   isrow                    ! index for ORCA1 starting row
+      INTEGER  ::   ie1e2u_v                 ! fag for u- & v-surface read in coordinate file or not
       !!----------------------------------------------------------------------
       !
@@ -134,7 +139,8 @@
          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
+         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
@@ -273,52 +279,34 @@
          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
+         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()
+            !SF CALL usr_def_hgr( nbench , jp_cfg  ,                    &
+            !SF &                 ff,                         &
+            !SF &                 glamt  , glamu   , glamv   , glamf   ,         &
+            !SF &                 gphit  , gphiu   , gphiv   , gphif   ,         &
+            !SF &                 e1t    , e1u     , e1v     , e1f     ,         &
+            !SF &                 e2t    , e2u     , e2v     , e2f     ,         & 
+            !SF &                 e1e2t  , e1e2u   , e1e2v   , e1e2f   ,         & 
+            !SF &                 e2_e1u , e1_e2v     )
+            !
+         !
+         ENDIF
          !
       CASE DEFAULT
@@ -420,4 +408,6 @@
          zf0   = 2. * omega * SIN( rad * zphi0 )                            ! compute f0 1st point south
          !
+         iff = 0 
+         CALL hgr_read( ie1e2u_v, iff )   
          ff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
          !
@@ -453,14 +443,15 @@
 
 
-   SUBROUTINE hgr_read( ke1e2u_v )
+   SUBROUTINE hgr_read( ke1e2u_v, kff )
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE hgr_read  ***
       !!
-      !! ** Purpose :   Read a coordinate file in NetCDF format using IOM
+      !! ** Purpose :   Read a mesh_mask 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, 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 ::   inum   ! temporary logical unit
@@ -469,9 +460,10 @@
       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 )
+      !SF CALL iom_open( 'coordinates', inum )
+      CALL iom_open( 'mesh_mask', inum )
       !
       CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr )
@@ -494,7 +486,9 @@
       CALL iom_get( inum, jpdom_data, 'e2v'  , e2v  , lrowattr=ln_use_jattr )
       CALL iom_get( inum, jpdom_data, 'e2f'  , e2f  , lrowattr=ln_use_jattr )
+      !SF add read coriolis in mesh_mask file
+      CALL iom_get( inum, jpdom_data, 'ff'   , ff   , lrowattr=ln_use_jattr )
       !
       IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in coordinates file'
+         IF(lwp) WRITE(numout,*) 'hgr_read : e1e2u & e1e2v read in mesh_mask 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 )
@@ -502,4 +496,12 @@
       ELSE
          ke1e2u_v = 0
+      ENDIF
+      !
+      IF( iom_varid( inum, 'kff', ldstop = .FALSE. ) > 0 ) THEN
+         IF(lwp) WRITE(numout,*) 'hgr_read : ff read in mesh_mask file'
+         CALL iom_get( inum, jpdom_data, 'ff'  , ff  , lrowattr=ln_use_jattr )
+         kff = 1
+      ELSE
+         kff = 0
       ENDIF
       !
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6410)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6422)
@@ -242,5 +242,5 @@
          &             nn_isplt, nn_jsplt, nn_jctls, nn_jctle,   &
          &             nn_bench, nn_timing, nn_diacfl
-      NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
+      NAMELIST/namcfg/ ln_read_cfg, cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio, ln_use_jattr
       !!----------------------------------------------------------------------
@@ -548,4 +548,5 @@
          WRITE(numout,*) '~~~~~~~ '
          WRITE(numout,*) '   Namelist namcfg'
+         WRITE(numout,*) '      read configuration definition files          ln_read_cfg = ', ln_read_cfg
          WRITE(numout,*) '      configuration name                               cp_cfg  = ', TRIM(cp_cfg)
          WRITE(numout,*) '      configuration zoom name                          cp_cfz  = ', TRIM(cp_cfz)
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 6410)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 6422)
@@ -28,4 +28,5 @@
    !!                   namcfg namelist parameters
    !!----------------------------------------------------------------------
+   LOGICAL       ::   ln_read_cfg      !: logical to read all needed files for a configuration
    CHARACTER(lc) ::   cp_cfg           !: name of the configuration
    CHARACTER(lc) ::   cp_cfz           !: name of the zoom of configuration
Index: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90
===================================================================
--- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90	(revision 6422)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90	(revision 6422)
@@ -0,0 +1,243 @@
+MODULE usrdef
+   !!==============================================================================
+   !!                       ***  MODULE usrdef   ***
+   !! User defined module: used like example to define init, sbc, bathy, ...
+   !!==============================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni) 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr       : initialize the horizontal mesh 
+   !!   usr_def_ini       : initial state 
+   !!   usr_def_sbc       : initialize the surface bounday conditions
+   !!   usr_def_xxx       : initialize the xxx
+   !!----------------------------------------------------------------------
+   !SF  first attempt to define a user_defined_module 
+   !SF  with of NEMO's routines
+   ! 
+   !SF    USE dom_oce        ! ocean space and time domain
+   !SF    USE par_oce        ! ocean space and time domain
+   !SF    USE phycst         ! physical constants
+   !SF    USE mppini         ! shared/distributed memory setting (mpp_init routine)
+   !SF    USE in_out_manager ! I/O manager
+   !SF    USE lib_mpp        ! MPP library
+   !SF    USE timing         ! Timing
+   USE step_oce       ! module used in the ocean time stepping module (step.F90)
+   IMPLICIT NONE
+   PRIVATE
+
+
+   PUBLIC usr_def_hgr
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id: usrdef.F90 6140 2016-03-21 11:35:23Z flavoni $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+SUBROUTINE usr_def_hgr()
+!SF without USE : pass all inputs and outputs 
+!SF   SUBROUTINE usr_def_hgr( nbench, jp_cfg, pff,    &
+!SF            &              pglamt, pglamu, pglamv, pglamf,     &
+!SF            &              pgphit, pgphiu, pgphiv, pgphif,     &          
+!SF            &              pe1t  , pe1u  , pe1v  , pe1f  ,     & 
+!SF            &              pe2t  , pe2u  , pe2v  , pe2f  ,     &
+!SF            &              pe1e2t, pe1e2u, pe1e2v, pe1e2f,     &
+!SF            &              pe2_e1u , pe1_e2v                   )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE usr_def_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, gphif: 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- i and f-points.
+      !!              - define ff: coriolis factor at f-point
+      !!
+      !! References :   Marti, Madec and Delecluse, 1992, JGR
+      !!                Madec, Imbard, 1996, Clim. Dyn.
+      !!----------------------------------------------------------------------
+!SF   all  varibales needed without using "USE"
+!SF       !!----------------------------------------------------------------------
+!SF       INTEGER                 , INTENT(in   ) ::   nbench, jp_cfg   ! parameter of namelist for benchmark, and dimension of GYRE
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff                             !  coriolis factor at f-point 
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs 
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  horizontal scale factors 
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  horizontal scale factors
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1e2t, pe1e2u, pe1e2v, pe1e2f   !  horizontal scale factors
+!SF       REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2_e1u, pe1_e2v   !  horizontal scale factors
+!SF       !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1, iff   ! temporary integers
+      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
+      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
+      REAL(wp) ::   zphi0, zbeta, znorme !
+      REAL(wp) ::   glam0, gphi0
+      REAL(wp) ::   zarg, zf0, zminff, zmaxff
+      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
+      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('usr_def_hgr')
+      !
+      !==  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)'
+      !
+      ! 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
+            !
+            !glamt(i,j) longitude at T-point
+            !gphit(i,j) latitude at T-point  
+            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(i,j) longitude at U-point
+            !gphiu(i,j) latitude at U-point
+            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(i,j) longitude at V-point
+            !gphiv(i,j) latitude at V-point
+            glamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            gphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+            !glamf(i,j) longitude at F-point
+            !gphif(i,j) latitude at F-point 
+            glamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            gphif(ji,jj) = gphi0 - zim1  * 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
+      !
+      e1e2u (:,:) = e1u(:,:) * e2u(:,:)   
+      e1e2v (:,:) = e1v(:,:) * e2v(:,:) 
+      !
+      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  !
+      ! ================= !
+      ! beta-plane and rotated domain (gyre configuration)
+      !
+      !SF old zbeta = 2. * omega * COS( rad * ppgphi0 ) / ra                     ! beta at latitude ppgphi0
+      zbeta = 2. * omega * COS( rad * gphi0 ) / ra                       ! beta at latitude gphi0
+      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)
+      iff = 1
+      !
+      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
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr')
+      !
+   END SUBROUTINE usr_def_hgr
+
+   !!======================================================================
+END MODULE usrdef
