Index: trunk/NEMO/OPA_SRC/cla_div.F90
===================================================================
--- trunk/NEMO/OPA_SRC/cla_div.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/cla_div.F90	(revision 467)
@@ -104,5 +104,4 @@
       !!          surface                        depth
       !!      The now divergence is given by :
-      !!       * z-coordinate (default key) and partial steps (key_partial_steps)
       !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
       !!
Index: trunk/NEMO/OPA_SRC/eosbn2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/eosbn2.F90	(revision 467)
@@ -489,5 +489,5 @@
          DO jj = 1, jpjm1
 !CDIR NOVERRCHK
-#if defined key_autotasking
+#if defined key_mpp_omp
             DO ji = 1, jpim1
 #else
@@ -501,5 +501,5 @@
          DO jj = 1, jpjm1                                 ! Horizontal slab
             !                                             ! ===============
-#if defined key_autotasking
+#if defined key_mpp_omp
             DO ji = 1, jpim1
 #else
@@ -556,5 +556,5 @@
          DO jj = 1, jpjm1                                 ! Horizontal slab
             !                                             ! ===============
-#if defined key_autotasking
+#if defined key_mpp_omp
             DO ji = 1, jpim1
 #else
@@ -573,5 +573,5 @@
          DO jj = 1, jpjm1                                 ! Horizontal slab
             !                                             ! ===============
-#if defined key_autotasking
+#if defined key_mpp_omp
             DO ji = 1, jpim1
 #else
Index: trunk/NEMO/OPA_SRC/istate.F90
===================================================================
--- trunk/NEMO/OPA_SRC/istate.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/istate.F90	(revision 467)
@@ -264,6 +264,6 @@
             IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
-            zh1 = gdept(  1  )
-            zh2 = gdept(jpkm1)
+            zh1 = gdept_0(  1  )
+            zh2 = gdept_0(jpkm1)
 
             zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )
@@ -534,6 +534,6 @@
          WRITE(numout,*)
          WRITE(numout,*) '              Initial temperature and salinity profiles:'
-         WRITE(numout, "(9x,' level   gdept   temperature   salinity   ')" )
-         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept(jk), tn(2,2,jk), sn(2,2,jk), jk = 1, jpk )
+         WRITE(numout, "(9x,' level   gdept_0   temperature   salinity   ')" )
+         WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_0(jk), tn(2,2,jk), sn(2,2,jk), jk = 1, jpk )
       ENDIF
 
Index: trunk/NEMO/OPA_SRC/mppini_2.h90
===================================================================
--- trunk/NEMO/OPA_SRC/mppini_2.h90	(revision 466)
+++ trunk/NEMO/OPA_SRC/mppini_2.h90	(revision 467)
@@ -115,10 +115,10 @@
 
       ! open the file
-         IF ( lk_zps ) THEN 
+         IF ( ln_zco ) THEN 
+            clname = 'bathy_level.nc'         ! Level bathymetry
+            clvar = 'Bathy_level'
+         ELSE
             clname = 'bathy_meter.nc'         ! Meter bathy in case of partial steps
             clvar = 'Bathymetry'
-         ELSE
-            clname = 'bathy_level.nc'                       ! Level bathymetry
-            clvar = 'Bathy_level'
          ENDIF
 #if defined key_agrif
@@ -152,5 +152,5 @@
       ! land/sea mask over the global/zoom domain
 
-      imask(:,:)=1
+      imask(:,:) = 1
       WHERE ( zdta(jpizoom:(jpizoom+jpiglo-1),jpjzoom:(jpjglo+jpjzoom-1)) <= 0. ) imask = 0
 
Index: trunk/NEMO/OPA_SRC/oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/oce.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/oce.F90	(revision 467)
@@ -42,7 +42,7 @@
       !! -----------------------
       CHARACTER(len=3), PUBLIC  ::   l_adv   !: 'ce2' centre scheme used
-              !                              !: 'tvd' TVD scheme used
-              !                              !: 'mus' MUSCL scheme used
-              !                              !: 'mu2' MUSCL2 scheme used
+         !                                   !: 'tvd' TVD scheme used
+         !                                   !: 'mus' MUSCL scheme used
+         !                                   !: 'mu2' MUSCL2 scheme used
 
    !! surface pressure gradient
@@ -51,29 +51,14 @@
       spgu, spgv             !: horizontal surface pressure gradient
 
-#if defined key_partial_steps     ||   defined key_esopa
-   !! interpolated gradient
+   !! interpolated gradient (only used in zps case)
    !! ---------------------
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
       gtu, gsu, gru,      &  !: t-, s- and rd horizontal gradient at u- and 
       gtv, gsv, grv          !: v-points at bottom ocean level 
-#else
-   REAL(wp), PUBLIC ::   &   !:
-      gtu, gsu, gru,      &  !: dummy scalars
-      gtv, gsv, grv          !:
-#endif
 
    !! free surface
    !! ------------
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      sshb, sshn,         &  !: before, now sea surface height (meters)
-      hu  , hv  ,         &  !: depth at u- and v-points (meters)
-      hur , hvr              !: inverse of u and v-points ocean depth (1/m)
-#if defined key_obc
-   REAL(wp), PUBLIC ::    &  !:
-      obcsurftot       !: Total lateral surface of open boundaries
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &   !:
-      obcumask, obcvmask     !: u-, v- Force filtering mask for the open 
-      !                      !  boundary condition on grad D
-#endif
+      sshb, sshn             !: before, now sea surface height (meters)
 
 #if defined key_dynspg_rl   ||   defined key_esopa
Index: trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- trunk/NEMO/OPA_SRC/opa.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/opa.F90	(revision 467)
@@ -7,6 +7,30 @@
    !!----------------------------------------------------------------------
    !!   opa_model      : solve ocean dynamics, tracer and/or sea-ice
+   !!   opa_init       : initialization of the opa model
    !!   opa_flg        : initialisation of algorithm flag 
    !!   opa_closefile  : close remaining files
+   !!----------------------------------------------------------------------
+   !! History :
+   !!   4.0  !  90-10  (C. Levy, G. Madec)  Original code
+   !!   7.0  !  91-11  (M. Imbard, C. Levy, G. Madec)
+   !!   7.1  !  93-03  (M. Imbard, C. Levy, G. Madec, O. Marti,
+   !!                   M. Guyon, A. Lazar, P. Delecluse, C. Perigaud,
+   !!                   G. Caniaux, B. Colot, C. Maes ) release 7.1 
+   !!        !  92-06  (L.Terray) coupling implementation
+   !!        !  93-11  (M.A. Filiberti) IGLOO sea-ice 
+   !!   8.0  !  96-03  (M. Imbard, C. Levy, G. Madec, O. Marti,
+   !!                   M. Guyon, A. Lazar, P. Delecluse, L.Terray,
+   !!                   M.A. Filiberti, J. Vialar, A.M. Treguier,
+   !!                   M. Levy)  release 8.0
+   !!   8.1  !  97-06  (M. Imbard, G. Madec)
+   !!   8.2  !  99-11  (M. Imbard, H. Goosse)  LIM sea-ice model 
+   !!        !  99-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP 
+   !!        !  00-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
+   !!   9.0  !  02-08  (G. Madec)  F90: Free form and modules
+   !!    "   !  04-08  (C. Talandier) New trends organization
+   !!    "   !  05-06  (C. Ethe) Add the 1D configuration possibility
+   !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!    "   !  06-03  (L. Debreu, C. Mazauric)  Agrif implementation
+   !!    "   !  06-04  (G. Madec, R. Benshila)  Step reorganization
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -30,5 +54,4 @@
 
    ! ocean physics
-   USE traqsr          ! solar radiation penetration   (tra_qsr_init routine)
    USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
    USE ldftra          ! lateral diffusivity setting    (ldftra_init routine)
@@ -59,4 +82,8 @@
    PRIVATE
 
+   !! * Module variables
+   CHARACTER (len=64) ::        &
+      cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
+
    !! * Routine accessibility
    PUBLIC opa_model      ! called by model.F90
@@ -65,5 +92,5 @@
    !!  OPA 9.0 , LOCEAN-IPSL (2005) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -83,40 +110,15 @@
       !!      Madec, Delecluse,Imbard, and Levy, 1997: reference manual.
       !!              internal report, IPSL.
-      !!
-      !! History :
-      !!   4.0  !  90-10  (C. Levy, G. Madec)  Original code
-      !!   7.0  !  91-11  (M. Imbard, C. Levy, G. Madec)
-      !!   7.1  !  93-03  (M. Imbard, C. Levy, G. Madec, O. Marti,
-      !!                   M. Guyon, A. Lazar, P. Delecluse, C. Perigaud,
-      !!                   G. Caniaux, B. Colot, C. Maes ) release 7.1 
-      !!        !  92-06  (L.Terray) coupling implementation
-      !!        !  93-11  (M.A. Filiberti) IGLOO sea-ice 
-      !!   8.0  !  96-03  (M. Imbard, C. Levy, G. Madec, O. Marti,
-      !!                   M. Guyon, A. Lazar, P. Delecluse, L.Terray,
-      !!                   M.A. Filiberti, J. Vialar, A.M. Treguier,
-      !!                   M. Levy)  release 8.0
-      !!   8.1  !  97-06  (M. Imbard, G. Madec)
-      !!   8.2  !  99-11  (M. Imbard, H. Goosse)  LIM sea-ice model 
-      !!        !  99-12  (V. Thierry, A-M. Treguier, M. Imbard, M-A. Foujols)  OPEN-MP 
-      !!        !  00-07  (J-M Molines, M. Imbard)  Open Boundary Conditions  (CLIPPER)
-      !!   9.0  !  02-08  (G. Madec)  F90: Free form and modules
-      !!    "   !  04-08  (C. Talandier) New trends organization
-      !!    "   !  05-06  (C. Ethe) Add the 1D configuration possibility
-      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER ::   istp       ! time step index
-      CHARACTER (len=64) ::        &
-         cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
       !!----------------------------------------------------------------------
 
 #if defined key_agrif
-
-      Call Agrif_Init_Grids()
+      CALL Agrif_Init_Grids()
 #endif
       
-      Call opa_init  ! Initializations
-
-      IF( lk_cfg_1d  )  THEN 
+      CALL opa_init  ! Initializations
+
+      IF( lk_cfg_1d ) THEN 
          istp = nit000
          DO WHILE ( istp <= nitend .AND. nstop == 0 )
@@ -156,4 +158,149 @@
 
 
+   SUBROUTINE opa_init
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE opa_init  ***
+      !!
+      !! ** Purpose :   initialization of the opa model
+      !!
+      !!----------------------------------------------------------------------
+#if defined key_coupled
+      INTEGER ::   itro, istp0        ! ???
+#endif
+      CHARACTER (len=20) ::   namelistname
+      CHARACTER (len=28) ::   file_out
+      !!----------------------------------------------------------------------
+
+      ! Initializations
+      ! ===============
+
+      file_out = 'ocean.output'
+      
+      ! open listing and namelist units
+      IF ( numout /= 0 .AND. numout /= 6 ) THEN 
+         CALL ctlopn( numout, file_out, 'UNKNOWN', 'FORMATTED',   &
+            &         'SEQUENTIAL', 1, numout, .FALSE., 1 )
+      ENDIF
+
+      namelistname = 'namelist'
+      CALL ctlopn( numnam, namelistname, 'OLD', 'FORMATTED', 'SEQUENTIAL',   &
+         &           1, numout, .FALSE., 1 )
+
+      WRITE(numout,*)
+      WRITE(numout,*) '                 L O D Y C - I P S L'
+      WRITE(numout,*) '                     O P A model'
+      WRITE(numout,*) '            Ocean General Circulation Model'
+      WRITE(numout,*) '               version OPA 9.0  (2005) '
+      WRITE(numout,*)
+      WRITE(numout,*)
+
+      ! Nodes selection
+      narea = mynode()
+      narea = narea + 1    ! mynode return the rank of proc (0 --> jpnij -1 )
+      lwp   = narea == 1
+
+      !                                     ! ============================== !
+      !                                     !  Model general initialization  !
+      !                                     ! ============================== !
+
+      IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
+
+                                            ! Domain decomposition
+      IF( jpni*jpnj == jpnij ) THEN
+         CALL mpp_init                          ! standard cutting out
+      ELSE
+         CALL mpp_init2                         ! eliminate land processors
+      ENDIF
+      
+      CALL phy_cst                          ! Physical constants
+
+      CALL dom_cfg                          ! Domain configuration
+      
+      CALL dom_init                         ! Domain
+
+      IF( ln_ctl )      CALL prt_ctl_init   ! Print control
+
+      IF( lk_cfg_1d )   CALL fcorio_1d      ! redefine Coriolis at T-point
+
+      IF( lk_obc    )   CALL obc_init       ! Open boundaries 
+
+      CALL day( nit000 )                    ! Calendar
+
+      CALL istate_init                      ! ocean initial state (Dynamics and tracers)
+
+      IF( lk_dynspg_flt .OR. lk_dynspg_rl ) THEN
+         CALL solver_init( nit000 )         ! Elliptic solver
+      ENDIF
+
+!!add
+                       CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
+
+                       CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency
+
+      IF( ln_zps .AND. .NOT. lk_cfg_1d )   &
+         &             CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before Horizontal DErivative
+                                            gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
+                                            gtv, gsv, grv )
+!!add
+
+      CALL oc_fz_pt                         ! Surface freezing point
+
+#if defined key_ice_lim
+      CALL ice_init                         ! Sea ice model
+#endif
+
+      !                                     ! Ocean scheme
+
+      CALL opa_flg                              ! Choice of algorithms
+
+      !                                     ! Ocean physics
+
+      CALL ldf_dyn_init                         ! Lateral ocean momentum physics
+
+      CALL ldf_tra_init                         ! Lateral ocean tracer physics
+
+      CALL zdf_init                             ! Vertical ocean physics
+
+      !                                     ! Ocean trends
+      ! Control parameters 
+      IF( lk_trdtra .OR. lk_trdmld )   l_trdtra = .TRUE.
+      IF( lk_trddyn .OR. lk_trdvor )   l_trddyn = .TRUE.
+
+      IF( lk_trddyn .OR. lk_trdtra )   &
+         &            CALL trd_icp_init         ! active tracers and/or momentum
+
+      IF( lk_trdmld ) CALL trd_mld_init         ! mixed layer
+
+      IF( lk_trdvor ) CALL trd_vor_init         ! vorticity
+
+#if defined key_passivetrc
+      CALL ini_trc                           ! Passive tracers
+#endif
+
+#if defined key_coupled
+      itro  = nitend - nit000 + 1           ! Coupled
+      istp0 = NINT( rdt )
+      CALL cpl_init( itro, nexco, istp0 )   ! Signal processing and process id exchange
+#endif
+
+      CALL flx_init                         ! Thermohaline forcing initialization
+
+      CALL flx_fwb_init                     ! FreshWater Budget correction
+
+      CALL dia_ptr_init                     ! Poleward TRansports initialization
+
+      !                                     ! =============== !
+      !                                     !  time stepping  !
+      !                                     ! =============== !
+
+      IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
+
+      IF( lk_cfg_1d  )  THEN 
+         CALL init_1d
+      ENDIF
+
+   END SUBROUTINE opa_init
+
+
    SUBROUTINE opa_flg
       !!----------------------------------------------------------------------
@@ -173,5 +320,5 @@
       !!----------------------------------------------------------------------
 
-      ! Read Namelist namflg : algorithm FLaG 
+      ! Read Namelist namflg : algorithm FLaG
       ! --------------------
       REWIND ( numnam )
@@ -193,4 +340,5 @@
    END SUBROUTINE opa_flg
 
+
    SUBROUTINE opa_closefile
       !!----------------------------------------------------------------------
@@ -198,6 +346,6 @@
       !!
       !! ** Purpose :   Close the files
-      !!           
-      !! ** Method  : 
+      !!
+      !! ** Method  :
       !!
       !! History :
@@ -236,157 +384,3 @@
 
    !!======================================================================
-   SUBROUTINE opa_init
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE opa_init  ***
-      !!
-      !! ** Purpose :   initialization of the opa model
-      !!
-      !! ** Method  : 
-      !!
-      !! References :
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-
-#if defined key_coupled
-      INTEGER ::   itro, istp0        ! ???
-#endif
-      CHARACTER (len=64) ::        &
-         cform_aaa="( /, 'AAAAAAAA', / ) "     ! flag for output listing
-      CHARACTER (len=20) :: namelistname
-      CHARACTER (len=28) :: file_out
-      !!----------------------------------------------------------------------
-
-      ! Initializations
-      ! ===============
-
-      file_out = 'ocean.output'
-      
-      ! open listing and namelist units
-      IF ( numout /= 0 .AND. numout /= 6 ) THEN 
-         CALL ctlopn(numout,file_out,'UNKNOWN', 'FORMATTED',   &
-                      'SEQUENTIAL',1,numout,.FALSE.,1)
-!         OPEN( UNIT=numout, FILE=TRIM(file_out), FORM='FORMATTED' )
-      ENDIF
-
-      namelistname = 'namelist'
-      CALL ctlopn(numnam,namelistname,'OLD', 'FORMATTED', 'SEQUENTIAL',   &
-                     1,numout,.FALSE.,1)
-!!!!      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
-
-      WRITE(numout,*)
-      WRITE(numout,*) '                 L O D Y C - I P S L'
-      WRITE(numout,*) '                     O P A model'
-      WRITE(numout,*) '            Ocean General Circulation Model'
-      WRITE(numout,*) '               version OPA 9.0  (2005) '
-      WRITE(numout,*)
-      WRITE(numout,*)
-
-      ! Nodes selection
-      narea = mynode()
-      narea = narea + 1    ! mynode return the rank of proc (0 --> jpnij -1 )
-      lwp   = narea == 1
-
-      !                                     ! ============================== !
-      !                                     !  Model general initialization  !
-      !                                     ! ============================== !
-
-      IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
-
-                                            ! Domain decomposition
-      IF( jpni*jpnj == jpnij ) THEN
-         CALL mpp_init                          ! standard cutting out
-      ELSE
-         CALL mpp_init2                         ! eliminate land processors
-      ENDIF
-      
-      CALL phy_cst                          ! Physical constants
-
-      CALL dom_cfg                          ! Domain configuration
-      
-      CALL dom_init                         ! Domain
-      IF( ln_ctl )      CALL prt_ctl_init   ! Print control
-
-      IF( lk_cfg_1d )   CALL fcorio_1d      ! redefine Coriolis at T-point
-
-      IF( lk_obc    )   CALL obc_init       ! Open boundaries 
-
-      CALL day( nit000 )                    ! Calendar
-
-      CALL istate_init                      ! ocean initial state (Dynamics and tracers)
-
-      IF( lk_dynspg_flt .OR. lk_dynspg_rl ) THEN
-         CALL solver_init( nit000 )         ! Elliptic solver
-      ENDIF
-
-!!add
-                       CALL eos( tb, sb, rhd, rhop )        ! before potential and in situ densities
-                       
-                       CALL bn2( tb, sb, rn2 )              ! before Brunt-Vaisala frequency
-
-      IF( lk_zps .AND. .NOT. lk_cfg_1d )   &
-         &             CALL zps_hde( nit000, tb, sb, rhd,  &  ! Partial steps: before Horizontal DErivative
-                                            gtu, gsu, gru, &  ! of t, s, rd at the bottom ocean level
-                                            gtv, gsv, grv )
-
-!!add
-
-      CALL oc_fz_pt                         ! Surface freezing point
-
-#if defined key_ice_lim
-      CALL ice_init                         ! Sea ice model
-#endif
-
-      !                                     ! Ocean scheme
-
-      CALL opa_flg                              ! Choice of algorithms
-
-      !                                     ! Ocean physics
-
-      CALL tra_qsr_init                         ! Solar radiation penetration
-
-      CALL ldf_dyn_init                         ! Lateral ocean momentum physics
-
-      CALL ldf_tra_init                         ! Lateral ocean tracer physics
-
-      CALL zdf_init                             ! Vertical ocean physics
-
-      !                                     ! Ocean trends
-      ! Control parameters 
-      IF( lk_trdtra .OR. lk_trdmld )   l_trdtra = .TRUE.
-      IF( lk_trddyn .OR. lk_trdvor )   l_trddyn = .TRUE.
-
-      IF( lk_trddyn .OR. lk_trdtra )   &
-         &            CALL trd_icp_init         ! active tracers and/or momentum
-
-      IF( lk_trdmld ) CALL trd_mld_init         ! mixed layer
-
-      IF( lk_trdvor ) CALL trd_vor_init         ! vorticity
-
-#if defined key_passivetrc
-      CALL ini_trc                           ! Passive tracers
-#endif
-
-#if defined key_coupled
-      itro  = nitend - nit000 + 1           ! Coupled
-      istp0 = NINT( rdt )
-      CALL cpl_init( itro, nexco, istp0 )   ! Signal processing and process id exchange
-#endif
-
-      CALL flx_init                         ! Thermohaline forcing initialization
-
-      CALL flx_fwb_init                     ! FreshWater Budget correction
-
-      CALL dia_ptr_init                     ! Poleward TRansports initialization
-
-      !                                     ! =============== !
-      !                                     !  time stepping  !
-      !                                     ! =============== !
-
-      IF(lwp) WRITE(numout,cform_aaa)       ! Flag AAAAAAA
-
-      IF( lk_cfg_1d  )  THEN 
-         CALL init_1d
-      ENDIF
-   END SUBROUTINE opa_init
-   !!======================================================================
 END MODULE opa
Index: trunk/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/par_oce.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/par_oce.F90	(revision 467)
@@ -209,5 +209,5 @@
 #endif
 
-#if defined key_autotasking
+#if defined key_mpp_omp
    LOGICAL, PUBLIC, PARAMETER ::   lk_jki = .TRUE.   !: j-k-i loop flag
 #else
Index: trunk/NEMO/OPA_SRC/restart.F90
===================================================================
--- trunk/NEMO/OPA_SRC/restart.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/restart.F90	(revision 467)
@@ -22,5 +22,5 @@
    USE flx_oce         ! sea-ice/ocean forcings variables
    USE dynspg_oce      ! free surface time splitting scheme variables
-   USE cpl_oce,         ONLY : lk_cpl              !
+   USE cpl_oce, ONLY : lk_cpl              !
 
    IMPLICIT NONE
@@ -151,5 +151,5 @@
          itime = 0
          CALL ymds2ju( nyear, nmonth, nday, 0.e0, zdate0 )
-         CALL restini( 'NONE', jpi, jpj, glamt, gphit, jpk, gdept, clname,   &
+         CALL restini( 'NONE', jpi, jpj, glamt, gphit, jpk, gdept_0, clname,   &
                         itime, zdate0, rdt*nstock ,numwrs, domain_id=nidom )
 
Index: trunk/NEMO/OPA_SRC/restart_dimg.h90
===================================================================
--- trunk/NEMO/OPA_SRC/restart_dimg.h90	(revision 466)
+++ trunk/NEMO/OPA_SRC/restart_dimg.h90	(revision 467)
@@ -108,17 +108,17 @@
        ! 'before' fields
        DO jk = 1, jpk
-          WRITE(inum,REC=irec) ub(:,:,jk)   ;    irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) vb(:,:,jk)   ;    irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) tb(:,:,jk)   ;    irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) sb(:,:,jk)   ;    irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) rotb(:,:,jk)   ;    irec = irec +1
+          WRITE(inum,REC=irec) ub   (:,:,jk)   ;    irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) vb   (:,:,jk)   ;    irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) tb   (:,:,jk)   ;    irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) sb   (:,:,jk)   ;    irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) rotb (:,:,jk)   ;    irec = irec +1
        END DO
        DO jk = 1, jpk
@@ -128,17 +128,17 @@
        ! 'now' fields
        DO jk = 1, jpk
-          WRITE(inum,REC=irec) un(:,:,jk)   ;   irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) vn(:,:,jk)   ;   irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) tn(:,:,jk)   ;   irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) sn(:,:,jk)   ;   irec = irec +1
-       END DO
-       DO jk = 1, jpk
-          WRITE(inum,REC=irec) rotn(:,:,jk)   ;   irec = irec +1
+          WRITE(inum,REC=irec) un   (:,:,jk)   ;   irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) vn   (:,:,jk)   ;   irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) tn   (:,:,jk)   ;   irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) sn   (:,:,jk)   ;   irec = irec +1
+       END DO
+       DO jk = 1, jpk
+          WRITE(inum,REC=irec) rotn (:,:,jk)   ;   irec = irec +1
        END DO
        DO jk = 1, jpk
@@ -147,5 +147,5 @@
 
        ! elliptic solver arrays
-       WRITE(inum,REC=irec ) gcx(1:jpi,1:jpj)   ;   irec = irec +1
+       WRITE(inum,REC=irec ) gcx (1:jpi,1:jpj)   ;   irec = irec +1
        WRITE(inum,REC=irec ) gcxb(1:jpi,1:jpj)   ;   irec = irec +1
 #if defined key_dynspg_rl
@@ -173,5 +173,5 @@
 #if defined key_zdftke
          DO jk = 1, jpk
-            WRITE(inum,REC=irec) en(:,:,jk) ; irec = irec + 1 
+            WRITE(inum,REC=irec) en(:,:,jk)   ;   irec = irec + 1 
          END DO
 #endif
@@ -179,17 +179,17 @@
 #if defined key_ice_lim
           zfice(1) = FLOAT( nfice )                                      ! Louvain La Neuve Sea Ice Model
-          WRITE(inum,REC=irec) zfice(:)     ; irec = irec + 1
-          WRITE(inum,REC=irec) sst_io(:,:)  ; irec = irec + 1
-          WRITE(inum,REC=irec) sss_io(:,:)  ; irec = irec + 1
-          WRITE(inum,REC=irec) u_io  (:,:)  ; irec = irec + 1
-          WRITE(inum,REC=irec) v_io  (:,:)  ; irec = irec + 1
+          WRITE(inum,REC=irec) zfice(:)     ;   irec = irec + 1
+          WRITE(inum,REC=irec) sst_io(:,:)  ;   irec = irec + 1
+          WRITE(inum,REC=irec) sss_io(:,:)  ;   irec = irec + 1
+          WRITE(inum,REC=irec) u_io  (:,:)  ;   irec = irec + 1
+          WRITE(inum,REC=irec) v_io  (:,:)  ;   irec = irec + 1
 #    if defined key_coupled
-          WRITE(inum,REC=irec) alb_ice(:,:)  ; irec = irec + 1
+          WRITE(inum,REC=irec) alb_ice(:,:)   ;   irec = irec + 1
 #    endif
 #endif
 # if defined key_flx_bulk_monthly || defined key_flx_bulk_daily
           zfblk(1) = FLOAT( nfbulk )                                 ! Bulk
-          WRITE(inum,REC=irec) zfblk(:)   ; irec = irec + 1
-          WRITE(inum,REC=irec) gsst(:,:)  ; irec = irec + 1
+          WRITE(inum,REC=irec) zfblk(:)     ;   irec = irec + 1
+          WRITE(inum,REC=irec) gsst (:,:)   ;   irec = irec + 1
 # endif
 
@@ -288,6 +288,6 @@
 
     READ(inum,REC=1) irecl8, ino1, it1, isor1, ipcg1, itke1, &
-     &  iice1, ibulk1, ios1, ios2, ios3, ios4, &
-     &  idast1, adatrj0,  ipi,ipj,ipk,ipni,ipnj,ipnij,iarea
+       &             iice1, ibulk1, ios1, ios2, ios3, ios4, &
+       &             idast1, adatrj0,  ipi,ipj,ipk,ipni,ipnj,ipnij,iarea
 
     ! Performs checks on the file
@@ -364,17 +364,17 @@
     ! 'before' fields
     DO jk = 1, jpk
-       READ(inum,REC=irec) ub(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) vb(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) tb(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) sb(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) rotb(:,:,jk)   ;   irec = irec +1
+       READ(inum,REC=irec) ub   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) vb   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) tb   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) sb   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) rotb (:,:,jk)   ;   irec = irec +1
     END DO
     DO jk = 1, jpk
@@ -384,17 +384,17 @@
     ! 'now' fields
     DO jk = 1, jpk
-       READ(inum,REC=irec) un(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) vn(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) tn(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) sn(:,:,jk)   ;   irec = irec +1
-    END DO
-    DO jk = 1, jpk
-       READ(inum,REC=irec) rotn(:,:,jk)   ;   irec = irec +1
+       READ(inum,REC=irec) un   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) vn   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) tn   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) sn   (:,:,jk)   ;   irec = irec +1
+    END DO
+    DO jk = 1, jpk
+       READ(inum,REC=irec) rotn (:,:,jk)   ;   irec = irec +1
     END DO
     DO jk = 1, jpk
@@ -403,5 +403,5 @@
 
     ! elliptic solver arrays
-    READ(inum,REC=irec ) gcx(1:jpi,1:jpj)   ;   irec = irec +1
+    READ(inum,REC=irec ) gcx (1:jpi,1:jpj)   ;   irec = irec +1
     READ(inum,REC=irec ) gcxb(1:jpi,1:jpj)   ;   irec = irec +1
 #if defined key_dynspg_rl
@@ -443,11 +443,11 @@
     ! check if it was in the previous run
     IF ( ios1 == 1 ) THEN
-       READ(inum,REC=irec) zfice(:)    ; irec = irec + 1
-       READ(inum,REC=irec) sst_io(:,:) ; irec = irec + 1
-       READ(inum,REC=irec) sss_io(:,:) ; irec = irec + 1
-       READ(inum,REC=irec) u_io  (:,:) ; irec = irec + 1
-       READ(inum,REC=irec) v_io  (:,:) ; irec = irec + 1
+       READ(inum,REC=irec) zfice (:)     ;   irec = irec + 1
+       READ(inum,REC=irec) sst_io(:,:)   ;   irec = irec + 1
+       READ(inum,REC=irec) sss_io(:,:)   ;   irec = irec + 1
+       READ(inum,REC=irec) u_io  (:,:)   ;   irec = irec + 1
+       READ(inum,REC=irec) v_io  (:,:)   ;   irec = irec + 1
 #  if defined key_coupled
-       READ(inum,REC=irec) alb_ice(:,:) ; irec = irec + 1
+       READ(inum,REC=irec) alb_ice(:,:)  ;   irec = irec + 1
 #  endif
     ENDIF
@@ -472,6 +472,6 @@
       ! bulk forcing 
       IF( ios2 == 1 ) THEN
-         READ(inum,REC=irec) zfblk(:)   ; irec = irec + 1
-         READ(inum,REC=irec) gsst (:,:) ; irec = irec + 1
+         READ(inum,REC=irec) zfblk(:)     ;   irec = irec + 1
+         READ(inum,REC=irec) gsst (:,:)   ;   irec = irec + 1
       ENDIF
       IF( zfblk(1) /= FLOAT(nfbulk)  .OR. ios2 == 0 ) THEN
@@ -480,23 +480,22 @@
          IF(lwp) WRITE(numout,*)
          gsst(:,:) = 0.e0
-         gsst(:,:) = gsst(:,:) + ( nfbulk-1 )*( tn(:,:,1) + rt0 )
+         gsst(:,:) = gsst(:,:) + ( nfbulk-1 ) * ( tn(:,:,1) + rt0 )
       ENDIF
 #endif
     CLOSE(inum)
   ! In case of restart with neuler = 0 then put all before fields = to now fields
-    IF ( neuler == 0 ) THEN
-      tb(:,:,:)=tn(:,:,:)
-      sb(:,:,:)=sn(:,:,:)
-      ub(:,:,:)=un(:,:,:)
-      vb(:,:,:)=vn(:,:,:)
-      rotb(:,:,:)=rotn(:,:,:)
-      hdivb(:,:,:)=hdivn(:,:,:)
+    IF( neuler == 0 ) THEN
+       tb(:,:,:) = tn(:,:,:)
+       sb(:,:,:) = sn(:,:,:)
+       ub(:,:,:) = un(:,:,:)
+       vb(:,:,:) = vn(:,:,:)
+       rotb (:,:,:) = rotn (:,:,:)
+       hdivb(:,:,:) = hdivn(:,:,:)
 #if defined key_dynspg_rl
-      bsfb(:,:)=bsfn(:,:)      ! rigid lid
+       bsfb(:,:)=bsfn(:,:)      ! rigid lid
 #else
-      sshb(:,:)=sshn(:,:)      ! free surface formulation (eta)
-#endif
-    ENDIF
-
-
-  END SUBROUTINE rst_read
+       sshb(:,:)=sshn(:,:)      ! free surface formulation (eta)
+#endif
+   ENDIF
+
+END SUBROUTINE rst_read
Index: trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- trunk/NEMO/OPA_SRC/step.F90	(revision 466)
+++ trunk/NEMO/OPA_SRC/step.F90	(revision 467)
@@ -4,5 +4,19 @@
    !! Time-stepping    : manager of the ocean, tracer and ice time stepping
    !!======================================================================
-
+   !! History :
+   !!        !  91-03  ()  Original code
+   !!        !  91-11  (G. Madec)
+   !!        !  92-06  (M. Imbard)  add a first output record
+   !!        !  96-04  (G. Madec)  introduction of dynspg
+   !!        !  96-04  (M.A. Foujols)  introduction of passive tracer
+   !!   8.0  !  97-06  (G. Madec)  new architecture of call
+   !!   8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
+   !!   8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
+   !!   8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
+   !!   9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
+   !!    "   !  04-08  (C. Talandier) New trends organization
+   !!    "   !  05-01  (C. Ethe) Add the KPP closure scheme
+   !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!    "   !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
    !!----------------------------------------------------------------------
    !!   stp            : OPA system time-stepping
@@ -12,6 +26,6 @@
    USE dom_oce         ! ocean space and time domain variables 
    USE zdf_oce         ! ocean vertical physics variables
-   USE ldftra_oce
-   USE ldfdyn_oce
+   USE ldftra_oce      ! ocean tracer   - trends
+   USE ldfdyn_oce      ! ocean dynamics - trends
    USE cpl_oce         ! coupled ocean-atmosphere variables
    USE in_out_manager  ! I/O manager
@@ -34,45 +48,28 @@
    USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
 
-   USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine)
-   USE dynhpg_atsk     ! hydrostatic pressure grad.  (dyn_hpg_atsk routine)
-   USE dynspg_oce      ! surface pressure gradient        (dyn_spg routine)
-   USE dynspg          ! surface pressure gradient        (dyn_spg routine)
-   USE dynkeg          ! kinetic energy gradient          (dyn_keg routine)
-   USE dynvor          ! vorticity term              (dyn_vor_... routines)
-   USE dynzad          ! vertical advection               (dyn_adv routine)
-   USE dynldf_bilapg   ! lateral mixing            (dyn_ldf_bilapg routine)
-   USE dynldf_bilap    ! lateral mixing             (dyn_ldf_bilap routine)
-   USE dynldf_iso      ! lateral mixing               (dyn_ldf_iso routine)
-   USE dynldf_lap      ! lateral mixing               (dyn_ldf_lap routine)
-   USE dynzdf_imp      ! vertical diffusion: implicit     (dyn_zdf routine)
-   USE dynzdf_imp_atsk ! vertical diffusion: implicit     (dyn_zdf routine)
-   USE dynzdf_iso      ! vertical diffusion: isopycnal    (dyn_zdf routine)
-   USE dynzdf_exp      ! vertical diffusion: explicit (dyn_zdf_exp routine)
-   USE dynnxt          ! time-stepping                    (dyn_nxt routine)
-
+   USE traqsr          ! solar radiation penetration      (tra_qsr routine)
+   USE trasbc          ! surface boundary condition       (tra_sbc routine)
    USE trabbc          ! bottom boundary condition        (tra_bbc routine)
    USE trabbl          ! bottom boundary layer            (tra_bbl routine)
    USE tradmp          ! internal damping                 (tra_dmp routine)
-   USE traldf_bilapg   ! lateral mixing            (tra_ldf_bilapg routine)
-   USE traldf_bilap    ! lateral mixing             (tra_ldf_bilap routine)
-   USE traldf_iso      ! lateral mixing               (tra_ldf_iso routine)
-   USE traldf_iso_zps  ! lateral mixing           (tra_ldf_iso_zps routine)
-   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine)
-   USE traqsr          ! solar radiation penetration      (tra_qsr routine)
+   USE traadv          ! advection scheme control     (tra_adv_ctl routine)
+   USE traldf          ! lateral mixing                   (tra_ldf routine)
+   USE cla             ! cross land advection             (tra_cla routine)
+   !   zdfkpp          ! KPP non-local tracer fluxes      (tra_kpp routine)
+   USE trazdf          ! vertical mixing                  (tra_zdf routine)
+   USE tranxt          ! time-stepping                    (tra_nxt routine)
    USE tranpc          ! non-penetrative convection       (tra_npc routine)
-   USE tranxt          ! time-stepping                    (tra_nxt routine)
-   USE traadv_ctl      ! advection scheme control     (tra_adv_ctl routine)
-   USE traadv_cen2     ! 2nd order centered scheme   (tra_adv_cen2 routine)
-   USE traadv_tvd      ! TVD scheme                (tra_adv_tvd    routine)
-   USE traadv_muscl    ! MUSCL scheme              (tra_adv_muscl  routine)
-   USE traadv_muscl2   ! MUSCL2 scheme             (tra_adv_muscl2 routine)
-   USE cla             ! cross land advection             (tra_cla routine)
-   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp routine)
-   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp routine)
-   USE trazdf_iso      ! vertical diffusion           (tra_zdf_exp routine)
-   USE trazdf_iso_vopt ! vertical diffusion           (tra_zdf_exp routine)
-   USE trasbc          ! surface boundary condition       (tra_sbc routine)
 
    USE eosbn2          ! equation of state                (eos_bn2 routine)
+
+   USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine)
+   USE dynkeg          ! kinetic energy gradient          (dyn_keg routine)
+   USE dynvor          ! vorticity term                   (dyn_vor routine)
+   USE dynzad          ! vertical advection               (dyn_adv routine)
+   USE dynldf          ! lateral momentum diffusion       (dyn_ldf routine)
+   USE dynzdf          ! vertical diffusion               (dyn_zdf routine)
+   USE dynspg_oce      ! surface pressure gradient        (dyn_spg routine)
+   USE dynspg          ! surface pressure gradient        (dyn_spg routine)
+   USE dynnxt          ! time-stepping                    (dyn_nxt routine)
 
    USE obc_par         ! open boundary condition variables
@@ -91,4 +88,5 @@
    USE zdfbfr          ! bottom friction                  (zdf_bfr routine)
    USE zdftke          ! TKE vertical mixing              (zdf_tke routine)
+   USE zdftke_jki      ! TKE vertical mixing              (zdf_tke routine)
    USE zdfkpp          ! KPP vertical mixing              (zdf_kpp routine)
    USE zdfddm          ! double diffusion mixing          (zdf_ddm routine)
@@ -140,9 +138,10 @@
 CONTAINS
 
-   SUBROUTINE stp( &
 #if !defined key_agrif
-   kstp &
-#endif   
-   )      !!----------------------------------------------------------------------
+   SUBROUTINE stp( kstp )
+#else
+   SUBROUTINE stp( )
+#endif
+      !!----------------------------------------------------------------------
       !!                     ***  ROUTINE stp  ***
       !!                      
@@ -160,25 +159,11 @@
       !!              -8- Outputs and diagnostics
       !!
-      !! History :
-      !!        !  91-03  ()  Original code
-      !!        !  91-11  (G. Madec)
-      !!        !  92-06  (M. Imbard)  add a first output record
-      !!        !  96-04  (G. Madec)  introduction of dynspg
-      !!        !  96-04  (M.A. Foujols)  introduction of passive tracer
-      !!   8.0  !  97-06  (G. Madec)  new architecture of call
-      !!   8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
-      !!   8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
-      !!   8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
-      !!   9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
-      !!    "   !  04-08  (C. Talandier) New trends organization
-      !!    "   !  05-01  (C. Ethe) Add the KPP closure scheme
-      !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
       !!----------------------------------------------------------------------
       !! * Arguments
-      INTEGER &
 #if !defined key_agrif   
-      , INTENT( in ) &
+      INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
+#else
+      INTEGER               :: kstp   ! ocean time-step index
 #endif      
-      ::   kstp   ! ocean time-step index
 
       !! * local declarations
@@ -206,5 +191,5 @@
       IF( lk_dtasst  )   CALL dta_sst( kstp )         ! Sea Surface Temperature data
 
-      IF( lk_dtasss  )   CALL dta_sss( kstp )         ! Sea Surface salinity data
+      IF( lk_dtasss  )   CALL dta_sss( kstp )         ! Sea Surface Salinity data
 
       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
@@ -258,6 +243,11 @@
       !                                                     ! Vertical eddy viscosity and diffusivity coefficients
       IF( lk_zdfric )   CALL zdf_ric( kstp )                       ! Richardson number dependent Kz
+#if defined key_mpp_omp
+      IF( lk_zdftke )   CALL zdf_tke_jki( kstp )                   ! TKE closure scheme for Kz - j-k-i loops
+#else
       IF( lk_zdftke )   CALL zdf_tke( kstp )                       ! TKE closure scheme for Kz
+#endif
       IF( lk_zdfkpp )   CALL zdf_kpp( kstp )                       ! KPP closure scheme for Kz
+
       IF( lk_zdfcst )   avt (:,:,:) = avt0 * tmask(:,:,:)          ! Constant Kz (reset avt to the background value)
 
@@ -266,4 +256,6 @@
             CASE ( 05 )                         ! ORCA R2 configuration
                avt  (:,:,2) = avt  (:,:,2) + 1.e-3 * upsrnfh(:,:)   ! increase diffusivity of rivers mouths
+            CASE ( 025 )                         ! ORCA R025 configuration
+               avt  (:,:,2) = avt  (:,:,2) + 2.e-3 * upsrnfh(:,:)   ! increase diffusivity of rivers mouths
          END SELECT
       ENDIF
@@ -310,58 +302,45 @@
       !-----------------------------------------------------------------------
 
-                               ta(:,:,:) = 0.e0               ! set tracer trends to zero
-                               sa(:,:,:) = 0.e0
-
-                               CALL tra_sbc( kstp )           ! surface boundary condition
-
-      IF( ln_traqsr        )   CALL tra_qsr( kstp )           ! penetrative solar radiation qsr
-
-      IF( lk_trabbc        )   CALL tra_bbc( kstp )           ! bottom heat flux
-
-      IF( lk_trabbl_dif    )   CALL tra_bbl_dif( kstp )           ! diffusive bottom boundary layer scheme
-      IF( lk_trabbl_adv    )   CALL tra_bbl_adv( kstp )           ! advective (and/or diffusive) bottom boundary layer scheme
-
-      IF( lk_tradmp        )   CALL tra_dmp( kstp )           ! internal damping trends
-
-      !                                                       ! horizontal & vertical advection
-      IF( kstp == nit000   )   CALL tra_adv_ctl                    ! chose/control the scheme used
-      IF( ln_traadv_cen2   )   CALL tra_adv_cen2  ( kstp )         ! 2nd order centered scheme
-      IF( ln_traadv_tvd    )   CALL tra_adv_tvd   ( kstp )         ! TVD scheme
-      IF( ln_traadv_muscl  )   CALL tra_adv_muscl ( kstp )         ! MUSCL scheme
-      IF( ln_traadv_muscl2 )   CALL tra_adv_muscl2( kstp )         ! MUSCL2 scheme
-
-      IF( n_cla == 1       )   CALL tra_cla( kstp )           ! Cross Land Advection (Update Hor. advection)
-
-      !                                                       ! lateral mixing 
-      IF( l_traldf_lap     )   CALL tra_ldf_lap    ( kstp )           ! iso-level laplacian
-      IF( l_traldf_bilap   )   CALL tra_ldf_bilap  ( kstp )           ! iso-level bilaplacian 
-      IF( l_traldf_bilapg  )   CALL tra_ldf_bilapg ( kstp )           ! s-coord. horizontal bilaplacian
-      IF( l_traldf_iso     )   CALL tra_ldf_iso    ( kstp )           ! iso-neutral/geopot. laplacian 
-      IF( l_traldf_iso_zps )   CALL tra_ldf_iso_zps( kstp )           ! partial step iso-neutral/geopot. laplacian
-
+                             ta(:,:,:) = 0.e0               ! set tracer trends to zero
+                             sa(:,:,:) = 0.e0
+
+                             CALL tra_sbc    ( kstp )       ! surface boundary condition
+
+      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr
+
+      IF( lk_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux
+
+      IF( lk_trabbl_dif  )   CALL tra_bbl_dif( kstp )       ! diffusive bottom boundary layer scheme
+      IF( lk_trabbl_adv  )   CALL tra_bbl_adv( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
+
+      IF( lk_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
+
+                             CALL tra_adv    ( kstp )       ! horizontal & vertical advection
+
+      IF( n_cla == 1     )   CALL tra_cla    ( kstp )       ! Cross Land Advection (Update Hor. advection)
+
+      IF( lk_zdfkpp )        CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
+
+                             CALL tra_ldf    ( kstp )       ! lateral mixing
 #if defined key_agrif
       IF (.NOT. Agrif_Root())  CALL Agrif_Sponge_tra( kstp )          ! tracers sponge
 #endif
-      !                                                       ! vertical diffusion
-      IF( l_trazdf_exp     )   CALL tra_zdf_exp     ( kstp )          ! explicit time stepping (time splitting scheme)
-      IF( l_trazdf_imp     )   CALL tra_zdf_imp     ( kstp )          ! implicit time stepping (euler backward)
-      IF( l_trazdf_iso     )   CALL tra_zdf_iso     ( kstp )          ! isopycnal
-      IF( l_trazdf_iso_vo  )   CALL tra_zdf_iso_vopt( kstp )          ! vector opt. isopycnal
-
-                               CALL tra_nxt( kstp )           ! tracer fields at next time step
-
-      IF( ln_zdfnpc        )   CALL tra_npc( kstp )           ! update the new (t,s) fields by non
-      !                                                       ! penetrative convective adjustment
-
-      IF( ln_dynhpg_imp    ) THEN                             ! semi-implicit hpg 
-                                  CALL eos( ta, sa, rhd, rhop )   ! Time-filtered in situ density used in dynhpg module
-         IF( lk_zps    )          CALL zps_hde( kstp, ta, sa, rhd,  & ! Partial steps: time filtered hor. gradient 
-            &                                        gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level
-            &                                        gtv, gsv, grv )  
-      ELSE                                                    ! centered hpg (default case)
-                                  CALL eos( tb, sb, rhd, rhop )       ! now (swap=before) in situ density for dynhpg module
-         IF( lk_zps    )          CALL zps_hde( kstp, tb, sb, rhd,  & ! Partial steps: now horizontal gradient
-            &                                        gtu, gsu, gru, & ! of t, s, rd at the bottom ocean level
-            &                                        gtv, gsv, grv )  
+                             CALL tra_zdf    ( kstp )       ! vertical mixing
+
+                             CALL tra_nxt( kstp )           ! tracer fields at next time step
+
+      IF( ln_zdfnpc      )   CALL tra_npc( kstp )           ! update the new (t,s) fields by non
+      !                                                     ! penetrative convective adjustment
+
+      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
+                               CALL eos( ta, sa, rhd, rhop )          ! Time-filtered in situ density used in dynhpg module
+         IF( ln_zps    )       CALL zps_hde( kstp, ta, sa, rhd,   &   ! Partial steps: time filtered hor. gradient
+            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
+            &                                     gtv, gsv, grv )
+      ELSE                                                  ! centered hpg (default case)
+                               CALL eos( tb, sb, rhd, rhop )          ! now (swap=before) in situ density for dynhpg module
+         IF( ln_zps    )       CALL zps_hde( kstp, tb, sb, rhd,   &   ! Partial steps: now horizontal gradient
+            &                                     gtu, gsu, gru,  &   ! of t, s, rd at the bottom ocean level
+            &                                     gtv, gsv, grv )
       ENDIF
 
@@ -371,4 +350,5 @@
       ! N.B. ta, sa arrays are used as workspace in this section 
       !-----------------------------------------------------------------------
+
 
                                ua(:,:,:) = 0.e0               ! set dynamics trends to zero
@@ -377,36 +357,17 @@
                                CALL dyn_keg( kstp )           ! horizontal gradient of kinetic energy
 
-      !                                                       ! vorticity term including Coriolis
-      IF( kstp == nit000   )   CALL dyn_vor_ctl                      ! chose/control the scheme used
-      IF( ln_dynvor_ens    )   CALL dyn_vor_enstrophy( kstp )        ! enstrophy conserving scheme
-      IF( ln_dynvor_ene    )   CALL dyn_vor_energy   ( kstp )        ! energy conserving scheme
-      IF( ln_dynvor_mix    )   CALL dyn_vor_mixed    ( kstp )        ! mixed energy/enstrophy conserving scheme
-      IF( ln_dynvor_een    )   CALL dyn_vor_ene_ens  ( kstp )        ! combined energy/enstrophy conserving scheme
-      
-      !                                                       ! lateral mixing 
-      IF( l_dynldf_lap     )   CALL dyn_ldf_lap    ( kstp )          ! iso-level laplacian
-      IF( l_dynldf_bilap   )   CALL dyn_ldf_bilap  ( kstp )          ! iso-level bilaplacian 
-      IF( l_dynldf_bilapg  )   CALL dyn_ldf_bilapg ( kstp )          ! s-coord. horizontal bilaplacian
-      IF( l_dynldf_iso     )   CALL dyn_ldf_iso    ( kstp )          ! iso-neutral laplacian 
-
+                               CALL dyn_vor( kstp )           ! vorticity term including Coriolis
+
+                               CALL dyn_ldf( kstp )           ! lateral mixing
 #if defined key_agrif
       IF (.NOT. Agrif_Root())  CALL Agrif_Sponge_dyn( kstp )         ! momemtum sponge
 #endif
-      !                                                       ! horizontal gradient of Hydrostatic pressure 
-      IF ( lk_jki ) THEN
-                               CALL dyn_hpg_atsk( kstp )             ! autotask case (j-k-i loop)
-      ELSE
-                               CALL dyn_hpg     ( kstp )             ! default case  (k-j-i loop)
-      ENDIF
-
-                               CALL dyn_zad    ( kstp )       ! vertical advection       
-
-      !                                                       ! vertical diffusion
-      IF( l_dynzdf_exp     )   CALL dyn_zdf_exp    ( kstp )          ! explicit time stepping (time splitting scheme)
-      IF( l_dynzdf_imp     )   CALL dyn_zdf_imp    ( kstp )          ! implicit time stepping (euler backward)
-      IF( l_dynzdf_imp_tsk )   CALL dyn_zdf_imp_tsk( kstp )          ! autotask implicit time stepping (euler backward)
-      IF( l_dynzdf_iso     )   CALL dyn_zdf_iso    ( kstp )          ! iso-neutral case
-
-      IF( lk_dynspg_rl ) THEN 
+                               CALL dyn_hpg( kstp )           ! horizontal gradient of Hydrostatic pressure
+
+                               CALL dyn_zad( kstp )           ! vertical advection
+
+                               CALL dyn_zdf( kstp )           ! vertical diffusion
+
+      IF( lk_dynspg_rl ) THEN
          IF( lk_obc    )       CALL obc_spg( kstp )           ! surface pressure gradient at open boundaries
       ENDIF
@@ -417,5 +378,5 @@
                                CALL dyn_spg( kstp, indic )    ! surface pressure gradient
 
-                               CALL dyn_nxt( kstp )           ! velocity at next time step 
+                               CALL dyn_nxt( kstp )           ! lateral velocity at next time step
 
 
