Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 6717)
@@ -79,12 +79,4 @@
    ln_use_jattr = .false.  !  use (T) the file attribute: open_ocean_jstart, if present
    !                       !  in netcdf input files, as the start j-row for reading
-/
-!-----------------------------------------------------------------------
-&namzgr        !   vertical coordinate                                  (default: NO selection)
-!-----------------------------------------------------------------------
-   ln_zco      = .false.   !  z-coordinate - full    steps
-   ln_zps      = .false.   !  z-coordinate - partial steps
-   ln_sco      = .false.   !  s- or hybrid z-s-coordinate
-   ln_isfcav   = .false.   !  ice shelf cavity
 /
 !-----------------------------------------------------------------------
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 6717)
@@ -144,6 +144,4 @@
       
          resto_ice(:,:,:) = 0._wp
-         !      Re-calculate the North and South boundary restoring term
-         !      because those boundaries may change with the prescribed zoom area.
          !
          irelax  = 16                     ! width of buffer zone with respect to close boundary
@@ -156,12 +154,12 @@
          ! REM: if there is no ice in the model and in the data, 
          !      no restoring even with non zero resto_ice
-         DO jj = mj0(jpjzoom - 1 + 1), mj1(jpjzoom -1 + irelax)
-            zreltim = zdmpmin + zfactor * ( mjg(jj) - jpjzoom + 1 )
+         DO jj = mj0(1), mj1( irelax)
+            zreltim = zdmpmin + zfactor * mjg(jj)
             resto_ice(:,jj,:) = 1._wp / ( zreltim * 86400._wp )
          END DO
 
          ! North boundary restoring term
-         DO jj =  mj0(jpjzoom -1 + jpjglo - irelax), mj1(jpjzoom - 1 + jpjglo)
-            zreltim = zdmpmin + zfactor * (jpjglo - ( mjg(jj) - jpjzoom + 1 ))
+         DO jj =  mj0(jpjglo - irelax), mj1(jpjglo)
+            zreltim = zdmpmin + zfactor * (jpjglo - mjg(jj))
             resto_ice(:,jj,:) = 1.e0 / ( zreltim * 86400 )
          END DO
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 6717)
@@ -449,12 +449,10 @@
       soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
       sice_0(:,:) = sice
-      !
-      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea 
-         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
-            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
-            soce_0(:,:) = 4._wp
-            sice_0(:,:) = 2._wp
-         END WHERE
-      ENDIF
+      !                                      ! decrease ocean & ice reference salinities in the Baltic sea 
+      WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+         &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+         soce_0(:,:) = 4._wp
+         sice_0(:,:) = 2._wp
+      END WHERE
       !                                      ! embedded sea ice
       IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
@@ -473,9 +471,8 @@
 !!gm
          IF( .NOT.ln_linssh ) THEN
-
-            do jk = 1,jpkm1                     ! adjust initial vertical scale factors
+            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
                e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
                e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
-            end do
+            END DO
             e3t_a(:,:,:) = e3t_b(:,:,:)
             ! Reconstruction of all vertical scale factors at now and before time steps
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6717)
@@ -316,12 +316,10 @@
       soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
       sice_0(:,:) = sice
-      !
-      IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea 
-         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
-            &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
-            soce_0(:,:) = 4._wp
-            sice_0(:,:) = 2._wp
-         END WHERE
-      ENDIF
+      !                                      ! decrease ocean & ice reference salinities in the Baltic Sea area
+      WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+         &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+         soce_0(:,:) = 4._wp
+         sice_0(:,:) = 2._wp
+      END WHERE
       !
       IF( .NOT. ln_rstart ) THEN
@@ -331,6 +329,6 @@
             snwice_mass_b(:,:) = snwice_mass(:,:)
          ELSE
-            snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
-            snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
+            snwice_mass  (:,:) = 0._wp          ! no mass exchanges
+            snwice_mass_b(:,:) = 0._wp          ! no mass exchanges
          ENDIF
          IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6717)
@@ -5,4 +5,8 @@
    !!=====================================================================
    !! History :  3.0  !  2002-11  (C. Ethe)  F90: Free form and module
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
    !!----------------------------------------------------------------------
    USE in_out_manager ! I/O manager
@@ -175,4 +179,13 @@
    END FUNCTION thd_ice_alloc
    
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE thd_ice_alloc          ! Empty routine
+   END SUBROUTINE thd_ice_alloc
+#endif
+ 
    !!======================================================================
 END MODULE thd_ice
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 6717)
@@ -32,7 +32,7 @@
 ! JC: change to allow for different vertical levels
 !     jpk is already set
-!     keep it jpk possibly different from jpkdta which 
+!     keep it jpk possibly different from jpkglo which 
 !     hold parent grid vertical levels number (set earlier)
-!      jpk     = jpkdta 
+!      jpk     = jpkglo 
       jpim1   = jpi-1 
       jpjm1   = jpj-1 
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 6717)
@@ -100,6 +100,6 @@
 
       DO ib_bdy = 1, nb_bdy
-         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
-
+         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
+            !
             td => tides(ib_bdy)
             nblen => idx_bdy(ib_bdy)%nblen
@@ -134,8 +134,6 @@
             ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
             ! relaxation area      
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
-               ilen0(:)=nblen(:)
-            ELSE
-               ilen0(:)=nblenrim(:)
+            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
+            ELSE                                   ;   ilen0(:) = nblenrim(:)
             ENDIF
 
@@ -156,17 +154,17 @@
             td%v   (:,:,:) = 0._wp
 
-            IF (ln_bdytide_2ddta) THEN
+            IF( ln_bdytide_2ddta ) THEN
                ! It is assumed that each data file contains all complex harmonic amplitudes
-               ! given on the data domain (ie global, jpidta x jpjdta)
-               !
-               CALL wrk_alloc( jpi, jpj, zti, ztr )
+               ! given on the global domain (ie global, jpiglo x jpjglo)
+               !
+               CALL wrk_alloc( jpi,jpj,   zti, ztr )
                !
                ! SSH fields
                clfile = TRIM(filtide)//'_grid_T.nc'
-               CALL iom_open (clfile , inum ) 
+               CALL iom_open( clfile , inum ) 
                igrd = 1                       ! Everything is at T-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -180,9 +178,9 @@
                ! U fields
                clfile = TRIM(filtide)//'_grid_U.nc'
-               CALL iom_open (clfile , inum ) 
+               CALL iom_open( clfile , inum ) 
                igrd = 2                       ! Everything is at U-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -196,9 +194,9 @@
                ! V fields
                clfile = TRIM(filtide)//'_grid_V.nc'
-               CALL iom_open (clfile , inum ) 
+               CALL iom_open( clfile , inum ) 
                igrd = 3                       ! Everything is at V-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -210,5 +208,5 @@
                CALL iom_close( inum )
                !
-               CALL wrk_dealloc( jpi, jpj, ztr, zti ) 
+               CALL wrk_dealloc( jpi,jpj,   ztr, zti ) 
                !
             ELSE            
@@ -219,10 +217,7 @@
                !
                ! Set map structure
-               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)
-               ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)
-               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)
-               ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)
-               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)
-               ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)
+               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)   ;   ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)
+               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)   ;   ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)
+               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)   ;   ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)
 
                ! Open files and read in tidal forcing data
@@ -258,4 +253,5 @@
                !
                DEALLOCATE( dta_read )
+               !
             ENDIF ! ln_bdytide_2ddta=.true.
             !
@@ -275,5 +271,5 @@
             dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
             !
-         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
+         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
          !
       END DO ! loop on ib_bdy
@@ -376,4 +372,5 @@
    END SUBROUTINE bdytide_update
 
+
    SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
       !!----------------------------------------------------------------------
@@ -422,16 +419,14 @@
 
       DO ib_bdy = 1,nb_bdy
-
-         IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
-
+         !
+         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
+            !
             nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
             nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
-
-            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
-               ilen0(:)=nblen(:)
-            ELSE
-               ilen0(:)=nblenrim(:)
+            !
+            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
+            ELSE                                   ;   ilen0(:) = nblenrim(:)
             ENDIF     
-
+            !
             ! We refresh nodal factors every day below
             ! This should be done somewhere else
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 6717)
@@ -392,5 +392,5 @@
         ENDIF               
 
-        IF( iptglo .NE. 0 )THEN
+        IF( iptglo /= 0 )THEN
              
            !read points'coordinates and directions 
@@ -399,9 +399,9 @@
            directemp(:) = 0                  !value of directions of each points
            DO jpt=1,iptglo
-              READ(numdct_in)i1,i2
+              READ(numdct_in) i1, i2
               coordtemp(jpt)%I = i1 
               coordtemp(jpt)%J = i2
            ENDDO
-           READ(numdct_in)directemp(1:iptglo)
+           READ(numdct_in) directemp(1:iptglo)
     
            !debug
@@ -416,11 +416,11 @@
            !Now each proc selects only points that are in its domain:
            !--------------------------------------------------------
-           iptloc = 0                    !initialize number of points selected
-           DO jpt=1,iptglo               !loop on listpoint read in the file
-                    
+           iptloc = 0                    ! initialize number of points selected
+           DO jpt = 1, iptglo            ! loop on listpoint read in the file
+              !      
               iiglo=coordtemp(jpt)%I          ! global coordinates of the point
               ijglo=coordtemp(jpt)%J          !  " 
 
-              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
+              IF( iiglo==jpiglo .AND. nimpp==1 )   iiglo = 2         !!gm BUG: Hard coded periodicity !
 
               iiloc=iiglo-nimpp+1   ! local coordinates of the point
@@ -428,12 +428,12 @@
 
               !verify if the point is on the local domain:(1,nlei)*(1,nlej)
-              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
-                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
+              IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
+                  ijloc >= 1 .AND. ijloc <= nlej       )THEN
                  iptloc = iptloc + 1                                                 ! count local points
                  secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
                  secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
               ENDIF
-
-           ENDDO
+              !
+           END DO
      
            secs(jsec)%nb_point=iptloc !store number of section's points
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90	(revision 6716)
+++ 	(revision )
@@ -1,472 +1,0 @@
-MODULE diafwb
-   !!======================================================================
-   !!                       ***  MODULE  diafwb  ***
-   !! Ocean diagnostics: freshwater budget
-   !!======================================================================
-   !! History :  8.2  !  01-02  (E. Durand)  Original code
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!---------------------------------------------------------------------- 
-   !!----------------------------------------------------------------------
-   !!   Only for ORCA2 ORCA1 and ORCA025
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   dia_fwb     : freshwater budget for global ocean configurations
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE sbc_oce         ! ???
-   USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distributed memory computing library
-   USE timing          ! preformance summary
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC dia_fwb    ! routine called by step.F90
-
-   REAL(wp)               ::   a_fwf ,          &
-      &                        a_sshb, a_sshn, a_salb, a_saln
-   REAL(wp), DIMENSION(4) ::   a_flxi, a_flxo, a_temi, a_temo, a_sali, a_salo
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dia_fwb( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dia_fwb  ***
-      !!     
-      !! ** Purpose :
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !!
-      INTEGER  :: inum             ! temporary logical unit
-      INTEGER  :: ji, jj, jk, jt   ! dummy loop indices
-      INTEGER  :: ii0, ii1, ij0, ij1
-      INTEGER  :: isrow         ! index for ORCA1 starting row
-      REAL(wp) :: zarea, zvol, zwei
-      REAL(wp) :: ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4)
-      REAL(wp) :: zt, zs, zu  
-      REAL(wp) :: zsm0, zfwfnew
-      IF( cp_cfg == "orca" .AND. jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN
-      !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )   CALL timing_start('dia_fwb')
-
-      ! Mean global salinity
-      zsm0 = 34.72654
-
-      ! To compute fwf mean value mean fwf
-
-      IF( kt == nit000 ) THEN
-
-         a_fwf    = 0.e0
-         a_sshb   = 0.e0 ! valeur de ssh au debut de la simulation
-         a_salb   = 0.e0 ! valeur de sal au debut de la simulation
-         ! sshb used because diafwb called after tranxt (i.e. after the swap)
-         a_sshb = SUM( e1e2t(:,:) * sshb(:,:) * tmask_i(:,:) )
-         IF( lk_mpp )   CALL mpp_sum( a_sshb )      ! sum over the global domain
-
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwei  = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-                  a_salb = a_salb + ( tsb(ji,jj,jk,jp_sal) - zsm0 ) * zwei
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( a_salb )      ! sum over the global domain
-      ENDIF
-      
-      a_fwf    = SUM( e1e2t(:,:) * ( emp(:,:)-rnf(:,:) ) * tmask_i(:,:) ) 
-      IF( lk_mpp )   CALL mpp_sum( a_fwf    )       ! sum over the global domain
-
-      IF( kt == nitend ) THEN
-         a_sshn = 0.e0
-         a_saln = 0.e0
-         zarea = 0.e0
-         zvol  = 0.e0
-         zfwfnew = 0.e0
-         ! Mean sea level at nitend
-         a_sshn = SUM( e1e2t(:,:) * sshn(:,:) * tmask_i(:,:) )
-         IF( lk_mpp )   CALL mpp_sum( a_sshn )      ! sum over the global domain
-         zarea  = SUM( e1e2t(:,:) *             tmask_i(:,:) )
-         IF( lk_mpp )   CALL mpp_sum( zarea  )      ! sum over the global domain
-         
-         DO jk = 1, jpkm1   
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zwei  = e1e2t(ji,jj) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
-                  a_saln = a_saln + ( tsn(ji,jj,jk,jp_sal) - zsm0 ) * zwei
-                  zvol  = zvol  + zwei
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( a_saln )      ! sum over the global domain
-         IF( lk_mpp )   CALL mpp_sum( zvol )      ! sum over the global domain
-         
-         ! Conversion in m3
-         a_fwf    = a_fwf * rdt * 1.e-3 
-         
-         ! fwf correction to bring back the mean ssh to zero
-         zfwfnew = a_sshn / ( ( nitend - nit000 + 1 ) * rdt ) * 1.e3 / zarea
-
-      ENDIF
-
-
-      ! Calcul des termes de transport 
-      ! ------------------------------
-      
-      ! 1 --> Gibraltar 
-      ! 2 --> Cadiz 
-      ! 3 --> Red Sea
-      ! 4 --> Baltic Sea
-
-      IF( kt == nit000 ) THEN
-         a_flxi(:) = 0.e0
-         a_flxo(:) = 0.e0
-         a_temi(:) = 0.e0
-         a_temo(:) = 0.e0
-         a_sali(:) = 0.e0
-         a_salo(:) = 0.e0
-      ENDIF
-
-      zflxi(:) = 0.e0
-      zflxo(:) = 0.e0
-      ztemi(:) = 0.e0
-      ztemo(:) = 0.e0
-      zsali(:) = 0.e0
-      zsalo(:) = 0.e0
-
-      ! Mean flow at Gibraltar
-
-      IF( cp_cfg == "orca" ) THEN  
-               
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration
-            !                                        ! =======================
-            ii0 = 70   ;   ii1 = 70
-            ij0 = 52   ;   ij1 = 52
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 140   ;   ii1 = 140
-            ij0 = 102   ;   ij1 = 102
-            !                                        ! =======================
-         CASE ( 1 )                                  !  ORCA_R1 configurations
-            !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            !
-            ii0 = 283           ;   ii1 = 283
-            ij0 = 241 - isrow   ;   ij1 = 241 - isrow
-            !                                        ! =======================
-         CASE DEFAULT                                !    ORCA R05 or R025
-            !                                        ! =======================
-            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
-            ! 
-         END SELECT
-         ! 
-         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
-            DO jj = mj0(ij0), mj1(ij1)
-               DO jk = 1, jpk 
-                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
-
-                  IF( un(ji,jj,jk) > 0.e0 ) THEN 
-                     zflxi(1) = zflxi(1) +    zu
-                     ztemi(1) = ztemi(1) + zt*zu
-                     zsali(1) = zsali(1) + zs*zu
-                  ELSE
-                     zflxo(1) = zflxo(1) +    zu
-                     ztemo(1) = ztemo(1) + zt*zu
-                     zsalo(1) = zsalo(1) + zs*zu
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-      
-      ! Mean flow at Cadiz
-      IF( cp_cfg == "orca" ) THEN
-               
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration
-            !                                        ! =======================
-            ii0 = 69   ;   ii1 = 69
-            ij0 = 52   ;   ij1 = 52
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 137   ;   ii1 = 137
-            ij0 = 101   ;   ij1 = 102
-            !                                        ! =======================
-         CASE ( 1 )                                  !  ORCA_R1 configurations
-            !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            ii0 = 282           ;   ii1 = 282
-            ij0 = 240 - isrow   ;   ij1 = 240 - isrow
-            !                                        ! =======================
-         CASE DEFAULT                                !    ORCA R05 or R025
-            !                                        ! =======================
-            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
-            ! 
-         END SELECT
-         ! 
-         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
-            DO jj = mj0(ij0), mj1(ij1)
-               DO jk = 1, jpk 
-                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
-                  
-                  IF( un(ji,jj,jk) > 0.e0 ) THEN 
-                     zflxi(2) = zflxi(2) +    zu
-                     ztemi(2) = ztemi(2) + zt*zu
-                     zsali(2) = zsali(2) + zs*zu
-                  ELSE
-                     zflxo(2) = zflxo(2) +    zu
-                     ztemo(2) = ztemo(2) + zt*zu
-                     zsalo(2) = zsalo(2) + zs*zu
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! Mean flow at Red Sea entrance
-      IF( cp_cfg == "orca" ) THEN
-               
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration
-            !                                        ! =======================
-            ii0 = 83   ;   ii1 = 83
-            ij0 = 45   ;   ij1 = 45
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 160   ;   ii1 = 160
-            ij0 = 88    ;   ij1 = 88 
-            !                                        ! =======================
-         CASE ( 1 )                                  !  ORCA_R1 configurations
-            !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            ii0 = 331           ;   ii1 = 331
-            ij0 = 215 - isrow   ;   ij1 = 215 - isrow
-            !                                        ! =======================
-         CASE DEFAULT                                !    ORCA R05 or R025
-            !                                        ! =======================
-            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
-            ! 
-         END SELECT
-         ! 
-         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
-            DO jj = mj0(ij0), mj1(ij1)
-               DO jk = 1, jpk 
-                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
-                  
-                  IF( un(ji,jj,jk) > 0.e0 ) THEN 
-                     zflxi(3) = zflxi(3) +    zu
-                     ztemi(3) = ztemi(3) + zt*zu
-                     zsali(3) = zsali(3) + zs*zu
-                  ELSE
-                     zflxo(3) = zflxo(3) +    zu
-                     ztemo(3) = ztemo(3) + zt*zu
-                     zsalo(3) = zsalo(3) + zs*zu
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! Mean flow at Baltic Sea entrance
-      IF( cp_cfg == "orca" ) THEN
-               
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration
-            !                                        ! =======================
-            ii0 = 1     ;   ii1 = 1  
-            ij0 = 1     ;   ij1 = 1  
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            ii0 = 146   ;   ii1 = 146 
-            ij0 = 116   ;   ij1 = 116
-            !                                        ! =======================
-         CASE ( 1 )                                  !  ORCA_R1 configurations
-            !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            ii0 = 297           ;   ii1 = 297
-            ij0 = 269 - isrow   ;   ij1 = 269 - isrow
-            !                                        ! =======================
-         CASE DEFAULT                                !    ORCA R05 or R025
-            !                                        ! =======================
-            CALL ctl_stop( ' dia_fwb Not yet implemented in ORCA_R05 or R025' )
-            ! 
-         END SELECT
-         ! 
-         DO ji = mi0(ii0), MIN(mi1(ii1),jpim1)
-            DO jj = mj0(ij0), mj1(ij1)
-               DO jk = 1, jpk
-                  zt = 0.5 * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
-                  zs = 0.5 * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
-                  zu = un(ji,jj,jk) * e3t_n(ji,jj,jk) * e2u(ji,jj) * tmask_i(ji,jj)
-                  
-                  IF( un(ji,jj,jk) > 0.e0 ) THEN 
-                     zflxi(4) = zflxi(4) +    zu
-                     ztemi(4) = ztemi(4) + zt*zu
-                     zsali(4) = zsali(4) + zs*zu
-                  ELSE
-                     zflxo(4) = zflxo(4) +    zu
-                     ztemo(4) = ztemo(4) + zt*zu
-                     zsalo(4) = zsalo(4) + zs*zu
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! Sum at each time-step
-      DO jt = 1, 4 
-         !
-         IF( zflxi(jt) /= 0.e0 ) THEN
-            a_flxi(jt) = a_flxi(jt) + zflxi(jt)
-            a_temi(jt) = a_temi(jt) + ztemi(jt)/zflxi(jt)
-            a_sali(jt) = a_sali(jt) + zsali(jt)/zflxi(jt)
-         ENDIF
-         !
-         IF( zflxo(jt) /= 0.e0 ) THEN
-            a_flxo(jt) = a_flxo(jt) + zflxo(jt)
-            a_temo(jt) = a_temo(jt) + ztemo(jt)/zflxo(jt)
-            a_salo(jt) = a_salo(jt) + zsalo(jt)/zflxo(jt)
-         ENDIF
-         !
-      END DO
-
-      IF( kt == nitend ) THEN
-         DO jt = 1, 4 
-            a_flxi(jt) = a_flxi(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
-            a_temi(jt) = a_temi(jt) /   FLOAT( nitend - nit000 + 1 )
-            a_sali(jt) = a_sali(jt) /   FLOAT( nitend - nit000 + 1 )
-            a_flxo(jt) = a_flxo(jt) / ( FLOAT( nitend - nit000 + 1 ) * 1.e6 )
-            a_temo(jt) = a_temo(jt) /   FLOAT( nitend - nit000 + 1 )
-            a_salo(jt) = a_salo(jt) /   FLOAT( nitend - nit000 + 1 )
-         END DO
-         IF( lk_mpp ) THEN
-            CALL mpp_sum( a_flxi, 4 )      ! sum over the global domain
-            CALL mpp_sum( a_temi, 4 )      ! sum over the global domain
-            CALL mpp_sum( a_sali, 4 )      ! sum over the global domain
-
-            CALL mpp_sum( a_flxo, 4 )      ! sum over the global domain
-            CALL mpp_sum( a_temo, 4 )      ! sum over the global domain
-            CALL mpp_sum( a_salo, 4 )      ! sum over the global domain
-         ENDIF
-      ENDIF
-
-
-      ! Ecriture des diagnostiques 
-      ! --------------------------
-
-      IF ( kt == nitend .AND. cp_cfg == "orca" .AND. lwp ) THEN
-
-         CALL ctl_opn( inum, 'STRAIT.dat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
-         WRITE(inum,*)
-         WRITE(inum,*)    'Net freshwater budget '
-         WRITE(inum,9010) '  fwf    = ',a_fwf,   ' m3 =', a_fwf   /(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
-         WRITE(inum,*)
-         WRITE(inum,9010) '  zarea =',zarea
-         WRITE(inum,9010) '  zvol  =',zvol
-         WRITE(inum,*)
-         WRITE(inum,*)    'Mean sea level : '
-         WRITE(inum,9010) '  at nit000 = ',a_sshb        ,' m3 '
-         WRITE(inum,9010) '  at nitend = ',a_sshn        ,' m3 '
-         WRITE(inum,9010) '  diff      = ',(a_sshn-a_sshb),' m3 =', (a_sshn-a_sshb)/(FLOAT(nitend-nit000+1)*rdt) * 1.e-6,' Sv'
-         WRITE(inum,9020) '  mean sea level elevation    =', a_sshn/zarea,' m'
-         WRITE(inum,*)
-         WRITE(inum,*)    'Anomaly of salinity content : '
-         WRITE(inum,9010) '  at nit000 = ',a_salb        ,' psu.m3 '
-         WRITE(inum,9010) '  at nitend = ',a_saln        ,' psu.m3 '
-         WRITE(inum,9010) '  diff      = ',(a_saln-a_salb),' psu.m3'
-         WRITE(inum,*)
-         WRITE(inum,*)    'Mean salinity : '
-         WRITE(inum,9020) '  at nit000 =',a_salb/zvol+zsm0   ,' psu '
-         WRITE(inum,9020) '  at nitend =',a_saln/zvol+zsm0   ,' psu '
-         WRITE(inum,9020) '  diff      =',(a_saln-a_salb)/zvol,' psu'
-         WRITE(inum,9020) '  S-SLevitus=',a_saln/zvol,' psu'
-         WRITE(inum,*)
-         WRITE(inum,*)    'Gibraltar : '
-         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(1)
-         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(1)
-         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(1)
-         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(1)
-         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(1)
-         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(1)
-         WRITE(inum,*)
-         WRITE(inum,*)    'Cadiz : '
-         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(2)
-         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(2)
-         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(2)
-         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(2)
-         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(2)
-         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(2)
-         WRITE(inum,*)
-         WRITE(inum,*)    'Bab el Mandeb : '
-         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(3)
-         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(3)
-         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(3)
-         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(3)
-         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(3)
-         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(3)
-         WRITE(inum,*)
-         WRITE(inum,*)    'Baltic : '
-         WRITE(inum,9030) '  Flux entrant (Sv) :', a_flxi(4)
-         WRITE(inum,9030) '  Flux sortant (Sv) :', a_flxo(4)
-         WRITE(inum,9030) '  T entrant (deg)   :', a_temi(4)
-         WRITE(inum,9030) '  T sortant (deg)   :', a_temo(4)
-         WRITE(inum,9030) '  S entrant (psu)   :', a_sali(4)
-         WRITE(inum,9030) '  S sortant (psu)   :', a_salo(4)
-         CLOSE(inum)
-      ENDIF
-
-      IF( nn_timing == 1 )   CALL timing_stop('dia_fwb')
-
- 9005 FORMAT(1X,A,ES24.16)
- 9010 FORMAT(1X,A,ES12.5,A,F10.5,A)
- 9020 FORMAT(1X,A,F10.5,A)
- 9030 FORMAT(1X,A,F9.4,A)
-  
-      ENDIF 
-
-   END SUBROUTINE dia_fwb
-
-   !!======================================================================
-END MODULE diafwb
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 6717)
@@ -666,5 +666,5 @@
          CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
+         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
 #endif
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90	(revision 6716)
+++ 	(revision )
@@ -1,460 +1,0 @@
-MODULE closea
-   !!======================================================================
-   !!                       ***  MODULE  closea  ***
-   !! Closed Seas  : specific treatments associated with closed seas
-   !!======================================================================
-   !! History :   8.2  !  00-05  (O. Marti)  Original code
-   !!             8.5  !  02-06  (E. Durand, G. Madec)  F90
-   !!             9.0  !  06-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
-   !!        NEMO 3.4  !  03-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dom_clo    : modification of the ocean domain for closed seas cases
-   !!   sbc_clo    : Special handling of closed seas
-   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
-   !!   clo_ups    : set mixed centered/upstream scheme in closed sea (see traadv_cen2)
-   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
-   !!----------------------------------------------------------------------
-   USE oce             ! dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE sbc_oce         ! ocean surface boundary conditions
-   USE lib_fortran,    ONLY: glob_sum, DDPDD
-   USE lbclnk          ! lateral boundary condition - MPP exchanges
-   USE lib_mpp         ! MPP library
-   USE timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC dom_clo      ! routine called by domain module
-   PUBLIC sbc_clo      ! routine called by step module
-   PUBLIC clo_rnf      ! routine called by sbcrnf module
-   PUBLIC clo_ups      ! routine called in traadv_cen2(_jki) module
-   PUBLIC clo_bat      ! routine called in domzgr module
-
-   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
-   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
-   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
-
-   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dom_clo
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_clo  ***
-      !!        
-      !! ** Purpose :   Closed sea domain initialization
-      !!
-      !! ** Method  :   if a closed sea is located only in a model grid point
-      !!                just the thermodynamic processes are applied.
-      !!
-      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
-      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
-      !!                ncsir(), ncsjr() : Location of runoff
-      !!                ncsnr            : number of point where run-off pours
-      !!                ncstt            : Type of closed sea
-      !!                                   =0 spread over the world ocean
-      !!                                   =2 put at location runoff
-      !!----------------------------------------------------------------------
-      INTEGER ::   jc      ! dummy loop indices
-      INTEGER ::   isrow   ! local index
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
-      IF(lwp) WRITE(numout,*)'~~~~~~~'
-      !
-      ! initial values
-      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1
-      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1
-      !
-      ! set the closed seas (in data domain indices)
-      ! -------------------
-      !
-      IF( cp_cfg == "orca" ) THEN
-         !
-         SELECT CASE ( jp_cfg )
-         !                                           ! =======================
-         CASE ( 1 )                                  ! ORCA_R1 configuration
-            !                                        ! =======================
-            ! This dirty section will be suppressed by simplification process:
-            ! all this will come back in input files
-            ! Currently these hard-wired indices relate to configuration with
-            ! extend grid (jpjglo=332)
-            isrow = 332 - jpjglo
-            !
-            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea
-            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow
-            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow
-            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
-            !                                        
-            !                                        ! =======================
-         CASE ( 2 )                                  !  ORCA_R2 configuration
-            !                                        ! =======================
-            !                                            ! Caspian Sea
-            ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
-            ncsi1(1)   =  11  ;  ncsj1(1)   = 103
-            ncsi2(1)   =  17  ;  ncsj2(1)   = 112
-            ncsir(1,1) =   1  ;  ncsjr(1,1) =   1 
-            !                                            ! Great North American Lakes
-            ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth
-            ncsi1(2)   =  97  ;  ncsj1(2)   = 107
-            ncsi2(2)   = 103  ;  ncsj2(2)   = 111
-            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111           
-            !                                            ! Black Sea (crossed by the cyclic boundary condition)
-            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea)
-            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         !
-            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106 
-            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105 
-            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105 
-            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea      
-            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.)
-            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea 
-            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.)
-             
-          
-
-            !                                        ! =======================
-         CASE ( 4 )                                  !  ORCA_R4 configuration
-            !                                        ! =======================
-            !                                            ! Caspian Sea
-            ncsnr(1)   =  1  ;  ncstt(1)   =  0  
-            ncsi1(1)   =  4  ;  ncsj1(1)   = 53 
-            ncsi2(1)   =  4  ;  ncsj2(1)   = 56
-            ncsir(1,1) =  1  ;  ncsjr(1,1) =  1
-            !                                            ! Great North American Lakes
-            ncsnr(2)   =  1  ;  ncstt(2)   =  2 
-            ncsi1(2)   = 49  ;  ncsj1(2)   = 55
-            ncsi2(2)   = 51  ;  ncsj2(2)   = 56
-            ncsir(2,1) = 57  ;  ncsjr(2,1) = 55
-            !                                            ! Black Sea
-            ncsnr(3)   =  4  ;  ncstt(3)   =  2  
-            ncsi1(3)   = 88  ;  ncsj1(3)   = 55 
-            ncsi2(3)   = 91  ;  ncsj2(3)   = 56
-            ncsir(3,1) = 86  ;  ncsjr(3,1) = 53
-            ncsir(3,2) = 87  ;  ncsjr(3,2) = 53 
-            ncsir(3,3) = 86  ;  ncsjr(3,3) = 52 
-            ncsir(3,4) = 87  ;  ncsjr(3,4) = 52
-            !                                            ! Baltic Sea
-            ncsnr(4)   =  1  ;  ncstt(4)   =  2
-            ncsi1(4)   = 75  ;  ncsj1(4)   = 59
-            ncsi2(4)   = 76  ;  ncsj2(4)   = 61
-            ncsir(4,1) = 84  ;  ncsjr(4,1) = 59 
-            !                                        ! =======================
-         CASE ( 025 )                                ! ORCA_R025 configuration
-            !                                        ! =======================
-            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
-            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
-            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
-            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
-            !                                        
-            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea 
-            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
-            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
-            ncsir(2,1) = 1    ; ncsjr(2,1) = 1
-            !
-         END SELECT
-         !
-      ENDIF
-
-      ! convert the position in local domain indices
-      ! --------------------------------------------
-      DO jc = 1, jpncs
-         ncsi1(jc)   = mi0( ncsi1(jc) )
-         ncsj1(jc)   = mj0( ncsj1(jc) )
-
-         ncsi2(jc)   = mi1( ncsi2(jc) )   
-         ncsj2(jc)   = mj1( ncsj2(jc) )  
-      END DO
-      !
-   END SUBROUTINE dom_clo
-
-
-   SUBROUTINE sbc_clo( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_clo  ***
-      !!                    
-      !! ** Purpose :   Special handling of closed seas
-      !!
-      !! ** Method  :   Water flux is forced to zero over closed sea
-      !!      Excess is shared between remaining ocean, or
-      !!      put as run-off in open ocean.
-      !!
-      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean model time step
-      !
-      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
-      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
-      REAL(wp)            ::   zze2, ztmp, zcorr     ! 
-      REAL(wp)            ::   zcoef, zcoef1         ! 
-      COMPLEX(wp)         ::   ctmp 
-      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
-      !                                                   !------------------!
-      IF( kt == nit000 ) THEN                             !  Initialisation  !
-         !                                                !------------------!
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
-         IF(lwp) WRITE(numout,*)'~~~~~~~'
-
-         surf(:) = 0.e0_wp
-         !
-         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
-         !
-         !                                        ! surface of closed seas 
-         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
-            DO jc = 1, jpncs
-               ctmp = CMPLX( 0.e0, 0.e0, wp )
-               DO jj = ncsj1(jc), ncsj2(jc)
-                  DO ji = ncsi1(jc), ncsi2(jc)
-                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
-                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-                  END DO 
-               END DO 
-               IF( lk_mpp )   CALL mpp_sum( ctmp )
-               surf(jc) = REAL(ctmp,wp)
-            END DO
-         ELSE                                          ! Standard calculation           
-            DO jc = 1, jpncs
-               DO jj = ncsj1(jc), ncsj2(jc)
-                  DO ji = ncsi1(jc), ncsi2(jc)
-                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
-                  END DO 
-               END DO 
-            END DO 
-            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain
-         ENDIF
-
-         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
-         DO jc = 1, jpncs
-            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
-         END DO
-
-         ! jpncs+1 : surface of sea, closed seas excluded
-         DO jc = 1, jpncs
-            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
-         END DO           
-         !
-      ENDIF
-      !                                                   !--------------------!
-      !                                                   !  update emp        !
-      zfwf = 0.e0_wp                                      !--------------------!
-      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
-         DO jc = 1, jpncs
-            ctmp = CMPLX( 0.e0, 0.e0, wp )
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
-                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
-               END DO  
-            END DO 
-            IF( lk_mpp )   CALL mpp_sum( ctmp )
-            zfwf(jc) = REAL(ctmp,wp)
-         END DO
-      ELSE                                          ! Standard calculation           
-         DO jc = 1, jpncs
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
-               END DO  
-            END DO 
-         END DO
-         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain
-      ENDIF
-
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
-         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
-         zfwf(3) = zze2
-         zfwf(4) = zze2
-      ENDIF
-
-      zcorr = 0._wp
-
-      DO jc = 1, jpncs
-         !
-         ! The following if avoids the redistribution of the round off
-         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
-            !
-            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
-               zcoef    = zfwf(jc) / surf(jpncs+1)
-               zcoef1   = rcp * zcoef
-               emp(:,:) = emp(:,:) + zcoef
-               qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-               ! accumulate closed seas correction
-               zcorr    = zcorr    + zcoef
-               !
-            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
-               IF ( zfwf(jc) <= 0.e0_wp ) THEN 
-                   DO jn = 1, ncsnr(jc)
-                     ji = mi0(ncsir(jc,jn))
-                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                     IF (      ji > 1 .AND. ji < jpi   &
-                         .AND. jj > 1 .AND. jj < jpj ) THEN 
-                         zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
-                         zcoef1     = rcp * zcoef
-                         emp(ji,jj) = emp(ji,jj) + zcoef
-                         qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                     ENDIF 
-                   END DO 
-               ELSE 
-                   zcoef    = zfwf(jc) / surf(jpncs+1)
-                   zcoef1   = rcp * zcoef
-                   emp(:,:) = emp(:,:) + zcoef
-                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
-                   ! accumulate closed seas correction
-                   zcorr    = zcorr    + zcoef
-               ENDIF
-            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
-               DO jn = 1, ncsnr(jc)
-                  ji = mi0(ncsir(jc,jn))
-                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
-                  IF(      ji > 1 .AND. ji < jpi    &
-                     .AND. jj > 1 .AND. jj < jpj ) THEN 
-                     zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
-                     zcoef1     = rcp * zcoef
-                     emp(ji,jj) = emp(ji,jj) + zcoef
-                     qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
-                  ENDIF 
-               END DO 
-            ENDIF 
-            !
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  zcoef      = zfwf(jc) / surf(jc)
-                  zcoef1     = rcp * zcoef
-                  emp(ji,jj) = emp(ji,jj) - zcoef
-                  qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
-               END DO  
-            END DO 
-            !
-         END IF
-      END DO 
-
-      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
-         DO jc = 1, jpncs                  ! only if it is large enough
-            DO jj = ncsj1(jc), ncsj2(jc)
-               DO ji = ncsi1(jc), ncsi2(jc)
-                  emp(ji,jj) = emp(ji,jj) - zcorr
-                  qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)
-               END DO  
-             END DO 
-          END DO
-      ENDIF
-      !
-      emp (:,:) = emp (:,:) * tmask(:,:,1)
-      !
-      CALL lbc_lnk( emp , 'T', 1._wp )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
-      !
-   END SUBROUTINE sbc_clo
-
-
-   SUBROUTINE clo_rnf( p_rnfmsk )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_rnf  ***
-      !!                    
-      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
-      !!                to be the same as river mouth grid-points
-      !!
-      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
-      !!                at the closed sea outflow grid-point.
-      !!
-      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
-      !
-      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
-             DO jn = 1, 4
-                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )
-                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )
-                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )
-                   END DO
-                END DO
-            END DO 
-         ENDIF 
-      END DO 
-      !
-   END SUBROUTINE clo_rnf
-
-   
-   SUBROUTINE clo_ups( p_upsmsk )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_rnf  ***
-      !!                    
-      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
-      !!                to be the same as river mouth grid-points
-      !!
-      !! ** Method  :   set to 0.5 the upstream mask (upsmsk, see traadv_cen2 
-      !!                module) over the closed seas.
-      !!
-      !! ** Action  :   update (p_)upsmsk (set 0.5 over closed seas)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_upsmsk   ! upstream mask (upsmsk array)
-      !
-      INTEGER  ::   jc, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         DO jj = ncsj1(jc), ncsj2(jc)
-            DO ji = ncsi1(jc), ncsi2(jc)
-               p_upsmsk(ji,jj) = 0.5_wp         ! mixed upstream/centered scheme over closed seas
-            END DO 
-         END DO 
-       END DO 
-       !
-   END SUBROUTINE clo_ups
-   
-      
-   SUBROUTINE clo_bat( pbat, kbat )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE clo_bat  ***
-      !!                    
-      !! ** Purpose :   suppress closed sea from the domain
-      !!
-      !! ** Method  :   set to 0 the meter and level bathymetry (given in 
-      !!                arguments) over the closed seas.
-      !!
-      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pbat   ! bathymetry in meters (bathy array)
-      INTEGER , DIMENSION(jpi,jpj), INTENT(inout) ::   kbat   ! bathymetry in levels (mbathy array)
-      !
-      INTEGER  ::   jc, ji, jj      ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      DO jc = 1, jpncs
-         DO jj = ncsj1(jc), ncsj2(jc)
-            DO ji = ncsi1(jc), ncsi2(jc)
-               pbat(ji,jj) = 0._wp   
-               kbat(ji,jj) = 0   
-            END DO 
-         END DO 
-       END DO 
-       !
-   END SUBROUTINE clo_bat
-
-   !!======================================================================
-END MODULE closea
-
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 6717)
@@ -166,9 +166,7 @@
 
    !!----------------------------------------------------------------------
-   !! masks, bathymetry
+   !! masks, top and bottom ocean point position
    !! ---------------------------------------------------------------------
-!!gm   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt, mbku, mbkv   !: bottom last wet T-, U- and V-level
-!!gm   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book)
@@ -178,5 +176,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   risfdep                 !: Iceshelf draft                    (ISF)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask, ssfmask    !: surface mask at T-,U-, V- and F-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssmask, ssumask, ssvmask             !: surface mask at T-,U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts
@@ -258,5 +256,5 @@
          &      njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) ,     &
          &                                      nleit(jpnij) , nlejt(jpnij) ,     &
-         &      mi0(jpidta)   , mi1 (jpidta),  mj0(jpjdta)   , mj1 (jpjdta) ,     &
+         &      mi0(jpiglo)   , mi1 (jpiglo),  mj0(jpjglo)   , mj1 (jpjglo) ,     &
          &      tpol(jpiglo)  , fpol(jpiglo)                                , STAT=ierr(2) )
          !
@@ -295,6 +293,6 @@
          &      e3t_1d  (jpk) , e3w_1d  (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , STAT=ierr(7) )
          !
-      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                       & 
-         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , &
+      ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                        & 
+         &      ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) ,     &
          &      mbkt   (jpi,jpj) , mbku   (jpi,jpj) , mbkv   (jpi,jpj) , STAT=ierr(9) )
          !
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 6717)
@@ -20,25 +20,25 @@
    !!   dom_nam       : read and contral domain namelists
    !!   dom_ctl       : control print for the ocean domain
-   !!   cfg_wri       : create the "domain_cfg.nc" file containing all required configuration information   
+   !!   cfg_write     : create the "domain_cfg.nc" file containing all required configuration information   
    !!----------------------------------------------------------------------
-   USE oce             ! ocean variables
-   USE dom_oce         ! domain: ocean
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE phycst          ! physical constants
-   USE closea          ! closed seas
-   USE domhgr          ! domain: set the horizontal mesh
-   USE domzgr          ! domain: set the vertical mesh
-   USE dommsk          ! domain: set the mask system
-   USE domwri          ! domain: write the meshmask file
-   USE domvvl          ! variable volume
-   USE c1d             ! 1D vertical configuration
-   USE dyncor_c1d      ! Coriolis term (c1d case)         (cor_c1d routine)
+   USE oce            ! ocean variables
+   USE dom_oce        ! domain: ocean
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE usrdef_closea  ! closed seas
+   USE domhgr         ! domain: set the horizontal mesh
+   USE domzgr         ! domain: set the vertical mesh
+   USE dommsk         ! domain: set the mask system
+   USE domwri         ! domain: write the meshmask file
+   USE domvvl         ! variable volume
+   USE c1d            ! 1D vertical configuration
+   USE dyncor_c1d     ! Coriolis term (c1d case)         (cor_c1d routine)
    !
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O library
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp         ! distributed memory computing library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
+   USE lib_mpp        ! distributed memory computing library
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -86,7 +86,7 @@
          WRITE(numout,*)     '      dimension of model'
          WRITE(numout,*)     '             Local domain      Global domain       Data domain '
-         WRITE(numout,cform) '        ','   jpi     : ', jpi, '   jpiglo  : ', jpiglo, '   jpidta  : ', jpidta
-         WRITE(numout,cform) '        ','   jpj     : ', jpj, '   jpjglo  : ', jpjglo, '   jpjdta  : ', jpjdta
-         WRITE(numout,cform) '        ','   jpk     : ', jpk, '   jpkglo  : ', jpkglo, '   jpkdta  : ', jpkdta
+         WRITE(numout,cform) '        ','   jpi     : ', jpi, '   jpiglo  : ', jpiglo
+         WRITE(numout,cform) '        ','   jpj     : ', jpj, '   jpjglo  : ', jpjglo
+         WRITE(numout,cform) '        ','   jpk     : ', jpk, '   jpkglo  : ', jpkglo
          WRITE(numout,cform) '       ' ,'   jpij    : ', jpij
          WRITE(numout,*)     '      mpp local domain info (mpp)'
@@ -100,7 +100,8 @@
       !      
       CALL dom_nam                     ! read namelist ( namrun, namdom )
-      CALL dom_clo                     ! Closed seas and lake
+      CALL dom_clo( cp_cfg, jp_cfg )   ! Closed seas and lake
       CALL dom_hgr                     ! Horizontal mesh
       CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry
+      IF( nn_closea == 0 )   CALL clo_bat( ik_top, ik_bot )    !==  remove closed seas or lakes  ==!
       CALL dom_msk( ik_top, ik_bot )   ! Masks
       !
@@ -171,5 +172,5 @@
       ENDIF
       !
-      IF( ln_write_cfg )   CALL cfg_wri           ! create the configuration file
+      IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file
       !
       IF( nn_timing == 1 )   CALL timing_stop('dom_init')
@@ -405,31 +406,17 @@
 
 
-   SUBROUTINE cfg_wri
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE cfg_wri  ***
+   SUBROUTINE cfg_write
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE cfg_write  ***
       !!                   
-      !! ** Purpose :   Create the NetCDF file(s) which contain(s) all the
-      !!      ocean domain informations (mesh and mask arrays). This (these)
-      !!      file(s) is (are) used for visualisation (SAXO software) and
-      !!      diagnostic computation.
+      !! ** Purpose :   Create the "domain_cfg" file, a NetCDF file which 
+      !!              contains all the ocean domain informations required to 
+      !!              define an ocean configuration.
       !!
-      !! ** Method  :   Write in a file all the arrays generated in routines
-      !!      domhgr, domzgr, and dommsk. Note: the file contain depends on
-      !!      the vertical coord. used (z-coord, partial steps, s-coord)
-      !!            MOD(nn_msh, 3) = 1  :   'mesh_mask.nc' file
-      !!                         = 2  :   'mesh.nc' and mask.nc' files
-      !!                         = 0  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
-      !!                                  'mask.nc' files
-      !!      For huge size domain, use option 2 or 3 depending on your 
-      !!      vertical coordinate.
+      !! ** Method  :   Write in a file all the arrays required to set up an
+      !!              ocean configuration.
       !!
-      !!      if     nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]
-      !!      if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 
-      !!                        corresponding to the depth of the bottom t- and w-points
-      !!      if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the
-      !!                        thickness (e3[tw]_ps) of the bottom points 
-      !!
-      !! ** output file :   meshmask.nc  : domain size, horizontal grid-point position,
-      !!                                   masks, depth and vertical scale factors
+      !! ** output file :   domain_cfg.nc : domain size, characteristics, horizontal mesh,
+      !!                              Coriolis parameter, depth and vertical scale factors
       !!----------------------------------------------------------------------
       INTEGER           ::   ji, jj, jk   ! dummy loop indices
@@ -441,6 +428,6 @@
       !
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'cfg_wri : create the "domain_cfg.nc" file containing all required configuration information'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
+      IF(lwp) WRITE(numout,*) 'cfg_write : create the "domain_cfg.nc" file containing all required configuration information'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
       !
       !                       ! ============================= !
@@ -532,5 +519,5 @@
       CALL iom_close( inum )
       !
-   END SUBROUTINE cfg_wri
+   END SUBROUTINE cfg_write
 
    !!======================================================================
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 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 6717)
@@ -26,5 +26,5 @@
    USE par_oce        ! ocean space and time domain
    USE phycst         ! physical constants
-   USE usrdef         ! User defined routine
+   USE usrdef_hgr     ! User defined routine
    !
    USE in_out_manager ! I/O manager
@@ -121,5 +121,5 @@
       ELSE
          IF( ln_read_cfg ) THEN
-            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been read in "mesh_mask" file'
+            IF(lwp) WRITE(numout,*) '          Coriolis parameter have been read in "domain_cfg" file'
          ELSE
             IF(lwp) WRITE(numout,*) '          Coriolis parameter have been set in usr_def_hgr routine'
@@ -210,5 +210,5 @@
       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 mesh_mask file'
+         IF(lwp) WRITE(numout,*) '           Coriolis factor at f- and t-points read in domain_cfg file'
          CALL iom_get( inum, jpdom_data, 'ff_f'  , ff_f  , lrowattr=ln_use_jattr )
          CALL iom_get( inum, jpdom_data, 'ff_t'  , ff_t  , lrowattr=ln_use_jattr )
@@ -219,5 +219,5 @@
       !
       IF( iom_varid( inum, 'e1e2u', ldstop = .FALSE. ) > 0 ) THEN
-         IF(lwp) WRITE(numout,*) '           e1e2u & e1e2v read in mesh_mask file'
+         IF(lwp) WRITE(numout,*) '           e1e2u & e1e2v read in domain_cfg 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 )
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 6717)
@@ -17,17 +17,19 @@
    !!            3.2  ! 2009-07  (R. Benshila) Suppression of rigid-lid option
    !!            3.6  ! 2015-05  (P. Mathiot) ISF: add wmask,wumask and wvmask
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dom_msk        : compute land/ocean mask
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
+   !!            4.0  ! 2016-06  (G. Madec, S. Flavoni)  domain configuration / user defined interface
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_msk       : compute land/ocean mask
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE usrdef_fmask   ! user defined fmask
    !
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         !
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! Massively Parallel Processing library
+   USE wrk_nemo       ! Memory allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -73,10 +75,5 @@
       !!                as MPP halos.
       !!      tmask_h : halo mask at t-point, i.e. excluding duplicated rows/lines
-      !!                due to cyclic or North Fold boundaries as well
-      !!                as MPP halos.
-      !!
-      !!      In case of open boundaries (lk_bdy=T):
-      !!        - tmask is set to 1 on the points to be computed by the open
-      !!          boundaries routines.
+      !!                due to cyclic or North Fold boundaries as well as MPP halos.
       !!
       !! ** Action :   tmask, umask, vmask, wmask, wumask, wvmask : land/ocean mask 
@@ -90,11 +87,10 @@
       INTEGER, DIMENSION(:,:), INTENT(in) ::   k_top, k_bot   ! first and last ocean level
       !
-      INTEGER  ::   ji, jj, jk               ! dummy loop indices
-      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
-      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
-      INTEGER  ::   iktop, ikbot             !   -       -
+      INTEGER  ::   ji, jj, jk     ! dummy loop indices
+      INTEGER  ::   iif, iil       ! local integers
+      INTEGER  ::   ijf, ijl       !   -       -
+      INTEGER  ::   iktop, ikbot   !   -       -
       INTEGER  ::   ios
-      INTEGER  ::   isrow                    ! index for ORCA1 starting row
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zwf   ! 2D workspace
       !!
       NAMELIST/namlbc/ rn_shlat, ln_vorlat
@@ -144,10 +140,44 @@
       END DO  
 
-      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
-      WHERE( k_bot(:,:) > 0 )   ;   ssmask(:,:) = 1._wp
-      ELSEWHERE                 ;   ssmask(:,:) = 0._wp
-      END WHERE
       
-      
+      !  Ocean/land mask at u-, v-, and f-points   (computed from tmask)
+      ! ----------------------------------------
+      ! NB: at this point, fmask is designed for free slip lateral boundary condition
+      DO jk = 1, jpk
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector loop
+               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
+               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
+            END DO
+            DO ji = 1, jpim1      ! NO vector opt.
+               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
+                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask  , 'V', 1._wp )
+      CALL lbc_lnk( fmask  , 'F', 1._wp )
+
+ 
+      ! Ocean/land mask at wu-, wv- and w points    (computed from tmask)
+      !-----------------------------------------
+      wmask (:,:,1) = tmask(:,:,1)     ! surface
+      wumask(:,:,1) = umask(:,:,1)
+      wvmask(:,:,1) = vmask(:,:,1)
+      DO jk = 2, jpk                   ! interior values
+         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
+      END DO
+
+
+      ! Ocean/land column mask at t-, u-, and v-points   (i.e. at least 1 wet cell in the vertical)
+      ! ----------------------------------------------
+      ssmask (:,:) = MAXVAL( tmask(:,:,:), DIM=3 )
+      ssumask(:,:) = MAXVAL( umask(:,:,:), DIM=3 )
+      ssvmask(:,:) = MAXVAL( vmask(:,:,:), DIM=3 )
+
+
       ! Interior domain mask  (used for global sum)
       ! --------------------
@@ -185,168 +215,51 @@
 
 
-      !  Ocean/land mask at u-, v-, and z-points (computed from tmask)
-      ! ----------------------------------------
-      DO jk = 1, jpk
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector loop
-               umask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)
-               vmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji  ,jj+1,jk)
-            END DO
-            DO ji = 1, jpim1      ! NO vector opt.
-               fmask(ji,jj,jk) = tmask(ji,jj  ,jk) * tmask(ji+1,jj  ,jk)   &
-                  &            * tmask(ji,jj+1,jk) * tmask(ji+1,jj+1,jk)
+      ! Lateral boundary conditions on velocity (modify fmask)
+      ! ---------------------------------------  
+      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition
+         !
+         CALL wrk_alloc( jpi,jpj,   zwf )
+         !
+         DO jk = 1, jpk
+            zwf(:,:) = fmask(:,:,jk)         
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  IF( fmask(ji,jj,jk) == 0._wp ) THEN
+                     fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
+                        &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
+                  ENDIF
+               END DO
+            END DO
+            DO jj = 2, jpjm1
+               IF( fmask(1,jj,jk) == 0._wp ) THEN
+                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+               ENDIF
+               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
+                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+               ENDIF
+            END DO         
+            DO ji = 2, jpim1
+               IF( fmask(ji,1,jk) == 0._wp ) THEN
+                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+               ENDIF
+               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
+                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+               ENDIF
             END DO
          END DO
-      END DO
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet cell at u point
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector loop
-!!gm  simpler :
-!            ssumask(ji,jj)  = MIN(  1._wp , SUM( umask(ji,jj,:) )  )
-!            ssvmask(ji,jj)  = MIN(  1._wp , SUM( vmask(ji,jj,:) )  )
-!!gm
-!!gm  faster :
-!         ssumask(ji,jj) = ssmask(ji,jj) * tmask(ji+1,jj  )
-!         ssvmask(ji,jj) = ssmask(ji,jj) * tmask(ji  ,jj+1)
-!!gm
-            ssumask(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            ssvmask(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
-!!end
-         END DO
-         DO ji = 1, jpim1      ! NO vector opt.
-!!gm faster
-!            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
-!               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1)
-!!gm 
-            ssfmask(ji,jj) =  ssmask(ji,jj  ) * ssmask(ji+1,jj  )   &
-               &            * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:)))
-!!gm
-         END DO
-      END DO
-      CALL lbc_lnk( umask  , 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask  , 'V', 1._wp )
-!      CALL lbc_lnk( fmask  , 'F', 1._wp )         ! applied after the specification of lateral b.c.
-      CALL lbc_lnk( ssumask, 'U', 1._wp )
-      CALL lbc_lnk( ssvmask, 'V', 1._wp )
-      CALL lbc_lnk( ssfmask, 'F', 1._wp )
-
-
-      ! Ocean/land mask at wu-, wv- and w points 
-      !----------------------------------------------
-      wmask (:,:,1) = tmask(:,:,1)     ! surface
-      wumask(:,:,1) = umask(:,:,1)
-      wvmask(:,:,1) = vmask(:,:,1)
-      DO jk = 2, jpk                   ! interior values
-         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
-      END DO
-
-
-      ! Lateral boundary conditions on velocity (modify fmask)
-      ! ---------------------------------------     
-      CALL wrk_alloc( jpi,jpj,   zwf  )
-      !
-      DO jk = 1, jpk
-         zwf(:,:) = fmask(:,:,jk)         
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( fmask(ji,jj,jk) == 0._wp ) THEN
-                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1),   &
-                     &                                           zwf(ji-1,jj), zwf(ji,jj-1)  )  )
-               ENDIF
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            IF( fmask(1,jj,jk) == 0._wp ) THEN
-               fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
-            ENDIF
-            IF( fmask(jpi,jj,jk) == 0._wp ) THEN
-               fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
-            ENDIF
-         END DO         
-         DO ji = 2, jpim1
-            IF( fmask(ji,1,jk) == 0._wp ) THEN
-               fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
-            ENDIF
-            IF( fmask(ji,jpj,jk) == 0._wp ) THEN
-               fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
-            ENDIF
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj,   zwf  )
-      !
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
-         !                                                 ! Increased lateral friction near of some straits
-         !                                ! Gibraltar strait  : partial slip (fmask=0.5)
-         ij0 = 101   ;   ij1 = 101
-         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
-         ij0 = 102   ;   ij1 = 102
-         ii0 = 139   ;   ii1 = 140   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
-         !
-         !                                ! Bab el Mandeb : partial slip (fmask=1)
-         ij0 =  87   ;   ij1 =  88
-         ii0 = 160   ;   ii1 = 160   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
-         ij0 =  88   ;   ij1 =  88
-         ii0 = 159   ;   ii1 = 159   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
-         !
-         !                                ! Danish straits  : strong slip (fmask > 2)
-! We keep this as an example but it is instable in this case 
-!         ij0 = 115   ;   ij1 = 115
-!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
-!         ij0 = 116   ;   ij1 = 116
-!         ii0 = 145   ;   ii1 = 146   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
-         !
-      ENDIF
-      !
-      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
-         !                                                 ! Increased lateral friction near of some straits
-         ! This dirty section will be suppressed by simplification process:
-         ! all this will come back in input files
-         ! Currently these hard-wired indices relate to configuration with
-         ! extend grid (jpjglo=332)
-         !
-         isrow = 332 - jpjglo
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
-         IF(lwp) WRITE(numout,*) '      Gibraltar '
-         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait 
-         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
-
-         IF(lwp) WRITE(numout,*) '      Bhosporus '
-         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait 
-         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
-
-         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
-         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top) 
-         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
-
-         IF(lwp) WRITE(numout,*) '      Lombok '
-         ii0 =  44           ;   ii1 =  44        ! Lombok Strait 
-         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
-
-         IF(lwp) WRITE(numout,*) '      Ombai '
-         ii0 =  53           ;   ii1 =  53        ! Ombai Strait 
-         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
-
-         IF(lwp) WRITE(numout,*) '      Timor Passage '
-         ii0 =  56           ;   ii1 =  56        ! Timor Passage 
-         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
-
-         IF(lwp) WRITE(numout,*) '      West Halmahera '
-         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait 
-         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
-
-         IF(lwp) WRITE(numout,*) '      East Halmahera '
-         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait 
-         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
-         !
-      ENDIF
-      !
-      CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
-      !
-      ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 )
+         !
+         CALL wrk_dealloc( jpi,jpj,   zwf )
+         !
+         CALL lbc_lnk( fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
+         !
+         ! CAUTION : The fmask may be further modified in dyn_vor_init ( dynvor.F90 ) depending on ln_vorlat
+         !
+      ENDIF
+      
+      ! User defined alteration of fmask (use to reduce ocean transport in specified straits)
+      ! -------------------------------- 
+      !
+      CALL usr_def_fmask( cp_cfg, jp_cfg, fmask )
+      !
       !
       IF( nn_timing == 1 )  CALL timing_stop('dom_msk')
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90	(revision 6717)
@@ -8,4 +8,5 @@
    !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90 and several file
    !!            3.0  ! 2008-01  (S. Masson)  add dom_uniq 
+   !!            4.0  ! 2016-01  (G. Madec)  simplified mesh_mask.nc file
    !!----------------------------------------------------------------------
 
@@ -27,5 +28,4 @@
 
    PUBLIC   dom_wri              ! routine called by inidom.F90
-   PUBLIC   dom_wri_coordinate   ! routine called by domhgr.F90
    PUBLIC   dom_stiff            ! routine called by inidom.F90
 
@@ -33,74 +33,9 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE dom_wri_coordinate
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_wri_coordinate  ***
-      !!                   
-      !! ** Purpose :   Create the NetCDF file which contains all the
-      !!              standard coordinate information plus the surface,
-      !!              e1e2u and e1e2v. By doing so, those surface will
-      !!              not be changed by the reduction of e1u or e2v scale 
-      !!              factors in some straits. 
-      !!                 NB: call just after the read of standard coordinate
-      !!              and the reduction of scale factors in some straits
-      !!
-      !! ** output file :   coordinate_e1e2u_v.nc
-      !!----------------------------------------------------------------------
-      INTEGER           ::   inum0    ! temprary units for 'coordinate_e1e2u_v.nc' file
-      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
-      !                                   !  workspaces
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dom_wri_coordinate')
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dom_wri_coordinate : create NetCDF coordinate file'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~'
-      
-      clnam0 = 'coordinate_e1e2u_v'  ! filename (mesh and mask informations)
-      
-      !  create 'coordinate_e1e2u_v.nc' file
-      ! ============================
-      !
-      CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
-      !
-      !                                                         ! horizontal mesh (inum3)
-      CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
-      CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
-      CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
-      CALL iom_rstput( 0, 0, inum0, 'e1u', e1u, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'e1v', e1v, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'e1f', e1f, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum0, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
-      CALL iom_rstput( 0, 0, inum0, 'e2u', e2u, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'e2v', e2v, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'e2f', e2f, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum0, 'e1e2u', e1e2u, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum0, 'e1e2v', e1e2v, ktype = jp_r8 )
-
-      CALL iom_close( inum0 )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dom_wri_coordinate')
-      !
-   END SUBROUTINE dom_wri_coordinate
-
 
    SUBROUTINE dom_wri
@@ -132,19 +67,11 @@
       !!                                   masks, depth and vertical scale factors
       !!----------------------------------------------------------------------
-      INTEGER           ::   inum0    ! temprary units for 'mesh_mask.nc' file
-      INTEGER           ::   inum1    ! temprary units for 'mesh.nc'      file
-      INTEGER           ::   inum2    ! temprary units for 'mask.nc'      file
-      INTEGER           ::   inum3    ! temprary units for 'mesh_hgr.nc'  file
-      INTEGER           ::   inum4    ! temprary units for 'mesh_zgr.nc'  file
-      CHARACTER(len=21) ::   clnam0   ! filename (mesh and mask informations)
-      CHARACTER(len=21) ::   clnam1   ! filename (mesh informations)
-      CHARACTER(len=21) ::   clnam2   ! filename (mask informations)
-      CHARACTER(len=21) ::   clnam3   ! filename (horizontal mesh informations)
-      CHARACTER(len=21) ::   clnam4   ! filename (vertical   mesh informations)
+      INTEGER           ::   inum    ! temprary units for 'mesh_mask.nc' file
+      CHARACTER(len=21) ::   clnam   ! filename (mesh and mask informations)
       INTEGER           ::   ji, jj, jk   ! dummy loop indices
       INTEGER           ::   izco, izps, isco, icav
-      !                                   !  workspaces
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zprt, zprw 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepu, zdepv
+      !                               
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zprt, zprw     ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdepu, zdepv   ! 3D workspace
       !!----------------------------------------------------------------------
       !
@@ -158,212 +85,115 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~'
       
-      clnam0 = 'mesh_mask'  ! filename (mesh and mask informations)
-      clnam1 = 'mesh'       ! filename (mesh informations)
-      clnam2 = 'mask'       ! filename (mask informations)
-      clnam3 = 'mesh_hgr'   ! filename (horizontal mesh informations)
-      clnam4 = 'mesh_zgr'   ! filename (vertical   mesh informations)
-      
-      SELECT CASE ( MOD(nn_msh, 3) )
-         !                                  ! ============================
-      CASE ( 1 )                            !  create 'mesh_mask.nc' file
-         !                                  ! ============================
-         CALL iom_open( TRIM(clnam0), inum0, ldwrt = .TRUE., kiolib = jprstlib )
-         inum2 = inum0                                            ! put all the informations
-         inum3 = inum0                                            ! in unit inum0
-         inum4 = inum0
-         
-         !                                  ! ============================
-      CASE ( 2 )                            !  create 'mesh.nc' and 
-         !                                  !         'mask.nc' files
-         !                                  ! ============================
-         CALL iom_open( TRIM(clnam1), inum1, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
-         inum3 = inum1                                            ! put mesh informations 
-         inum4 = inum1                                            ! in unit inum1 
-         !                                  ! ============================
-      CASE ( 0 )                            !  create 'mesh_hgr.nc'
-         !                                  !         'mesh_zgr.nc' and
-         !                                  !         'mask.nc'     files
-         !                                  ! ============================
-         CALL iom_open( TRIM(clnam2), inum2, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_open( TRIM(clnam3), inum3, ldwrt = .TRUE., kiolib = jprstlib )
-         CALL iom_open( TRIM(clnam4), inum4, ldwrt = .TRUE., kiolib = jprstlib )
-         !
-      END SELECT
-
+      clnam = 'mesh_mask'  ! filename (mesh and mask informations)
+      
+      !                                  ! ============================
+      !                                  !  create 'mesh_mask.nc' file
+      !                                  ! ============================
+      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
+      !
       !                                                         ! global domain size
-      CALL iom_rstput( 0, 0, inum2, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
-      CALL iom_rstput( 0, 0, inum2, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
-      CALL iom_rstput( 0, 0, inum2, 'jpkglo', REAL( jpk   , wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpkglo, wp), ktype = jp_i4 )
 
       !                                                         ! domain characteristics
-      CALL iom_rstput( 0, 0, inum2, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 )
       !                                                         ! type of vertical coordinate
       IF( ln_zco    ) THEN   ;   izco = 1   ;   ELSE   ;   izco = 0   ;   ENDIF
       IF( ln_zps    ) THEN   ;   izps = 1   ;   ELSE   ;   izps = 0   ;   ENDIF
       IF( ln_sco    ) THEN   ;   isco = 1   ;   ELSE   ;   isco = 0   ;   ENDIF
-      CALL iom_rstput( 0, 0, inum2, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
-      CALL iom_rstput( 0, 0, inum2, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
-      CALL iom_rstput( 0, 0, inum2, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_zco'   , REAL( izco, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_zps'   , REAL( izps, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_sco'   , REAL( isco, wp), ktype = jp_i4 )
       !                                                         ! ocean cavities under iceshelves
       IF( ln_isfcav ) THEN   ;   icav = 1   ;   ELSE   ;   icav = 0   ;   ENDIF
-      CALL iom_rstput( 0, 0, inum2, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
+      CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
   
-      !                                                         ! masks (inum2) 
-      CALL iom_rstput( 0, 0, inum2, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
-      CALL iom_rstput( 0, 0, inum2, 'umask', umask, ktype = jp_i1 )
-      CALL iom_rstput( 0, 0, inum2, 'vmask', vmask, ktype = jp_i1 )
-      CALL iom_rstput( 0, 0, inum2, 'fmask', fmask, ktype = jp_i1 )
+      !                                                         ! masks
+      CALL iom_rstput( 0, 0, inum, 'tmask', tmask, ktype = jp_i1 )     !    ! land-sea mask
+      CALL iom_rstput( 0, 0, inum, 'umask', umask, ktype = jp_i1 )
+      CALL iom_rstput( 0, 0, inum, 'vmask', vmask, ktype = jp_i1 )
+      CALL iom_rstput( 0, 0, inum, 'fmask', fmask, ktype = jp_i1 )
       
       CALL dom_uniq( zprw, 'T' )
       DO jj = 1, jpj
          DO ji = 1, jpi
-            jk=mikt(ji,jj) 
-            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+            zprt(ji,jj) = ssmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
          END DO
       END DO                             !    ! unique point mask
-      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )  
+      CALL iom_rstput( 0, 0, inum, 'tmaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'U' )
       DO jj = 1, jpj
          DO ji = 1, jpi
-            jk=miku(ji,jj) 
-            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+            zprt(ji,jj) = ssumask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
          END DO
       END DO
-      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )  
+      CALL iom_rstput( 0, 0, inum, 'umaskutil', zprt, ktype = jp_i1 )  
       CALL dom_uniq( zprw, 'V' )
       DO jj = 1, jpj
          DO ji = 1, jpi
-            jk=mikv(ji,jj) 
-            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
+            zprt(ji,jj) = ssvmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
          END DO
       END DO
-      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )  
-      CALL dom_uniq( zprw, 'F' )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            jk=mikf(ji,jj) 
-            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask
-         END DO
-      END DO
-      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )  
+      CALL iom_rstput( 0, 0, inum, 'vmaskutil', zprt, ktype = jp_i1 )  
+!!gm  ssfmask has been removed  ==>> find another solution to defined fmaskutil
+!!    Here we just remove the output of fmaskutil.
+!      CALL dom_uniq( zprw, 'F' )
+!      DO jj = 1, jpj
+!         DO ji = 1, jpi
+!            zprt(ji,jj) = ssfmask(ji,jj) * zprw(ji,jj)                        !    ! unique point mask
+!         END DO
+!      END DO
+!      CALL iom_rstput( 0, 0, inum, 'fmaskutil', zprt, ktype = jp_i1 )  
+!!gm
 
       !                                                         ! horizontal mesh (inum3)
-      CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
-      CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
-      CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
-      CALL iom_rstput( 0, 0, inum3, 'e1u', e1u, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'e1v', e1v, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'e1f', e1f, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum3, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
-      CALL iom_rstput( 0, 0, inum3, 'e2u', e2u, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'e2v', e2v, ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum3, 'e2f', e2f, ktype = jp_r8 )
-      
-      CALL iom_rstput( 0, 0, inum3, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor
-      CALL iom_rstput( 0, 0, inum3, 'ff_t', ff_t, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'glamt', glamt, ktype = jp_r8 )     !    ! latitude
+      CALL iom_rstput( 0, 0, inum, 'glamu', glamu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'glamv', glamv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'glamf', glamf, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum, 'gphit', gphit, ktype = jp_r8 )     !    ! longitude
+      CALL iom_rstput( 0, 0, inum, 'gphiu', gphiu, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gphiv', gphiv, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gphif', gphif, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum, 'e1t', e1t, ktype = jp_r8 )         !    ! e1 scale factors
+      CALL iom_rstput( 0, 0, inum, 'e1u', e1u, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e1v', e1v, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e1f', e1f, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum, 'e2t', e2t, ktype = jp_r8 )         !    ! e2 scale factors
+      CALL iom_rstput( 0, 0, inum, 'e2u', e2u, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e2v', e2v, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'e2f', e2f, ktype = jp_r8 )
+      
+      CALL iom_rstput( 0, 0, inum, 'ff_f', ff_f, ktype = jp_r8 )           !    ! coriolis factor
+      CALL iom_rstput( 0, 0, inum, 'ff_t', ff_t, ktype = jp_r8 )
       
       ! note that mbkt is set to 1 over land ==> use surface tmask
       zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum, 'mbathy', zprt, ktype = jp_i4 )     !    ! nb of ocean T-points
       zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum, 'misf', zprt, ktype = jp_i4 )       !    ! nb of ocean T-points
       zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp )
-      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points
+      CALL iom_rstput( 0, 0, inum, 'isfdraft', zprt, ktype = jp_r8 )       !    ! nb of ocean T-points
             
-      IF( ln_sco ) THEN                                         ! s-coordinate
-         CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         !    ! scale factors
-         CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
-         CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
-         CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
-         !
-         CALL dom_stiff( zprt )
-         CALL iom_rstput( 0, 0, inum4, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
-         !
-         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d )
-         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
-      ENDIF
-      
-      IF( ln_zps ) THEN                                         ! z-coordinate - partial steps
-         !
-         IF( nn_msh <= 6 ) THEN                                   !    ! 3D vertical scale factors
-            CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 )         
-            CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 )
-            CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 )
-            CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 )
-         ELSE                                                   !    ! 2D masked bottom ocean scale factors
-            DO jj = 1,jpj   
-               DO ji = 1,jpi
-                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
-                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * ssmask(ji,jj)
-               END DO
-            END DO
-            CALL iom_rstput( 0, 0, inum4, 'e3t_ps', e3tp )      
-            CALL iom_rstput( 0, 0, inum4, 'e3w_ps', e3wp )
-         END IF
-         !
-         IF( nn_msh <= 3 ) THEN                                   !    ! 3D depth
-            CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 )
-            DO jk = 1,jpk   
-               DO jj = 1, jpjm1   
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj  ,jk) )
-                     zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji  ,jj+1,jk) )
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( zdepu, 'U', 1. )   ;   CALL lbc_lnk( zdepv, 'V', 1. )
-            CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 )
-            CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 )
-            CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 )
-         ELSE                                                   !    ! 2D bottom depth
-            DO jj = 1,jpj   
-               DO ji = 1,jpi
-                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * ssmask(ji,jj)
-                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * ssmask(ji,jj)
-               END DO
-            END DO
-            CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 )
-            CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 )
-         ENDIF
-         !
-         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! reference z-coord.
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
-         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )
-         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
-      ENDIF
-      
-      IF( ln_zco ) THEN
-         !                                                      ! z-coordinate - full steps
-         CALL iom_rstput( 0, 0, inum4, 'gdept_1d', gdept_1d )   !    ! depth
-         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d', gdepw_1d )
-         CALL iom_rstput( 0, 0, inum4, 'e3t_1d'  , e3t_1d   )   !    ! scale factors
-         CALL iom_rstput( 0, 0, inum4, 'e3w_1d'  , e3w_1d   )
-      ENDIF
+      CALL iom_rstput( 0, 0, inum, 'e3t_0', e3t_0 )         !    ! scale factors
+      CALL iom_rstput( 0, 0, inum, 'e3u_0', e3u_0 )
+      CALL iom_rstput( 0, 0, inum, 'e3v_0', e3v_0 )
+      CALL iom_rstput( 0, 0, inum, 'e3w_0', e3w_0 )
+      !
+      CALL dom_stiff( zprt )
+      CALL iom_rstput( 0, 0, inum, 'stiffness', zprt )      !    ! Max. grid stiffness ratio
+      !
+      CALL iom_rstput( 0, 0, inum, 'gdept_1d' , gdept_1d )  !    ! stretched system
+      CALL iom_rstput( 0, 0, inum, 'gdepw_1d' , gdepw_1d )
+      CALL iom_rstput( 0, 0, inum, 'gdept_0', gdept_0, ktype = jp_r8 )
+      CALL iom_rstput( 0, 0, inum, 'gdepw_0', gdepw_0, ktype = jp_r8 )
+      
       !                                     ! ============================
-      !                                     !        close the files 
+      CALL iom_close( inum )                !        close the files 
       !                                     ! ============================
-      SELECT CASE ( MOD(nn_msh, 3) )
-      CASE ( 1 )                
-         CALL iom_close( inum0 )
-      CASE ( 2 )
-         CALL iom_close( inum1 )
-         CALL iom_close( inum2 )
-      CASE ( 0 )
-         CALL iom_close( inum2 )
-         CALL iom_close( inum3 )
-         CALL iom_close( inum4 )
-      END SELECT
       !
       CALL wrk_dealloc( jpi, jpj, zprt, zprw )
@@ -384,5 +214,4 @@
       !!                2) check which elements have been changed
       !!----------------------------------------------------------------------
-      !
       CHARACTER(len=1)        , INTENT(in   ) ::   cdgrd   ! 
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   puniq   ! 
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 6717)
@@ -24,5 +24,5 @@
    !!   dom_zgr       : read or set the ocean vertical coordinate system
    !!   zgr_read      : read the vertical domain coordinate and mask in domain_cfg file
-   !!   zgr_tb_level  : ocean top and bottom level for t-, u, and v-points with 1 as minimum value
+   !!   zgr_top_bot   : ocean top and bottom level for t-, u, and v-points with 1 as minimum value
    !!---------------------------------------------------------------------
    USE oce            ! ocean variables
@@ -73,19 +73,8 @@
       INTEGER  ::   ioptio, ibat, ios   ! local integer
       REAL(wp) ::   zrefdep             ! depth of the reference level (~10m)
-      !!
-!      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('dom_zgr')
       !
-!      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
-!      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
-!901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
-!
-!      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
-!      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
-!902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
-!      IF(lwm) WRITE ( numond, namzgr )
-
       IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
@@ -110,5 +99,5 @@
       ELSE                          !==  User defined configuration  ==!
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '          User defined horizontal mesh (usr_def_hgr)'
+         IF(lwp) WRITE(numout,*) '          User defined vertical mesh (usr_def_zgr)'
          !
          CALL usr_def_zgr( ln_zco  , ln_zps  , ln_sco, ln_isfcav,   & 
@@ -130,5 +119,5 @@
       IF(lwp) THEN                     ! Control print
          WRITE(numout,*)
-         WRITE(numout,*) '   Read in domain_cfg.nc or user defined type of vertical coordinate:'
+         WRITE(numout,*) '   Type of vertical coordinate (read in domain_cfg.nc or set through user defined routines) :'
          WRITE(numout,*) '      z-coordinate - full steps      ln_zco    = ', ln_zco
          WRITE(numout,*) '      z-coordinate - partial steps   ln_zps    = ', ln_zps
@@ -145,5 +134,5 @@
 
       !                                ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
-      CALL zgr_tb_level( k_top, k_bot )      ! with a minimum value set to 1
+      CALL zgr_top_bot( k_top, k_bot )      ! with a minimum value set to 1
       
 
@@ -207,6 +196,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'hgr_read : read the vertical coordinates in "domain_cfg.nc" file'
-         WRITE(numout,*) '~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
+         WRITE(numout,*) '   zgr_read : read the vertical coordinates in "domain_cfg.nc" file'
+         WRITE(numout,*) '   ~~~~~~~~      jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
       ENDIF
       !
@@ -256,7 +245,7 @@
 
 
-   SUBROUTINE zgr_tb_level( k_top, k_bot )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_tb_level  ***
+   SUBROUTINE zgr_top_bot( k_top, k_bot )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_top_bot  ***
       !!
       !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
@@ -282,6 +271,6 @@
       !
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_tb_level : ocean top and bottom k-index of T-, U-, V- and W-levels '
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) '    zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels '
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
       !
       mikt(:,:) = MAX( k_top(:,:) , 1 )    ! top    ocean k-index of T-level (=1 over land)
@@ -313,5 +302,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level')
       !
-   END SUBROUTINE zgr_tb_level
+   END SUBROUTINE zgr_top_bot
 
    !!======================================================================
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 6717)
@@ -14,6 +14,5 @@
    !!            3.3  !  2010-10  (C. Ethe) merge TRC-TRA
    !!            3.4  !  2011-04  (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 
-   !!            3.7  !  2016-04  (S. Flavoni) change configuration's interface:
-   !!                              read file or CALL usr_def module to compute initial state (example given for GYRE)
+   !!            3.7  !  2016-04  (S. Flavoni) introduce user defined initial state 
    !!----------------------------------------------------------------------
 
@@ -22,17 +21,13 @@
    !!   istate_uvg    : initial velocity in geostropic balance
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers 
-   USE dom_oce         ! ocean space and time domain 
-   USE c1d             ! 1D vertical configuration
-   USE daymod          ! calendar
-   USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
-   USE ldftra          ! lateral physics: ocean active tracers
-   USE zdf_oce         ! ocean vertical physics
-   USE phycst          ! physical constants
-   USE dtatsd          ! data temperature and salinity   (dta_tsd routine)
-   USE dtauvd          ! data: U & V current             (dta_uvd routine)
+   USE oce            ! ocean dynamics and active tracers 
+   USE dom_oce        ! ocean space and time domain 
+   USE daymod         ! calendar
+   USE divhor         ! horizontal divergence            (div_hor routine)
+   USE dtatsd         ! data temperature and salinity   (dta_tsd routine)
+   USE dtauvd         ! data: U & V current             (dta_uvd routine)
    USE domvvl          ! varying vertical mesh
    USE iscplrst        ! ice sheet coupling
-   USE usrdef          ! User defined routine
+   USE usrdef_istate   ! User defined initial state
    !
    USE in_out_manager  ! I/O manager
@@ -47,5 +42,4 @@
 
    PUBLIC   istate_init   ! routine called by step.F90
-   !SF PUBLIC   ini_read      ! subroutine ini_read 
 
    !! * Substitutions
@@ -70,12 +64,14 @@
       IF( nn_timing == 1 )   CALL timing_start('istate_init')
       !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'istate_init : Initialization of the dynamics and tracers'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
-
-     !SF initialisation of T & S with data file
+!!gm  Why not include in the first call of dta_tsd ?  
+!!gm  probably associated with the use of internal damping...
                      CALL dta_tsd_init        ! Initialisation of T & S input data
-      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
+!!gm to be moved in usrdef of C1D case
+!      IF( lk_c1d )   CALL dta_uvd_init        ! Initialization of U & V input data
+!!gm
 
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
@@ -87,8 +83,8 @@
          !                                    ! -------------------
          CALL rst_read                        ! Read the restart file
-         IF (ln_iscpl)       CALL iscpl_stp   ! extraloate restart to wet and dry
+         IF (ln_iscpl)       CALL iscpl_stp   ! extrapolate restart to wet and dry
          CALL day_init                        ! model calendar (using both namelist and restart infos)
-      ELSE
-         !                                    ! Start from rest
+         !
+      ELSE                                    ! Start from rest
          !                                    ! ---------------
          numror = 0                           ! define numror = 0 -> no restart file to read
@@ -96,32 +92,46 @@
          CALL day_init                        ! model calendar (using both namelist and restart infos)
          !                                    ! Initialization of ocean to zero
-         !   before fields      !       now fields     
-         sshb (:,:)   = 0._wp   ;   sshn (:,:)   = 0._wp
-         ub   (:,:,:) = 0._wp   ;   un   (:,:,:) = 0._wp
-         vb   (:,:,:) = 0._wp   ;   vn   (:,:,:) = 0._wp  
-                                    hdivn(:,:,:) = 0._wp
          !
-         IF( ln_tsd_init ) THEN               ! read 3D T and S data at nit000
-            CALL dta_tsd( nit000, tsb )  
+         IF( ln_tsd_init ) THEN               
+            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
+            !
+            sshb(:,:)   = 0._wp               ! set the ocean at rest
+            ub  (:,:,:) = 0._wp
+            vb  (:,:,:) = 0._wp  
+            !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_ini( tsb )                
+            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
          ENDIF
-         tsn(:,:,:,:) = tsb(:,:,:,:)          ! set now to before values
+         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
+         sshn (:,:)     = sshb(:,:)   
+         un   (:,:,:)   = ub  (:,:,:)
+         vn   (:,:,:)   = vb  (:,:,:)
+         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
+         CALL div_hor( 0 )                    ! compute interior hdivn value  
+!!gm                                    hdivn(:,:,:) = 0._wp
+
+!!gm POTENTIAL BUG :
+!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
+!!             as well as gdept and gdepw....   !!!!! 
+!!      ===>>>>   probably a call to domvvl initialisation here....
+
+
          !
-         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
-            CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
-            CALL dta_uvd( nit000, zuvd )
-            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
-            CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
-         ENDIF
+!!gm to be moved in usrdef of C1D case
+!         IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000
+!            CALL wrk_alloc( jpi,jpj,jpk,2,   zuvd )
+!            CALL dta_uvd( nit000, zuvd )
+!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
+!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+!            CALL wrk_dealloc( jpi,jpj,jpk,2,   zuvd )
+!         ENDIF
          !
 !!gm This is to be changed !!!!
-         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
-         IF( .NOT.ln_linssh ) THEN
-            DO jk = 1, jpk
-               e3t_b(:,:,jk) = e3t_n(:,:,jk)
-            END DO
-         ENDIF
+!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+!         IF( .NOT.ln_linssh ) THEN
+!            DO jk = 1, jpk
+!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+!            END DO
+!         ENDIF
 !!gm 
          ! 
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6717)
@@ -789,5 +789,5 @@
                   ENDIF
                   IF( PRESENT(pv_r3d) ) THEN
-                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkdta
+                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkglo
                      ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3)
                      ELSE                                                           ; icnt(3) = jpk
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 6717)
@@ -18,5 +18,5 @@
    PRIVATE
 
-   INTEGER, PARAMETER, PUBLIC ::   jpdom_data          = 1   !: ( 1  :jpidta, 1  :jpjdta)
+   INTEGER, PARAMETER, PUBLIC ::   jpdom_data          = 1   !: ( 1  :jpiglo, 1  :jpjglo)    !!gm to be suppressed
    INTEGER, PARAMETER, PUBLIC ::   jpdom_global        = 2   !: ( 1  :jpiglo, 1  :jpjglo)
    INTEGER, PARAMETER, PUBLIC ::   jpdom_local         = 3   !: One of the 3 following cases
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 6717)
@@ -6,5 +6,5 @@
       !!     FOR USING THIS VERSION, A PREPROCESSING TRAITMENT IS RECOMMENDED
       !!     FOR DEFINING BETTER CUTTING OUT.
-      !!       This routine is used with a the bathymetry file.
+      !!       This routine requires the presence of the "domain_cfg.nc" file.
       !!       In this version, the land processors are avoided and the adress
       !!     processor (nproc, narea,noea, ...) are calculated again.
@@ -32,10 +32,10 @@
       !!                    nono      : number for local neighboring processor
       !!
-      !! History :
-      !!        !  94-11  (M. Guyon)  Original code
-      !!        !  95-04  (J. Escobar, M. Imbard)
-      !!        !  98-02  (M. Guyon)  FETI method
-      !!        !  98-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
-      !!   9.0  !  04-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
+      !! History :       !  1994-11  (M. Guyon)  Original code
+      !!  OPA            !  1995-04  (J. Escobar, M. Imbard)
+      !!                 !  1998-02  (M. Guyon)  FETI method
+      !!                 !  1998-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
+      !!  NEMO      1.0  !  2004-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
+      !!            4.0  !  2016-06  (G. Madec)  use domain_cfg file instead of bathymetry file
       !!----------------------------------------------------------------------
       USE in_out_manager  ! I/O Manager
@@ -65,75 +65,34 @@
          ione  , ionw  , iose  , iosw  ,   &  !    "           "
          ibne  , ibnw  , ibse  , ibsw         !    "           "
-      INTEGER,  DIMENSION(jpiglo,jpjglo) ::   &
-         imask                                ! temporary global workspace
-      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   &
-         zdta, zdtaisf                     ! temporary data workspace
-      REAL(wp) ::   zidom , zjdom          ! temporary scalars
-
-      ! read namelist for ln_zco
-      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
-
+      INTEGER,  DIMENSION(jpiglo,jpjglo) ::   imask        ! global workspace
+      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zbot, ztop   ! global workspace
+      REAL(wp) ::   zidom , zjdom          ! local scalars
       !!----------------------------------------------------------------------
-      !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+      !! NEMO/OPA 4.0 , NEMO Consortium (2016)
       !! $Id$
-      !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+      !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
       !!----------------------------------------------------------------------
 
-      REWIND( numnam_ref )              ! Namelist namzgr in reference namelist : Vertical coordinate
-      READ  ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
-
-      REWIND( numnam_cfg )              ! Namelist namzgr in configuration namelist : Vertical coordinate
-      READ  ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namzgr )
-
       IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'mpp_init : Message Passing MPI'
-      IF(lwp)WRITE(numout,*) '~~~~~~~~'
+      IF(lwp)WRITE(numout,*) 'mpp_init_2 : Message Passing MPI'
+      IF(lwp)WRITE(numout,*) '~~~~~~~~~~'
       IF(lwp)WRITE(numout,*) ' '
 
-      IF( jpni*jpnj < jpnij ) CALL ctl_stop( ' jpnij > jpni x jpnj impossible' )
+      IF( jpni*jpnj < jpnij )   CALL ctl_stop( ' jpnij > jpni x jpnj impossible' )
 
       ! 0. initialisation
       ! -----------------
-
-      ! open the file
-      ! Remember that at this level in the code, mpp is not yet initialized, so
-      ! the file must be open with jpdom_unknown, and kstart and kcount forced 
-      jstartrow = 1
-      IF ( ln_zco ) THEN 
-         CALL iom_open ( 'bathy_level.nc', inum )   ! Level bathymetry
-          ! Optionally use a file attribute (open_ocean_jstart) to set a start row for reading from the global file
-          ! This allows the unextended grid bathymetry to be stored in the same file as the under ice-shelf extended bathymetry
-         CALL iom_getatt(inum, 'open_ocean_jstart', jstartrow ) ! -999 is returned if the attribute is not found
-         jstartrow = MAX(1,jstartrow)
-         CALL iom_get( inum, jpdom_unknown, 'Bathy_level', zdta, kstart=(/1,1+jstartrow-1/), kcount=(/jpiglo,jpjglo/) )
-      ELSE
-         CALL iom_open ( 'bathy_meter.nc', inum )   ! Meter bathy in case of partial steps
-         IF ( ln_isfcav ) THEN
-             CALL iom_get( inum, jpdom_unknown, 'Bathymetry_isf' , zdta, kstart=(/1,1/), kcount=(/jpiglo,jpjglo/) )
-         ELSE
-             ! Optionally use a file attribute (open_ocean_jstart) to set a start row for reading from the global file
-             ! This allows the unextended grid bathymetry to be stored in the same file as the under ice-shelf extended bathymetry
-             CALL iom_getatt(inum, 'open_ocean_jstart', jstartrow ) ! -999 is returned if the attribute is not found
-             jstartrow = MAX(1,jstartrow)
-             CALL iom_get( inum, jpdom_unknown, 'Bathymetry' , zdta, kstart=(/1,1+jstartrow-1/), kcount=(/jpiglo,jpjglo/) )
-         ENDIF
-      ENDIF
-      CALL iom_close (inum)
-      
-      ! used to compute the land processor in case of not masked bathy file.
-      zdtaisf(:,:) = 0.0_wp
-      IF ( ln_isfcav ) THEN
-         CALL iom_open ( 'bathy_meter.nc', inum )   ! Meter bathy in case of partial steps
-         CALL iom_get ( inum, jpdom_unknown, 'isf_draft' , zdtaisf, kstart=(/1,1/), kcount=(/jpiglo,jpjglo/) )
-      END IF
-      CALL iom_close (inum)
-
-      ! land/sea mask over the global domain
-
-      imask(:,:)=1
-      WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0
+      CALL iom_open( 'domain_cfg', inum )
+      !
+      !                                   ! ocean top and bottom level
+      CALL iom_get( inum, jpdom_data, 'bottom level' , zbot    )  ! nb of ocean T-points
+      CALL iom_get( inum, jpdom_data, 'top    level' , ztop     )   ! nb of ocean T-points (ISF)
+      !
+      CALL iom_close( inum )
+      !
+      ! 2D ocean mask (=1 if at least one level of the water column is ocean, =0 otherwise)
+      WHERE( zbot(:,:) - ztop(:,:) + 1 > 0 )   ;   imask(:,:) = 1
+      ELSEWHERE                                ;   imask(:,:) = 0
+      END WHERE
 
       !  1. Dimension arrays for subdomains
@@ -320,5 +279,5 @@
          DO jj = 1+jprecj, ilj-jprecj
             DO  ji = 1+jpreci, ili-jpreci
-               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1) isurf = isurf+1
+               IF( imask(ji+iimppt(ii,ij)-1, jj+ijmppt(ii,ij)-1) == 1)   isurf = isurf+1
             END DO
          END DO
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90	(revision 6717)
@@ -13,7 +13,5 @@
    !!   obs_sor       : Sort the observation arrays
    !!---------------------------------------------------------------------
-   !! * Modules used
-   USE par_kind, ONLY : & ! Precision variables
-      & wp   
+   USE par_kind, ONLY : wp ! Precision variables
    USE in_out_manager     ! I/O manager
    USE obs_profiles_def   ! Definitions for storage arrays for profiles
@@ -24,16 +22,12 @@
    USE obs_inter_sup      ! Interpolation support
    USE obs_oper           ! Observation operators
-   USE lib_mpp, ONLY : &
-      & ctl_warn, ctl_stop
+   USE lib_mpp, ONLY :   ctl_warn, ctl_stop
 
    IMPLICIT NONE
-
-   !! * Routine accessibility
    PRIVATE
 
-   PUBLIC &
-      & obs_pre_prof, &    ! First level check and screening of profile obs
-      & obs_pre_surf, &    ! First level check and screening of surface obs
-      & calc_month_len     ! Calculate the number of days in the months of a year
+   PUBLIC   obs_pre_prof     ! First level check and screening of profile obs
+   PUBLIC   obs_pre_surf     ! First level check and screening of surface obs
+   PUBLIC   calc_month_len   ! Calculate the number of days in the months of a year
 
    !!----------------------------------------------------------------------
@@ -63,17 +57,11 @@
       !!        !  2015-02  (M. Martin) Combined routine for surface types.
       !!----------------------------------------------------------------------
-      !! * Modules used
-      USE domstp              ! Domain: set the time-step
       USE par_oce             ! Ocean parameters
-      USE dom_oce, ONLY : &   ! Geographical information
-         & glamt,   &
-         & gphit,   &
-         & tmask,   &
-         & nproc
+      USE dom_oce, ONLY       :   glamt, gphit, tmask, nproc   ! Geographical information
       !! * Arguments
       TYPE(obs_surf), INTENT(INOUT) :: surfdata    ! Full set of surface data
       TYPE(obs_surf), INTENT(INOUT) :: surfdataqc   ! Subset of surface data not failing screening
       LOGICAL, INTENT(IN) :: ld_nea         ! Switch for rejecting observation near land
-      !! * Local declarations
+      !
       INTEGER :: iyea0        ! Initial date
       INTEGER :: imon0        !  - (year, month, day, hour, minute)
@@ -94,12 +82,12 @@
       INTEGER :: inlasobsmpp    !  - close to land
       INTEGER :: igrdobsmpp     !  - fail the grid search
-      LOGICAL, DIMENSION(:), ALLOCATABLE :: & 
-         & llvalid            ! SLA data selection
+      LOGICAL, DIMENSION(:), ALLOCATABLE ::   llvalid            ! SLA data selection
       INTEGER :: jobs         ! Obs. loop variable
       INTEGER :: jstp         ! Time loop variable
       INTEGER :: inrc         ! Time index variable
-
-      IF(lwp) WRITE(numout,*)'obs_pre_surf : Preparing the surface observations...'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*) 'obs_pre_surf : Preparing the surface observations...'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
       
       ! Initial date initialization (year, month, day, hour, minute)
@@ -253,10 +241,6 @@
       !!
       !!----------------------------------------------------------------------
-      !! * Modules used
-      USE domstp              ! Domain: set the time-step
-      USE par_oce             ! Ocean parameters
-      USE dom_oce, ONLY : &   ! Geographical information
-         & gdept_1d,             &
-         & nproc
+      USE par_oce           ! Ocean parameters
+      USE dom_oce, ONLY :   gdept_1d, nproc   ! Geographical information
 
       !! * Arguments
@@ -314,4 +298,5 @@
       INTEGER :: jstp         ! Time loop variable
       INTEGER :: inrc         ! Time index variable
+      !!----------------------------------------------------------------------
 
       IF(lwp) WRITE(numout,*)'obs_pre_prof: Preparing the profile data...'
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 6717)
@@ -43,9 +43,9 @@
    USE sbcisf         ! surface boundary condition: ice shelf
    USE sbcfwb         ! surface boundary condition: freshwater budget
-   USE closea         ! closed sea
    USE icbstp         ! Icebergs
    USE traqsr         ! active tracers: light penetration
    USE sbcwave        ! Wave module
    USE bdy_par        ! Require lk_bdy
+   USE usrdef_closea  ! closed sea
    !
    USE prtctl         ! Print control                    (prt_ctl routine)
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 6717)
@@ -20,6 +20,6 @@
    USE sbc_oce        ! surface boundary condition variables
    USE sbcisf         ! PM we could remove it I think
-   USE closea         ! closed seas
    USE eosbn2         ! Equation Of State
+   USE usrdef_closea  ! closed seas
    !
    USE in_out_manager ! I/O manager
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 6717)
@@ -581,9 +581,4 @@
             ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
             !
-         CASE ( 4 )                          ! ORCA_R4
-            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
-            ii0 =  70   ;   ii1 =  71
-            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
-            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
          END SELECT
          !
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_closea.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_closea.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_closea.F90	(revision 6717)
@@ -0,0 +1,418 @@
+MODULE usrdef_closea
+   !!======================================================================
+   !!                   ***  MODULE  usrdef_closea  ***
+   !! User define : specific treatments associated with closed seas
+   !!======================================================================
+   !! History :   8.2  !  2000-05  (O. Marti)  Original code
+   !!   NEMO      1.0  !  2002-06  (E. Durand, G. Madec)  F90
+   !!             3.0  !  2006-07  (G. Madec)  add clo_rnf, clo_ups, clo_bat
+   !!             3.4  !  2014-12  (P.G. Fogli) sbc_clo bug fix & mpp reproducibility
+   !!             4.0  !  2016-06  (G. Madec)  move to usrdef_closea, remove clo_ups
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dom_clo    : modification of the ocean domain for closed seas cases
+   !!   sbc_clo    : Special handling of closed seas
+   !!   clo_rnf    : set close sea outflows as river mouths (see sbcrnf)
+   !!   clo_bat    : set to zero a field over closed sea (see domzrg)
+   !!----------------------------------------------------------------------
+   USE oce             ! dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE sbc_oce         ! ocean surface boundary conditions
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_fortran,    ONLY: glob_sum, DDPDD
+   USE lbclnk          ! lateral boundary condition - MPP exchanges
+   USE lib_mpp         ! MPP library
+   USE timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dom_clo      ! routine called by domain module
+   PUBLIC sbc_clo      ! routine called by step module
+   PUBLIC clo_rnf      ! routine called by sbcrnf module
+   PUBLIC clo_bat      ! routine called in domzgr module
+
+   INTEGER, PUBLIC, PARAMETER          ::   jpncs   = 4      !: number of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncstt            !: Type of closed sea
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi1, ncsj1     !: south-west closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsi2, ncsj2     !: north-east closed sea limits (i,j)
+   INTEGER, PUBLIC, DIMENSION(jpncs)   ::   ncsnr            !: number of point where run-off pours
+   INTEGER, PUBLIC, DIMENSION(jpncs,4) ::   ncsir, ncsjr     !: Location of runoff
+
+   REAL(wp), DIMENSION (jpncs+1)       ::   surf             ! closed sea surface
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dom_clo( cd_cfg, kcfg )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_clo  ***
+      !!        
+      !! ** Purpose :   Closed sea domain initialization
+      !!
+      !! ** Method  :   if a closed sea is located only in a model grid point
+      !!                just the thermodynamic processes are applied.
+      !!
+      !! ** Action  :   ncsi1(), ncsj1() : south-west closed sea limits (i,j)
+      !!                ncsi2(), ncsj2() : north-east Closed sea limits (i,j)
+      !!                ncsir(), ncsjr() : Location of runoff
+      !!                ncsnr            : number of point where run-off pours
+      !!                ncstt            : Type of closed sea
+      !!                                   =0 spread over the world ocean
+      !!                                   =2 put at location runoff
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_cfg   ! configuration name
+      INTEGER                   , INTENT(in   ) ::   kcfg     ! configuration identifier 
+      !
+      INTEGER ::   jc      ! dummy loop indices
+      INTEGER ::   isrow   ! local index
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*)'dom_clo : closed seas '
+      IF(lwp) WRITE(numout,*)'~~~~~~~'
+      !
+      ! initial values
+      ncsnr(:) = 1  ;  ncsi1(:) = 1  ;  ncsi2(:) = 1  ;  ncsir(:,:) = 1
+      ncstt(:) = 0  ;  ncsj1(:) = 1  ;  ncsj2(:) = 1  ;  ncsjr(:,:) = 1
+      !
+      ! set the closed seas (in data domain indices)
+      ! -------------------
+      !
+      IF( cd_cfg == "orca" ) THEN      !==  ORCA configuration  ==!
+         !
+         SELECT CASE ( kcfg )
+         !                                           ! =======================
+         CASE ( 1 )                                  !  ORCA_R1 configuration
+            !                                        ! =======================
+            IF(lwp) WRITE(numout,*)'   ORCA_R1 closed seas :  only the Caspian Sea'
+            ! This dirty section will be suppressed by simplification process:
+            ! all this will come back in input files
+            ! Currently these hard-wired indices relate to configuration with
+            ! extend grid (jpjglo=332)
+            isrow = 332 - jpjglo
+            !
+            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea  (spread over the globe)
+            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow
+            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow
+            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
+            !                                        
+            !                                        ! =======================
+         CASE ( 2 )                                  !  ORCA_R2 configuration
+            !                                        ! =======================
+            IF(lwp) WRITE(numout,*)'   ORCA_R2 closed seas and lakes : '
+            !                                            ! Caspian Sea
+            IF(lwp) WRITE(numout,*)'      Caspian Sea  '
+            ncsnr(1)   =   1  ;  ncstt(1)   =   0           ! spread over the globe
+            ncsi1(1)   =  11  ;  ncsj1(1)   = 103
+            ncsi2(1)   =  17  ;  ncsj2(1)   = 112
+            ncsir(1,1) =   1  ;  ncsjr(1,1) =   1 
+            !                                            ! Great North American Lakes
+            IF(lwp) WRITE(numout,*)'      Great North American Lakes  '
+            ncsnr(2)   =   1  ;  ncstt(2)   =   2           ! put at St Laurent mouth
+            ncsi1(2)   =  97  ;  ncsj1(2)   = 107
+            ncsi2(2)   = 103  ;  ncsj2(2)   = 111
+            ncsir(2,1) = 110  ;  ncsjr(2,1) = 111           
+            !                                            ! Black Sea (crossed by the cyclic boundary condition)
+            IF(lwp) WRITE(numout,*)'      Black Sea  '
+            ncsnr(3:4) =   4  ;  ncstt(3:4) =   2           ! put in Med Sea (north of Aegean Sea)
+            ncsir(3:4,1) = 171;  ncsjr(3:4,1) = 106         !
+            ncsir(3:4,2) = 170;  ncsjr(3:4,2) = 106 
+            ncsir(3:4,3) = 171;  ncsjr(3:4,3) = 105 
+            ncsir(3:4,4) = 170;  ncsjr(3:4,4) = 105 
+            ncsi1(3)   = 174  ;  ncsj1(3)   = 107           ! 1 : west part of the Black Sea      
+            ncsi2(3)   = 181  ;  ncsj2(3)   = 112           !            (ie west of the cyclic b.c.)
+            ncsi1(4)   =   2  ;  ncsj1(4)   = 107           ! 2 : east part of the Black Sea 
+            ncsi2(4)   =   6  ;  ncsj2(4)   = 112           !           (ie east of the cyclic b.c.)
+            !
+            !                                        ! =========================
+         CASE ( 025 )                                !  ORCA_R025 configuration
+            !                                        ! =========================
+            IF(lwp) WRITE(numout,*)'   ORCA_R025 closed seas : '
+            !                                            ! Caspian Sea
+            IF(lwp) WRITE(numout,*)'      Caspian Sea  '
+            ncsnr(1)   = 1    ; ncstt(1)   = 0               ! Caspian + Aral sea
+            ncsi1(1)   = 1330 ; ncsj1(1)   = 645
+            ncsi2(1)   = 1400 ; ncsj2(1)   = 795
+            ncsir(1,1) = 1    ; ncsjr(1,1) = 1
+            !                                        
+            IF(lwp) WRITE(numout,*)'      Azov Sea  '
+            ncsnr(2)   = 1    ; ncstt(2)   = 0               ! Azov Sea 
+            ncsi1(2)   = 1284 ; ncsj1(2)   = 722
+            ncsi2(2)   = 1304 ; ncsj2(2)   = 747
+            ncsir(2,1) = 1    ; ncsjr(2,1) = 1
+            !
+         END SELECT
+         !
+      ELSE                             !==  No closed sea in the configuration  ==!
+         !
+         IF(lwp) WRITE(numout,*)'   No closed seas or lakes in the configuration '
+         !
+      ENDIF
+
+      ! convert the position in local domain indices
+      ! --------------------------------------------
+      DO jc = 1, jpncs
+         ncsi1(jc)   = mi0( ncsi1(jc) )
+         ncsj1(jc)   = mj0( ncsj1(jc) )
+         !
+         ncsi2(jc)   = mi1( ncsi2(jc) )   
+         ncsj2(jc)   = mj1( ncsj2(jc) )  
+      END DO
+      !
+   END SUBROUTINE dom_clo
+
+
+   SUBROUTINE sbc_clo( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_clo  ***
+      !!                    
+      !! ** Purpose :   Special handling of closed seas
+      !!
+      !! ** Method  :   Water flux is forced to zero over closed sea
+      !!      Excess is shared between remaining ocean, or
+      !!      put as run-off in open ocean.
+      !!
+      !! ** Action  :   emp updated surface freshwater fluxes and associated heat content at kt
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean model time step
+      !
+      INTEGER             ::   ji, jj, jc, jn   ! dummy loop indices
+      REAL(wp), PARAMETER ::   rsmall = 1.e-20_wp    ! Closed sea correction epsilon
+      REAL(wp)            ::   zze2, ztmp, zcorr     ! 
+      REAL(wp)            ::   zcoef, zcoef1         ! 
+      COMPLEX(wp)         ::   ctmp 
+      REAL(wp), DIMENSION(jpncs) ::   zfwf   ! 1D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('sbc_clo')
+      !                                                   !------------------!
+      IF( kt == nit000 ) THEN                             !  Initialisation  !
+         !                                                !------------------!
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'sbc_clo : closed seas '
+         IF(lwp) WRITE(numout,*)'~~~~~~~'
+
+         surf(:) = 0._wp
+         !
+         surf(jpncs+1) = glob_sum( e1e2t(:,:) )   ! surface of the global ocean
+         !
+         !                                        ! surface of closed seas 
+         IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
+            DO jc = 1, jpncs
+               ctmp = CMPLX( 0.e0, 0.e0, wp )
+               DO jj = ncsj1(jc), ncsj2(jc)
+                  DO ji = ncsi1(jc), ncsi2(jc)
+                     ztmp = e1e2t(ji,jj) * tmask_i(ji,jj)
+                     CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+                  END DO 
+               END DO 
+               IF( lk_mpp )   CALL mpp_sum( ctmp )
+               surf(jc) = REAL(ctmp,wp)
+            END DO
+         ELSE                                          ! Standard calculation           
+            DO jc = 1, jpncs
+               DO jj = ncsj1(jc), ncsj2(jc)
+                  DO ji = ncsi1(jc), ncsi2(jc)
+                     surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj)      ! surface of closed seas
+                  END DO 
+               END DO 
+            END DO 
+            IF( lk_mpp )   CALL mpp_sum ( surf, jpncs )       ! mpp: sum over all the global domain
+         ENDIF
+
+         IF(lwp) WRITE(numout,*)'     Closed sea surfaces'
+         DO jc = 1, jpncs
+            IF(lwp)WRITE(numout,FMT='(1I3,4I4,5X,F16.2)') jc, ncsi1(jc), ncsi2(jc), ncsj1(jc), ncsj2(jc), surf(jc)
+         END DO
+
+         ! jpncs+1 : surface of sea, closed seas excluded
+         DO jc = 1, jpncs
+            surf(jpncs+1) = surf(jpncs+1) - surf(jc)
+         END DO           
+         !
+      ENDIF
+      !                                                   !--------------------!
+      !                                                   !  update emp        !
+      zfwf = 0.e0_wp                                      !--------------------!
+      IF( lk_mpp_rep ) THEN                         ! MPP reproductible calculation
+         DO jc = 1, jpncs
+            ctmp = CMPLX( 0.e0, 0.e0, wp )
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj)
+                  CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp )
+               END DO  
+            END DO 
+            IF( lk_mpp )   CALL mpp_sum( ctmp )
+            zfwf(jc) = REAL(ctmp,wp)
+         END DO
+      ELSE                                          ! Standard calculation           
+         DO jc = 1, jpncs
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 
+               END DO  
+            END DO 
+         END DO
+         IF( lk_mpp )   CALL mpp_sum ( zfwf(:) , jpncs )       ! mpp: sum over all the global domain
+      ENDIF
+
+      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      ! Black Sea case for ORCA_R2 configuration
+         zze2    = ( zfwf(3) + zfwf(4) ) * 0.5_wp
+         zfwf(3) = zze2
+         zfwf(4) = zze2
+      ENDIF
+
+      zcorr = 0._wp
+
+      DO jc = 1, jpncs
+         !
+         ! The following if avoids the redistribution of the round off
+         IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN
+            !
+            IF( ncstt(jc) == 0 ) THEN           ! water/evap excess is shared by all open ocean
+               zcoef    = zfwf(jc) / surf(jpncs+1)
+               zcoef1   = rcp * zcoef
+               emp(:,:) = emp(:,:) + zcoef
+               qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+               ! accumulate closed seas correction
+               zcorr    = zcorr    + zcoef
+               !
+            ELSEIF( ncstt(jc) == 1 ) THEN       ! Excess water in open sea, at outflow location, excess evap shared
+               IF ( zfwf(jc) <= 0.e0_wp ) THEN 
+                   DO jn = 1, ncsnr(jc)
+                     ji = mi0(ncsir(jc,jn))
+                     jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
+                     IF (      ji > 1 .AND. ji < jpi   &
+                         .AND. jj > 1 .AND. jj < jpj ) THEN 
+                         zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) )
+                         zcoef1     = rcp * zcoef
+                         emp(ji,jj) = emp(ji,jj) + zcoef
+                         qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
+                     ENDIF 
+                   END DO 
+               ELSE 
+                   zcoef    = zfwf(jc) / surf(jpncs+1)
+                   zcoef1   = rcp * zcoef
+                   emp(:,:) = emp(:,:) + zcoef
+                   qns(:,:) = qns(:,:) - zcoef1 * sst_m(:,:)
+                   ! accumulate closed seas correction
+                   zcorr    = zcorr    + zcoef
+               ENDIF
+            ELSEIF( ncstt(jc) == 2 ) THEN       ! Excess e-p-r (either sign) goes to open ocean, at outflow location
+               DO jn = 1, ncsnr(jc)
+                  ji = mi0(ncsir(jc,jn))
+                  jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean
+                  IF(      ji > 1 .AND. ji < jpi    &
+                     .AND. jj > 1 .AND. jj < jpj ) THEN 
+                     zcoef      = zfwf(jc) / ( REAL(ncsnr(jc)) *  e1e2t(ji,jj) )
+                     zcoef1     = rcp * zcoef
+                     emp(ji,jj) = emp(ji,jj) + zcoef
+                     qns(ji,jj) = qns(ji,jj) - zcoef1 * sst_m(ji,jj)
+                  ENDIF 
+               END DO 
+            ENDIF 
+            !
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  zcoef      = zfwf(jc) / surf(jc)
+                  zcoef1     = rcp * zcoef
+                  emp(ji,jj) = emp(ji,jj) - zcoef
+                  qns(ji,jj) = qns(ji,jj) + zcoef1 * sst_m(ji,jj)
+               END DO  
+            END DO 
+            !
+         END IF
+      END DO 
+
+      IF ( ABS(zcorr) > rsmall ) THEN      ! remove the global correction from the closed seas
+         DO jc = 1, jpncs                  ! only if it is large enough
+            DO jj = ncsj1(jc), ncsj2(jc)
+               DO ji = ncsi1(jc), ncsi2(jc)
+                  emp(ji,jj) = emp(ji,jj) - zcorr
+                  qns(ji,jj) = qns(ji,jj) + rcp * zcorr * sst_m(ji,jj)
+               END DO  
+             END DO 
+          END DO
+      ENDIF
+      !
+      emp (:,:) = emp (:,:) * tmask(:,:,1)
+      !
+      CALL lbc_lnk( emp , 'T', 1._wp )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('sbc_clo')
+      !
+   END SUBROUTINE sbc_clo
+
+
+   SUBROUTINE clo_rnf( p_rnfmsk )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!                    
+      !! ** Purpose :   allow the treatment of closed sea outflow grid-points
+      !!                to be the same as river mouth grid-points
+      !!
+      !! ** Method  :   set to 1 the runoff mask (mskrnf, see sbcrnf module)
+      !!                at the closed sea outflow grid-point.
+      !!
+      !! ** Action  :   update (p_)mskrnf (set 1 at closed sea outflow)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   p_rnfmsk   ! river runoff mask (rnfmsk array)
+      !
+      INTEGER  ::   jc, jn, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         IF( ncstt(jc) >= 1 ) THEN            ! runoff mask set to 1 at closed sea outflows
+             DO jn = 1, 4
+                DO jj =    mj0( ncsjr(jc,jn) ), mj1( ncsjr(jc,jn) )
+                   DO ji = mi0( ncsir(jc,jn) ), mi1( ncsir(jc,jn) )
+                      p_rnfmsk(ji,jj) = MAX( p_rnfmsk(ji,jj), 1.0_wp )
+                   END DO
+                END DO
+            END DO 
+         ENDIF 
+      END DO 
+      !
+   END SUBROUTINE clo_rnf
+   
+      
+   SUBROUTINE clo_bat( k_top, k_bot )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE clo_bat  ***
+      !!                    
+      !! ** Purpose :   suppress closed sea from the domain
+      !!
+      !! ** Method  :   set first and last ocean level to 0 over the closed seas.
+      !!
+      !! ** Action  :   set pbat=0 and kbat=0 over closed seas
+      !!----------------------------------------------------------------------
+      INTEGER, DIMENSION(:,:), INTENT(inout) ::   k_top, k_bot   ! ocean first and last level indices
+      !
+      INTEGER  ::   jc, ji, jj      ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jc = 1, jpncs
+         DO jj = ncsj1(jc), ncsj2(jc)
+            DO ji = ncsi1(jc), ncsi2(jc)
+               k_top(ji,jj) = 0   
+               k_bot(ji,jj) = 0   
+            END DO 
+         END DO 
+       END DO 
+       !
+   END SUBROUTINE clo_bat
+
+   !!======================================================================
+END MODULE usrdef_closea
+
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_fmask.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_fmask.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_fmask.F90	(revision 6717)
@@ -0,0 +1,146 @@
+MODULE usrdef_fmask
+   !!======================================================================
+   !!                       ***  MODULE usrdef_fmask   ***
+   !! User Defined : alteration of land/sea f-point mask in some straits
+   !!======================================================================
+   !! History :  4.0  ! 2016-06  (G. Madec, S. Flavoni)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_fmask  : alteration of f-point land/ocean mask in some straits
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   !
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! Massively Parallel Processing library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_fmask    ! routine called by dommsk.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , LODYC-IPSL  (2016)
+   !! $Id: $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_fmask( cd_cfg, kcfg, pfmsk )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !!
+      !! ** Purpose :   User defined alteration of the lateral boundary 
+      !!              condition on velocity.
+      !!
+      !! ** Method  :   Local change of the value of fmask at lateral ocean/land 
+      !!              boundary in straits in order to increase the viscous 
+      !!              boundary layer and thus reduce the transport through the 
+      !!              corresponding straits.
+      !!                Here only alterations in ORCA R2 and R1 cases
+      !!
+      !! ** Action :   fmask : land/ocean mask at f-point with increased value 
+      !!                       in some user defined straits
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_cfg   ! configuration name
+      INTEGER                   , INTENT(in   ) ::   kcfg     ! configuration identifier 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pfmsk    ! Ocean/Land f-point mask including lateral boundary cond.
+      !
+      INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
+      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
+      INTEGER  ::   isrow                    ! index for ORCA1 starting row
+      !!----------------------------------------------------------------------
+      !
+      IF( cd_cfg == "orca" ) THEN      !==  ORCA Configurations  ==!
+         !
+         SELECT CASE ( kcfg )
+         !
+         CASE( 2 )                           !  R2 case
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R2: increase lateral friction near the following straits:'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+            !
+            IF(lwp) WRITE(numout,*) '      Gibraltar '
+            ij0 = 101   ;   ij1 = 101           ! Gibraltar strait  : partial slip (pfmsk=0.5)
+            ii0 = 139   ;   ii1 = 140   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
+            ij0 = 102   ;   ij1 = 102
+            ii0 = 139   ;   ii1 = 140   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  0.5_wp
+            !
+            IF(lwp) WRITE(numout,*) '      Bab el Mandeb '
+            ij0 =  87   ;   ij1 =  88           ! Bab el Mandeb : partial slip (pfmsk=1)
+            ii0 = 160   ;   ii1 = 160   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
+            ij0 =  88   ;   ij1 =  88
+            ii0 = 159   ;   ii1 = 159   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) =  1._wp
+            !
+            ! We keep this as an example but it is instable in this case 
+            !IF(lwp) WRITE(numout,*) '      Danish straits '
+            !         ij0 = 115   ;   ij1 = 115 ! Danish straits  : strong slip (pfmsk > 2)
+            !         ii0 = 145   ;   ii1 = 146   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
+            !         ij0 = 116   ;   ij1 = 116
+            !         ii0 = 145   ;   ii1 = 146   ;   pfmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 4._wp
+            !
+         CASE( 1 )                           ! R1 case
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R1: increase lateral friction near the following straits:'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'   
+!!gm    ! This dirty section will be suppressed by simplification process:
+!!gm    ! all this will come back in input files
+!!gm    ! Currently these hard-wired indices relate to configuration with extend grid (jpjglo=332)
+            !
+            isrow = 332 - jpjglo
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
+            IF(lwp) WRITE(numout,*) '      Gibraltar '
+            ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait 
+            ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      Bhosporus '
+            ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait 
+            ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      Makassar (Top) '
+            ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top) 
+            ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      Lombok '
+            ii0 =  44           ;   ii1 =  44        ! Lombok Strait 
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      Ombai '
+            ii0 =  53           ;   ii1 =  53        ! Ombai Strait 
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      Timor Passage '
+            ii0 =  56           ;   ii1 =  56        ! Timor Passage 
+            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      West Halmahera '
+            ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait 
+            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+            !
+            IF(lwp) WRITE(numout,*) '      East Halmahera '
+            ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait 
+            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   pfmsk( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp  
+            !
+         CASE DEFAULT
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'usr_def_fmask : ORCA_R', kcfg,' : NO alteration of fmask in specific straits '
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'   
+         END SELECT
+      ELSE
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'usr_def_fmask : NO alteration of fmask in specific straits '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
+      CALL lbc_lnk( pfmsk, 'F', 1._wp )      ! Lateral boundary conditions on fmask
+      !
+   END SUBROUTINE usr_def_fmask
+   
+   !!======================================================================
+END MODULE usrdef_fmask
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_hgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_hgr.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_hgr.F90	(revision 6717)
@@ -0,0 +1,160 @@
+MODULE usrdef_hgr
+   !!==============================================================================
+   !!                       ***  MODULE usrdef_hgr   ***
+   !! User defined module: used like example to define domain, init, sbc, ...
+   !!==============================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni) 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE usrdef_nam     !
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_hgr   ! called in domhgr.F90  module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
+      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
+      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain 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 is used in strait(s))
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE usr_def_hgr  ***
+      !!
+      !! ** Purpose :   user defined mesh and Coriolis parameter
+      !!
+      !! ** Method  :   set all intent(out) argument to a proper value
+      !!
+      !!                Here GYRE configuration :
+      !!          Rectangular mid-latitude domain 
+      !!          - with axes rotated by 45 degrees
+      !!          - a constant horizontal resolution of 106 km 
+      !!          - on a beta-plane
+      !!
+      !! ** 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)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
+      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
+      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 computed here, =0 otherwise 
+      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
+      !
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
+      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      !   -      -
+      !!-------------------------------------------------------------------------------
+      !
+      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,*) 'usr_def_hgr : GYRE configuration (beta-plane with rotated regular grid-spacing)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+      !
+      !                       !==  grid point position  ==!
+      !
+      zlam1 = -85._wp                           ! position of gridpoint (i,j) = (1,jpjglo)
+      zphi1 =  29._wp
+      !
+      ze1 = 106000._wp / REAL( nn_GYRE , wp )   ! gridspacing in meters
+      !
+      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp     ! angle: 45 degrees
+      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
+      ze1deg = ze1 / (ra * rad)
+      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
+      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
+      !   
+      IF( ln_bench ) THEN     ! benchmark: forced the resolution to be 106 km 
+         ze1 = 106000._wp     ! but keep (lat,lon) at the right nn_GYRE resolution
+         CALL ctl_warn( ' GYRE used as Benchmark: e1=e2=106km, no need to adjust rdt, ahm,aht ' )
+      ENDIF
+      IF( nprint==1 .AND. lwp )   THEN
+         WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
+         WRITE(numout,*) 'ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
+      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  
+            plamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+            pphit(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
+            plamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
+            pphiu(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
+            plamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            pphiv(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 
+            plamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
+            pphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
+         END DO
+      END DO
+      !
+      !                       !== 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
+      !
+      !
+      !                       !==  Coriolis parameter  ==!
+      kff = 1                                            !  indicate not to compute ff afterward
+      !
+      zbeta = 2. * omega * COS( rad * zphi1 ) / ra       ! beta at latitude zphi1
+      !SF we overwrite zphi0 (south point in latitude) used just above to define pphif (value of zphi0=15.5190567531966)
+      !SF for computation of Coriolis we keep the parameter of Hazeleger, W., and S. S. Drijfhout, JPO 1998.
+      zphi0 = 15._wp                                     !  latitude of the most southern grid point  
+      zf0   = 2. * omega * SIN( rad * zphi0 )            !  compute f0 1st point south
+      !
+      pff_f(:,:) = ( zf0 + zbeta * ABS( pphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
+      pff_t(:,:) = ( zf0 + zbeta * ABS( pphit(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
+      !
+      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
+
+   !!======================================================================
+END MODULE usrdef_hgr
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_istate.F90	(revision 6717)
@@ -0,0 +1,82 @@
+MODULE usrdef_istate
+   !!==============================================================================
+   !!                       ***  MODULE usrdef_istate   ***
+   !! User defined : set the initial state of a user configuration
+   !!==============================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni) Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!  usr_def_istate : initial state in Temperature and salinity
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_istate   ! called in istate.F90  module
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+  
+   SUBROUTINE usr_def_istate( pdept, ptmask, pts, pu, pv, pssh )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE usr_def_istate  ***
+      !! 
+      !! ** Purpose :   Initialization of the dynamics and tracers
+      !!                Here GYRE configuration example : (double gyre with rotated domain)
+      !!
+      !! ** Method  : - set temprature field
+      !!              - set salinity   field
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pdept   ! depth of t-point               [m]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   ptmask  ! t-point ocean mask             [m]
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts     ! T & S fields      [Celcius ; g/kg]
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pu      ! i-component of the velocity  [m/s] 
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   pv      ! j-component of the velocity  [m/s] 
+      REAL(wp), DIMENSION(jpi,jpj)         , INTENT(  out) ::   pssh    ! sea-surface height
+      !
+      INTEGER :: ji, jj, jk  ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_istate : analytical definition of initial state '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   Ocean at rest, with an horizontally uniform T and S profiles'
+      !
+      pu  (:,:,:) = 0._wp        ! ocean at rest
+      pu  (:,:,:) = 0._wp
+      pssh(:,:)   = 0._wp
+      !
+      DO jk = 1, jpk             ! horizontally uniform T & S profiles
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   &
+                    &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             &
+                    &           + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) )            &
+                    &           - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.)                        &
+                    &           + 7.  * (1500. - pdept(ji,jj,jk) ) / 1500.)                        &
+                    &           * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * ptmask(ji,jj,jk)
+
+               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) )  &
+                    &         * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2                  &
+                    &         + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000.                 &
+                    &         - 1.62 * TANH( (pdept(ji,jj,jk) - 60.  ) / 650. )                    &
+                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 35.  ) / 100. )                    &
+                    &         + 0.2  * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) )                  &
+                    &         * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * ptmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !   
+   END SUBROUTINE usr_def_istate
+
+   !!======================================================================
+END MODULE usrdef_istate
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_nam.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_nam.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_nam.F90	(revision 6717)
@@ -0,0 +1,95 @@
+MODULE usrdef_nam
+   !!==============================================================================
+   !!                       ***  MODULE usrdef_nam   ***
+   !! User defined : set the domain characteristics
+   !!
+   !!                ===     Here  GYRE configuration      ===
+   !!
+   !!==============================================================================
+   !! History :  NEMO ! 2016-03  (S. Flavoni, G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_nam   : read user defined namelist and set global domain size
+   !!   usr_def_hgr   : initialize the horizontal mesh 
+   !!----------------------------------------------------------------------
+   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain
+   USE par_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+   
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_nam   ! called in nemogcm.F90 module
+
+   !                              !!* namusr_def namelist *!!
+   LOGICAL, PUBLIC ::   ln_bench   ! =T benchmark test with gyre: the gridsize is constant (no need to adjust timestep or viscosity)
+   INTEGER, PUBLIC ::   nn_GYRE    ! 1/nn_GYRE = the resolution chosen in degrees and thus defining the horizontal domain size
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_nam( ldtxt, ldnam, kpi, kpj, kpk, kperio )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_nam  ***
+      !!                    
+      !! ** Purpose :   read user defined namelist and define the domain size
+      !!
+      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
+      !!
+      !!                Here GYRE configuration
+      !!
+      !! ** input   : - namusr_def namelist found in namelist_cfg
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information
+      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
+      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
+      !
+      INTEGER ::   ios, ii   ! Local integer
+      !!
+      NAMELIST/namusr_def/ nn_GYRE, ln_bench, jpkglo
+      !!----------------------------------------------------------------------
+      !
+      ii = 1
+      !
+      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
+      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. )
+      !
+!!gm  This does not work...  I don't know how to write namusr_def in "output.namelist.dyn"
+      WRITE( ldnam(ii), namusr_def )
+!!gm
+      !
+      kpi = 30 * nn_GYRE + 2        ! Global Domain size
+      kpj = 20 * nn_GYRE + 2
+      kpk = jpkglo
+      !
+      !                             ! control print
+      WRITE(ldtxt(ii),*) '   '                                                                            ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'     ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                   ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '   Namelist namusr_def : GYRE case'                                             ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      GYRE used as Benchmark (=T)                      ln_bench  = ', ln_bench  ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      inverse resolution & implied domain size         nn_GYRE   = ', nn_GYRE   ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         jpiglo = 30*nn_GYRE+2                            jpiglo = ', kpi       ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '         jpjglo = 20*nn_GYRE+2                            jpjglo = ', kpj       ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '      number of model levels                           jpkglo    = ', kpk       ;   ii = ii + 1
+
+      !                             ! Set the lateral boundary condition of the global domain
+      kperio = 0                    ! GYRE configuration : closed domain
+      !
+      WRITE(ldtxt(ii),*) '   '                                                                        ;   ii = ii + 1
+      WRITE(ldtxt(ii),*) '   Lateral b.c. of the global domain set to closed     jperio = ', kperio   ;   ii = ii + 1
+      !
+   END SUBROUTINE usr_def_nam
+
+   !!======================================================================
+END MODULE usrdef_nam
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_sbc.F90	(revision 6717)
@@ -0,0 +1,233 @@
+MODULE usrdef_sbc
+   !!======================================================================
+   !!                       ***  MODULE  usrdef_sbc  ***
+   !! Ocean forcing:  user defined momentum, heat and freshwater forcings
+   !!
+   !!                ===     Here  GYRE configuration      ===
+   !!
+   !!=====================================================================
+   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_sbc    : user defined surface bounday conditions in GYRE case
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
+   !
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_fortran     !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_sbc    ! routine called in sbcmod module
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE usr_def_sbc( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE usr_def_sbc  ***
+      !!              
+      !! ** Purpose :   provide at each time-step the GYRE surface boundary
+      !!              condition, i.e. the momentum, heat and freshwater fluxes.
+      !!
+      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
+      !!                CAUTION : never mask the surface stress field !
+      !!
+      !! ** Action  : - set the ocean surface boundary condition, i.e.   
+      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
+      !!
+      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!
+      INTEGER  ::   ji, jj                 ! dummy loop indices
+      INTEGER  ::   zyear0                 ! initial year 
+      INTEGER  ::   zmonth0                ! initial month
+      INTEGER  ::   zday0                  ! initial day
+      INTEGER  ::   zday_year0             ! initial day since january 1st
+      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
+      REAL(wp) ::   ztime                  ! time in hour
+      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
+      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
+      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
+      REAL(wp) ::   ztaun                  ! intensity
+      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
+      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
+      REAL(wp) ::   zsumemp, zsurf
+      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
+      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
+      REAL(wp) ::   zyydd                 ! number of days in one year
+      !!---------------------------------------------------------------------
+      zyydd = REAL(nyear_len(1),wp)
+
+      ! ---------------------------- !
+      !  heat and freshwater fluxes  !
+      ! ---------------------------- !
+      !same temperature, E-P as in HAZELEGER 2000
+
+      zyear0     =   ndate0 / 10000                             ! initial year
+      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
+      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
+      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
+
+      ! current day (in hours) since january the 1st of the current year
+      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
+         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
+
+      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
+      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
+      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
+      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
+      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
+
+      ! amplitudes
+      zemp_S    = 0.7       ! intensity of COS in the South
+      zemp_N    = 0.8       ! intensity of COS in the North
+      zemp_sais = 0.1
+      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
+
+      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
+      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
+      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
+
+      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
+      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
+            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
+            ! 64.5 in summer, 42.5 in winter
+            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
+               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
+               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
+            ! 23.5 deg : tropics
+            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
+            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
+            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
+               emp  (ji,jj) =   zemp_S * zconv   &
+                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
+                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
+            ELSE
+               emp (ji,jj) =  - zemp_N * zconv   &
+                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
+                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
+            ENDIF
+         END DO
+      END DO
+
+!!gm  ===>>> The following lines are to be de-commented
+!!    Commented here as they change the GYRE results 
+!#if !defined key_mpp_rep
+!      CALL ctl_stop( 'usr_def_sbc: key_mpp_rep is required to run GYRE configuration')
+!#endif
+!!gm
+
+      zsumemp = GLOB_SUM( emp  (:,:)   ) 
+      zsurf   = GLOB_SUM( tmask(:,:,1) ) 
+      zsumemp = zsumemp / zsurf         ! Default GYRE configuration
+
+      ! freshwater (mass flux) and update of qns with heat content of emp
+      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
+      sfx (:,:) = 0.0_wp                                   ! no salt flux
+      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
+
+
+      ! ---------------------------- !
+      !       momentum fluxes        !
+      ! ---------------------------- !
+      ! same wind as in Wico
+      !test date0 : ndate0 = 010203
+      zyear0  =   ndate0 / 10000
+      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
+      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
+      !Calculates nday_year, day since january 1st
+      zday_year0 = (zmonth0-1)*30.+zday0
+
+      !accumulates days of previous months of this year
+      ! day (in hours) since january the 1st
+      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
+         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
+      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
+      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
+      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
+
+      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
+      ztau = 0.105 / SQRT( 2. )
+      ! seasonal oscillation intensity
+      ztau_sais = 0.015
+      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+           ! domain from 15deg to 50deg and 1/2 period along 14deg
+           ! so 5/4 of half period with seasonal cycle
+           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
+           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
+         END DO
+      END DO
+
+      ! module of wind stress and wind speed at T-point
+      zcoef = 1. / ( zrhoa * zcdrag ) 
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vect. opt.
+            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
+            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
+            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
+            taum(ji,jj) = zmod
+            wndm(ji,jj) = SQRT( zmod * zcoef )
+         END DO
+      END DO
+      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
+
+      ! ---------------------------------- !
+      !  control print at first time-step  !
+      ! ---------------------------------- !
+      IF( kt == nit000 .AND. lwp ) THEN 
+         WRITE(numout,*)
+         WRITE(numout,*)'usr_def_sbc : analytical surface fluxes for GYRE configuration'               
+         WRITE(numout,*)'~~~~~~~~~~~ ' 
+         WRITE(numout,*)'           nyear      = ', nyear
+         WRITE(numout,*)'           nmonth     = ', nmonth
+         WRITE(numout,*)'           nday       = ', nday
+         WRITE(numout,*)'           nday_year  = ', nday_year
+         WRITE(numout,*)'           ztime      = ', ztime
+         WRITE(numout,*)'           ztimemax   = ', ztimemax
+         WRITE(numout,*)'           ztimemin   = ', ztimemin
+         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
+         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
+         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
+         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
+         WRITE(numout,*)'           zyear0     = ', zyear0
+         WRITE(numout,*)'           zmonth0    = ', zmonth0
+         WRITE(numout,*)'           zday0      = ', zday0
+         WRITE(numout,*)'           zday_year0 = ', zday_year0
+         WRITE(numout,*)'           zyydd      = ', zyydd
+         WRITE(numout,*)'           zemp_S     = ', zemp_S
+         WRITE(numout,*)'           zemp_N     = ', zemp_N
+         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
+         WRITE(numout,*)'           zTstar     = ', zTstar
+         WRITE(numout,*)'           zsumemp    = ', zsumemp
+         WRITE(numout,*)'           zsurf      = ', zsurf
+         WRITE(numout,*)'           ztrp       = ', ztrp
+         WRITE(numout,*)'           zconv      = ', zconv
+         WRITE(numout,*)'           ndastp     = ', ndastp
+         WRITE(numout,*)'           adatrj     = ', adatrj
+      ENDIF
+      !
+   END SUBROUTINE usr_def_sbc
+
+   !!======================================================================
+END MODULE usrdef_sbc
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90	(revision 6717)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/USR/usrdef_zgr.F90	(revision 6717)
@@ -0,0 +1,270 @@
+MODULE usrdef_zgr
+   !!==============================================================================
+   !!                       ***  MODULE usrdef_zgr  ***
+   !! Ocean domain : user defined vertical coordinate system 
+   !!
+   !!                       ===      GYRE case      ===
+   !!
+   !!==============================================================================
+   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   usr_def_zgr      : user defined vertical coordinate system
+   !!       zgr_z        : reference 1D z-coordinate 
+   !!       zgr_top_bot  : ocean top and bottom level indices
+   !!       zgr_zco      : 3D verticl coordinate in pure z-coordinate case
+   !!---------------------------------------------------------------------
+   USE oce               ! ocean variables
+   USE dom_oce           ! ocean domain
+   !
+   USE in_out_manager    ! I/O manager
+   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp           ! distributed memory computing library
+   USE wrk_nemo          ! Memory allocation
+   USE timing            ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   usr_def_zgr        ! called by domzgr.F90
+
+  !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id: domzgr.F90 6624 2016-05-26 08:59:48Z gm $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS             
+
+   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
+      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
+      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
+      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
+      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
+      &                    k_top  , k_bot    )                             ! top & bottom ocean level
+      !!---------------------------------------------------------------------
+      !!              ***  ROUTINE usr_def_zgr  ***
+      !!
+      !! ** Purpose :   User defined the vertical coordinates
+      !!
+      !!----------------------------------------------------------------------
+      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
+      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors 
+      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
+      !
+      INTEGER  ::   inum   ! local logical unit
+      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'usr_def_zgr : GYRE configuration (z-coordinate closed flat box ocean without cavities)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+      !
+      !
+      ! type of vertical coordinate
+      ! ---------------------------
+      ld_zco    = .TRUE.         ! GYRE case:  z-coordinate & no ocean cavities
+      ld_zps    = .FALSE.
+      ld_sco    = .FALSE.
+      ld_isfcav = .FALSE.
+      !
+      !
+      ! Build the vertical coordinate system
+      ! ------------------------------------
+      CALL zgr_z  ( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )  ! Reference z-coordinate system
+      !
+      CALL zgr_top_bot( k_top , k_bot )                      ! top and bottom ocean level indices
+      !
+      !                                                      ! z-coordinate (3D arrays) from the 1D z-coord.
+      CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate
+         &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth
+         &          pe3t    , pe3u    , pe3v   , pe3f   ,   &   !       vertical scale factors
+         &          pe3w    , pe3uw   , pe3vw             )     !           -      -      -
+      !
+   END SUBROUTINE usr_def_zgr
+
+
+   SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zgr_z  ***
+      !!
+      !! ** Purpose :   set the depth of model levels and the resulting 
+      !!      vertical scale factors.
+      !!
+      !! ** Method  :   z-coordinate system (use in all type of coordinate)
+      !!      The depth of model levels is defined from an analytical
+      !!      function the derivative of which gives the scale factors.
+      !!      both depth and scale factors only depend on k (1d arrays).
+      !!              w-level: pdepw_1d  = pdep(k)
+      !!                       pe3w_1d(k) = dk(pdep)(k)     = e3(k)
+      !!              t-level: pdept_1d  = pdep(k+0.5)
+      !!                       pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5)
+      !!
+      !!      Here the Madec & Imbard (1996) function is used
+      !!
+      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
+      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
+      !!
+      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
+      !!             Madec and Imbard, 1996, Clim. Dyn.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
+      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
+      !
+      INTEGER  ::   jk       ! dummy loop indices
+      REAL(wp) ::   zt, zw   ! local scalars
+      REAL(wp) ::   zsur, za0, za1, zkth, zacr   ! Values for the Madec & Imbard (1996) function  
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
+      !
+      ! Set variables from parameters
+      ! ------------------------------
+      zsur = -2033.194295283385_wp       
+      za0  =   155.8325369664153_wp 
+      za1  =   146.3615918601890_wp
+      zkth =    17.28520372419791_wp   
+      zacr =     5.0_wp       
+      !
+      IF(lwp) THEN                         ! Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
+         WRITE(numout,*) '    ~~~~~~~'
+         WRITE(numout,*) '       GYRE case : MI96 function with the following coefficients :'
+         WRITE(numout,*) '                 zsur = ', zsur
+         WRITE(numout,*) '                 za0  = ', za0
+         WRITE(numout,*) '                 za1  = ', za1
+         WRITE(numout,*) '                 zkth = ', zkth
+         WRITE(numout,*) '                 zacr = ', zacr
+      ENDIF
+
+
+      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
+      ! ----------------------
+      DO jk = 1, jpk
+         zw = REAL( jk , wp )
+         zt = REAL( jk , wp ) + 0.5_wp
+         pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  )
+         pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  )
+         pe3w_1d (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
+         pe3t_1d (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
+      END DO
+      pdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
+
+
+!!gm   This should become the reference !
+!      IF ( ln_isfcav ) THEN
+! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
+! define pe3t_0 and pe3w_0 as the differences between pdept and pdepw respectively
+!         DO jk = 1, jpkm1
+!            pe3t_1d(jk) = pdepw_1d(jk+1)-pdepw_1d(jk) 
+!         END DO
+!         pe3t_1d(jpk) = pe3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
+!
+!         DO jk = 2, jpk
+!            pe3w_1d(jk) = pdept_1d(jk) - pdept_1d(jk-1) 
+!         END DO
+!         pe3w_1d(1  ) = 2._wp * (pdept_1d(1) - pdepw_1d(1)) 
+!      END IF
+!!gm end
+
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
+         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
+         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
+      ENDIF
+      DO jk = 1, jpk                      ! control positivity
+         IF( pe3w_1d (jk) <= 0._wp .OR. pe3t_1d (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
+         IF( pdepw_1d(jk) <  0._wp .OR. pdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
+      !
+   END SUBROUTINE zgr_z
+
+
+   SUBROUTINE zgr_top_bot( k_top , k_bot )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_top_bot  ***
+      !!
+      !! ** Purpose :   set the top and bottom ocean levels
+      !!
+      !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities
+      !!                   k_top = 1     except along north, south, east and west boundaries
+      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
+      !!
+      !! ** Action  : - k_top : first ocean level index
+      !!              - k_bot : last  ocean level index
+      !!----------------------------------------------------------------------
+      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last ocean level
+      !
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1   ! local indices
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom ocean levels.'
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
+      !
+      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
+      !
+      CALL lbc_lnk( z2d, 'T', 1. )           ! applied the lateral boundary condition (here jperio=0 ==>> closed)
+      !
+      k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
+      !
+      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! =1     over the ocean point, =0 elsewhere
+      !
+   END SUBROUTINE zgr_top_bot
+   
+
+   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
+      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
+      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
+      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zgr_zco  ***
+      !!
+      !! ** Purpose :   define the reference z-coordinate system
+      !!
+      !! ** Method  :   set 3D coord. arrays to reference 1D array 
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
+      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
+      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
+      !
+      INTEGER  ::   jk
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
+      !
+      DO jk = 1, jpk
+         pdept(:,:,jk) = pdept_1d(jk)
+         pdepw(:,:,jk) = pdepw_1d(jk)
+         pe3t (:,:,jk) = pe3t_1d (jk)
+         pe3u (:,:,jk) = pe3t_1d (jk)
+         pe3v (:,:,jk) = pe3t_1d (jk)
+         pe3f (:,:,jk) = pe3t_1d (jk)
+         pe3w (:,:,jk) = pe3w_1d (jk)
+         pe3uw(:,:,jk) = pe3w_1d (jk)
+         pe3vw(:,:,jk) = pe3w_1d (jk)
+      END DO
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
+      !
+   END SUBROUTINE zgr_zco
+
+   !!======================================================================
+END MODULE usrdef_zgr
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 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 6717)
@@ -93,5 +93,5 @@
    USE dia25h         ! 25h mean output
    USE wet_dry        ! Wetting and drying setting   (wad_init routine)
-   USE usrdef         ! user defined configuration
+   USE usrdef_nam     ! user defined configuration
 
    IMPLICIT NONE
@@ -105,5 +105,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -300,7 +300,4 @@
       ENDIF
       jpk    = jpkglo
-      jpidta = jpiglo         !!gm  jpidta, jpjdta : to be suppressed
-      jpjdta = jpjglo         !!gm   
-      jpkdta = jpkglo
       !
 #if defined key_agrif
@@ -310,6 +307,4 @@
          jpi     = ( jpiglo-2*jpreci + (jpni-1+0) ) / jpni + 2*jpreci
          jpj     = ( jpjglo-2*jprecj + (jpnj-1+0) ) / jpnj + 2*jprecj
-         jpidta  = jpiglo
-         jpjdta  = jpjglo
          nperio  = 0
          jperio  = 0
@@ -387,12 +382,14 @@
 #endif
       ENDIF
-               
-      jpk = jpkdta                                             ! third dim
-      
+
+!!gm ???    why here  it has already been done in line 301 !
+      jpk = jpkglo                                             ! third dim
+!!gm end
+
 #if defined key_agrif
       ! simple trick to use same vertical grid as parent but different number of levels: 
-      ! Save maximum number of levels in jpkdta, then define all vertical grids with this number.
+      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
       ! Suppress once vertical online interpolation is ok
-      IF(.NOT.Agrif_Root())   jpkdta = Agrif_Parent( jpkdta )
+      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo )
 #endif
       jpim1 = jpi-1                                            ! inner domain indices
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 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 6717)
@@ -31,11 +31,4 @@
    !! Domain Matrix size 
    !!---------------------------------------------------------------------
-
-!!gm  TO BE SUPRESSED
-   ! data size                        !!! * size of all input files *
-   INTEGER       ::   jpidta           !: 1st lateral dimension 
-   INTEGER       ::   jpjdta           !: 2nd    "         "   
-   INTEGER       ::   jpkdta           !: number of levels     
-!!gm  END
 
    ! global domain size               !!! * total computational domain *
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6717)
@@ -237,5 +237,4 @@
       IF( nn_diacfl == 1 )   CALL dia_cfl( kstp )         ! Courant number diagnostics
       IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
-      IF(.NOT.ln_cpl )   CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
       IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
       IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
Index: /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6716)
+++ /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 6717)
@@ -84,5 +84,4 @@
    USE diaar5           ! AR5 diagnosics                   (dia_ar5 routine)
    USE diahth           ! thermocline depth                (dia_hth routine)
-   USE diafwb           ! freshwater budget                (dia_fwb routine)
    USE diahsb           ! heat, salt and volume budgets    (dia_hsb routine)
    USE diaharm
Index: anches/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 6716)
+++ 	(revision )
@@ -1,265 +1,0 @@
-MODULE usrdef
-   !!==============================================================================
-   !!                       ***  MODULE usrdef   ***
-   !! User defined module: used like example to define domain, init, sbc, ...
-   !!==============================================================================
-   !! History :  NEMO ! 2016-03  (S. Flavoni) 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   usr_def_nam       : read user defined namelist and set global domain size
-   !!   usr_def_hgr       : initialize the horizontal mesh 
-   !!   usr_def_ini       : initial state 
-   !!----------------------------------------------------------------------
-   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain
-   USE par_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE timing         ! Timing
-   
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   usr_def_nam   ! called in nemogcm.F90 module
-   PUBLIC   usr_def_hgr   ! called in domhgr.F90  module
-   PUBLIC   usr_def_ini   ! called in istate.F90  module
-
-   !                      !!* namusr_def namelist *!!
-   LOGICAL ::   ln_bench   ! =T benchmark test with gyre: the gridsize is constant (no need to adjust timestep or viscosity)
-   INTEGER ::   nn_GYRE    ! 1/nn_GYRE = the resolution chosen in degrees and thus defining the horizontal domain size
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
-   !! $Id:$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE usr_def_nam( ldtxt, ldnam, kpi, kpj, kpk, kperio )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_nam  ***
-      !!                    
-      !! ** Purpose :   read user defined namelist and define the domain size
-      !!
-      !! ** Method  :   read in namusr_def containing all the user specific namelist parameter
-      !!
-      !!                Here GYRE configuration
-      !!
-      !! ** input   : - namusr_def namelist found in namelist_cfg
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt, ldnam    ! stored print information
-      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
-      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
-      !
-      INTEGER ::   ios, ii   ! Local integer
-      !!
-      NAMELIST/namusr_def/ nn_GYRE, ln_bench, jpkglo
-      !!----------------------------------------------------------------------
-      !
-      ii = 1
-      !
-      REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
-      READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', .TRUE. )
-      !
-!!gm  This does not work...  I don't know how to write namusr_def in "output.namelist.dyn"
-      WRITE( ldnam(ii), namusr_def )
-      !
-      kpi = 30 * nn_GYRE + 2        ! Global Domain size
-      kpj = 20 * nn_GYRE + 2
-      kpk = jpkglo
-      !
-      !                             ! control print
-      WRITE(ldtxt(ii),*) '   '                                                                            ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) 'usr_def_nam  : read the user defined namelist (namusr_def) in namelist_cfg'     ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '~~~~~~~~~~~ '                                                                   ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '   Namelist namusr_def : GYRE case'                                             ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '      GYRE used as Benchmark (=T)                      ln_bench  = ', ln_bench  ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '      inverse resolution & implied domain size         nn_GYRE   = ', nn_GYRE   ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '         jpiglo = 30*nn_GYRE+2                            jpiglo = ', kpi       ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '         jpjglo = 20*nn_GYRE+2                            jpjglo = ', kpj       ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '      number of model levels                           jpkglo    = ', kpk       ;   ii = ii + 1
-
-      !                             ! Set the lateral boundary condition of the global domain
-      kperio = 0                    ! GYRE configuration : closed domain
-      !
-      WRITE(ldtxt(ii),*) '   '                                                                        ;   ii = ii + 1
-      WRITE(ldtxt(ii),*) '   Lateral b.c. of the global domain set to closed     jperio = ', kperio   ;   ii = ii + 1
-      !
-   END SUBROUTINE usr_def_nam
-
-
-   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
-      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
-      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain 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 is used in strait(s))
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE usr_def_hgr  ***
-      !!
-      !! ** Purpose :   user defined mesh and Coriolis parameter
-      !!
-      !! ** Method  :   set all intent(out) argument to a proper value
-      !!
-      !!                Here GYRE configuration :
-      !!          Rectangular mid-latitude domain 
-      !!          - with axes rotated by 45 degrees
-      !!          - a constant horizontal resolution of 106 km 
-      !!          - on a beta-plane
-      !!
-      !! ** 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)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
-      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
-      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 computed here, =0 otherwise 
-      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
-      !
-      INTEGER  ::   ji, jj               ! dummy loop indices
-      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
-      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      !   -      -
-      !!-------------------------------------------------------------------------------
-      !
-      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,*) 'usr_def_hgr : GYRE configuration (beta-plane with rotated regular grid-spacing)'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      !
-      !
-      !                       !==  grid point position  ==!
-      !
-      zlam1 = -85._wp                           ! position of gridpoint (i,j) = (1,jpjglo)
-      zphi1 =  29._wp
-      !
-      ze1 = 106000._wp / REAL( nn_GYRE , wp )   ! gridspacing in meters
-      !
-      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp     ! angle: 45 degrees
-      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
-      ze1deg = ze1 / (ra * rad)
-      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
-      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
-      !   
-      IF( ln_bench ) THEN     ! benchmark: forced the resolution to be 106 km 
-         ze1 = 106000._wp     ! but keep (lat,lon) at the right nn_GYRE resolution
-         CALL ctl_warn( ' GYRE used as Benchmark: e1=e2=106km, no need to adjust rdt, ahm,aht ' )
-      ENDIF
-      IF( nprint==1 .AND. lwp )   THEN
-         WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
-         WRITE(numout,*) 'ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
-      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  
-            plamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pphit(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
-            plamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
-            pphiu(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
-            plamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pphiv(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 
-            plamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
-            pphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
-         END DO
-      END DO
-      !
-      !                       !== 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
-      !
-      !
-      !                       !==  Coriolis parameter  ==!
-      kff = 1                                            !  indicate not to compute ff afterward
-      !
-      zbeta = 2. * omega * COS( rad * zphi1 ) / ra       ! beta at latitude zphi1
-      !SF we overwrite zphi0 (south point in latitude) used just above to define pphif (value of zphi0=15.5190567531966)
-      !SF for computation of Coriolis we keep the parameter of Hazeleger, W., and S. S. Drijfhout, JPO 1998.
-      zphi0 = 15._wp                                     !  latitude of the most southern grid point  
-      zf0   = 2. * omega * SIN( rad * zphi0 )            !  compute f0 1st point south
-      !
-      pff_f(:,:) = ( zf0 + zbeta * ABS( pphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
-      pff_t(:,:) = ( zf0 + zbeta * ABS( pphit(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
-      !
-      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_ini( pts )
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE usr_def_ini  ***
-      !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers
-      !!                Here GYRE configuration example : (double gyre with rotated domain)
-      !!
-      !! ** Method  : - set temprature field
-      !!              - set salinity   field
-      !!----------------------------------------------------------------------
-      USE dom_oce, ONLY : gdept_0, tmask
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   pts   ! T & S data
-      !
-      INTEGER :: ji, jj, jk  ! dummy loop indices
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'usr_def_ini : analytical definition of initial state '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   T and S profiles deduced from LEVITUS '
-      !
-      DO jk = 1, jpk
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (gdept_0(ji,jj,jk) - 400) / 700 ) )   &
-                    &           * (-TANH( (500-gdept_0(ji,jj,jk)) / 150. ) + 1.) / 2.                &
-                    &           + ( 15. * ( 1. - TANH( (gdept_0(ji,jj,jk)-50.) / 1500.) )            &
-                    &           - 1.4 * TANH((gdept_0(ji,jj,jk)-100.) / 100.)                        &
-                    &           + 7.  * (1500. - gdept_0(ji,jj,jk) ) / 1500.)                        &
-                    &           * (-TANH( (gdept_0(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * tmask(ji,jj,jk)
-
-               pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (gdept_0(ji,jj,jk) - 305) / 460 ) )  &
-                    &         * (-TANH((500. - gdept_0(ji,jj,jk)) / 150.) + 1.) / 2                  &
-                    &         + ( 35.55 + 1.25 * (5000. - gdept_0(ji,jj,jk)) / 5000.                 &
-                    &         - 1.62 * TANH( (gdept_0(ji,jj,jk) - 60.  ) / 650. )                    &
-                    &         + 0.2  * TANH( (gdept_0(ji,jj,jk) - 35.  ) / 100. )                    &
-                    &         + 0.2  * TANH( (gdept_0(ji,jj,jk) - 1000.) / 5000.) )                  &
-                    &         * (-TANH( (gdept_0(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * tmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !   
-   END SUBROUTINE usr_def_ini
-
-   !!======================================================================
-END MODULE usrdef
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef_sbc.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef_sbc.F90	(revision 6716)
+++ 	(revision )
@@ -1,231 +1,0 @@
-MODULE usrdef_sbc
-   !!======================================================================
-   !!                       ***  MODULE  usrdef_sbc  ***
-   !! Ocean forcing:  analytical momentum, heat and freshwater forcings
-   !!=====================================================================
-   !! History :  3.0   ! 2006-06  (G. Madec)  Original code
-   !!            3.2   ! 2009-07  (G. Madec)  Style only
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   usr_def_sbc    : user defined surface bounday conditions
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE phycst          ! physical constants
-   !
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_fortran     !
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   usr_def_sbc    ! routine called in sbcmod module
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id:$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE usr_def_sbc( kt )
-      !!---------------------------------------------------------------------
-      !!                    ***  ROUTINE usr_def_sbc  ***
-      !!              
-      !! ** Purpose :   provide at each time-step the GYRE surface boundary
-      !!              condition, i.e. the momentum, heat and freshwater fluxes.
-      !!
-      !! ** Method  :   analytical seasonal cycle for GYRE configuration.
-      !!                CAUTION : never mask the surface stress field !
-      !!
-      !! ** Action  : - set the ocean surface boundary condition, i.e.   
-      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx
-      !!
-      !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000.
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      !!
-      INTEGER  ::   ji, jj                 ! dummy loop indices
-      INTEGER  ::   zyear0                 ! initial year 
-      INTEGER  ::   zmonth0                ! initial month
-      INTEGER  ::   zday0                  ! initial day
-      INTEGER  ::   zday_year0             ! initial day since january 1st
-      REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle
-      REAL(wp) ::   ztime                  ! time in hour
-      REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january
-      REAL(wp) ::   ztaun                  ! intensity
-      REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar
-      REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star
-      REAL(wp) ::   zsumemp, zsurf
-      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables
-      REAL(wp) ::   zyydd                 ! number of days in one year
-      !!---------------------------------------------------------------------
-      zyydd = REAL(nyear_len(1),wp)
-
-      ! ---------------------------- !
-      !  heat and freshwater fluxes  !
-      ! ---------------------------- !
-      !same temperature, E-P as in HAZELEGER 2000
-
-      zyear0     =   ndate0 / 10000                             ! initial year
-      zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month
-      zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30
-      zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360
-
-      ! current day (in hours) since january the 1st of the current year
-      ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours)
-         &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours)
-
-      ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours
-      ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours
-      ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours
-      ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours
-      !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours
-
-      ! amplitudes
-      zemp_S    = 0.7       ! intensity of COS in the South
-      zemp_N    = 0.8       ! intensity of COS in the North
-      zemp_sais = 0.1
-      zTstar    = 28.3      ! intemsity from 28.3 a -5 deg
-
-      ! 1/2 period between 21th June and 21th December and between 21th July and 21th January
-      zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi ) 
-      zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
-
-      ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K)
-      zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3
-            ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
-            ! 64.5 in summer, 42.5 in winter
-            t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                &
-               &                    * COS( rpi * (gphit(ji,jj) - 5.)               &
-               &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
-            ! 23.5 deg : tropics
-            qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
-            qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj)
-            IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg
-               emp  (ji,jj) =   zemp_S * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_S * zcos_sais1)
-            ELSE
-               emp (ji,jj) =  - zemp_N * zconv   &
-                  &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  &
-                  &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
-            ENDIF
-         END DO
-      END DO
-
-!!gm  ===>>> The following lines are to be de-commented
-!!    Commented here as they change the GYRE results 
-!#if !defined key_mpp_rep
-!      CALL ctl_stop( 'usr_def_sbc: key_mpp_rep is required to run GYRE configuration')
-!#endif
-!!gm
-
-      zsumemp = GLOB_SUM( emp  (:,:)   ) 
-      zsurf   = GLOB_SUM( tmask(:,:,1) ) 
-      zsumemp = zsumemp / zsurf         ! Default GYRE configuration
-
-      ! freshwater (mass flux) and update of qns with heat content of emp
-      emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average)
-      sfx (:,:) = 0.0_wp                                   ! no salt flux
-      qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST
-
-
-      ! ---------------------------- !
-      !       momentum fluxes        !
-      ! ---------------------------- !
-      ! same wind as in Wico
-      !test date0 : ndate0 = 010203
-      zyear0  =   ndate0 / 10000
-      zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100
-      zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100
-      !Calculates nday_year, day since january 1st
-      zday_year0 = (zmonth0-1)*30.+zday0
-
-      !accumulates days of previous months of this year
-      ! day (in hours) since january the 1st
-      ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour
-         &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years
-      ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours
-      ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours
-      !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours
-
-      ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle
-      ztau = 0.105 / SQRT( 2. )
-      ! seasonal oscillation intensity
-      ztau_sais = 0.015
-      ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-           ! domain from 15deg to 50deg and 1/2 period along 14deg
-           ! so 5/4 of half period with seasonal cycle
-           utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) )
-           vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) )
-         END DO
-      END DO
-
-      ! module of wind stress and wind speed at T-point
-      zcoef = 1. / ( zrhoa * zcdrag ) 
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            ztx = utau(ji-1,jj  ) + utau(ji,jj) 
-            zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
-            zmod = 0.5 * SQRT( ztx * ztx + zty * zty )
-            taum(ji,jj) = zmod
-            wndm(ji,jj) = SQRT( zmod * zcoef )
-         END DO
-      END DO
-      CALL lbc_lnk( taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( wndm(:,:), 'T', 1. )
-
-      ! ---------------------------------- !
-      !  control print at first time-step  !
-      ! ---------------------------------- !
-      IF( kt == nit000 .AND. lwp ) THEN 
-         WRITE(numout,*)
-         WRITE(numout,*)'usr_def_sbc : analytical surface fluxes for GYRE configuration'               
-         WRITE(numout,*)'~~~~~~~~~~~ ' 
-         WRITE(numout,*)'           nyear      = ', nyear
-         WRITE(numout,*)'           nmonth     = ', nmonth
-         WRITE(numout,*)'           nday       = ', nday
-         WRITE(numout,*)'           nday_year  = ', nday_year
-         WRITE(numout,*)'           ztime      = ', ztime
-         WRITE(numout,*)'           ztimemax   = ', ztimemax
-         WRITE(numout,*)'           ztimemin   = ', ztimemin
-         WRITE(numout,*)'           ztimemax1  = ', ztimemax1
-         WRITE(numout,*)'           ztimemin1  = ', ztimemin1
-         WRITE(numout,*)'           ztimemax2  = ', ztimemax2
-         WRITE(numout,*)'           ztimemin2  = ', ztimemin2
-         WRITE(numout,*)'           zyear0     = ', zyear0
-         WRITE(numout,*)'           zmonth0    = ', zmonth0
-         WRITE(numout,*)'           zday0      = ', zday0
-         WRITE(numout,*)'           zday_year0 = ', zday_year0
-         WRITE(numout,*)'           zyydd      = ', zyydd
-         WRITE(numout,*)'           zemp_S     = ', zemp_S
-         WRITE(numout,*)'           zemp_N     = ', zemp_N
-         WRITE(numout,*)'           zemp_sais  = ', zemp_sais
-         WRITE(numout,*)'           zTstar     = ', zTstar
-         WRITE(numout,*)'           zsumemp    = ', zsumemp
-         WRITE(numout,*)'           zsurf      = ', zsurf
-         WRITE(numout,*)'           ztrp       = ', ztrp
-         WRITE(numout,*)'           zconv      = ', zconv
-         WRITE(numout,*)'           ndastp     = ', ndastp
-         WRITE(numout,*)'           adatrj     = ', adatrj
-      ENDIF
-      !
-   END SUBROUTINE usr_def_sbc
-
-   !!======================================================================
-END MODULE usrdef_sbc
Index: anches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef_zgr.F90
===================================================================
--- /branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef_zgr.F90	(revision 6716)
+++ 	(revision )
@@ -1,275 +1,0 @@
-MODULE usrdef_zgr
-   !!==============================================================================
-   !!                       ***  MODULE usrdef_zgr  ***
-   !! Ocean domain : user defined vertical coordinate system 
-   !!
-   !!                       ===      GYRE case      ===
-   !!
-   !!==============================================================================
-   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   usr_def_zgr      : user defined vertical coordinate system
-   !!       zgr_z        : reference 1D z-coordinate 
-   !!       zgr_top_bot  : ocean top and bottom level indices
-   !!       zgr_zco      : 3D verticl coordinate in pure z-coordinate case
-   !!---------------------------------------------------------------------
-   USE oce               ! ocean variables
-   USE dom_oce           ! ocean domain
-   USE dommsk            ! domain: set the mask system
-   USE wet_dry           ! wetting and drying
-   USE closea            ! closed seas
-   USE c1d               ! 1D vertical configuration
-   !
-   USE in_out_manager    ! I/O manager
-   USE iom               ! I/O library
-   USE lbclnk            ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp           ! distributed memory computing library
-   USE wrk_nemo          ! Memory allocation
-   USE timing            ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   usr_def_zgr        ! called by domzgr.F90
-
-  !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id: domzgr.F90 6624 2016-05-26 08:59:48Z gm $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS             
-
-   SUBROUTINE usr_def_zgr( ld_zco  , ld_zps  , ld_sco  , ld_isfcav,    &   ! type of vertical coordinate
-      &                    pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d  ,    &   ! 1D reference vertical coordinate
-      &                    pdept , pdepw ,                             &   ! 3D t & w-points depth
-      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
-      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
-      &                    k_top  , k_bot    )                             ! top & bottom ocean level
-      !!---------------------------------------------------------------------
-      !!              ***  ROUTINE usr_def_zgr  ***
-      !!
-      !! ** Purpose :   User defined the vertical coordinates
-      !!
-      !!----------------------------------------------------------------------
-      LOGICAL                   , INTENT(out) ::   ld_zco, ld_zps, ld_sco      ! vertical coordinate flags
-      LOGICAL                   , INTENT(out) ::   ld_isfcav                   ! under iceshelf cavity flag
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth     [m]
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d           ! 1D grid-point depth     [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pdept, pdepw                ! grid-point depth        [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors  [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pe3w , pe3uw, pe3vw         ! i-scale factors 
-      INTEGER , DIMENSION(:,:)  , INTENT(out) ::   k_top, k_bot                ! first & last ocean level
-      !
-      INTEGER  ::   inum   ! local logical unit
-      REAL(WP) ::   z_zco, z_zps, z_sco, z_cav
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'usr_def_zgr : GYRE configuration (z-coordinate closed flat box ocean without cavities)'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      !
-      !
-      ! type of vertical coordinate
-      ! ---------------------------
-      ld_zco    = .TRUE.         ! GYRE case:  z-coordinate & no ocean cavities
-      ld_zps    = .FALSE.
-      ld_sco    = .FALSE.
-      ld_isfcav = .FALSE.
-      !
-      !
-      ! Build the vertical coordinate system
-      ! ------------------------------------
-      CALL zgr_z  ( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )  ! Reference z-coordinate system
-      !
-      CALL zgr_top_bot( k_top , k_bot )                      ! top and bottom ocean level indices
-      !
-      !                                                      ! z-coordinate (3D arrays) from the 1D z-coord.
-      CALL zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in  : 1D reference vertical coordinate
-         &          pdept   , pdepw   ,                     &   ! out : 3D t & w-points depth
-         &          pe3t    , pe3u    , pe3v   , pe3f   ,   &   !       vertical scale factors
-         &          pe3w    , pe3uw   , pe3vw             )     !           -      -      -
-      !
-   END SUBROUTINE usr_def_zgr
-
-
-   SUBROUTINE zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE zgr_z  ***
-      !!
-      !! ** Purpose :   set the depth of model levels and the resulting 
-      !!      vertical scale factors.
-      !!
-      !! ** Method  :   z-coordinate system (use in all type of coordinate)
-      !!      The depth of model levels is defined from an analytical
-      !!      function the derivative of which gives the scale factors.
-      !!      both depth and scale factors only depend on k (1d arrays).
-      !!              w-level: pdepw_1d  = pdep(k)
-      !!                       pe3w_1d(k) = dk(pdep)(k)     = e3(k)
-      !!              t-level: pdept_1d  = pdep(k+0.5)
-      !!                       pe3t_1d(k) = dk(pdep)(k+0.5) = e3(k+0.5)
-      !!
-      !!      Here the Madec & Imbard (1996) function is used
-      !!
-      !! ** Action  : - pdept_1d, pdepw_1d : depth of T- and W-point (m)
-      !!              - pe3t_1d , pe3w_1d  : scale factors at T- and W-levels (m)
-      !!
-      !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
-      !!             Madec and Imbard, 1996, Clim. Dyn.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pdept_1d, pdepw_1d   ! 1D grid-point depth        [m]
-      REAL(wp), DIMENSION(:)    , INTENT(out) ::   pe3t_1d , pe3w_1d    ! 1D vertical scale factors  [m]
-      !
-      INTEGER  ::   jk       ! dummy loop indices
-      REAL(wp) ::   zt, zw   ! local scalars
-      REAL(wp) ::   zsur, za0, za1, zkth, zacr   ! Values for the Madec & Imbard (1996) function  
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_z')
-      !
-      ! Set variables from parameters
-      ! ------------------------------
-      zsur = -2033.194295283385_wp       
-      za0  =   155.8325369664153_wp 
-      za1  =   146.3615918601890_wp
-      zkth =    17.28520372419791_wp   
-      zacr =     5.0_wp       
-      !
-      IF(lwp) THEN                         ! Parameter print
-         WRITE(numout,*)
-         WRITE(numout,*) '    zgr_z   : Reference vertical z-coordinates '
-         WRITE(numout,*) '    ~~~~~~~'
-         WRITE(numout,*) '       GYRE case : MI96 function with the following coefficients :'
-         WRITE(numout,*) '                 zsur = ', zsur
-         WRITE(numout,*) '                 za0  = ', za0
-         WRITE(numout,*) '                 za1  = ', za1
-         WRITE(numout,*) '                 zkth = ', zkth
-         WRITE(numout,*) '                 zacr = ', zacr
-      ENDIF
-
-
-      ! Reference z-coordinate (depth - scale factor at T- and W-points)   ! Madec & Imbard 1996 function
-      ! ----------------------
-      DO jk = 1, jpk
-         zw = REAL( jk , wp )
-         zt = REAL( jk , wp ) + 0.5_wp
-         pdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr *  LOG( COSH( (zw-zkth) / zacr ) )  )
-         pdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr *  LOG( COSH( (zt-zkth) / zacr ) )  )
-         pe3w_1d (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
-         pe3t_1d (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
-      END DO
-      pdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero
-
-
-!!gm   This should become the reference !
-!      IF ( ln_isfcav ) THEN
-! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
-! define pe3t_0 and pe3w_0 as the differences between pdept and pdepw respectively
-!         DO jk = 1, jpkm1
-!            pe3t_1d(jk) = pdepw_1d(jk+1)-pdepw_1d(jk) 
-!         END DO
-!         pe3t_1d(jpk) = pe3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO
-!
-!         DO jk = 2, jpk
-!            pe3w_1d(jk) = pdept_1d(jk) - pdept_1d(jk-1) 
-!         END DO
-!         pe3w_1d(1  ) = 2._wp * (pdept_1d(1) - pdepw_1d(1)) 
-!      END IF
-!!gm end
-
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) '              Reference 1D z-coordinate depth and scale factors:'
-         WRITE(numout, "(9x,' level  gdept_1d  gdepw_1d  e3t_1d   e3w_1d  ')" )
-         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
-      ENDIF
-      DO jk = 1, jpk                      ! control positivity
-         IF( pe3w_1d (jk) <= 0._wp .OR. pe3t_1d (jk) <= 0._wp )   CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 '    )
-         IF( pdepw_1d(jk) <  0._wp .OR. pdept_1d(jk) <  0._wp )   CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
-      END DO
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_z')
-      !
-   END SUBROUTINE zgr_z
-
-
-   SUBROUTINE zgr_top_bot( k_top , k_bot )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_top_bot  ***
-      !!
-      !! ** Purpose :   set the top and bottom ocean levels
-      !!
-      !! ** Method  :   GYRE case = closed flat box ocean without ocean cavities
-      !!                   k_top = 1     except along north, south, east and west boundaries
-      !!                   k_bot = jpk-1 except along north, south, east and west boundaries
-      !!
-      !! ** Action  : - k_top : first ocean level index
-      !!              - k_bot : last  ocean level index
-      !!----------------------------------------------------------------------
-      INTEGER , DIMENSION(:,:), INTENT(out) ::   k_top , k_bot   ! first & last ocean level
-      !
-      INTEGER  ::   ji, jj               ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! local indices
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! 2D local workspace
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_top_bot : defines the top and bottom ocean levels.'
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~'
-      IF(lwp) WRITE(numout,*) '       GYRE case : closed flat box ocean without ocean cavities'
-      !
-      z2d(:,:) = REAL( jpkm1 , wp )          ! flat bottom
-      !
-      CALL lbc_lnk( z2d, 'T', 1. )           ! applied the lateral boundary condition (here jperio=0 ==>> closed)
-      !
-      k_bot(:,:) = INT( z2d(:,:) )           ! =jpkm1 over the ocean point, =0 elsewhere
-      !
-      k_top(:,:) = MIN( 1 , k_bot(:,:) )     ! =1     over the ocean point, =0 elsewhere
-      !
-   END SUBROUTINE zgr_top_bot
-   
-
-   SUBROUTINE zgr_zco( pdept_1d, pdepw_1d, pe3t_1d, pe3w_1d,   &   ! in : 1D reference vertical coordinate
-      &                pdept   , pdepw   ,                     &   ! out: 3D t & w-points depth
-      &                pe3t    , pe3u    , pe3v   , pe3f   ,   &   !      vertical scale factors
-      &                pe3w    , pe3uw   , pe3vw             )     !          -      -      -
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE zgr_zco  ***
-      !!
-      !! ** Purpose :   define the reference z-coordinate system
-      !!
-      !! ** Method  :   set 3D coord. arrays to reference 1D array 
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pdept_1d, pdepw_1d          ! 1D grid-point depth       [m]
-      REAL(wp), DIMENSION(:)    , INTENT(in   ) ::   pe3t_1d , pe3w_1d           ! 1D vertical scale factors [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pdept, pdepw                ! grid-point depth          [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3t , pe3u , pe3v , pe3f   ! vertical scale factors    [m]
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pe3w , pe3uw, pe3vw         !    -       -      -
-      !
-      INTEGER  ::   jk
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('zgr_zco')
-      !
-      DO jk = 1, jpk
-         pdept(:,:,jk) = pdept_1d(jk)
-         pdepw(:,:,jk) = pdepw_1d(jk)
-         pe3t (:,:,jk) = pe3t_1d (jk)
-         pe3u (:,:,jk) = pe3t_1d (jk)
-         pe3v (:,:,jk) = pe3t_1d (jk)
-         pe3f (:,:,jk) = pe3t_1d (jk)
-         pe3w (:,:,jk) = pe3w_1d (jk)
-         pe3uw(:,:,jk) = pe3w_1d (jk)
-         pe3vw(:,:,jk) = pe3w_1d (jk)
-      END DO
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('zgr_zco')
-      !
-   END SUBROUTINE zgr_zco
-
-   !!======================================================================
-END MODULE usrdef_zgr
