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 6482)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6484)
@@ -297,11 +297,10 @@
          ELSE
             IF(lwp) WRITE(numout,*) '          ln_read_cfg CALL user_defined module'
-            CALL usr_def_hgr( nbench , jp_cfg  ,                    &
-            &                 ff,                         &
-            &                 glamt  , glamu   , glamv   , glamf   ,         &
-            &                 gphit  , gphiu   , gphiv   , gphif   ,         &
-            &                 e1t    , e1u     , e1v     , e1f     ,         &
-            &                 e2t    , e2u     , e2v     , e2f     ,         & 
-            &                 e1e2u  , e1e2v   , e2_e1u  , e1_e2v     )
+            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
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 6482)
+++ branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90	(revision 6484)
@@ -2,11 +2,20 @@
    !!==============================================================================
    !!                       ***  MODULE usrdef   ***
-   !! User defined module: used like example to define init, sbc, bathy, ...
+   !! User defined module:  set analytically the minimum information needed to
+   !!       define a configuration (domain, initial state, forcing)
+   !!
+   !!              ===  GYRE configuration given as an example  ===
+   !! 
+   !!        This module is given as an example, it can be modified
+   !!      following the user's desiderata. 
+   !!      second order accuracy schemes.
    !!==============================================================================
-   !! History :  NEMO ! 2016-03  (S. Flavoni) 
+   !! History :  NEMO ! 2016-03  (S. Flavoni)  Original code 
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   usr_def_hgr       : compute the horizontal mesh 
+   !!   usr_def_zgr       : compute the vertical mesh 
+   !!  to be defined:
    !!   usr_def_ini       : initial state 
    !!   usr_def_sbc       : compute the surface bounday conditions
@@ -21,91 +30,53 @@
 
    PUBLIC usr_def_hgr
+   PUBLIC usr_def_zgr
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id: usrdef.F90 6140 2016-03-21 11:35:23Z flavoni $ 
+   !! NEMO/OPA 3.7 , NEMO Consortium (2016)
+   !! $Id: $ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE usr_def_hgr( nbench, jp_cfg, pff,    &
-           &              pglamt, pglamu, pglamv , pglamf,     &
-           &              pgphit, pgphiu, pgphiv , pgphif,     &          
-           &              pe1t  , pe1u  , pe1v   , pe1f  ,     & 
-           &              pe2t  , pe2u  , pe2v   , pe2f  ,     &
-           &              pe1e2u, pe1e2v, pe2_e1u, pe1_e2v      )
-      !!----------------------------------------------------------------------
+
+   SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff, pff,       & ! Coriolis parameter  (if domain not on the sphere)
+           &               pglamt, pglamu, pglamv , pglamf, & ! gridpoints position (required)
+           &               pgphit, pgphiu, pgphiv , pgphif, & !         
+           &               pe1t  , pe1u  , pe1v   , pe1f  , & ! scale factors (required)
+           &               pe2t  , pe2u  , pe2v   , pe2f  , & !
+           &               ke1e2u_v                       )! u- & v-surfaces (if gridsize reduction is used in strait(s))
+      !!-------------------------------------------------------------------------------------------------
       !!                  ***  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).
+      !! ** Purpose :   user defined mesh and Coriolis parameter
       !!
-      !! ** 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)
+      !! ** Method  :   set all intent(out) argument to a proper value
       !!
-      !!      The coriolis factor is given at z-point by:
-      !!         ff = 2.*omega*sin(gphif)      (in s-1)
+      !!                Here GYRE configuration :
+      !!          Rectangular mid-latitude domain 
+      !!          - with axes rotated by 45 degrees
+      !!          - a constant horizontal resolution of 106 km 
+      !!          - on a beta-plane
       !!
-      !!      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.
+      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees) 
+      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
+      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
+      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
+      !!-------------------------------------------------------------------------------------------------
+
       !!----------------------------------------------------------------------
       INTEGER                 , INTENT(in   ) ::   nbench, jp_cfg   ! parameter of namelist for benchmark, and dimension of GYRE
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff              !  coriolis factor at f-point 
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs 
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  horizontal scale factors 
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  horizontal scale factors
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1e2u , pe1e2v    !  horizontal scale factors
-      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2_e1u, pe1_e2v   !  horizontal scale factors
-       !!----------------------------------------------------------------------
+      INTEGER                 , INTENT(  out) ::   kff              ! =1 Coriolis parameter computed here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pff              ! Coriolis factor at f-point                  [1/s]
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pglamt, pglamu, pglamv, pglamf !  longitude outputs        [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgphit, pgphiu, pgphiv, pgphif !  latitude outputs         [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe1t, pe1u, pe1v, pe1f   !  i-scale factors                      [m]
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pe2t, pe2u, pe2v, pe2f   !  j-scale factors                      [m]
+      INTEGER                 , INTENT(  out) ::   ke1e2u_v                 ! =1 u- & v-surfaces read here, =0 otherwise 
+      !!-------------------------------------------------------------------------------
       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
-      !!----------------------------------------------------------------------
+      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0   ! local scalars
+      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme        ! local scalars
+      !!-------------------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('usr_def_hgr')
@@ -128,10 +99,10 @@
       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 )
+      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
+      zphi0 = 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
+         WRITE(numout,*) '          ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
       ENDIF
       !
@@ -143,87 +114,56 @@
             !glamt(i,j) longitude at T-point
             !gphit(i,j) latitude at T-point  
-            pglamt(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pgphit(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
+            pglamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+            pgphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
             !
             !glamu(i,j) longitude at U-point
             !gphiu(i,j) latitude at U-point
-            pglamu(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pgphiu(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
+            pglamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+            pgphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
             !
             !glamv(i,j) longitude at V-point
             !gphiv(i,j) latitude at V-point
-            pglamv(ji,jj) = glam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pgphiv(ji,jj) = gphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+            pglamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            pgphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+            !
             !glamf(i,j) longitude at F-point
             !gphif(i,j) latitude at F-point 
-            pglamf(ji,jj) = glam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pgphif(ji,jj) = gphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+            pglamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            pgphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
          END DO
       END DO
       !
-      ! Horizontal scale factors (in meters)
-      !                              ======
+      !       !== Horizontal scale factors ==! (in meters)
+      !                     
+      !                             ! constant grid spacing
       pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1
       pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1
       pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1
       pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1
+      !                             ! NO reduction of grid size in some straits 
+      ke1e2u_v = 0                  !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
       !
-      pe1e2u (:,:) = pe1u(:,:) * pe2u(:,:)   
-      pe1e2v (:,:) = pe1v(:,:) * pe2v(:,:) 
       !
-      pe2_e1u(:,:) = pe2u(:,:) / pe1u(:,:)
-      pe1_e2v(:,:) = pe1v(:,:) / pe2v(:,:)
 
-      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, pglamt(ji,1), pglamu(ji,1),   &
-            pglamv(ji,1), pglamf(ji,1),   &
-            pe1t(ji,1), pe1u(ji,1),   &
-            pe1v(ji,1), pe1f(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, pgphit(1,jj), pgphiu(1,jj),   &
-            &                     pgphiv(1,jj), pgphif(1,jj),   &
-            &                     pe2t  (1,jj), pe2u  (1,jj),   &
-            &                     pe2v  (1,jj), pe2f  (1,jj), jj = 1, jpj, 10 )
-      ENDIF
-
-
-      ! ================= !
-      !  Coriolis factor  !
-      ! ================= !
-      ! beta-plane and rotated domain (gyre configuration)
+      !                      !==  Coriolis parameter  ==!
+      kff = 1                                                           !  indicate not to compute ff afterward
       !
-      !SF old ppsphi0: not more necessary?? 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
+      zbeta = 2. * omega * COS( rad * zphi1 ) / ra                      ! beta at latitude zphi1
+      !SF we overwrite zphi0 (south point in latitude) used just above to define pgphif (value of zphi0=15.5190567531966)
+      !SF for computation of coriolis we keep the parameter of XXXX
+      zphi0 = 15._wp                                                     !  latitude of the most southern grid point  
+      zf0   = 2. * omega * SIN( rad * zphi0 )                            !  compute f0 1st point south
       !
-      pff(:,:) = ( zf0 + zbeta * ABS( pgphif(:,:) - zphi0 ) * rad * ra )   ! f = f0 +beta* y ( y=0 at south)
-      iff = 1
+      pff(:,:) = ( 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 ', pff(nldi,nldj),' to ', pff(nldi,nlej)
-      ENDIF
-      !
-      IF( lk_mpp ) THEN 
-         zminff=pff(nldi,nldj)
-         zmaxff=pff(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(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)'
       !
       IF( nn_timing == 1 )  CALL timing_stop('usr_def_hgr')
       !
    END SUBROUTINE usr_def_hgr
+  
+   SUBROUTINE usr_def_zgr()           ! Coriolis parameter  (if domain not on the sphere)
+       ! subroutine for vertical grid
+   END SUBROUTINE usr_def_zgr
 
    !!======================================================================
