Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/ARCH/arch-macport_osx.fcm
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/ARCH/arch-macport_osx.fcm	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/ARCH/arch-macport_osx.fcm	(revision 7522)
@@ -54,5 +54,5 @@
 %CPP	             cpp-mp-4.8
 %FC	             mpif90 
-%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer 
+%FCFLAGS             -fdefault-real-8 -O3 -funroll-all-loops -fcray-pointer -ffree-line-length-none
 %FFLAGS              %FCFLAGS
 %LD                  %FC
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist_cfg	(revision 7522)
@@ -101,29 +101,23 @@
 !          !  file name  ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
 !          !             !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-   sn_tem  = 'dyna_grid_T' ,    120            , 'votemper' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_sal  = 'dyna_grid_T' ,    120            , 'vosaline' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_mld  = 'dyna_grid_T' ,    120            , 'somixhgt' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflup' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_fmf  = 'dyna_grid_T' ,    120            , 'iowaflup' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_ice  = 'dyna_grid_T' ,    120            , 'soicecov' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_qsr  = 'dyna_grid_T' ,    120            , 'soshfldo' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_wnd  = 'dyna_grid_T' ,    120            , 'sowindsp' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_uwd  = 'dyna_grid_U' ,    120            , 'vozocrtx' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_vwd  = 'dyna_grid_V' ,    120            , 'vomecrty' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_wwd  = 'dyna_grid_W' ,    120            , 'vovecrtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_avt  = 'dyna_grid_W' ,    120            , 'voddmavs' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_ubl  = 'dyna_grid_U' ,    120            , 'sobblcox' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_vbl  = 'dyna_grid_V' ,    120            , 'sobblcoy' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_ahu  = 'dyna_grid_U' ,    120            , 'vozoahtu' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_ahv  = 'dyna_grid_V' ,    120            , 'vomeahtv' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_ahw  = 'dyna_grid_W' ,    120            , 'voveahtz' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_eiu  = 'dyna_grid_U' ,    120            , 'vozoaeiu' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_eiv  = 'dyna_grid_V' ,    120            , 'vomeaeiv' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
-   sn_eiw  = 'dyna_grid_W' ,    120            , 'soleaeiw' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_tem  = 'dyna_grid_T' ,    120            , 'votemper'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_sal  = 'dyna_grid_T' ,    120            , 'vosaline'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_mld  = 'dyna_grid_T' ,    120            , 'somixhgt'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_emp  = 'dyna_grid_T' ,    120            , 'sowaflup'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_fmf  = 'dyna_grid_T' ,    120            , 'iowaflup'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_ice  = 'dyna_grid_T' ,    120            , 'soicecov'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_qsr  = 'dyna_grid_T' ,    120            , 'soshfldo'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_wnd  = 'dyna_grid_T' ,    120            , 'sowindsp'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_uwd  = 'dyna_grid_U' ,    120            , 'uocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_vwd  = 'dyna_grid_V' ,    120            , 'vocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_wwd  = 'dyna_grid_W' ,    120            , 'wocetr_eff' ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_avt  = 'dyna_grid_W' ,    120            , 'voddmavs'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_ubl  = 'dyna_grid_U' ,    120            , 'sobblcox'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
+   sn_vbl  = 'dyna_grid_V' ,    120            , 'sobblcoy'   ,  .true.    , .true. ,   'yearly'  , ''       , ''    , ''
 !
-   cn_dir      = './'      !  root directory for the location of the dynamical files
-   ln_degrad   =  .false.  !  flag for degradation -                requires ("key_degrad")
-   ln_dynwzv   =  .true.   !  computation of vertical velocity instead of using the one read in file
-   ln_dynbbl   =  .true.   !  bbl coef are in files, so read them - requires ("key_trabbl")
+   cn_dir          = './'       !  root directory for the location of the dynamical files
+   ln_dynrnf       =  .false.   !  runoffs option enabled (T) or not (F)
+   ln_dynrnf_depth =  .false.   ! runoffs is spread in vertical (T) or not (F)
+!   fwbcorr      = 3.786e-06    ! annual global mean of empmr for ssh correction
 /
 !-----------------------------------------------------------------------
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 7522)
@@ -17,4 +17,5 @@
    USE lib_mpp         ! distributed memory computing library
 
+   USE iom
    USE domstp          ! domain: set the time-step
 
@@ -73,10 +74,8 @@
 
       CALL dom_nam      ! read namelist ( namrun, namdom, namcla )
+      CALL dom_msk      ! Masks
+      CALL dom_hgr      ! Horizontal grid
       CALL dom_zgr      ! Vertical mesh and bathymetry option
-      CALL dom_grd      ! Create a domain file
-
-     !
-      ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines
-      !        but could be usefull in many other routines
+      !
       e12t    (:,:) = e1t(:,:) * e2t(:,:)
       e1e2t   (:,:) = e1t(:,:) * e2t(:,:)
@@ -91,16 +90,5 @@
       re1v_e2v(:,:) = e1v(:,:) / e2v(:,:)
       !
-      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points
-      hv(:,:) = 0._wp
-      DO jk = 1, jpk
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
-
       CALL dom_stp      ! Time step
-      CALL dom_msk      ! Masks
       CALL dom_ctl      ! Domain control
 
@@ -178,6 +166,4 @@
       nstocklist = nn_stocklist
       nwrite = nn_write
-
-
       !                             ! control of output frequency
       IF ( nstock == 0 .OR. nstock > nitend ) THEN
@@ -222,4 +208,5 @@
 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namdom )
+
 
       IF(lwp) THEN
@@ -321,10 +308,259 @@
    END SUBROUTINE dom_nam
 
+   SUBROUTINE dom_msk
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk  ***
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean mask informations and defines the interior domain T-mask.
+      !!
+      !! ** Method  :  Read in a file all the arrays generated in routines
+      !!               dommsk:   'mask.nc' file
+      !!              The interior ocean/land mask is computed from tmask
+      !!              setting to zero the duplicated row and lines due to
+      !!              MPP exchange halos, est-west cyclic and north fold
+      !!              boundary conditions.
+      !!
+      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
+      !!               tpol     : ???
+      !!----------------------------------------------------------------------
+      !
+      INTEGER  ::  inum   ! local integers
+      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
+      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
+      !
+      !!---------------------------------------------------------------------
+      
+
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      CALL wrk_alloc( jpi, jpj, zmbk )
+      zmbk(:,:) = 0._wp
+
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mask', inum )
+
+         !                                                         ! masks (inum2) 
+      CALL iom_get( inum, jpdom_data, 'tmask', tmask )
+      CALL iom_get( inum, jpdom_data, 'umask', umask )
+      CALL iom_get( inum, jpdom_data, 'vmask', vmask )
+      CALL iom_get( inum, jpdom_data, 'fmask', fmask )
+
+      CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
+      CALL lbc_lnk( umask, 'U', 1._wp )      
+      CALL lbc_lnk( vmask, 'V', 1._wp )
+      CALL lbc_lnk( fmask, 'F', 1._wp )
+
+#if defined key_c1d
+      ! set umask and vmask equal tmask in 1D configuration
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
+      IF(lwp) WRITE(numout,*) '**********                                                     ********'
+
+      umask(:,:,:) = tmask(:,:,:)
+      vmask(:,:,:) = tmask(:,:,:)
+#endif
+
+#if defined key_degrad
+      CALL iom_get( inum, jpdom_data, 'facvolt', facvol )
+#endif
+
+      CALL iom_get( inum, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
+      mbathy (:,:) = INT( zmbk(:,:) )
+      misfdep(:,:) = 1                                               ! ice shelf case not yet done
+      
+      CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
+
+      !                                     ! ============================
+      !                                     !        close the files 
+      !                                     ! ============================
+
+      !
+      ! Interior domain mask (used for global sum)
+      ! --------------------
+      ssmask(:,:)  = tmask(:,:,1)
+      tmask_i(:,:) = tmask(:,:,1)
+      iif = jpreci                        ! thickness of exchange halos in i-axis
+      iil = nlci - jpreci + 1
+      ijf = jprecj                        ! thickness of exchange halos in j-axis
+      ijl = nlcj - jprecj + 1
+      !
+      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
+      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
+      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
+      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
+      !
+      !                                   ! north fold mask
+      tpol(1:jpiglo) = 1._wp
+      !                                
+      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
+      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
+      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
+         IF( mjg(ijl-1) == jpjglo-1 ) THEN
+            DO ji = iif+1, iil-1
+               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
+            END DO
+         ENDIF
+      ENDIF 
+      !
+      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
+      ! least 1 wet u point
+      DO jj = 1, jpjm1
+         DO ji = 1, fs_jpim1   ! vector loop
+            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
+            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
+         END DO
+         DO ji = 1, jpim1      ! NO vector opt.
+            fmask_i(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,:)))
+         END DO
+      END DO
+      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
+      CALL lbc_lnk( vmask_i, 'V', 1._wp )
+      CALL lbc_lnk( fmask_i, 'F', 1._wp )
+
+      ! 3. Ocean/land mask at wu-, wv- and w points 
+      !----------------------------------------------
+      wmask (:,:,1) = tmask(:,:,1) ! ????????
+      wumask(:,:,1) = umask(:,:,1) ! ????????
+      wvmask(:,:,1) = vmask(:,:,1) ! ????????
+      DO jk = 2, jpk
+         wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1)
+         wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1)   
+         wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1)
+      END DO
+      !
+      CALL wrk_dealloc( jpi, jpj, zmbk )
+      !
+      CALL iom_close( inum )
+      !
+   END SUBROUTINE dom_msk
+
+   SUBROUTINE zgr_bot_level
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE zgr_bot_level  ***
+      !!
+      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
+      !!
+      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
+      !!
+      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
+      !!                                     ocean level at t-, u- & v-points
+      !!                                     (min value = 1 over land)
+      !!----------------------------------------------------------------------
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
+      !!----------------------------------------------------------------------
+
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
+      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
+      !
+      CALL wrk_alloc( jpi, jpj, zmbk )
+      !
+      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
+      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
+      !                                     ! bottom k-index of W-level = mbkt+1
+      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
+         DO ji = 1, jpim1
+            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
+            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
+         END DO
+      END DO
+      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
+      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
+      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
+      !
+      CALL wrk_dealloc( jpi, jpj, zmbk )
+      !
+   END SUBROUTINE zgr_bot_level
+
+   SUBROUTINE dom_hgr
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dom_hgr  ***
+      !!                   
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean horizontal mesh informations 
+      !!
+      !! ** Method  :   Read in a file all the arrays generated in routines
+      !!                domhgr:   'mesh_hgr.nc' file
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER ::   ji, jj   ! dummy loop indices
+      INTEGER  ::  inum    ! local integers
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)'
+      IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mesh_hgr', inum )
+
+      !                                                         ! horizontal mesh (inum3)
+      CALL iom_get( inum, jpdom_data, 'glamt', glamt )
+      CALL iom_get( inum, jpdom_data, 'glamu', glamu )
+      CALL iom_get( inum, jpdom_data, 'glamv', glamv )
+      CALL iom_get( inum, jpdom_data, 'glamf', glamf )
+
+      CALL iom_get( inum, jpdom_data, 'gphit', gphit )
+      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu )
+      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv )
+      CALL iom_get( inum, jpdom_data, 'gphif', gphif )
+
+      CALL iom_get( inum, jpdom_data, 'e1t', e1t )
+      CALL iom_get( inum, jpdom_data, 'e1u', e1u )
+      CALL iom_get( inum, jpdom_data, 'e1v', e1v )
+      
+      CALL iom_get( inum, jpdom_data, 'e2t', e2t )
+      CALL iom_get( inum, jpdom_data, 'e2u', e2u )
+      CALL iom_get( inum, jpdom_data, 'e2v', e2v )
+
+      CALL iom_get( inum, jpdom_data, 'ff', ff )
+
+
+      ! Control printing : Grid informations (if not restart)
+      ! ----------------
+
+      IF(lwp .AND. .NOT.ln_rstart ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '          longitude and e1 scale factors'
+         WRITE(numout,*) '          ------------------------------'
+         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
+            glamv(ji,1), glamf(ji,1),   &
+            e1t(ji,1), e1u(ji,1),   &
+            e1v(ji,1), ji = 1, jpi,10)
+
+         WRITE(numout,*)
+         WRITE(numout,*) '          latitude and e2 scale factors'
+         WRITE(numout,*) '          -----------------------------'
+         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
+            &                     gphiv(1,jj), gphif(1,jj),   &
+            &                     e2t  (1,jj), e2u  (1,jj),   &
+            &                     e2v  (1,jj), jj = 1, jpj, 10 )
+      ENDIF
+
+      !                                     ! ============================
+      !                                     !        close the files 
+      !                                     ! ============================
+      CALL iom_close( inum )
+      !
+9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
+            f19.10, 1x, f19.10, 1x, f19.10 )
+   END SUBROUTINE dom_hgr
+
+
    SUBROUTINE dom_zgr
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_zgr  ***
       !!                   
-      !! ** Purpose :  set the depth of model levels and the resulting 
-      !!      vertical scale factors.
+      !! ** Purpose :  Read the NetCDF file(s) which contain(s) all the
+      !!      ocean horizontal mesh informations and/or set the depth of model levels 
+      !!      and the resulting vertical scale factors.
       !!
       !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d)
@@ -338,8 +574,11 @@
       !! ** Action  :   define gdep., e3., mbathy and bathy
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio = 0   ! temporary integer
-      INTEGER ::   ios
+      INTEGER  ::  ioptio = 0   ! temporary integer
+      INTEGER  ::  inum, ios
+      INTEGER  ::  ji, jj, jk, ik
+      REAL(wp) ::  zrefdep
       !!
       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
+      REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw
       !!----------------------------------------------------------------------
 
@@ -372,234 +611,74 @@
       IF ( ioptio == 33 )   CALL ctl_stop( ' isf cavity with off line module not yet done    ' )
 
-   END SUBROUTINE dom_zgr
-
-   SUBROUTINE dom_ctl
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE dom_ctl  ***
-      !!
-      !! ** Purpose :   Domain control.
-      !!
-      !! ** Method  :   compute and print extrema of masked scale factors
-      !!
-      !! History :
-      !!   8.5  !  02-08  (G. Madec)    Original code
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
-      INTEGER, DIMENSION(2) ::   iloc      ! 
-      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
-      !!----------------------------------------------------------------------
-
-      ! Extrema of the scale factors
-
-      IF(lwp)WRITE(numout,*)
-      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
-      IF(lwp)WRITE(numout,*) '~~~~~~~'
-
-      IF (lk_mpp) THEN
-         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
-         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
-         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
-         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
+      CALL iom_open( 'mesh_zgr', inum )
+
+      CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
+      CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d )
+      IF( ln_zco .OR. ln_zps ) THEN
+         CALL iom_get( inum, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
+         CALL iom_get( inum, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
+      ENDIF
+
+!!gm BUG in s-coordinate this does not work!
+      ! deepest/shallowest W level Above/Below ~10m
+      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
+      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
+      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
+!!gm end bug
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
+         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
+         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
+      ENDIF
+
+      DO jk = 1, jpk
+         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
+         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
+      END DO
+
+      IF( lk_vvl ) THEN
+          CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) )
+          CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) )
+          ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at  T-points
+          DO jk = 1, jpk
+             ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
+          END DO
       ELSE
-         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
-
-         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi1 = iloc(1) + nimpp - 1
-         ijmi1 = iloc(2) + njmpp - 1
-         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iimi2 = iloc(1) + nimpp - 1
-         ijmi2 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima1 = iloc(1) + nimpp - 1
-         ijma1 = iloc(2) + njmpp - 1
-         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
-         iima2 = iloc(1) + nimpp - 1
-         ijma2 = iloc(2) + njmpp - 1
-      ENDIF
-
-      IF(lwp) THEN
-         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 dom_grd
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_grd  ***
-      !!                   
-      !! ** Purpose :  Read 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.
-      !!
-      !! ** Method  :   Read 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)
-      !!                    nmsh = 1  :   'mesh_mask.nc' file
-      !!                         = 2  :   'mesh.nc' and mask.nc' files
-      !!                         = 3  :   'mesh_hgr.nc', 'mesh_zgr.nc' and
-      !!                                  'mask.nc' files
-      !!      For huge size domain, use option 2 or 3 depending on your 
-      !!      vertical coordinate.
-      !!
-      !! ** input file : 
-      !!      meshmask.nc  : domain size, horizontal grid-point position,
-      !!                     masks, depth and vertical scale factors
-      !!----------------------------------------------------------------------
-      USE iom
-      !!
-      INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
-      REAL(wp) ::   zrefdep         ! local real
-      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw
-      !!----------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
-      IF(lwp) WRITE(numout,*) '~~~~~~~'
-
-      CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw )
-
-      zmbk(:,:) = 0._wp
-
-      SELECT CASE (nmsh)
-         !                                     ! ============================
-         CASE ( 1 )                            !  create 'mesh_mask.nc' file
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          one file in "mesh_mask.nc" '
-            CALL iom_open( 'mesh_mask', inum0 )
-
-            inum2 = inum0                                            ! put all the informations
-            inum3 = inum0                                            ! in unit inum0
-            inum4 = inum0
-
-            !                                  ! ============================
-         CASE ( 2 )                            !  create 'mesh.nc' and 
-            !                                  !         'mask.nc' files
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          two files in "mesh.nc" and "mask.nc" '
-            CALL iom_open( 'mesh', inum1 )
-            CALL iom_open( 'mask', inum2 )
-
-            inum3 = inum1                                            ! put mesh informations 
-            inum4 = inum1                                            ! in unit inum1 
-
-            !                                  ! ============================
-         CASE ( 3 )                            !  create 'mesh_hgr.nc'
-            !                                  !         'mesh_zgr.nc' and
-            !                                  !         'mask.nc'     files
-            !                                  ! ============================
-
-            IF(lwp) WRITE(numout,*) '          three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" '
-            CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc'
-            CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc'
-            CALL iom_open( 'mask'    , inum2 ) ! create 'mask.nc'
-
-            !                                  ! ===========================
-         CASE DEFAULT                          ! return error 
-            !                                  ! mesh has to be provided
-            !                                  ! ===========================
-            CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ',   &
-            &                                 ' Invalid nn_msh value in the namelist (0 is not allowed)' )
-
-         END SELECT
-
-         !                                                         ! masks (inum2) 
-         CALL iom_get( inum2, jpdom_data, 'tmask', tmask )
-         CALL iom_get( inum2, jpdom_data, 'umask', umask )
-         CALL iom_get( inum2, jpdom_data, 'vmask', vmask )
-         CALL iom_get( inum2, jpdom_data, 'fmask', fmask )
-
-         CALL lbc_lnk( tmask, 'T', 1._wp )    ! Lateral boundary conditions
-         CALL lbc_lnk( umask, 'U', 1._wp )      
-         CALL lbc_lnk( vmask, 'V', 1._wp )
-         CALL lbc_lnk( fmask, 'F', 1._wp )
-
-#if defined key_c1d
-         ! set umask and vmask equal tmask in 1D configuration
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '**********  1D configuration : set umask and vmask equal tmask ********'
-         IF(lwp) WRITE(numout,*) '**********                                                     ********'
-
-         umask(:,:,:) = tmask(:,:,:)
-         vmask(:,:,:) = tmask(:,:,:)
-#endif
-
-#if defined key_degrad
-         CALL iom_get( inum2, jpdom_data, 'facvolt', facvol )
-#endif
-
-         !                                                         ! horizontal mesh (inum3)
-         CALL iom_get( inum3, jpdom_data, 'glamt', glamt )
-         CALL iom_get( inum3, jpdom_data, 'glamu', glamu )
-         CALL iom_get( inum3, jpdom_data, 'glamv', glamv )
-         CALL iom_get( inum3, jpdom_data, 'glamf', glamf )
-
-         CALL iom_get( inum3, jpdom_data, 'gphit', gphit )
-         CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu )
-         CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv )
-         CALL iom_get( inum3, jpdom_data, 'gphif', gphif )
-
-         CALL iom_get( inum3, jpdom_data, 'e1t', e1t )
-         CALL iom_get( inum3, jpdom_data, 'e1u', e1u )
-         CALL iom_get( inum3, jpdom_data, 'e1v', e1v )
-         
-         CALL iom_get( inum3, jpdom_data, 'e2t', e2t )
-         CALL iom_get( inum3, jpdom_data, 'e2u', e2u )
-         CALL iom_get( inum3, jpdom_data, 'e2v', e2v )
-
-         CALL iom_get( inum3, jpdom_data, 'ff', ff )
-
-         CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk )              ! number of ocean t-points
-         mbathy (:,:) = INT( zmbk(:,:) )
-         misfdep(:,:) = 1                                               ! ice shelf case not yet done
-         
-         CALL zgr_bot_level                                             ! mbk. arrays (deepest ocean t-, u- & v-points
-         !
          IF( ln_sco ) THEN                                         ! s-coordinate
-            CALL iom_get( inum4, jpdom_data, 'hbatt', hbatt )
-            CALL iom_get( inum4, jpdom_data, 'hbatu', hbatu )
-            CALL iom_get( inum4, jpdom_data, 'hbatv', hbatv )
-            CALL iom_get( inum4, jpdom_data, 'hbatf', hbatf )
+            CALL iom_get( inum, jpdom_data, 'hbatt', hbatt )
+            CALL iom_get( inum, jpdom_data, 'hbatu', hbatu )
+            CALL iom_get( inum, jpdom_data, 'hbatv', hbatv )
+            CALL iom_get( inum, jpdom_data, 'hbatf', hbatf )
             
-            CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
-            CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw )
-            CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w ) 
-            CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt )
-            CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw )
-
-            CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
-            CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
-            CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
-            CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
-
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
+            CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef.
+            CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw )
+            CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w ) 
+            CALL iom_get( inum, jpdom_unknown, 'esigt', esigt )
+            CALL iom_get( inum, jpdom_unknown, 'esigw', esigw )
+
+            CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors
+            CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
+            CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
+            CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
          ENDIF
 
  
          IF( ln_zps ) THEN                                           ! z-coordinate - partial steps
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d )  ! reference depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
-            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )    ! reference scale factors
-            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
             !
-            IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
-               CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
+            IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )
             ELSE                                                        ! 2D bottom scale factors
-               CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp )
-               CALL iom_get( inum4, jpdom_data, 'e3w_ps', e3wp )
+               CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp )
+               CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp )
                !                                                        ! deduces the 3D scale factors
                DO jk = 1, jpk
@@ -633,10 +712,12 @@
             END IF
 
-            IF( iom_varid( inum4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
-               CALL iom_get( inum4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
-               CALL iom_get( inum4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
+            IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
+               CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )
+               CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )
             ELSE                                                           ! 2D bottom depth
-               CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
-               CALL iom_get( inum4, jpdom_data, 'hdepw', zprw )
+               CALL wrk_alloc( jpi, jpj, zprt, zprw )
+               !
+               CALL iom_get( inum, jpdom_data, 'hdept', zprt )
+               CALL iom_get( inum, jpdom_data, 'hdepw', zprw )
                !
                DO jk = 1, jpk                                              ! deduces the 3D depth
@@ -654,4 +735,5 @@
                   END DO
                END DO
+               CALL wrk_dealloc( jpi, jpj, zprt, zprw )
             ENDIF
             !
@@ -659,8 +741,4 @@
 
          IF( ln_zco ) THEN           ! Vertical coordinates and scales factors
-            CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
-            CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )
-            CALL iom_get( inum4, jpdom_unknown, 'e3t_1d'  , e3t_1d   )
-            CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
             DO jk = 1, jpk
                fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
@@ -672,212 +750,70 @@
             END DO
          ENDIF
-
-!!gm BUG in s-coordinate this does not work!
-      ! deepest/shallowest W level Above/Below ~10m
-      zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) )                 ! ref. depth with tolerance (10% of minimum layer thickness)
-      nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
-      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m
-!!gm end bug
-
-      ! Control printing : Grid informations (if not restart)
-      ! ----------------
-
-      IF(lwp .AND. .NOT.ln_rstart ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '          longitude and e1 scale factors'
-         WRITE(numout,*) '          ------------------------------'
-         WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1),   &
-            glamv(ji,1), glamf(ji,1),   &
-            e1t(ji,1), e1u(ji,1),   &
-            e1v(ji,1), ji = 1, jpi,10)
-9300     FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x,    &
-            f19.10, 1x, f19.10, 1x, f19.10 )
-
-         WRITE(numout,*)
-         WRITE(numout,*) '          latitude and e2 scale factors'
-         WRITE(numout,*) '          -----------------------------'
-         WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj),   &
-            &                     gphiv(1,jj), gphif(1,jj),   &
-            &                     e2t  (1,jj), e2u  (1,jj),   &
-            &                     e2v  (1,jj), jj = 1, jpj, 10 )
-      ENDIF
-
-
-      IF( nprint == 1 .AND. lwp ) THEN
-         WRITE(numout,*) '          e1u e2u '
-         CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         WRITE(numout,*) '          e1v e2v  '
-         CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-         CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout )
-      ENDIF
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:'
-         WRITE(numout, "(9x,' level   gdept    gdepw     e3t      e3w  ')" )
-         WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
-      ENDIF
-
-      DO jk = 1, jpk
-         IF( e3w_1d  (jk) <= 0._wp .OR. e3t_1d  (jk) <= 0._wp )   CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' )
-         IF( gdepw_1d(jk) <  0._wp .OR. gdept_1d(jk) <  0._wp )   CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' )
-      END DO
+         !
+      ENDIF
       !                                     ! ============================
       !                                     !        close the files 
       !                                     ! ============================
-      SELECT CASE ( nmsh )
-         CASE ( 1 )                
-            CALL iom_close( inum0 )
-         CASE ( 2 )
-            CALL iom_close( inum1 )
-            CALL iom_close( inum2 )
-         CASE ( 3 )
-            CALL iom_close( inum2 )
-            CALL iom_close( inum3 )
-            CALL iom_close( inum4 )
-      END SELECT
-      !
-      CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw )
-      !
-   END SUBROUTINE dom_grd
-
-
-   SUBROUTINE zgr_bot_level
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE zgr_bot_level  ***
-      !!
-      !! ** Purpose :   defines the vertical index of ocean bottom (mbk. arrays)
-      !!
-      !! ** Method  :   computes from mbathy with a minimum value of 1 over land
-      !!
-      !! ** Action  :   mbkt, mbku, mbkv :   vertical indices of the deeptest
-      !!                                     ocean level at t-, u- & v-points
-      !!                                     (min value = 1 over land)
-      !!----------------------------------------------------------------------
-      !
-      INTEGER ::   ji, jj   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
-      !!----------------------------------------------------------------------
-
-      !
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
-      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
-      !
-      CALL wrk_alloc( jpi, jpj, zmbk )
-      !
-      mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
-      mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf)
-      !                                     ! bottom k-index of W-level = mbkt+1
-      DO jj = 1, jpjm1                      ! bottom k-index of u- (v-) level
-         DO ji = 1, jpim1
-            mbku(ji,jj) = MIN(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )
-            mbkv(ji,jj) = MIN(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
-         END DO
-      END DO
-      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
-      zmbk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
-      zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
-      !
-      CALL wrk_dealloc( jpi, jpj, zmbk )
-      !
-   END SUBROUTINE zgr_bot_level
-
-   SUBROUTINE dom_msk
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE dom_msk  ***
-      !!
-      !! ** Purpose :   Off-line case: defines the interior domain T-mask.
-      !!
-      !! ** Method  :   The interior ocean/land mask is computed from tmask
-      !!              setting to zero the duplicated row and lines due to
-      !!              MPP exchange halos, est-west cyclic and north fold
-      !!              boundary conditions.
-      !!
-      !! ** Action :   tmask_i  : interiorland/ocean mask at t-point
-      !!               tpol     : ???
-      !!----------------------------------------------------------------------
-      !
-      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
-      INTEGER  ::   iif, iil, ijf, ijl       ! local integers
-      INTEGER, POINTER, DIMENSION(:,:) ::  imsk 
-      !
-      !!---------------------------------------------------------------------
-      
-      CALL wrk_alloc( jpi, jpj, imsk )
-      !
-      ! Interior domain mask (used for global sum)
-      ! --------------------
-      ssmask(:,:)  = tmask(:,:,1)
-      tmask_i(:,:) = tmask(:,:,1)
-      iif = jpreci                        ! thickness of exchange halos in i-axis
-      iil = nlci - jpreci + 1
-      ijf = jprecj                        ! thickness of exchange halos in j-axis
-      ijl = nlcj - jprecj + 1
-      !
-      tmask_i( 1 :iif,   :   ) = 0._wp    ! first columns
-      tmask_i(iil:jpi,   :   ) = 0._wp    ! last  columns (including mpp extra columns)
-      tmask_i(   :   , 1 :ijf) = 0._wp    ! first rows
-      tmask_i(   :   ,ijl:jpj) = 0._wp    ! last  rows (including mpp extra rows)
-      !
-      !                                   ! north fold mask
-      tpol(1:jpiglo) = 1._wp
-      !                                
-      IF( jperio == 3 .OR. jperio == 4 )   tpol(jpiglo/2+1:jpiglo) = 0._wp    ! T-point pivot
-      IF( jperio == 5 .OR. jperio == 6 )   tpol(     1    :jpiglo) = 0._wp    ! F-point pivot
-      IF( jperio == 3 .OR. jperio == 4 ) THEN      ! T-point pivot: only half of the nlcj-1 row
-         IF( mjg(ijl-1) == jpjglo-1 ) THEN
-            DO ji = iif+1, iil-1
-               tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji))
-            END DO
-         ENDIF
-      ENDIF 
-      !
-      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at
-      ! least 1 wet u point
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector loop
-            umask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:)))
-            vmask_i(ji,jj)  = ssmask(ji,jj) * ssmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:)))
-         END DO
-         DO ji = 1, jpim1      ! NO vector opt.
-            fmask_i(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,:)))
-         END DO
-      END DO
-      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions
-      CALL lbc_lnk( vmask_i, 'V', 1._wp )
-      CALL lbc_lnk( fmask_i, 'F', 1._wp )
-
-      ! 3. Ocean/land mask at wu-, wv- and w points 
-      !----------------------------------------------
-      wmask (:,:,1) = tmask(:,:,1) ! ????????
-      wumask(:,:,1) = umask(:,:,1) ! ????????
-      wvmask(:,:,1) = vmask(:,:,1) ! ????????
-      DO jk=2,jpk
-         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1)
-         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)   
-         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1)
-      END DO
-      !
-      IF( nprint == 1 .AND. lwp ) THEN    ! Control print
-         imsk(:,:) = INT( tmask_i(:,:) )
-         WRITE(numout,*) ' tmask_i : '
-         CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         WRITE (numout,*)
-         WRITE (numout,*) ' dommsk: tmask for each level'
-         WRITE (numout,*) ' ----------------------------'
-         DO jk = 1, jpk
-            imsk(:,:) = INT( tmask(:,:,jk) )
-            WRITE(numout,*)
-            WRITE(numout,*) ' level = ',jk
-            CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout)
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, imsk )
-      !
-   END SUBROUTINE dom_msk
+      CALL iom_close( inum )
+      !
+      !
+   END SUBROUTINE dom_zgr
+
+   SUBROUTINE dom_ctl
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dom_ctl  ***
+      !!
+      !! ** Purpose :   Domain control.
+      !!
+      !! ** Method  :   compute and print extrema of masked scale factors
+      !!
+      !! History :
+      !!   8.5  !  02-08  (G. Madec)    Original code
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2
+      INTEGER, DIMENSION(2) ::   iloc      ! 
+      REAL(wp) ::   ze1min, ze1max, ze2min, ze2max
+      !!----------------------------------------------------------------------
+
+      ! Extrema of the scale factors
+
+      IF(lwp)WRITE(numout,*)
+      IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors'
+      IF(lwp)WRITE(numout,*) '~~~~~~~'
+
+      IF (lk_mpp) THEN
+         CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 )
+         CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 )
+         CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 )
+         CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 )
+      ELSE
+         ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+         ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )    
+
+         iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi1 = iloc(1) + nimpp - 1
+         ijmi1 = iloc(2) + njmpp - 1
+         iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iimi2 = iloc(1) + nimpp - 1
+         ijmi2 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima1 = iloc(1) + nimpp - 1
+         ijma1 = iloc(2) + njmpp - 1
+         iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )
+         iima2 = iloc(1) + nimpp - 1
+         ijma2 = iloc(2) + njmpp - 1
+      ENDIF
+
+      IF(lwp) THEN
+         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
 
    !!======================================================================
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 7522)
@@ -22,4 +22,5 @@
    USE c1d             ! 1D configuration: lk_c1d
    USE dom_oce         ! ocean domain: variables
+   USE domvvl          ! variable volume
    USE zdf_oce         ! ocean vertical physics: variables
    USE sbc_oce         ! surface module: variables
@@ -28,4 +29,5 @@
    USE trabbl          ! active tracer: bottom boundary layer
    USE ldfslp          ! lateral diffusion: iso-neutral slopes
+   USE sbcrnf          ! river runoffs
    USE ldfeiv          ! eddy induced velocity coef. 
    USE ldftra_oce      ! ocean tracer   lateral physics
@@ -39,5 +41,7 @@
    USE prtctl          ! print control
    USE fldread         ! read input fields 
+   USE wrk_nemo        ! Memory allocation 
    USE timing          ! Timing
+   USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc
 
    IMPLICIT NONE
@@ -46,48 +50,40 @@
    PUBLIC   dta_dyn_init   ! called by opa.F90
    PUBLIC   dta_dyn        ! called by step.F90
-
-   CHARACTER(len=100) ::   cn_dir       !: Root directory for location of ssr files
-   LOGICAL            ::   ln_dynwzv    !: vertical velocity read in a file (T) or computed from u/v (F)
-   LOGICAL            ::   ln_dynbbl    !: bbl coef read in a file (T) or computed (F)
-   LOGICAL            ::   ln_degrad    !: degradation option enabled or not
-   LOGICAL            ::   ln_dynrnf    !: read runoff data in file (T) or set to zero (F)
-
-   INTEGER  , PARAMETER ::   jpfld = 21     ! maximum number of fields to read
+   PUBLIC   dta_dyn_swp   ! called by step.F90
+
+   CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files
+   LOGICAL            ::   ln_dynrnf       !: read runoff data in file (T) or set to zero (F)
+   LOGICAL            ::   ln_dynrnf_depth       !: read runoff data in file (T) or set to zero (F)
+   REAL(wp)           ::   fwbcorr
+
+
+   INTEGER  , PARAMETER ::   jpfld = 20     ! maximum number of fields to read
    INTEGER  , SAVE      ::   jf_tem         ! index of temperature
    INTEGER  , SAVE      ::   jf_sal         ! index of salinity
-   INTEGER  , SAVE      ::   jf_uwd         ! index of u-wind
-   INTEGER  , SAVE      ::   jf_vwd         ! index of v-wind
-   INTEGER  , SAVE      ::   jf_wwd         ! index of w-wind
+   INTEGER  , SAVE      ::   jf_uwd         ! index of u-transport
+   INTEGER  , SAVE      ::   jf_vwd         ! index of v-transport
+   INTEGER  , SAVE      ::   jf_wwd         ! index of v-transport
    INTEGER  , SAVE      ::   jf_avt         ! index of Kz
    INTEGER  , SAVE      ::   jf_mld         ! index of mixed layer deptht
    INTEGER  , SAVE      ::   jf_emp         ! index of water flux
+   INTEGER  , SAVE      ::   jf_empb        ! index of water flux
    INTEGER  , SAVE      ::   jf_qsr         ! index of solar radiation
    INTEGER  , SAVE      ::   jf_wnd         ! index of wind speed
    INTEGER  , SAVE      ::   jf_ice         ! index of sea ice cover
    INTEGER  , SAVE      ::   jf_rnf         ! index of river runoff
+   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
    INTEGER  , SAVE      ::   jf_ubl         ! index of u-bbl coef
    INTEGER  , SAVE      ::   jf_vbl         ! index of v-bbl coef
-   INTEGER  , SAVE      ::   jf_ahu         ! index of u-diffusivity coef
-   INTEGER  , SAVE      ::   jf_ahv         ! index of v-diffusivity coef 
-   INTEGER  , SAVE      ::   jf_ahw         ! index of w-diffusivity coef
-   INTEGER  , SAVE      ::   jf_eiu         ! index of u-eiv
-   INTEGER  , SAVE      ::   jf_eiv         ! index of v-eiv
-   INTEGER  , SAVE      ::   jf_eiw         ! index of w-eiv
-   INTEGER  , SAVE      ::   jf_fmf         ! index of downward salt flux
-
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
+   INTEGER  , SAVE      ::   jf_div         ! index of e3t
+
+
+   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn  ! structure of input fields (file informations, fields read)
    !                                               ! 
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at 2 time step
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:  ) :: wnow       ! vertical velocity at 2 time step
    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
    REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: uslpnow    ! zonal isopycnal slopes
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: vslpnow    ! meridional isopycnal slopes
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpinow   ! zonal diapycnal slopes
-   REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: wslpjnow   ! meridional diapycnal slopes
-
-   INTEGER :: nrecprev_tem , nrecprev_uwd
+
+   INTEGER, SAVE  :: nprevrec, nsecdyn
 
    !! * Substitutions
@@ -113,16 +109,9 @@
       !!----------------------------------------------------------------------
       !
-      USE oce, ONLY:  zts    => tsa 
-      USE oce, ONLY:  zuslp  => ua   , zvslp  => va
-      USE oce, ONLY:  zwslpi => rotb , zwslpj => rotn
-      USE oce, ONLY:  zu     => ub   , zv     => vb,  zw => hdivb
-      !
+      USE oce, ONLY:  zhdivtr => ua
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      !
-      INTEGER  ::   ji, jj     ! dummy loop indices
-      INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
-      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
-      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
-      INTEGER  ::   iswap_tem, iswap_uwd    ! 
+      INTEGER             ::   ji, jj, jk
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zemp
+      !
       !!----------------------------------------------------------------------
       
@@ -130,119 +119,49 @@
       IF( nn_timing == 1 )  CALL timing_start( 'dta_dyn')
       !
-      isecsbc = nsec_year + nsec1jan000 
-      !
-      IF( kt == nit000 ) THEN
-         nrecprev_tem = 0
-         nrecprev_uwd = 0
-         !
-         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
-         !
-         IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
-            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
-            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity 
-            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
-            CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
-            uslpdta (:,:,:,1) = zuslp (:,:,:) 
-            vslpdta (:,:,:,1) = zvslp (:,:,:) 
-            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
-            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
-         ENDIF
-         IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint )  THEN    ! compute vertical velocity from u/v
-            zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1)
-            zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1)
-            CALL dta_dyn_wzv( zu, zv, zw )
-            wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:)
-         ENDIF
-      ELSE
-         nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2)
-         nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2)
-         !
-         CALL fld_read( kt, 1, sf_dyn )      !==   read data at kt time step   ==!
-         !
-      ENDIF
-      ! 
-      IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN    ! Computes slopes (here avt is used as workspace)                       
-         iswap_tem = 0
-         IF(  kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 )  iswap_tem = 1
-         IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
-            IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
-            IF( sf_dyn(jf_tem)%ln_tint ) THEN                 ! time interpolation of data
-               IF( kt /= nit000 ) THEN
-                  uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
-                  vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)  
-                  wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
-                  wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
-               ENDIF
-               !
-               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
-               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity 
-               avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
-               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
-               !
-               uslpdta (:,:,:,2) = zuslp (:,:,:) 
-               vslpdta (:,:,:,2) = zvslp (:,:,:) 
-               wslpidta(:,:,:,2) = zwslpi(:,:,:) 
-               wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
-            ELSE
-               zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)
-               zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)
-               avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)
-               CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
-               uslpnow (:,:,:)   = zuslp (:,:,:) 
-               vslpnow (:,:,:)   = zvslp (:,:,:) 
-               wslpinow(:,:,:)   = zwslpi(:,:,:) 
-               wslpjnow(:,:,:)   = zwslpj(:,:,:) 
-            ENDIF
-         ENDIF
-         IF( sf_dyn(jf_tem)%ln_tint )  THEN
-            ztinta =  REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp )  &
-               &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
-            ztintb =  1. - ztinta
-            uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)  
-            vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)  
-            wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)  
-            wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)  
-         ELSE
-            uslp (:,:,:) = uslpnow (:,:,:)
-            vslp (:,:,:) = vslpnow (:,:,:)
-            wslpi(:,:,:) = wslpinow(:,:,:)
-            wslpj(:,:,:) = wslpjnow(:,:,:)
-         ENDIF
-      ENDIF
-      !
-      IF( ln_dynwzv )  THEN    ! compute vertical velocity from u/v
-         iswap_uwd = 0
-         IF(  kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 )  iswap_uwd = 1
-         IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 )  THEN    ! read/update the after data
-            IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt
-            IF(lwp) WRITE(numout,*)
-            IF( sf_dyn(jf_uwd)%ln_tint ) THEN                 ! time interpolation of data
-               IF( kt /= nit000 )  THEN
-                  wdta(:,:,:,1) =  wdta(:,:,:,2)     ! swap the data for initialisation
-               ENDIF
-               zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2)
-               zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2)
-               CALL dta_dyn_wzv( zu, zv, zw )
-               wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:)
-            ELSE
-               zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 
-               zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:)
-               CALL dta_dyn_wzv( zu, zv, zw )
-               wnow(:,:,:)  = zw(:,:,:) * tmask(:,:,:)
-            ENDIF
-         ENDIF
-         IF( sf_dyn(jf_uwd)%ln_tint )  THEN
-            ztinta =  REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp )  &
-               &    / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp )
-            ztintb =  1. - ztinta
-            wn(:,:,:) = ztintb * wdta(:,:,:,1)  + ztinta * wdta(:,:,:,2)  
-         ELSE
-            wn(:,:,:) = wnow(:,:,:)
-         ENDIF
-      ENDIF
+      !
+      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
+      !
+      IF( kt == nit000 ) THEN    ;    nprevrec = 0
+      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
+      ENDIF
+      !
+      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
+      !
+      IF( lk_ldfslp .AND. .NOT.lk_c1d )   CALL  dta_dyn_slp( kt )    ! Computation of slopes
       !
       tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
       tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
-      !
+      wndm(:,:)         = sf_dyn(jf_wnd)%fnow(:,:,1)  * tmask(:,:,1)    ! wind speed - needed for gas exchange
+      fmmflx(:,:)       = sf_dyn(jf_fmf)%fnow(:,:,1)  * tmask(:,:,1)    ! downward salt flux (v3.5+)
+      fr_i(:,:)         = sf_dyn(jf_ice)%fnow(:,:,1)  * tmask(:,:,1)    ! Sea-ice fraction
+      qsr (:,:)         = sf_dyn(jf_qsr)%fnow(:,:,1)  * tmask(:,:,1)    ! solar radiation
+      emp (:,:)         = sf_dyn(jf_emp)%fnow(:,:,1)  * tmask(:,:,1)    ! E-P
+      IF( ln_dynrnf ) THEN 
+         rnf (:,:)      = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
+         IF( ln_dynrnf_depth .AND. lk_vvl )    CALL  dta_dyn_hrnf
+      ENDIF
+      !
+      un(:,:,:)        = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:)    ! effective u-transport
+      vn(:,:,:)        = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:)    ! effective v-transport
+      wn(:,:,:)        = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    ! effective v-transport
+      !
+      IF( lk_vvl ) THEN
+         CALL wrk_alloc(jpi, jpj, zemp )
+         zhdivtr(:,:,:)    = sf_dyn(jf_div)%fnow(:,:,:)  * tmask(:,:,:)    ! effective u-transport
+         emp_b (:,:)       = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
+         zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1)
+         CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) )  !=  ssh, vertical scale factor & vertical transport
+         CALL wrk_dealloc(jpi, jpj, zemp )
+         !                                           Write in the tracer restart file
+         !                                          *******************************
+         IF( lrst_trc ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ',   &
+               &                    'at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*) '~~~~'
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha )
+            CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn )
+         ENDIF
+      ENDIF
       !
       CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
@@ -251,53 +170,16 @@
 
       rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl
-      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl 
-      !
-      avt(:,:,:)       = sf_dyn(jf_avt)%fnow(:,:,:)  * tmask(:,:,:)    ! vertical diffusive coefficient 
-      un (:,:,:)       = sf_dyn(jf_uwd)%fnow(:,:,:)  * umask(:,:,:)    ! u-velocity
-      vn (:,:,:)       = sf_dyn(jf_vwd)%fnow(:,:,:)  * vmask(:,:,:)    ! v-velocity 
-      IF( .NOT.ln_dynwzv ) &                                          ! w-velocity read in file 
-         wn (:,:,:)    = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:)    
-      hmld(:,:)        = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
-      wndm(:,:)        = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1)    ! wind speed - needed for gas exchange
-      emp (:,:)        = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1)    ! E-P
-      fmmflx(:,:)      = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1)    ! downward salt flux (v3.5+)
-      fr_i(:,:)        = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1)    ! Sea-ice fraction
-      qsr (:,:)        = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1)    ! solar radiation
-      IF( ln_dynrnf ) &
-      rnf (:,:)        = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1)    ! river runoffs 
-
-      !                                                      ! bbl diffusive coef
+      CALL zdf_mxl( kt )                                                   ! In any case, we need mxl
+      !
+      hmld(:,:)         = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1)    ! mixed layer depht
+      avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)    ! vertical diffusive coefficient 
+      !
 #if defined key_trabbl && ! defined key_c1d
-      IF( ln_dynbbl ) THEN                                        ! read in a file
-         ahu_bbl(:,:)  = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)
-         ahv_bbl(:,:)  = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
-      ELSE                                                        ! Compute bbl coefficients if needed
-         tsb(:,:,:,:) = tsn(:,:,:,:)
-         CALL bbl( kt, nit000, 'TRC')
-      END IF
+      ahu_bbl(:,:)      = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1)    ! bbl diffusive coef
+      ahv_bbl(:,:)      = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1)
 #endif
-#if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d 
-      aeiw(:,:)        = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1)    ! w-eiv
-      !                                                           ! Computes the horizontal values from the vertical value
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )  ! Average the diffusive coefficient at u- v- points
-            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )  ! at u- v- points
-         END DO
-      END DO
-      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
-#endif
-      
-#if defined key_degrad && ! defined key_c1d 
-      !                                          ! degrad option : diffusive and eiv coef are 3D
-      ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:)
-      ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:)
-      ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:)
-#  if defined key_traldf_eiv 
-      aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:)
-      aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:)
-      aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:)
-#  endif
-#endif
+      !
+      !
+      CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
       !
       IF(ln_ctl) THEN                  ! print control
@@ -308,10 +190,10 @@
          CALL prt_ctl(tab3d_1=wn               , clinfo1=' wn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
          CALL prt_ctl(tab3d_1=avt              , clinfo1=' kz      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
-         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
-         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask, ovlap=1 )
+!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask, ovlap=1 )
       ENDIF
       !
@@ -335,16 +217,24 @@
       INTEGER  :: inum, idv, idimv                   ! local integer
       INTEGER  :: ios                                ! Local integer output status for namelist read
-      !!
-      CHARACTER(len=100)            ::  cn_dir   !   Root directory for location of core files
-      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d    ! array of namelist informations on the fields to read
-      TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf  ! informations about the fields to be read
-      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl          !   "                                 "
-      TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf  !   "                                 "
-      !!----------------------------------------------------------------------
-      !
-      NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf,    &
-         &                sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf,  &
-         &                sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl,          &
-         &                sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf
+      INTEGER  :: ji, jj, jk
+      REAL(wp) :: zcoef
+      INTEGER  :: nkrnf_max
+      REAL(wp) :: hrnf_max
+      !!
+      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
+      TYPE(FLD_N), DIMENSION(jpfld) ::  slf_d         ! array of namelist informations on the fields to read
+      TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp  ! informations about the fields to be read
+      TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt   !   "                 "
+      TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf   !   "               "
+      TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf    !   "              "
+      TYPE(FLD_N) :: sn_div  ! informations about the fields to be read
+      !!----------------------------------------------------------------------
+
+      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
+         &                sn_uwd, sn_vwd, sn_wwd, sn_emp,    &
+         &                sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr ,   &
+         &                sn_wnd, sn_ice, sn_fmf,                    &
+         &                sn_ubl, sn_vbl, sn_rnf,                   &
+         &                sn_empb, sn_div 
       !
       REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
@@ -363,74 +253,49 @@
          WRITE(numout,*) '~~~~~~~ '
          WRITE(numout,*) '   Namelist namdta_dyn'
-         WRITE(numout,*) '      vertical velocity read from file (T) or computed (F) ln_dynwzv  = ', ln_dynwzv
-         WRITE(numout,*) '      bbl coef read from file (T) or computed (F)          ln_dynbbl  = ', ln_dynbbl
-         WRITE(numout,*) '      degradation option enabled (T) or not (F)            ln_degrad  = ', ln_degrad
-         WRITE(numout,*) '      river runoff option enabled (T) or not (F)           ln_dynrnf  = ', ln_dynrnf
+         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
+         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
+         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
          WRITE(numout,*)
       ENDIF
       ! 
-      IF( ln_degrad .AND. .NOT.lk_degrad ) THEN
-         CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' )
-         ln_degrad = .FALSE.
-      ENDIF
-      IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN
-         CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' )
-         ln_dynbbl = .FALSE.
-      ENDIF
-
-      jf_tem = 1   ;   jf_sal = 2   ;  jf_mld = 3   ;  jf_emp = 4   ;   jf_fmf  = 5   ;  jf_ice = 6   ;   jf_qsr = 7
-      jf_wnd = 8   ;   jf_uwd = 9   ;  jf_vwd = 10  ;  jf_wwd = 11  ;   jf_avt  = 12  ;  jfld  = jf_avt
-      !
-      slf_d(jf_tem) = sn_tem   ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
-      slf_d(jf_emp) = sn_emp   ;   slf_d(jf_fmf ) = sn_fmf   ;   slf_d(jf_ice) = sn_ice 
-      slf_d(jf_qsr) = sn_qsr   ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_avt) = sn_avt 
-      slf_d(jf_uwd) = sn_uwd   ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
-
+
+      jf_uwd  = 1     ;   jf_vwd  = 2    ;   jf_wwd = 3    ;   jf_emp = 4    ;   jf_avt = 5
+      jf_tem  = 6     ;   jf_sal  = 7    ;   jf_mld = 8    ;   jf_qsr = 9
+      jf_wnd  = 10    ;   jf_ice  = 11   ;   jf_fmf = 12   ;   jfld   = jf_fmf
+
+      !
+      slf_d(jf_uwd)  = sn_uwd    ;   slf_d(jf_vwd)  = sn_vwd   ;   slf_d(jf_wwd) = sn_wwd
+      slf_d(jf_emp)  = sn_emp    ;   slf_d(jf_avt)  = sn_avt
+      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal   ;   slf_d(jf_mld) = sn_mld
+      slf_d(jf_qsr)  = sn_qsr    ;   slf_d(jf_wnd)  = sn_wnd   ;   slf_d(jf_ice) = sn_ice
+      slf_d(jf_fmf)  = sn_fmf
+
+
+      !
+      IF( lk_vvl ) THEN
+                 jf_div  = jfld + 1    ;         jf_empb  = jfld + 2      ;      jfld = jf_empb
+           slf_d(jf_div) = sn_div      ;   slf_d(jf_empb) = sn_empb
+      ENDIF
+      !
+      IF( lk_trabbl ) THEN
+                 jf_ubl  = jfld + 1    ;         jf_vbl  = jfld + 2     ;      jfld = jf_vbl
+           slf_d(jf_ubl) = sn_ubl      ;   slf_d(jf_vbl) = sn_vbl
+      ENDIF
       !
       IF( ln_dynrnf ) THEN
-                jf_rnf = jfld + 1  ;  jfld  = jf_rnf
-         slf_d(jf_rnf) = sn_rnf
+                jf_rnf  = jfld + 1     ;     jfld  = jf_rnf
+          slf_d(jf_rnf) = sn_rnf    
       ELSE
-         rnf (:,:) = 0._wp
-      ENDIF
-
-      !
-      IF( .NOT.ln_degrad ) THEN     ! no degrad option
-         IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN        ! eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;        jf_eiw  = jfld + 3   ;   jfld = jf_eiw
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl  ;   slf_d(jf_eiw) = sn_eiw
-         ENDIF
-         IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN   ! no eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-         ENDIF
-         IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN   ! eiv & no bbl
-           jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw
-         ENDIF
-      ELSE
-              jf_ahu  = jfld + 1 ;        jf_ahv  = jfld + 2 ;        jf_ahw  = jfld + 3  ;  jfld = jf_ahw
-        slf_d(jf_ahu) = sn_ahu  ;   slf_d(jf_ahv) = sn_ahv  ;   slf_d(jf_ahw) = sn_ahw
-        IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN         ! eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-                 jf_eiu  = jfld + 3 ;        jf_eiv  = jfld + 4 ;    jf_eiw  = jfld + 5   ;  jfld = jf_eiw 
-           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;    slf_d(jf_eiw) = sn_eiw
-        ENDIF
-        IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN    ! no eiv & bbl
-                 jf_ubl  = jfld + 1 ;        jf_vbl  = jfld + 2 ;  jfld = jf_vbl
-           slf_d(jf_ubl) = sn_ubl  ;   slf_d(jf_vbl) = sn_vbl
-        ENDIF
-        IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN    ! eiv & no bbl
-                 jf_eiu  = jfld + 1 ;         jf_eiv  = jfld + 2 ;    jf_eiw  = jfld + 3   ; jfld = jf_eiw 
-           slf_d(jf_eiu) = sn_eiu  ;   slf_d(jf_eiv) = sn_eiv  ;   slf_d(jf_eiw) = sn_eiw
-        ENDIF
-      ENDIF
+         rnf(:,:) = 0._wp
+      ENDIF
+
   
       ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
-      IF( ierr > 0 ) THEN
+      IF( ierr > 0 )  THEN
          CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
       ENDIF
       !                                         ! fill sf with slf_i and control print
       CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
+      !
       ! Open file for each variable to get his number of dimension
       DO ifpr = 1, jfld
@@ -456,75 +321,353 @@
             ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
             &         wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 )
-         ELSE
-            ALLOCATE( uslpnow (jpi,jpj,jpk)  , vslpnow (jpi,jpj,jpk)  ,    &
-            &         wslpinow(jpi,jpj,jpk)  , wslpjnow(jpi,jpj,jpk)  , STAT=ierr2 )
-         ENDIF 
-         IF( ierr2 > 0 ) THEN
-            CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
+            !
+            IF( ierr2 > 0 )  THEN
+               CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' )   ;   RETURN
+            ENDIF
          ENDIF
       ENDIF
-      IF( ln_dynwzv ) THEN                  ! slopes 
-         IF( sf_dyn(jf_uwd)%ln_tint ) THEN      ! time interpolation
-            ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 )
-         ELSE
-            ALLOCATE( wnow(jpi,jpj,jpk)  , STAT=ierr3 )
-         ENDIF 
-         IF( ierr3 > 0 ) THEN
-            CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' )   ;   RETURN
-         ENDIF
-      ENDIF
-      !
-      CALL dta_dyn( nit000 )
-      !
-   END SUBROUTINE dta_dyn_init
-
-   SUBROUTINE dta_dyn_wzv( pu, pv, pw )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE wzv  ***
-      !!
-      !! ** Purpose :   Compute the now vertical velocity after the array swap
-      !!
-      !! ** Method  : - compute the now divergence given by :
-      !!         * z-coordinate ONLY !!!!
-      !!         hdiv = 1/(e1t*e2t) [ di(e2u  u) + dj(e1v  v) ]
-      !!     - Using the incompressibility hypothesis, the vertical
-      !!      velocity is computed by integrating the horizontal divergence
-      !!      from the bottom to the surface.
-      !!        The boundary conditions are w=0 at the bottom (no flux).
-      !!----------------------------------------------------------------------
-      USE oce, ONLY:  zhdiv => hdivn
-      !
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  vertical velocity
-      !!
-      INTEGER  ::  ji, jj, jk
-      REAL(wp) ::  zu, zu1, zv, zv1, zet
-      !!----------------------------------------------------------------------
-      !
-      ! Computation of vertical velocity using horizontal divergence
-      zhdiv(:,:,:) = 0._wp
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk)
-               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk)
-               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk)
-               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk)
-               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
+      !
+      IF( lk_vvl ) THEN
+        IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND.    &                     ! Restart: read in restart file
+           iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN
+           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
+           CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:)   )
+           CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:)   )
+        ELSE
+           IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation'
+           CALL iom_open( 'restart', inum )
+           CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:)   )
+           CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:)   )
+           CALL iom_close( inum )                                        ! close file
+        ENDIF
+        !
+        DO jk = 1, jpkm1
+           fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
+        ENDDO
+        fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
+
+        ! Horizontal scale factor interpolations
+        ! --------------------------------------
+        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+
+        ! Vertical scale factor interpolations
+        ! ------------------------------------
+        CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' )
+  
+        fse3t_b(:,:,:)  = fse3t_n(:,:,:)
+        fse3u_b(:,:,:)  = fse3u_n(:,:,:)
+        fse3v_b(:,:,:)  = fse3v_n(:,:,:)
+
+        ! t- and w- points depth
+        ! ----------------------
+        fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+        fsdepw_n(:,:,1) = 0.0_wp
+
+        DO jk = 2, jpk
+           DO jj = 1,jpj
+              DO ji = 1,jpi
+                !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere
+                !    tmask = wmask, ie everywhere expect at jk = mikt
+                                                                   ! 1 for jk =
+                                                                   ! mikt
+                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
+                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
+              END DO
+           END DO
+        END DO
+
+        fsdept_b(:,:,:) = fsdept_n(:,:,:)
+        fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+        !
+      ENDIF
+      !
+      IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN       ! read depht over which runoffs are distributed
+         IF(lwp) WRITE(numout,*) 
+         IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed'
+         CALL iom_open ( "runoffs", inum )                           ! open file
+         CALL iom_get  ( inum, jpdom_data, 'rodepth', h_rnf )   ! read the river mouth array
+         CALL iom_close( inum )                                        ! close file
+         !
+         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( h_rnf(ji,jj) > 0._wp ) THEN
+                  jk = 2
+                  DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ;  jk = jk + 1
+                  END DO
+                  nk_rnf(ji,jj) = jk
+               ELSEIF( h_rnf(ji,jj) == -1._wp   ) THEN   ;  nk_rnf(ji,jj) = 1
+               ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN   ;  nk_rnf(ji,jj) = mbkt(ji,jj)
+               ELSE
+                  CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999'  )
+                  WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj)
+               ENDIF
             END DO
          END DO
+         DO jj = 1, jpj                                ! set the associated depth
+            DO ji = 1, jpi
+               h_rnf(ji,jj) = 0._wp
+               DO jk = 1, nk_rnf(ji,jj)
+                  h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE                                       ! runoffs applied at the surface
+         nk_rnf(:,:) = 1
+         h_rnf (:,:) = fse3t(:,:,1)
+      ENDIF
+      nkrnf_max = MAXVAL( nk_rnf(:,:) )
+      hrnf_max = MAXVAL( h_rnf(:,:) )
+      IF( lk_mpp )  THEN
+         CALL mpp_max( nkrnf_max )                 ! max over the  global domain
+         CALL mpp_max( hrnf_max )                 ! max over the  global domain
+      ENDIF
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,'    max level  : ', nkrnf_max
+      IF(lwp) WRITE(numout,*) ' '
+      !
+      CALL dta_dyn( nit000 )
+      !
+   END SUBROUTINE dta_dyn_init
+
+   SUBROUTINE dta_dyn_swp( kt )
+     !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE dta_dyn_swp  ***
+      !!
+      !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point
+      !!              and the depht
+      !!
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt       ! time step
+      INTEGER             :: ji, jj, jk
+      REAL(wp)            :: zcoef
+      !
+      !!---------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+
+      sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:))  ! before <-- now filtered
+      sshn(:,:) = ssha(:,:)
+
+      fse3t_n(:,:,:) = fse3t_a(:,:,:)
+
+      ! Reconstruction of all vertical scale factors at now and before time steps
+      ! =============================================================================
+
+      ! Horizontal scale factor interpolations
+      ! --------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+
+      ! Vertical scale factor interpolations
+      ! ------------------------------------
+      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' )
+
+      fse3t_b(:,:,:)  = fse3t_n(:,:,:)
+      fse3u_b(:,:,:)  = fse3u_n(:,:,:)
+      fse3v_b(:,:,:)  = fse3v_n(:,:,:)
+
+      ! t- and w- points depth
+      ! ----------------------
+      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.0_wp
+
+      DO jk = 2, jpk
+         DO jj = 1,jpj
+            DO ji = 1,jpi
+                 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
+                 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1)
+                 fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  &
+                     &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk))
+              END DO
+           END DO
+        END DO
+
+      fsdept_b(:,:,:) = fsdept_n(:,:,:)
+      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
+
+      !
+   END SUBROUTINE dta_dyn_swp
+
+   SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb,  pemp, pssha, pe3ta )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dta_dyn_wzv  ***
+      !!                   
+      !! ** Purpose :   compute the after ssh (ssha) and the now vertical velocity
+      !!
+      !! ** Method  : Using the incompressibility hypothesis, 
+      !!        - the ssh increment is computed by integrating the horizontal divergence 
+      !!          and multiply by the time step.
+      !!
+      !!        - compute the after scale factor : repartition of ssh INCREMENT proportionnaly
+      !!                                           to the level thickness ( z-star case )
+      !!
+      !!        - the vertical velocity is computed by integrating the horizontal divergence  
+      !!          from the bottom to the surface minus the scale factor evolution.
+      !!          The boundary conditions are w=0 at the bottom (no flux)
+      !!
+      !! ** action  :   ssha / e3t_a / wn
+      !!
+      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER,                                   INTENT(in )    :: kt        !  time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk)          , INTENT(in )   :: phdivtr   ! horizontal divergence transport
+      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: psshb     ! now ssh
+      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(in )   :: pemp      ! evaporation minus precipitation
+      REAL(wp), DIMENSION(jpi,jpj)    , OPTIONAL, INTENT(inout) :: pssha     ! after ssh
+      REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out)   :: pe3ta     ! after vertical scale factor
+      !! * Local declarations
+      INTEGER                       :: jk
+      REAL(wp), DIMENSION(jpi,jpj)  :: zhdiv  
+      REAL(wp)  :: z2dt  
+      !!----------------------------------------------------------------------
+      
+      !
+      z2dt = 2._wp * rdt
+      !
+      zhdiv(:,:) = 0._wp
+      DO jk = 1, jpkm1
+         zhdiv(:,:) = zhdiv(:,:) +  phdivtr(:,:,jk) * tmask(:,:,jk)
       END DO
-      CALL lbc_lnk( zhdiv, 'T', 1. )      ! Lateral boundary conditions on zhdiv
-      !
-      ! computation of vertical velocity from the bottom
-      pw(:,:,jpk) = 0._wp
-      DO jk = jpkm1, 1, -1
-         pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)
+      !                                                ! Sea surface  elevation time-stepping
+      pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:)  + zhdiv(:,:) ) ) * ssmask(:,:)
+      !                                                 ! 
+      !                                                 ! After acale factors at t-points ( z_star coordinate )
+      DO jk = 1, jpkm1
+        pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) )
       END DO
       !
-   END SUBROUTINE dta_dyn_wzv
-
-   SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj )
+   END SUBROUTINE dta_dyn_ssh
+
+
+   SUBROUTINE dta_dyn_hrnf
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE sbc_rnf  ***
+      !!
+      !! ** Purpose :   update the horizontal divergence with the runoff inflow
+      !!
+      !! ** Method  :
+      !!                CAUTION : rnf is positive (inflow) decreasing the
+      !!                          divergence and expressed in m/s
+      !!
+      !! ** Action  :   phdivn   decreased by the runoff inflow
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      DO jj = 1, jpj                   ! update the depth over which runoffs are distributed
+         DO ji = 1, jpi
+            h_rnf(ji,jj) = 0._wp
+            DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres
+                h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk)   ! to the bottom of the relevant grid box
+            END DO
+        END DO
+      END DO
+      !
+   END SUBROUTINE dta_dyn_hrnf
+
+
+
+   SUBROUTINE dta_dyn_slp( kt )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE dta_dyn_slp  ***
+      !!
+      !! ** Purpose : Computation of slope
+      !!
+      !!---------------------------------------------------------------------
+      USE oce, ONLY:  zts => tsa 
+      !
+      INTEGER,  INTENT(in) :: kt       ! time step
+      !
+      INTEGER  ::   ji, jj     ! dummy loop indices
+      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation
+      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation
+      INTEGER  ::   iswap
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj
+      !!---------------------------------------------------------------------
+      !
+      CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
+      !
+      IF( sf_dyn(jf_tem)%ln_tint ) THEN    ! Computes slopes (here avt is used as workspace)                       
+         IF( kt == nit000 ) THEN
+            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:)   ! temperature
+            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:)   ! salinity 
+            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:)   ! vertical diffusive coef.
+            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+            uslpdta (:,:,:,1) = zuslp (:,:,:) 
+            vslpdta (:,:,:,1) = zvslp (:,:,:) 
+            wslpidta(:,:,:,1) = zwslpi(:,:,:) 
+            wslpjdta(:,:,:,1) = zwslpj(:,:,:) 
+            !
+            zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
+            zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity 
+            avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
+            CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+            uslpdta (:,:,:,2) = zuslp (:,:,:) 
+            vslpdta (:,:,:,2) = zvslp (:,:,:) 
+            wslpidta(:,:,:,2) = zwslpi(:,:,:) 
+            wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
+         ELSE
+           ! 
+           iswap = 0
+           IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 )  iswap = 1
+           IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 )  THEN    ! read/update the after data
+              IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt
+              uslpdta (:,:,:,1) =  uslpdta (:,:,:,2)         ! swap the data
+              vslpdta (:,:,:,1) =  vslpdta (:,:,:,2)  
+              wslpidta(:,:,:,1) =  wslpidta(:,:,:,2) 
+              wslpjdta(:,:,:,1) =  wslpjdta(:,:,:,2) 
+              !
+              zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:)   ! temperature
+              zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:)   ! salinity 
+              avt(:,:,:)        = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:)   ! vertical diffusive coef.
+              CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+              !
+              uslpdta (:,:,:,2) = zuslp (:,:,:) 
+              vslpdta (:,:,:,2) = zvslp (:,:,:) 
+              wslpidta(:,:,:,2) = zwslpi(:,:,:) 
+              wslpjdta(:,:,:,2) = zwslpj(:,:,:) 
+            ENDIF
+         ENDIF
+      ENDIF
+      !
+      IF( sf_dyn(jf_tem)%ln_tint )  THEN
+         ztinta =  REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp )  &
+            &    / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp )
+         ztintb =  1. - ztinta
+#if defined key_ldfslp && ! defined key_c1d
+         uslp (:,:,:) = ztintb * uslpdta (:,:,:,1)  + ztinta * uslpdta (:,:,:,2)  
+         vslp (:,:,:) = ztintb * vslpdta (:,:,:,1)  + ztinta * vslpdta (:,:,:,2)  
+         wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1)  + ztinta * wslpidta(:,:,:,2)  
+         wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1)  + ztinta * wslpjdta(:,:,:,2)  
+#endif
+      ELSE
+         zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:)   ! temperature
+         zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:)   ! salinity 
+         avt(:,:,:)        = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:)   ! vertical diffusive coef.
+         CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj )
+         !
+#if defined key_ldfslp && ! defined key_c1d
+         uslp (:,:,:) = zuslp (:,:,:)
+         vslp (:,:,:) = zvslp (:,:,:)
+         wslpi(:,:,:) = zwslpi(:,:,:)
+         wslpj(:,:,:) = zwslpj(:,:,:)
+#endif
+      ENDIF
+      !
+      CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj )
+      !
+   END SUBROUTINE dta_dyn_slp
+
+   SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE dta_dyn_slp  ***
@@ -568,5 +711,6 @@
 #endif
       !
-   END SUBROUTINE dta_dyn_slp
+   END SUBROUTINE compute_slopes
+
    !!======================================================================
 END MODULE dtadyn
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 7522)
@@ -34,5 +34,5 @@
    USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
    USE dtadyn          ! Lecture and interpolation of the dynamical fields
-   !              ! I/O & MPP
+   !                   ! I/O & MPP
    USE iom             ! I/O library
    USE in_out_manager  ! I/O manager
@@ -50,5 +50,4 @@
    USE trcnam
    USE trcrst
-   USE diaptr         ! Need to initialise this as some variables are used in if statements later
 
    IMPLICIT NONE
@@ -94,13 +93,14 @@
       istp = nit000
       ! 
-      CALL iom_init( cxios_context )            ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
-      ! 
       DO WHILE ( istp <= nitend .AND. nstop == 0 )    ! time stepping
          !
-         IF( istp /= nit000 )   CALL day      ( istp )         ! Calendar (day was already called at nit000 in day_init)
-                                CALL iom_setkt( istp - nit000 + 1, "nemo" )   ! say to iom that we are at time step kstp
-                                CALL dta_dyn  ( istp )         ! Interpolation of the dynamical fields
-                                CALL trc_stp  ( istp )         ! time-stepping
-                                CALL stp_ctl  ( istp, indic )  ! Time loop: control and print
+         IF( istp == nit000 )  CALL iom_init( cxios_context )            ! iom_put initialization
+         IF( istp /= nit000 )   CALL day        ( istp )         ! Calendar (day was already called at nit000 in day_init)
+                                CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp
+                                CALL trc_rst_opn( istp )         ! Open tracer                                !   restart file
+                                CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields
+                                CALL trc_stp    ( istp )         ! time-stepping
+         IF( lk_vvl )           CALL dta_dyn_swp( istp )         ! swap of sea surface height and vertical scale factors
+                                CALL stp_ctl    ( istp, indic )  ! Time loop: control and print
          istp = istp + 1
          IF( lk_mpp )   CALL mpp_max( nstop )
@@ -265,19 +265,20 @@
       IF( nn_timing == 1 )  CALL timing_start( 'nemo_init')
       !
-                            CALL     phy_cst    ! Physical constants
-                            CALL     eos_init   ! Equation of state
-      IF( lk_c1d        )   CALL     c1d_init   ! 1D column configuration
-                            CALL     dom_cfg    ! Domain configuration
+                            CALL  phy_cst    ! Physical constants
+                            CALL  eos_init   ! Equation of state
+      IF( lk_c1d        )   CALL  c1d_init   ! 1D column configuration
+                            CALL  dom_cfg    ! Domain configuration
+      !
       !
       INQUIRE( FILE='coordinates.nc', EXIST = llexist )   ! Check if coordinate file exist
       !
-      IF( llexist )  THEN  ;  CALL  dom_init   !  compute the grid from coordinates and bathymetry
-      ELSE                 ;  CALL  dom_rea    !  read grid from the meskmask
+      IF( llexist )  THEN ;  CALL  dom_init   !  compute the grid from coordinates and bathymetry
+      ELSE                ;  CALL  dom_rea    !  read grid from the meskmask
       ENDIF
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
-      IF( ln_nnogather )    CALL nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined)
-
-      IF( ln_ctl        )   CALL prt_ctl_init   ! Print control
+      IF( ln_nnogather )    CALL  nemo_northcomms   ! Initialise the northfold neighbour lists (must be done after the masks are defined)
+
+      IF( ln_ctl       )    CALL  prt_ctl_init   ! Print control
 
                             CALL     sbc_init   ! Forcings : surface module
@@ -289,12 +290,10 @@
 
                             CALL tra_qsr_init   ! penetrative solar radiation qsr
-      IF( lk_trabbl     )   CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( lk_trabbl   )     CALL tra_bbl_init   ! advective (and/or diffusive) bottom boundary layer scheme
 
                             CALL trc_nam_run    ! Needed to get restart parameters for passive tracers
                             CALL trc_rst_cal( nit000, 'READ' )   ! calendar
                             CALL dta_dyn_init   ! Initialization for the dynamics
-
                             CALL     trc_init   ! Passive tracers initialization
-                            CALL dia_ptr_init   ! Initialise diaptr as some variables are used 
       !                                         ! in various advection and diffusion routines
       IF(lwp) WRITE(numout,cform_aaa)           ! Flag AAAAAAA
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 7522)
@@ -23,4 +23,5 @@
    USE dom_oce         ! domain: ocean
    USE sbc_oce         ! surface boundary condition: ocean
+   USE trc_oce         ! shared ocean-passive tracers variables
    USE phycst          ! physical constants
    USE closea          ! closed seas
@@ -97,20 +98,23 @@
       END DO
       !
-      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh
-      !
-      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
-      !
-      !
-      hu(:,:) = 0._wp                          ! Ocean depth at U-points
-      hv(:,:) = 0._wp                          ! Ocean depth at V-points
-      ht(:,:) = 0._wp                          ! Ocean depth at T-points
-      DO jk = 1, jpkm1
-         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
-         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
-         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
-      END DO
-      !                                        ! Inverse of the local depth
-      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
-      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+      IF( lk_c1d )           CALL cor_c1d      ! 1D configuration: Coriolis set at T-point
+      !
+      IF( .NOT.lk_offline ) THEN
+        !
+        IF( lk_vvl )         CALL dom_vvl_init ! Vertical variable mesh
+        !
+        hu(:,:) = 0._wp                          ! Ocean depth at U-points
+        hv(:,:) = 0._wp                          ! Ocean depth at V-points
+        ht(:,:) = 0._wp                          ! Ocean depth at T-points
+        DO jk = 1, jpkm1
+           hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
+           hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
+           ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
+        END DO
+        !                                        ! Inverse of the local depth
+        hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:)
+        hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:)
+        !
+      ENDIF
 
                              CALL dom_stp      ! time step
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 7522)
@@ -263,10 +263,4 @@
       ENDIF
 
-      ! set the number of level over which river runoffs are applied 
-      ! online configuration : computed in sbcrnf
-      IF( lk_offline ) THEN
-        nk_rnf(:,:) = 1
-        h_rnf (:,:) = e3t_0(:,:,1)
-      ENDIF
 
       ! dust input from the atmosphere
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90	(revision 7522)
@@ -88,26 +88,35 @@
          r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog)
       ENDIF
-      !                                                   ! effective transport
-      DO jk = 1, jpkm1
-         !                                                ! eulerian transport only
-         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)
-         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
-         !
-      END DO
-      !
-      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
-         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
-         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
+      !  
+      IF( lk_offline ) THEN
+         zun(:,:,:) = un(:,:,:)     ! effective transport already in un/vn/wn
+         zvn(:,:,:) = vn(:,:,:)
+         zwn(:,:,:) = wn(:,:,:)
+      ELSE
+         !                                                         ! effective transport
+         DO jk = 1, jpkm1
+            !                                                ! eulerian transport only
+            zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)
+            zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+            zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
+            !
+         END DO
+         !
+         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
+            zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
+            zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
+         ENDIF
+         !
+         zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
+         zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
+         zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
+         !
+
+         IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary)
+            &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )
+         !
+         IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary)
+         !
       ENDIF
-      !
-      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
-
-      IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary)
-         &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )
-      !
-      IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary)
       !
       SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90	(revision 7522)
@@ -65,5 +65,5 @@
       IF( lk_c14b    )   CALL trc_rad_sms( kt, trb, trn, jp_c14b0, jp_c14b1              )  ! bomb C14
       IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model
-      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1               )  ! MY_TRC model
+      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1, cpreserv='Y' )  ! MY_TRC model
 
       !
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7522)
@@ -128,7 +128,7 @@
       ! Coupling offline : runoff are in emp which contains E-P-R
       !
-      IF( .NOT. lk_offline .AND. lk_vvl ) THEN  ! online coupling with vvl
+      IF( lk_vvl ) THEN                         ! linear free surface vvl
          zsfx(:,:) = 0._wp
-      ELSE                                      ! online coupling free surface or offline with free surface
+      ELSE                                      ! no vvl
          zsfx(:,:) = emp(:,:)
       ENDIF
Index: /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7521)
+++ /branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 7522)
@@ -86,5 +86,5 @@
          tra(:,:,:,:) = 0.e0
          !
-                                   CALL trc_rst_opn  ( kt )       ! Open tracer restart file 
+         IF( .NOT.lk_offline )     CALL trc_rst_opn  ( kt )       ! Open tracer restart file 
          IF( lrst_trc )            CALL trc_rst_cal  ( kt, 'WRITE' )   ! calendar
          IF( lk_iomput ) THEN  ;   CALL trc_wri      ( kt )       ! output of passive tracers with iom I/O manager
