Index: NEMO/trunk/tests/WAD/MY_SRC/bdyini.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/bdyini.F90	(revision 9659)
+++ NEMO/trunk/tests/WAD/MY_SRC/bdyini.F90	(revision 9925)
@@ -44,9 +44,8 @@
    INTEGER, DIMENSION(jp_nseg) ::   jpjnob, jpindt, jpinft, npckgn   !
    INTEGER, DIMENSION(jp_nseg) ::   jpjsob, jpisdt, jpisft, npckgs   !
-
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
-   !! $Id: bdyini.F90 7421 2016-12-01 17:10:41Z flavoni $ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: bdyini.F90 9807 2018-06-15 22:51:16Z lovato $ 
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -86,4 +85,6 @@
       IF(lwm) WRITE ( numond, nambdy )
 
+      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
+      
       ! -----------------------------------------
       ! unstructured open boundaries use control
@@ -112,5 +113,5 @@
    END SUBROUTINE bdy_init
 
-   
+
    SUBROUTINE bdy_segs
       !!----------------------------------------------------------------------
@@ -124,7 +125,4 @@
       !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
       !!----------------------------------------------------------------------      
-
-      ! local variables
-      !-------------------
       INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
       INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
@@ -146,5 +144,5 @@
       INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
-      REAL(wp), DIMENSION(jpi,jpj)       ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
+      REAL(wp), TARGET, DIMENSION(jpi,jpj) ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
       !!
       CHARACTER(LEN=1)                     ::   ctypebdy   !     -        - 
@@ -345,18 +343,18 @@
 
 #if defined key_si3
-        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
-        SELECT CASE( cn_ice(ib_bdy) )                  
-          CASE('none')
+         IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
+         SELECT CASE( cn_ice(ib_bdy) )                  
+         CASE('none')
              IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-             dta_bdy(ib_bdy)%ll_a_i  = .false.
-             dta_bdy(ib_bdy)%ll_ht_i = .false.
-             dta_bdy(ib_bdy)%ll_ht_s = .false.
-          CASE('frs')
+             dta_bdy(ib_bdy)%ll_a_i = .false.
+             dta_bdy(ib_bdy)%ll_h_i = .false.
+             dta_bdy(ib_bdy)%ll_h_s = .false.
+         CASE('frs')
              IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-             dta_bdy(ib_bdy)%ll_a_i  = .true.
-             dta_bdy(ib_bdy)%ll_ht_i = .true.
-             dta_bdy(ib_bdy)%ll_ht_s = .true.
-          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
-        END SELECT
+             dta_bdy(ib_bdy)%ll_a_i = .true.
+             dta_bdy(ib_bdy)%ll_h_i = .true.
+             dta_bdy(ib_bdy)%ll_h_s = .true.
+         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
+         END SELECT
         IF( cn_ice(ib_bdy) /= 'none' ) THEN 
            SELECT CASE( nn_ice_dta(ib_bdy) )                   ! 
@@ -374,8 +372,8 @@
         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
         IF(lwp) WRITE(numout,*)
-
-      ENDDO
-
-     IF (nb_bdy .gt. 0) THEN
+         !
+      END DO
+
+     IF( nb_bdy > 0 ) THEN
         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
           IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
@@ -417,7 +415,7 @@
 
       DO ib_bdy = 1, nb_bdy
-         !
+
          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
-            !
+ 
             icount = icount + 1
             ! No REWIND here because may need to read more than one nambdy_index namelist.
@@ -494,6 +492,4 @@
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
-               !clem nblendta(igrd,ib_bdy) = kdimsz(1)
-               !clem jpbdtau = MAX(jpbdtau, kdimsz(1))
                nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
                jpbdtau = MAX(jpbdtau, MAXVAL(kdimsz))
@@ -885,21 +881,21 @@
                   IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
                ENDIF
-            ENDDO
+            END DO
             idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
             idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc        
-         ENDDO  ! igrd
+         END DO  ! igrd
 
          ! Allocate index arrays for this boundary set
          !--------------------------------------------
          ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) ,   &
+            &      idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
 
          ! Dispatch mapping indices and discrete distances on each processor
@@ -1114,5 +1110,5 @@
          END DO 
 
-      ENDDO
+      END DO
 
       ! ------------------------------------------------------
@@ -1127,34 +1123,15 @@
       bdytmask(:,:) = ssmask(:,:)
 
-      IF( ln_mask_file ) THEN
-         CALL iom_open( cn_mask_file, inum )
-         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
-         CALL iom_close( inum )
-
-         ! Derive mask on U and V grid from mask on T grid
-         bdyumask(:,:) = 0._wp
-         bdyvmask(:,:) = 0._wp
-         DO ij=1, jpjm1
-            DO ii=1, jpim1
-               bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
-               bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
-            END DO
+      ! Derive mask on U and V grid from mask on T grid
+
+      bdyumask(:,:) = 0._wp
+      bdyvmask(:,:) = 0._wp
+      DO ij = 1, jpjm1
+         DO ii = 1, jpim1
+            bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
+            bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
          END DO
-         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
-
-      ENDIF ! ln_mask_file=.TRUE.
-      
-      IF( .NOT.ln_mask_file ) THEN
-         ! If .not. ln_mask_file then we need to derive mask on U and V grid from mask on T grid here.
-         bdyumask(:,:) = 0._wp
-         bdyvmask(:,:) = 0._wp
-         DO ij = 1, jpjm1
-            DO ii = 1, jpim1
-               bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1, ij )
-               bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
-            END DO
-         END DO
-         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
-      ENDIF
+      END DO
+      CALL lbc_lnk_multi( bdyumask, 'U', 1. , bdyvmask, 'V', 1. )   ! Lateral boundary cond. 
 
       ! bdy masks are now set to zero on boundary points:
@@ -1178,9 +1155,10 @@
         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
           bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0._wp
-        ENDDO
-      ENDDO
+        END DO
+      END DO
 
       ! For the flagu/flagv calculation below we require a version of fmask without
       ! the land boundary condition (shlat) included:
+      zfmask(:,:) = 0
       DO ij = 2, jpjm1
          DO ii = 2, jpim1
@@ -1191,8 +1169,6 @@
 
       ! Lateral boundary conditions
-      CALL lbc_lnk( zfmask       , 'F', 1. )
-      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
-      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
-
+      CALL lbc_lnk( zfmask, 'F', 1. ) 
+      CALL lbc_lnk_multi( bdyumask, 'U', 1. , bdyvmask, 'V', 1., bdytmask, 'T', 1. )
       DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
 
@@ -1206,5 +1182,5 @@
          ! flagu =  1 : u is normal to the boundary and is direction is inward
   
-         DO igrd = 1,jpbgrd 
+         DO igrd = 1, jpbgrd 
             SELECT CASE( igrd )
                CASE( 1 )   ;   pmask => umask   (:,:,1)   ;   i_offset = 0
@@ -1310,5 +1286,4 @@
       !--------
       IF( nb_bdy>0 )   DEALLOCATE( nbidta, nbjdta, nbrdta )
-      !
       !
    END SUBROUTINE bdy_segs
@@ -1690,4 +1665,5 @@
    END SUBROUTINE bdy_ctl_seg
 
+
    SUBROUTINE bdy_ctl_corn( ib1, ib2 )
       !!----------------------------------------------------------------------
Index: NEMO/trunk/tests/WAD/MY_SRC/domain.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/domain.F90	(revision 9659)
+++ 	(revision )
@@ -1,676 +1,0 @@
-MODULE domain
-   !!==============================================================================
-   !!                       ***  MODULE domain   ***
-   !! Ocean initialization : domain initialization
-   !!==============================================================================
-   !! History :  OPA  !  1990-10  (C. Levy - G. Madec)  Original code
-   !!                 !  1992-01  (M. Imbard) insert time step initialization
-   !!                 !  1996-06  (G. Madec) generalized vertical coordinate 
-   !!                 !  1997-02  (G. Madec) creation of domwri.F
-   !!                 !  2001-05  (E.Durand - G. Madec) insert closed sea
-   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
-   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
-   !!            3.3  !  2010-11  (G. Madec)  initialisation in C1D configuration
-   !!            3.6  !  2013     ( J. Simeon, C. Calone, G. Madec, C. Ethe ) Online coarsening of outputs
-   !!            3.7  !  2015-11  (G. Madec, A. Coward)  time varying zgr by default
-   !!            4.0  !  2016-10  (G. Madec, S. Flavoni)  domain configuration / user defined interface
-   !!----------------------------------------------------------------------
-   
-   !!----------------------------------------------------------------------
-   !!   dom_init      : initialize the space and time domain
-   !!   dom_glo       : initialize global domain <--> local domain indices
-   !!   dom_nam       : read and contral domain namelists
-   !!   dom_ctl       : control print for the ocean domain
-   !!   domain_cfg    : read the global domain size in domain configuration file
-   !!   cfg_write     : create the domain configuration file
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean variables
-   USE dom_oce        ! domain: ocean
-   USE sbc_oce        ! surface boundary condition: ocean
-   USE trc_oce        ! shared ocean & passive tracers variab
-   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 configuration
-   USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine)
-   USE wet_dry        ! wetting and drying
-   !
-   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
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dom_init     ! called by nemogcm.F90
-   PUBLIC   domain_cfg   ! called by nemogcm.F90
-
-   !!-------------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: domain.F90 7587 2017-01-20 15:17:56Z acc $
-   !! Software governed by the CeCILL licence        (NEMOGCM/NEMO_CeCILL.txt)
-   !!-------------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dom_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_init  ***
-      !!                    
-      !! ** Purpose :   Domain initialization. Call the routines that are 
-      !!              required to create the arrays which define the space 
-      !!              and time domain of the ocean model.
-      !!
-      !! ** Method  : - dom_msk: compute the masks from the bathymetry file
-      !!              - dom_hgr: compute or read the horizontal grid-point position
-      !!                         and scale factors, and the coriolis factor
-      !!              - dom_zgr: define the vertical coordinate and the bathymetry
-      !!              - dom_wri: create the meshmask file (ln_meshmask=T)
-      !!              - 1D configuration, move Coriolis, u and v at T-point
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj, jk, ik   ! dummy loop indices
-      INTEGER ::   iconf = 0    ! local integers
-      CHARACTER (len=64) ::   cform = "(A12, 3(A13, I7))" 
-      INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) THEN         ! Ocean domain Parameters (control print)
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_init : domain initialization'
-         WRITE(numout,*) '~~~~~~~~'
-         !
-         WRITE(numout,*)     '   Domain info'
-         WRITE(numout,*)     '      dimension of model:'
-         WRITE(numout,*)     '             Local domain      Global domain       Data domain '
-         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):'
-         WRITE(numout,*)     '              jpni    : ', jpni, '   nn_hls  : ', nn_hls
-         WRITE(numout,*)     '              jpnj    : ', jpnj, '   nn_hls  : ', nn_hls
-         WRITE(numout,*)     '              jpnij   : ', jpnij
-         WRITE(numout,*)     '      lateral boundary of the Global domain : jperio  = ', jperio
-         SELECT CASE ( jperio )
-         CASE( 0 )   ;   WRITE(numout,*) '         (i.e. closed)'
-         CASE( 1 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west)'
-         CASE( 2 )   ;   WRITE(numout,*) '         (i.e. equatorial symmetric)'
-         CASE( 3 )   ;   WRITE(numout,*) '         (i.e. north fold with T-point pivot)'
-         CASE( 4 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north fold with T-point pivot)'
-         CASE( 5 )   ;   WRITE(numout,*) '         (i.e. north fold with F-point pivot)'
-         CASE( 6 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north fold with F-point pivot)'
-         CASE( 7 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north-south)'
-         CASE DEFAULT
-            CALL ctl_stop( 'jperio is out of range' )
-         END SELECT
-         WRITE(numout,*)     '      Ocean model configuration used:'
-         WRITE(numout,*)     '         cn_cfg = ', TRIM( cn_cfg ), '   nn_cfg = ', nn_cfg
-      ENDIF
-      !
-      !           !==  Reference coordinate system  ==!
-      !
-      CALL dom_glo                     ! global domain versus local domain
-      CALL dom_nam                     ! read namelist ( namrun, namdom )
-      CALL dom_hgr                     ! Horizontal mesh
-      CALL dom_zgr( ik_top, ik_bot )   ! Vertical mesh and bathymetry
-      CALL dom_msk( ik_top, ik_bot )   ! Masks
-      IF( ln_closea ) CALL dom_clo     ! ln_closea=T : closed seas included in the simulation
-                                       ! Read in masks to define closed seas and lakes 
-      !
-      DO jj = 1, jpj                   ! depth of the iceshelves
-         DO ji = 1, jpi
-            ik = mikt(ji,jj)
-            risfdep(ji,jj) = gdepw_0(ji,jj,ik)
-         END DO
-      END DO
-      !
-      ht_0(:,:) = 0._wp  ! Reference ocean thickness
-      hu_0(:,:) = 0._wp
-      hv_0(:,:) = 0._wp
-      DO jk = 1, jpk
-         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
-         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
-         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !
-      !           !==  time varying part of coordinate system  ==!
-      !
-      IF( ln_linssh ) THEN       != Fix in time : set to the reference one for all
-      !
-         !       before        !          now          !       after         !
-            gdept_b = gdept_0  ;   gdept_n = gdept_0   !        ---          ! depth of grid-points
-            gdepw_b = gdepw_0  ;   gdepw_n = gdepw_0   !        ---          !
-                                   gde3w_n = gde3w_0   !        ---          !
-         !                                                                  
-              e3t_b =   e3t_0  ;     e3t_n =   e3t_0   ;   e3t_a =  e3t_0    ! scale factors
-              e3u_b =   e3u_0  ;     e3u_n =   e3u_0   ;   e3u_a =  e3u_0    !
-              e3v_b =   e3v_0  ;     e3v_n =   e3v_0   ;   e3v_a =  e3v_0    !
-                                     e3f_n =   e3f_0   !        ---          !
-              e3w_b =   e3w_0  ;     e3w_n =   e3w_0   !        ---          !
-             e3uw_b =  e3uw_0  ;    e3uw_n =  e3uw_0   !        ---          !
-             e3vw_b =  e3vw_0  ;    e3vw_n =  e3vw_0   !        ---          !
-         !
-         z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) )     ! _i mask due to ISF
-         z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
-         !
-         !        before       !          now          !       after         !
-                                      ht_n =    ht_0   !                     ! water column thickness
-               hu_b =    hu_0  ;      hu_n =    hu_0   ;    hu_a =    hu_0   ! 
-               hv_b =    hv_0  ;      hv_n =    hv_0   ;    hv_a =    hv_0   !
-            r1_hu_b = z1_hu_0  ;   r1_hu_n = z1_hu_0   ; r1_hu_a = z1_hu_0   ! inverse of water column thickness
-            r1_hv_b = z1_hv_0  ;   r1_hv_n = z1_hv_0   ; r1_hv_a = z1_hv_0   !
-         !
-         !
-      ELSE                       != time varying : initialize before/now/after variables
-         !
-         IF( .NOT.l_offline )  CALL dom_vvl_init 
-         !
-      ENDIF
-      !
-      IF( lk_c1d         )   CALL cor_c1d       ! 1D configuration: Coriolis set at T-point
-      !
-      IF( ln_meshmask .AND. .NOT.ln_iscpl )                        CALL dom_wri     ! Create a domain file
-      IF( ln_meshmask .AND.      ln_iscpl .AND. .NOT.ln_rstart )   CALL dom_wri     ! Create a domain file
-      IF(                                       .NOT.ln_rstart )   CALL dom_ctl     ! Domain control
-      !
-      IF( ln_write_cfg )   CALL cfg_write         ! create the configuration file
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_init :   ==>>>   END of domain initialization'
-         WRITE(numout,*) '~~~~~~~~'
-         WRITE(numout,*) 
-      ENDIF
-      !
-   END SUBROUTINE dom_init
-
-
-   SUBROUTINE dom_glo
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_glo  ***
-      !!
-      !! ** Purpose :   initialization of global domain <--> local domain indices
-      !!
-      !! ** Method  :   
-      !!
-      !! ** Action  : - mig , mjg : local  domain indices ==> global domain indices
-      !!              - mi0 , mi1 : global domain indices ==> local  domain indices
-      !!              - mj0,, mj1   (global point not in the local domain ==> mi0>mi1 and/or mj0>mj1)
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj   ! dummy loop argument
-      !!----------------------------------------------------------------------
-      !
-      DO ji = 1, jpi                 ! local domain indices ==> global domain indices
-        mig(ji) = ji + nimpp - 1
-      END DO
-      DO jj = 1, jpj
-        mjg(jj) = jj + njmpp - 1
-      END DO
-      !                              ! global domain indices ==> local domain indices
-      !                                   ! (return (m.0,m.1)=(1,0) if data domain gridpoint is to the west/south of the 
-      !                                   ! local domain, or (m.0,m.1)=(jp.+1,jp.) to the east/north of local domain. 
-      DO ji = 1, jpiglo
-        mi0(ji) = MAX( 1 , MIN( ji - nimpp + 1, jpi+1 ) )
-        mi1(ji) = MAX( 0 , MIN( ji - nimpp + 1, jpi   ) )
-      END DO
-      DO jj = 1, jpjglo
-        mj0(jj) = MAX( 1 , MIN( jj - njmpp + 1, jpj+1 ) )
-        mj1(jj) = MAX( 0 , MIN( jj - njmpp + 1, jpj   ) )
-      END DO
-      IF(lwp) THEN                   ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_glo : domain: global <<==>> local '
-         WRITE(numout,*) '~~~~~~~ '
-         WRITE(numout,*) '   global domain:   jpiglo = ', jpiglo, ' jpjglo = ', jpjglo, ' jpkglo = ', jpkglo
-         WRITE(numout,*) '   local  domain:   jpi    = ', jpi   , ' jpj    = ', jpj   , ' jpk    = ', jpk
-         WRITE(numout,*)
-         WRITE(numout,*) '   conversion from local to global domain indices (and vise versa) done'
-         IF( nn_print >= 1 ) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) '          conversion local  ==> global i-index domain (mig)'
-            WRITE(numout,25)              (mig(ji),ji = 1,jpi)
-            WRITE(numout,*)
-            WRITE(numout,*) '          conversion global ==> local  i-index domain'
-            WRITE(numout,*) '             starting index (mi0)'
-            WRITE(numout,25)              (mi0(ji),ji = 1,jpiglo)
-            WRITE(numout,*) '             ending index (mi1)'
-            WRITE(numout,25)              (mi1(ji),ji = 1,jpiglo)
-            WRITE(numout,*)
-            WRITE(numout,*) '          conversion local  ==> global j-index domain (mjg)'
-            WRITE(numout,25)              (mjg(jj),jj = 1,jpj)
-            WRITE(numout,*)
-            WRITE(numout,*) '          conversion global ==> local  j-index domain'
-            WRITE(numout,*) '             starting index (mj0)'
-            WRITE(numout,25)              (mj0(jj),jj = 1,jpjglo)
-            WRITE(numout,*) '             ending index (mj1)'
-            WRITE(numout,25)              (mj1(jj),jj = 1,jpjglo)
-         ENDIF
-      ENDIF
- 25   FORMAT( 100(10x,19i4,/) )
-      !
-   END SUBROUTINE dom_glo
-
-
-   SUBROUTINE dom_nam
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_nam  ***
-      !!                    
-      !! ** Purpose :   read domaine namelists and print the variables.
-      !!
-      !! ** input   : - namrun namelist
-      !!              - namdom namelist
-      !!              - namnc4 namelist   ! "key_netcdf4" only
-      !!----------------------------------------------------------------------
-      USE ioipsl
-      !!
-      INTEGER  ::   ios   ! Local integer
-      !
-      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,                 &
-         &             nn_no   , cn_exp   , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl ,     &
-         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     &
-         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     &
-         &             ln_cfmeta, ln_iscpl
-      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask
-#if defined key_netcdf4
-      NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip
-#endif
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_nam  : domain initialization through namelist read'
-         WRITE(numout,*) '~~~~~~~ '
-      ENDIF
-      !
-      REWIND( numnam_ref )              ! Namelist namrun in reference namelist : Parameters of the run
-      READ  ( numnam_ref, namrun, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namrun in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namrun in configuration namelist : Parameters of the run
-      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 )
-902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp )
-      IF(lwm) WRITE ( numond, namrun )
-      !
-      IF(lwp) THEN                  ! control print
-         WRITE(numout,*) '   Namelist : namrun'
-         WRITE(numout,*) '      job number                      nn_no           = ', nn_no
-         WRITE(numout,*) '      experiment name for output      cn_exp          = ', TRIM( cn_exp           )
-         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in    = ', TRIM( cn_ocerst_in     )
-         WRITE(numout,*) '      restart input directory         cn_ocerst_indir = ', TRIM( cn_ocerst_indir  )
-         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out   = ', TRIM( cn_ocerst_out    )
-         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', TRIM( cn_ocerst_outdir )
-         WRITE(numout,*) '      restart logical                 ln_rstart       = ', ln_rstart
-         WRITE(numout,*) '      start with forward time step    nn_euler        = ', nn_euler
-         WRITE(numout,*) '      control of time step            nn_rstctl       = ', nn_rstctl
-         WRITE(numout,*) '      number of the first time step   nn_it000        = ', nn_it000
-         WRITE(numout,*) '      number of the last time step    nn_itend        = ', nn_itend
-         WRITE(numout,*) '      initial calendar date aammjj    nn_date0        = ', nn_date0
-         WRITE(numout,*) '      initial time of day in hhmm     nn_time0        = ', nn_time0
-         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy        = ', nn_leapy
-         WRITE(numout,*) '      initial state output            nn_istate       = ', nn_istate
-         IF( ln_rst_list ) THEN
-            WRITE(numout,*) '      list of restart dump times      nn_stocklist    =', nn_stocklist
-         ELSE
-            WRITE(numout,*) '      frequency of restart file       nn_stock        = ', nn_stock
-         ENDIF
-         WRITE(numout,*) '      frequency of output file        nn_write        = ', nn_write
-         WRITE(numout,*) '      mask land points                ln_mskland      = ', ln_mskland
-         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta       = ', ln_cfmeta
-         WRITE(numout,*) '      overwrite an existing file      ln_clobber      = ', ln_clobber
-         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz      = ', nn_chunksz
-         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl        = ', ln_iscpl
-      ENDIF
-
-      no = nn_no                    ! conversion DOCTOR names into model names (this should disappear soon)
-      cexper = cn_exp
-      nrstdt = nn_rstctl
-      nit000 = nn_it000
-      nitend = nn_itend
-      ndate0 = nn_date0
-      nleapy = nn_leapy
-      ninist = nn_istate
-      nstock = nn_stock
-      nstocklist = nn_stocklist
-      nwrite = nn_write
-      neuler = nn_euler
-      IF( neuler == 1 .AND. .NOT. ln_rstart ) THEN
-         IF(lwp) WRITE(numout,*)  
-         IF(lwp) WRITE(numout,*)'   ==>>>   Start from rest (ln_rstart=F)'
-         IF(lwp) WRITE(numout,*)'           an Euler initial time step is used : nn_euler is forced to 0 '   
-         neuler = 0
-      ENDIF
-      !                             ! control of output frequency
-      IF( nstock == 0 .OR. nstock > nitend ) THEN
-         WRITE(ctmp1,*) 'nstock = ', nstock, ' it is forced to ', nitend
-         CALL ctl_warn( ctmp1 )
-         nstock = nitend
-      ENDIF
-      IF ( nwrite == 0 ) THEN
-         WRITE(ctmp1,*) 'nwrite = ', nwrite, ' it is forced to ', nitend
-         CALL ctl_warn( ctmp1 )
-         nwrite = nitend
-      ENDIF
-
-#if defined key_agrif
-      IF( Agrif_Root() ) THEN
-#endif
-      SELECT CASE ( nleapy )        ! Choose calendar for IOIPSL
-      CASE (  1 ) 
-         CALL ioconf_calendar('gregorian')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "gregorian", i.e. leap year'
-      CASE (  0 )
-         CALL ioconf_calendar('noleap')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "noleap", i.e. no leap year'
-      CASE ( 30 )
-         CALL ioconf_calendar('360d')
-         IF(lwp) WRITE(numout,*) '   The IOIPSL calendar is "360d", i.e. 360 days in a year'
-      END SELECT
-#if defined key_agrif
-      ENDIF
-#endif
-
-      REWIND( numnam_ref )              ! Namelist namdom in reference namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_ref, namdom, IOSTAT = ios, ERR = 903)
-903   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdom in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namdom in configuration namelist : space & time domain (bathymetry, mesh, timestep)
-      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 )
-904   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namdom )
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist : namdom   ---   space & time domain'
-         WRITE(numout,*) '      linear free surface (=T)                ln_linssh   = ', ln_linssh
-         WRITE(numout,*) '      create mesh/mask file                   ln_meshmask = ', ln_meshmask
-         WRITE(numout,*) '      treshold to open the isf cavity         rn_isfhmin  = ', rn_isfhmin, ' [m]'
-         WRITE(numout,*) '      ocean time step                         rn_rdt      = ', rn_rdt
-         WRITE(numout,*) '      asselin time filter parameter           rn_atfp     = ', rn_atfp
-         WRITE(numout,*) '      online coarsening of dynamical fields   ln_crs      = ', ln_crs
-      ENDIF
-      !
-      !          ! conversion DOCTOR names into model names (this should disappear soon)
-      atfp = rn_atfp
-      rdt  = rn_rdt
-
-#if defined key_netcdf4
-      !                             ! NetCDF 4 case   ("key_netcdf4" defined)
-      REWIND( numnam_ref )              ! Namelist namnc4 in reference namelist : NETCDF
-      READ  ( numnam_ref, namnc4, IOSTAT = ios, ERR = 907)
-907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namnc4 in reference namelist', lwp )
-      REWIND( numnam_cfg )              ! Namelist namnc4 in configuration namelist : NETCDF
-      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 )
-908   IF( ios >  0 )   CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp )
-      IF(lwm) WRITE( numond, namnc4 )
-
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) '   Namelist namnc4 - Netcdf4 chunking parameters'
-         WRITE(numout,*) '      number of chunks in i-dimension             nn_nchunks_i = ', nn_nchunks_i
-         WRITE(numout,*) '      number of chunks in j-dimension             nn_nchunks_j = ', nn_nchunks_j
-         WRITE(numout,*) '      number of chunks in k-dimension             nn_nchunks_k = ', nn_nchunks_k
-         WRITE(numout,*) '      apply netcdf4/hdf5 chunking & compression   ln_nc4zip    = ', ln_nc4zip
-      ENDIF
-
-      ! Put the netcdf4 settings into a simple structure (snc4set, defined in in_out_manager module)
-      ! Note the chunk size in the unlimited (time) dimension will be fixed at 1
-      snc4set%ni   = nn_nchunks_i
-      snc4set%nj   = nn_nchunks_j
-      snc4set%nk   = nn_nchunks_k
-      snc4set%luse = ln_nc4zip
-#else
-      snc4set%luse = .FALSE.        ! No NetCDF 4 case
-#endif
-      !
-   END SUBROUTINE dom_nam
-
-
-   SUBROUTINE dom_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_ctl  ***
-      !!
-      !! ** Purpose :   Domain control.
-      !!
-      !! ** Method  :   compute and print extrema of masked scale factors
-      !!----------------------------------------------------------------------
-      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
-      INTEGER, DIMENSION(2) ::   iloc   ! 
-      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
-      !!----------------------------------------------------------------------
-      !
-      IF(lk_mpp) THEN
-         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 )
-         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 )
-         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 )
-         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 )
-      ELSE
-         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
-         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
-         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )    
-         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )    
-         !
-         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
-         iimi1 = iloc(1) + nimpp - 1
-         ijmi1 = iloc(2) + njmpp - 1
-         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
-         iimi2 = iloc(1) + nimpp - 1
-         ijmi2 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp )
-         iima1 = iloc(1) + nimpp - 1
-         ijma1 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp )
-         iima2 = iloc(1) + nimpp - 1
-         ijma2 = iloc(2) + njmpp - 1
-      ENDIF
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
-         WRITE(numout,*) '~~~~~~~'
-         WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1
-         WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1
-         WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2
-         WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2
-      ENDIF
-      !
-   END SUBROUTINE dom_ctl
-
-
-   SUBROUTINE domain_cfg( ldtxt, cd_cfg, kk_cfg, kpi, kpj, kpk, kperio )
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_nam  ***
-      !!                    
-      !! ** Purpose :   read the domain size in domain configuration file
-      !!
-      !! ** Method  :   read the cn_domcfg NetCDF file
-      !!----------------------------------------------------------------------
-      CHARACTER(len=*), DIMENSION(:), INTENT(out) ::   ldtxt           ! stored print information
-      CHARACTER(len=*)              , INTENT(out) ::   cd_cfg          ! configuration name
-      INTEGER                       , INTENT(out) ::   kk_cfg          ! configuration resolution
-      INTEGER                       , INTENT(out) ::   kpi, kpj, kpk   ! global domain sizes 
-      INTEGER                       , INTENT(out) ::   kperio          ! lateral global domain b.c. 
-      !
-      INTEGER ::   inum, ii   ! local integer
-      REAL(wp) ::   zorca_res                     ! local scalars
-      REAL(wp) ::   ziglo, zjglo, zkglo, zperio   !   -      -
-      !!----------------------------------------------------------------------
-      !
-      ii = 1
-      WRITE(ldtxt(ii),*) '           '                                                    ;   ii = ii+1
-      WRITE(ldtxt(ii),*) 'domain_cfg : domain size read in ', TRIM( cn_domcfg ), ' file'  ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '~~~~~~~~~~ '                                                    ;   ii = ii+1
-      !
-      CALL iom_open( cn_domcfg, inum )
-      !
-      !                                   !- ORCA family specificity
-      IF(  iom_varid( inum, 'ORCA'       , ldstop = .FALSE. ) > 0  .AND.  &
-         & iom_varid( inum, 'ORCA_index' , ldstop = .FALSE. ) > 0    ) THEN
-         !
-         cd_cfg = 'ORCA'
-         CALL iom_get( inum, 'ORCA_index', zorca_res )   ;   kk_cfg = INT( zorca_res )
-         !
-         WRITE(ldtxt(ii),*) '       '                                                     ;   ii = ii+1
-         WRITE(ldtxt(ii),*) '       ==>>>   ORCA configuration '                          ;   ii = ii+1
-         WRITE(ldtxt(ii),*) '       '                                                     ;   ii = ii+1
-         !
-      ELSE                                !- cd_cfg & k_cfg are not used
-         cd_cfg = 'UNKNOWN'
-         kk_cfg = -9999999
-                                          !- or they may be present as global attributes 
-                                          !- (netcdf only)  
-         IF( iom_file(inum)%iolib == jpnf90 ) THEN
-            CALL iom_getatt( inum, 'cn_cfg', cd_cfg )  ! returns   !  if not found
-            CALL iom_getatt( inum, 'nn_cfg', kk_cfg )  ! returns -999 if not found
-            IF( TRIM(cd_cfg) .EQ. '!') cd_cfg = 'UNKNOWN'
-            IF( kk_cfg .EQ. -999     ) kk_cfg = -9999999
-         ENDIF
-         !
-      ENDIF
-      !
-      CALL iom_get( inum, 'jpiglo', ziglo  )   ;   kpi = INT( ziglo )
-      CALL iom_get( inum, 'jpjglo', zjglo  )   ;   kpj = INT( zjglo )
-      CALL iom_get( inum, 'jpkglo', zkglo  )   ;   kpk = INT( zkglo )
-      CALL iom_get( inum, 'jperio', zperio )   ;   kperio = INT( zperio )
-      CALL iom_close( inum )
-      !
-      WRITE(ldtxt(ii),*) '   cn_cfg = ', TRIM(cd_cfg), '   nn_cfg = ', kk_cfg             ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpiglo = ', kpi                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpjglo = ', kpj                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   jpkglo = ', kpk                                              ;   ii = ii+1
-      WRITE(ldtxt(ii),*) '   type of global domain lateral boundary   jperio = ', kperio  ;   ii = ii+1
-      !        
-   END SUBROUTINE domain_cfg
-   
-   
-   SUBROUTINE cfg_write
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE cfg_write  ***
-      !!                   
-      !! ** Purpose :   Create the "cn_domcfg_out" 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 required to set up an
-      !!              ocean configuration.
-      !!
-      !! ** output file :   domcfg_out.nc : domain size, characteristics, horizontal 
-      !!                       mesh, Coriolis parameter, and vertical scale factors
-      !!                    NB: also contain ORCA family information
-      !!----------------------------------------------------------------------
-      INTEGER           ::   ji, jj, jk   ! dummy loop indices
-      INTEGER           ::   izco, izps, isco, icav
-      INTEGER           ::   inum     ! local units
-      CHARACTER(len=21) ::   clnam    ! filename (mesh and mask informations)
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
-      !!----------------------------------------------------------------------
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'cfg_write : create the domain configuration file (', TRIM(cn_domcfg_out),'.nc)'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
-      !
-      !                       ! ============================= !
-      !                       !  create 'domcfg_out.nc' file  !
-      !                       ! ============================= !
-      !         
-      clnam = cn_domcfg_out  ! filename (configuration information)
-      CALL iom_open( TRIM(clnam), inum, ldwrt = .TRUE., kiolib = jprstlib )
-      
-      !
-      !                             !==  ORCA family specificities  ==!
-      IF( cn_cfg == "ORCA" ) THEN
-         CALL iom_rstput( 0, 0, inum, 'ORCA'      , 1._wp            , ktype = jp_i4 )
-         CALL iom_rstput( 0, 0, inum, 'ORCA_index', REAL( nn_cfg, wp), ktype = jp_i4 )         
-      ENDIF
-      !
-      !                             !==  global domain size  ==!
-      !
-      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( jpk   , wp), ktype = jp_i4 )
-      !
-      !                             !==  domain characteristics  ==!
-      !
-      !                                   ! lateral boundary of the global domain
-      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, 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, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )
-      !
-      !                             !==  horizontal mesh  !
-      !
-      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 )   ! i-scale factors (e1.)
-      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 )   ! j-scale factors (e2.)
-      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 )
-      !
-      !                             !==  vertical mesh  ==!
-      !                                                     
-      CALL iom_rstput( 0, 0, inum, 'e3t_1d'  , e3t_1d , ktype = jp_r8 )   ! reference 1D-coordinate
-      CALL iom_rstput( 0, 0, inum, 'e3w_1d'  , e3w_1d , ktype = jp_r8 )
-      !
-      CALL iom_rstput( 0, 0, inum, 'e3t_0'   , e3t_0  , ktype = jp_r8 )   ! vertical scale factors
-      CALL iom_rstput( 0, 0, inum, 'e3u_0'   , e3u_0  , ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'e3v_0'   , e3v_0  , ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'e3f_0'   , e3f_0  , ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'e3w_0'   , e3w_0  , ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'e3uw_0'  , e3uw_0 , ktype = jp_r8 )
-      CALL iom_rstput( 0, 0, inum, 'e3vw_0'  , e3vw_0 , ktype = jp_r8 )
-      !                                         
-      !                             !==  wet top and bottom level  ==!   (caution: multiplied by ssmask)
-      !
-      CALL iom_rstput( 0, 0, inum, 'top_level'    , REAL( mikt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points (ISF)
-      CALL iom_rstput( 0, 0, inum, 'bottom_level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 )   ! nb of ocean T-points
-      !
-      IF( ln_sco ) THEN             ! s-coordinate: store grid stiffness ratio  (Not required anyway)
-         CALL dom_stiff( z2d )
-         CALL iom_rstput( 0, 0, inum, 'stiffness', z2d )        !    ! Max. grid stiffness ratio
-      ENDIF
-      !
-      IF( ll_wd ) THEN              ! wetting and drying domain
-         CALL iom_rstput( 0, 0, inum, 'ht_0'   , ht_0   , ktype = jp_r8 )
-      ENDIF
-      !
-      ! Add some global attributes ( netcdf only )
-      IF( iom_file(inum)%iolib == jpnf90 ) THEN
-         CALL iom_putatt( inum, 'nn_cfg', nn_cfg )
-         CALL iom_putatt( inum, 'cn_cfg', TRIM(cn_cfg) )
-      ENDIF
-      !
-      !                                ! ============================
-      !                                !        close the files 
-      !                                ! ============================
-      CALL iom_close( inum )
-      !
-   END SUBROUTINE cfg_write
-
-   !!======================================================================
-END MODULE domain
Index: NEMO/trunk/tests/WAD/MY_SRC/usrdef_hgr.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/usrdef_hgr.F90	(revision 9659)
+++ NEMO/trunk/tests/WAD/MY_SRC/usrdef_hgr.F90	(revision 9925)
@@ -7,5 +7,5 @@
    !! User defined :   mesh and Coriolis parameter of a user configuration
    !!======================================================================
-   !! History :  NEMO  ! 2016-08  (S. Flavoni, G. Madec)  Original code
+   !! History :  4.0 ! 2016-03  (S. Flavoni) 
    !!----------------------------------------------------------------------
 
@@ -27,7 +27,7 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
Index: NEMO/trunk/tests/WAD/MY_SRC/usrdef_istate.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/usrdef_istate.F90	(revision 9659)
+++ NEMO/trunk/tests/WAD/MY_SRC/usrdef_istate.F90	(revision 9925)
@@ -7,5 +7,5 @@
    !! User defined : set the initial state of a user configuration
    !!======================================================================
-   !! History :  NEMO ! 2016-03  (S. Flavoni, G. Madec) Original code
+   !! History :  4.0 ! 2016-03  (S. Flavoni) Original code
    !!----------------------------------------------------------------------
 
@@ -24,10 +24,10 @@
    PRIVATE
 
-   PUBLIC   usr_def_istate   ! called by istate.F90
+   PUBLIC   usr_def_istate   ! called in istate.F90
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
Index: NEMO/trunk/tests/WAD/MY_SRC/usrdef_nam.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/usrdef_nam.F90	(revision 9659)
+++ NEMO/trunk/tests/WAD/MY_SRC/usrdef_nam.F90	(revision 9925)
@@ -7,5 +7,5 @@
    !! User defined : set the domain characteristics of a user configuration
    !!======================================================================
-   !! History :  NEMO ! 2016-03  (S. Flavoni, G. Madec)  Original code
+   !! History :  4.0  ! 2016-03  (S. Flavoni, G. Madec)  Original code
    !!----------------------------------------------------------------------
 
@@ -33,7 +33,7 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id:$ 
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS
Index: NEMO/trunk/tests/WAD/MY_SRC/usrdef_zgr.F90
===================================================================
--- NEMO/trunk/tests/WAD/MY_SRC/usrdef_zgr.F90	(revision 9659)
+++ NEMO/trunk/tests/WAD/MY_SRC/usrdef_zgr.F90	(revision 9925)
@@ -7,10 +7,10 @@
    !! Ocean domain : user defined vertical coordinate system 
    !!======================================================================
-   !! History :  4.0  ! 2016-08  (G. Madec, S. Flavoni)  Original code
+   !! History :  4.0  ! 2016-06  (G. Madec)  Original code
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   usr_def_zgr   : user defined vertical coordinate system (required)
-   !!       zgr_z1d   : reference 1D z-coordinate 
+   !!       zgr_z     : reference 1D z-coordinate 
    !!---------------------------------------------------------------------
    USE oce            ! ocean variables
@@ -27,12 +27,12 @@
    PRIVATE
 
-   PUBLIC   usr_def_zgr   ! called by domzgr.F90
-
-  !! * Substitutions
+   PUBLIC   usr_def_zgr        ! called by domzgr.F90
+
+   !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id:$
+   !! Software governed by the CeCILL licence     (./LICENSE)
    !!----------------------------------------------------------------------
 CONTAINS             
@@ -41,6 +41,6 @@
       &                    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,                      &   !     -      -      -
+      &                    pe3t  , pe3u  , pe3v   , pe3f ,             &   ! vertical scale factors
+      &                    pe3w  , pe3uw , pe3vw         ,             &   !     -      -      -
       &                    k_top  , k_bot    )                             ! top & bottom ocean level
       !!---------------------------------------------------------------------
@@ -255,5 +255,5 @@
       END DO
       !     
-      CALL zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
+      CALL zgr_z( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! Reference z-coordinate system
       !
       !
@@ -334,7 +334,7 @@
 
 
-   SUBROUTINE zgr_z1d( pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d )   ! 1D reference vertical coordinate
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE zgr_z1d  ***
+   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 
@@ -363,5 +363,5 @@
       IF(lwp) THEN                         ! Parameter print
          WRITE(numout,*)
-         WRITE(numout,*) '    zgr_z1d : Reference vertical z-coordinates: uniform dz = ', rn_dz
+         WRITE(numout,*) '    zgr_z : Reference vertical z-coordinates: uniform dz = ', rn_dz
          WRITE(numout,*) '    ~~~~~~~'
       ENDIF
@@ -385,5 +385,5 @@
       ENDIF
       !
-   END SUBROUTINE zgr_z1d
+   END SUBROUTINE zgr_z
    
    !!======================================================================
