Index: trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 2528)
@@ -2,6 +2,16 @@
    !!======================================================================
    !!                       ***  MODULE  dtadyn  ***
-   !! OFFLINE : interpolation of the physical fields
-   !!=====================================================================
+   !! Off-line : interpolation of the physical fields
+   !!======================================================================
+   !! History :   OPA  ! 1992-01 (M. Imbard) Original code
+   !!             8.0  ! 1998-04 (L.Bopp MA Foujols) slopes for isopyc. 
+   !!              -   ! 1998-05 (L. Bopp) read output of coupled run
+   !!             8.2  ! 2001-01 (M. Levy et M. Benjelloul) add netcdf FORMAT
+   !!   NEMO      1.0  ! 2005-03 (O. Aumont and A. El Moussaoui) F90
+   !!              -   ! 2005-12 (C. Ethe) Adapted for DEGINT
+   !!             3.0  ! 2007-06 (C. Ethe) use of iom module
+   !!              -   ! 2007-09  (C. Ethe)  add swap_dyn_data
+   !!             3.3  ! 2010-11 (C. Ethe) Full reorganization of the off-line: phasing with the on-line
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,104 +19,77 @@
    !!   dta_dyn      : Interpolation of the fields
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables
-   USE zdf_oce         ! ocean vertical physics
+   USE c1d             ! 1D configuration: lk_c1d
+   USE dom_oce         ! ocean domain: variables
+   USE zdf_oce         ! ocean vertical physics: variables
+   USE sbc_oce         ! surface module: variables
+   USE phycst          ! physical constants
+   USE trabbl          ! active tracer: bottom boundary layer
+   USE ldfslp          ! lateral diffusion: iso-neutral slopes
+   USE ldfeiv          ! eddy induced velocity coef. 
+   USE ldftra_oce      ! ocean tracer   lateral physics
+   USE zdfmxl          ! vertical physics: mixed layer depth
+   USE eosbn2          ! equation of state - Brunt Vaisala frequency
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE zpshde          ! z-coord. with partial steps: horizontal derivatives
    USE in_out_manager  ! I/O manager
-   USE phycst          ! physical constants
-   USE sbc_oce
-   USE ldfslp
-   USE ldfeiv          ! eddy induced velocity coef.      (ldf_eiv routine)
-   USE ldftra_oce      ! ocean tracer   lateral physics
-   USE zdfmxl
-   USE trabbl
-   USE eosbn2
-   USE zdfddm          ! vertical  physics: double diffusion
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE zpshde
+   USE iom             ! I/O library
    USE lib_mpp         ! distributed memory computing library
+   USE prtctl          !  print control
 
    IMPLICIT NONE
    PRIVATE
 
-   !! *  Routine accessibility
-   PUBLIC dta_dyn_init   ! called by opa.F90
-   PUBLIC dta_dyn        ! called by step.F90
-
-   LOGICAL , PUBLIC :: &
-      lperdyn = .TRUE. , & ! boolean for periodic fields or not
-      lfirdyn = .TRUE.     ! boolean for the first call or not
-
-   INTEGER , PUBLIC :: &
-      ndtadyn = 73 ,  & ! Number of dat in one year
-      ndtatot = 73 ,  & ! Number of data in the input field
-      nsptint = 1 ,   & ! type of spatial interpolation
-      nficdyn = 2       ! number of dynamical fields 
-
-   CHARACTER(len=45)  ::  &
-      cfile_grid_T = 'dyna_grid_T.nc', &   !: name of the grid_T file
-      cfile_grid_U = 'dyna_grid_U.nc', &   !: name of the grid_U file
-      cfile_grid_V = 'dyna_grid_V.nc', &   !: name of the grid_V file
-      cfile_grid_W = 'dyna_grid_W.nc'      !: name of the grid_W file
+   PUBLIC   dta_dyn_init   ! called by opa.F90
+   PUBLIC   dta_dyn        ! called by step.F90
+
+   LOGICAL, PUBLIC ::   lperdyn = .TRUE.   !: boolean for periodic fields or not
+   LOGICAL, PUBLIC ::   lfirdyn = .TRUE.   !: boolean for the first call or not
+
+   INTEGER, PUBLIC ::   ndtadyn = 73       !: Number of dat in one year
+   INTEGER, PUBLIC ::   ndtatot = 73       !: Number of data in the input field
+   INTEGER, PUBLIC ::   nsptint = 1        !: type of spatial interpolation
+
+   CHARACTER(len=45) ::   cfile_grid_T = 'dyna_grid_T.nc'   ! name of the grid_T file
+   CHARACTER(len=45) ::   cfile_grid_U = 'dyna_grid_U.nc'   ! name of the grid_U file
+   CHARACTER(len=45) ::   cfile_grid_V = 'dyna_grid_V.nc'   ! name of the grid_V file
+   CHARACTER(len=45) ::   cfile_grid_W = 'dyna_grid_W.nc'   ! name of the grid_W file
    
-   REAL(wp)      ::   &
-      rnspdta  ,       &  !: number of time step per 2 consecutives data
-      rnspdta2            !: rnspdta * 0.5
-
-   INTEGER ::     &
-      ndyn1, ndyn2 , &
-      nlecoff = 0  , & ! switch for the first read
-      numfl_t, numfl_u, &
-      numfl_v, numfl_w
-
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   &
-      tdta   ,   & ! temperature at two consecutive times
-      sdta   ,   & ! salinity at two consecutive times
-      udta   ,   & ! zonal velocity at two consecutive times
-      vdta   ,   & ! meridional velocity at two consecutive times
-      wdta   ,   & ! vertical velocity at two consecutive times
-#if defined key_trc_diatrd
-      hdivdta,   & ! horizontal divergence
-#endif
-      avtdta       ! vertical diffusivity coefficient
-
-   REAL(wp), DIMENSION(jpi,jpj,2) ::       &
-      hmlddta,   & ! mixed layer depth at two consecutive times
-      wspddta,   & ! wind speed at two consecutive times
-      frlddta,   & ! sea-ice fraction at two consecutive times
-      empdta ,   & ! E-P at two consecutive times
-      qsrdta       ! short wave heat flux at two consecutive times
-
+   REAL(wp) ::   rnspdta    ! number of time step per 2 consecutives data
+   REAL(wp) ::   rnspdta2   ! rnspdta * 0.5
+
+   INTEGER ::   ndyn1, ndyn2    !
+   INTEGER ::   nlecoff = 0     ! switch for the first read
+   INTEGER ::   numfl_t, numfl_u, numfl_v, numfl_w
+
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: tdta       ! temperature at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: sdta       ! salinity at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: udta       ! zonal velocity at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: vdta       ! meridional velocity at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wdta       ! vertical velocity at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: avtdta     ! vertical diffusivity coefficient
+
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: hmlddta    ! mixed layer depth at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: wspddta    ! wind speed at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: frlddta    ! sea-ice fraction at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: empdta     ! E-P at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: qsrdta     ! short wave heat flux at two consecutive times
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: bblxdta    ! frequency of bbl in the x direction at 2 consecutive times 
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: bblydta    ! frequency of bbl in the y direction at 2 consecutive times 
+   LOGICAL :: l_offbbl
 #if defined key_ldfslp
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   &
-      uslpdta ,  & ! zonal isopycnal slopes
-      vslpdta ,  & ! meridional isopycnal slopes
-      wslpidta , & ! zonal diapycnal slopes
-      wslpjdta     ! meridional diapycnal slopes
-#endif
-
-#if ! defined key_off_degrad &&  defined key_traldf_c2d
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &
-      ahtwdta      ! Lateral diffusivity
-# if defined key_trcldf_eiv 
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   &
-      aeiwdta      ! G&M coefficient
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: uslpdta    ! zonal isopycnal slopes
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: vslpdta    ! meridional isopycnal slopes
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wslpidta   ! zonal diapycnal slopes
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wslpjdta   ! meridional diapycnal slopes
+#endif
+#if ! defined key_degrad &&  defined key_traldf_c2d && defined key_traldf_eiv 
+   REAL(wp), DIMENSION(jpi,jpj    ,2) :: aeiwdta    ! G&M coefficient
+#endif
+#if defined key_degrad
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ahtudta, ahtvdta, ahtwdta   ! Lateral diffusivity
+# if defined key_traldf_eiv
+   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: aeiudta, aeivdta, aeiwdta   ! G&M coefficient
 # endif
-#endif
-
-#if defined key_off_degrad
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   &
-      ahtudta, ahtvdta, ahtwdta  !  Lateral diffusivity
-# if defined key_trcldf_eiv
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   &
-      aeiudta, aeivdta, aeiwdta  ! G&M coefficient
-# endif
-
-#endif
-
-#if defined key_trcbbl_dif   ||   defined key_trcbbl_adv
-   REAL(wp), DIMENSION(jpi,jpj,2) ::       &
-      bblxdta ,  & ! frequency of bbl in the x direction at 2 consecutive times
-      bblydta      ! frequency of bbl in the y direction at 2 consecutive times
 #endif
 
@@ -115,9 +98,8 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL  (2005)
-   !!   $Id$
-   !!   This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! NEMO/OFF 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -126,51 +108,23 @@
       !!                  ***  ROUTINE dta_dyn  ***
       !!
-      !! ** Purpose : Prepares dynamics and physics fields from an 
-      !!              OPA9 simulation  for an off-line simulation
-      !!               for passive tracer
+      !! ** Purpose :   Prepares dynamics and physics fields from an NEMO run
+      !!              for an off-line simulation of passive tracers
       !!
       !! ** Method : calculates the position of DATA to read READ DATA 
       !!             (example month changement) computes slopes IF needed
       !!             interpolates DATA IF needed
-      !!
-      !! ** History :
-      !!   ! original  : 92-01 (M. Imbard: sub domain)
-      !!   ! addition  : 98-04 (L.Bopp MA Foujols: slopes for isopyc.) 
-      !!   ! addition  : 98-05 (L. Bopp read output of coupled run)
-      !!   ! addition  : 05-03 (O. Aumont and A. El Moussaoui) F90
-      !!   ! addition  : 05-12 (C. Ethe) Adapted for DEGINT
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER ::   iper, iperm1, iswap, izt   
-
-      REAL(wp) :: zpdtan, zpdtpe, zdemi, zt
-      REAL(wp) :: zweigh
-
-      ! 0. Initialization
-      ! -----------------
-
-      IF( lfirdyn ) THEN
-         ! first time step MUST BE nit000
-         IF( kt /= nit000 ) THEN
-            IF (lwp) THEN 
-               WRITE (numout,*) ' kt MUST BE EQUAL to nit000. kt = ',kt ,' nit000 = ',nit000 
-              STOP 'dtadyn'
-            ENDIF
-          ENDIF 
-          ! Initialize the parameters of the interpolation
-          CALL dta_dyn_init
-      ENDIF
-
-      zt       = ( FLOAT (kt) + rnspdta2 ) / rnspdta
-      izt      = INT( zt )
-      zweigh   = zt - FLOAT( INT(zt) )
-
-      IF( lperdyn ) THEN
-         iperm1 = MOD( izt, ndtadyn )
-      ELSE
-         iperm1 = MOD( izt, ndtatot - 1 ) + 1
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   iper, iperm1, iswap, izt   ! local integers 
+      REAL(wp) ::   zt, zweigh                 ! local scalars
+      !!----------------------------------------------------------------------
+
+      zt     = ( REAL(kt,wp) + rnspdta2 ) / rnspdta
+      izt    = INT( zt )
+      zweigh = zt - REAL( INT(zt), wp )
+
+      IF( lperdyn ) THEN   ;   iperm1 = MOD( izt, ndtadyn )
+      ELSE                 ;   iperm1 = MOD( izt, ndtatot - 1 ) + 1
       ENDIF
 
@@ -181,7 +135,6 @@
           ELSE 
               IF( lfirdyn ) THEN 
-                  IF (lwp) WRITE (numout,*) & 
-                      &   ' dynamic file is not periodic with or without interpolation  &
-                      &   we take the first value for the previous period iperm1 = 0  '
+                  IF(lwp) WRITE (numout,*) 'dta_dyn:  dynamic file is not periodic with or without interpolation    &
+                     &                                we take the first value for the previous period iperm1 = 0  '
               END IF
           END IF 
@@ -194,73 +147,60 @@
 
       IF( lfirdyn ) THEN
-         ! store the information of the period read
-         ndyn1 = iperm1
+         ndyn1 = iperm1         ! store the information of the period read
          ndyn2 = iper
          
-         IF (lwp) THEN
-            WRITE (numout,*) ' dynamics data read for the period ndyn1 =',ndyn1, &
-               &             ' and for the period ndyn2 = ',ndyn2
+         IF(lwp) THEN
+            WRITE (numout,*) ' dynamics data read for the period ndyn1 =', ndyn1,   &
+               &             ' and for the period ndyn2 = ', ndyn2
             WRITE (numout,*) ' time step is : ', kt
-            WRITE (numout,*) ' we have ndtadyn = ',ndtadyn,' records in the dynamic file for one year'
+            WRITE (numout,*) ' we have ndtadyn = ', ndtadyn, ' records in the dynamic file for one year'
          END IF
          !
-         IF( iperm1 /= 0 ) THEN         ! data read for the iperm1 period
-            CALL dynrea( kt, iperm1 ) 
-         ELSE 
-            CALL dynrea( kt, 1 )
-         ENDIF
+         CALL dynrea( kt, MAX( 1, iperm1) )           ! data read for the iperm1 period
          
-#if defined key_ldfslp
-         ! Computes slopes
-         ! Caution : here tn, sn and avt are used as workspace
-         tn (:,:,:) = tdta  (:,:,:,2)
-         sn (:,:,:) = sdta  (:,:,:,2)
-         avt(:,:,:) = avtdta(:,:,:,2)
+         IF( lk_ldfslp .AND. .NOT. lk_c1d ) THEN      ! Computes slopes (here tsn and avt are used as workspace)
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
          
-         CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-         CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
-         IF( ln_zps )   &
-            &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-            &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-            &                 gtv, gsv, grv )
-         CALL zdf_mxl( kt )              ! mixed layer depth
-         CALL ldf_slp( kt, rhd, rn2 )
+            CALL eos( tsn, rhd, rhop )   ! Time-filtered in situ density 
+            CALL bn2( tsn, rn2 )         ! before Brunt-Vaisala frequency
+            IF( ln_zps )   &
+               &   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                           rhd, gru , grv   )    ! of t, s, rd at the bottom ocean level
+            CALL zdf_mxl( kt )           ! mixed layer depth
+            CALL ldf_slp( kt, rhd, rn2 )
          
-         uslpdta (:,:,:,2) = uslp (:,:,:)
-         vslpdta (:,:,:,2) = vslp (:,:,:)
-         wslpidta(:,:,:,2) = wslpi(:,:,:)
-         wslpjdta(:,:,:,2) = wslpj(:,:,:)
-#endif
-         
-         ! swap from record 2 to 1
-         CALL swap_dyn_data
-         
+            uslpdta (:,:,:,2) = uslp (:,:,:)
+            vslpdta (:,:,:,2) = vslp (:,:,:)
+            wslpidta(:,:,:,2) = wslpi(:,:,:)
+            wslpjdta(:,:,:,2) = wslpj(:,:,:)
+         END IF
+         !
+         CALL swap_dyn_data            ! swap from record 2 to 1
+         !
          iswap = 1        !  indicates swap
-         
-         CALL dynrea( kt, iper )    ! data read for the iper period
-         
-#if defined key_ldfslp
-         ! Computes slopes
-         ! Caution : here tn, sn and avt are used as workspace
-         tn (:,:,:) = tdta  (:,:,:,2)
-         sn (:,:,:) = sdta  (:,:,:,2)
-         avt(:,:,:) = avtdta(:,:,:,2)
-         
-         CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-         CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
-         IF( ln_zps )   &
-            &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-            &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-            &                 gtv, gsv, grv )
-         CALL zdf_mxl( kt )              ! mixed layer depth
-         CALL ldf_slp( kt, rhd, rn2 )
-         
-         uslpdta (:,:,:,2) = uslp (:,:,:)
-         vslpdta (:,:,:,2) = vslp (:,:,:)
-         wslpidta(:,:,:,2) = wslpi(:,:,:)
-         wslpjdta(:,:,:,2) = wslpj(:,:,:)
-#endif
-         !
-         lfirdyn=.FALSE.    ! trace the first call
+         !
+         CALL dynrea( kt, iper )       ! data read for the iper period
+         !
+         IF( lk_ldfslp .AND. .NOT. lk_c1d ) THEN      ! Computes slopes (here tsn and avt are used as workspace)
+            tsn (:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn (:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)         = avtdta(:,:,:,2)
+            !
+                           CALL eos( tsn, rhd, rhop )                   ! now in situ density 
+                           CALL bn2( tsn, rn2 )                         ! now Brunt-Vaisala frequency
+            IF( ln_zps )   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                                   rhd, gru , grv   )   ! of t, s, rd at the bottom ocean level
+                           CALL zdf_mxl( kt )                           ! mixed layer depth
+                           CALL ldf_slp( kt, rhd, rn2 )                 ! slope of iso-neutral surfaces
+            !
+            uslpdta (:,:,:,2) = uslp (:,:,:)
+            vslpdta (:,:,:,2) = vslp (:,:,:)
+            wslpidta(:,:,:,2) = wslpi(:,:,:)
+            wslpjdta(:,:,:,2) = wslpj(:,:,:)
+         END IF
+         !
+         lfirdyn = .FALSE.    ! trace the first call
       ENDIF
       !
@@ -269,7 +209,7 @@
       !
       IF( iperm1 /= ndyn1 ) THEN 
-
-         IF( iperm1 == 0. ) THEN
-            IF (lwp) THEN
+         !
+         IF( iperm1 == 0 ) THEN
+            IF(lwp) THEN
                WRITE (numout,*) ' dynamic file is not periodic with periodic interpolation'
                WRITE (numout,*) ' we take the last value for the last period '
@@ -280,41 +220,35 @@
          ENDIF
          !
-         ! We have to prepare a new read of data : swap from record 2 to 1
-         !
-         CALL swap_dyn_data
-
-         iswap = 1        !  indicates swap
-         
+         CALL swap_dyn_data         ! We have to prepare a new read of data : swap from record 2 to 1
+         !
+         iswap = 1                  !  indicates swap
+         !
          CALL dynrea( kt, iper )    ! data read for the iper period
-
-#if defined key_ldfslp
-         ! Computes slopes
-         ! Caution : here tn, sn and avt are used as workspace
-         tn (:,:,:) = tdta  (:,:,:,2)
-         sn (:,:,:) = sdta  (:,:,:,2)
-         avt(:,:,:) = avtdta(:,:,:,2)
-         
-         CALL eos( tn, sn, rhd, rhop )   ! Time-filtered in situ density 
-         CALL bn2( tn, sn, rn2 )         ! before Brunt-Vaisala frequency
-         IF( ln_zps )   &
-            &   CALL zps_hde( kt, tn , sn , rhd,  &  ! Partial steps: before Horizontal DErivative
-            &                 gtu, gsu, gru,  &  ! of t, s, rd at the bottom ocean level
-            &                 gtv, gsv, grv )
-         CALL zdf_mxl( kt )              ! mixed layer depth
-         CALL ldf_slp( kt, rhd, rn2 )
-         
-         uslpdta (:,:,:,2) = uslp (:,:,:)
-         vslpdta (:,:,:,2) = vslp (:,:,:)
-         wslpidta(:,:,:,2) = wslpi(:,:,:)
-         wslpjdta(:,:,:,2) = wslpj(:,:,:)
-#endif
-       
-         ! store the information of the period read
-         ndyn1 = ndyn2
+         !
+         IF( lk_ldfslp .AND. .NOT. lk_c1d ) THEN
+            ! Computes slopes. Caution : here tsn and avt are used as workspace
+            tsn(:,:,:,jp_tem) = tdta  (:,:,:,2)
+            tsn(:,:,:,jp_sal) = sdta  (:,:,:,2)
+            avt(:,:,:)        = avtdta(:,:,:,2)
+            !
+                           CALL eos( tsn, rhd, rhop )                   ! now in situ density 
+                           CALL bn2( tsn, rn2 )                         ! now Brunt-Vaisala frequency
+            IF( ln_zps )   CALL zps_hde( kt, jpts, tsn, gtsu, gtsv,  &  ! Partial steps: before Horizontal DErivative
+               &                                   rhd, gru , grv   )   ! of t, s, rd at the bottom ocean level
+            CALL zdf_mxl( kt )                                          ! mixed layer depth
+            CALL ldf_slp( kt, rhd, rn2 )                                ! slope of iso-neutral surfaces
+            !
+            uslpdta (:,:,:,2) = uslp (:,:,:)
+            vslpdta (:,:,:,2) = vslp (:,:,:)
+            wslpidta(:,:,:,2) = wslpi(:,:,:)
+            wslpjdta(:,:,:,2) = wslpj(:,:,:)
+         END IF
+         !
+         ndyn1 = ndyn2         ! store the information of the period read
          ndyn2 = iper
-
-         IF (lwp) THEN
-            WRITE (numout,*) ' dynamics data read for the period ndyn1 =',ndyn1, &
-               &             ' and for the period ndyn2 = ',ndyn2
+         !
+         IF(lwp) THEN
+            WRITE (numout,*) ' dynamics data read for the period ndyn1 =', ndyn1,   &
+               &             ' and for the period ndyn2 = ', ndyn2
             WRITE (numout,*) ' time step is : ', kt
          END IF
@@ -325,26 +259,41 @@
       !----------------------------------------     
 
-      IF( nsptint == 0 ) THEN   
-         ! No spatial interpolation, data are probably correct
-         ! We have to initialize data if we have changed the period         
-         CALL assign_dyn_data          
-      ELSE IF( nsptint == 1 ) THEN
-         ! linear interpolation
+      IF( nsptint == 0 ) THEN          ! No space interpolation, data are probably correct
+         !                             ! We have to initialize data if we have changed the period         
+         CALL assign_dyn_data
+      ELSEIF( nsptint == 1 ) THEN      ! linear interpolation
          CALL linear_interp_dyn_data( zweigh )
-      ELSE 
-         ! other interpolation
+      ELSE                             ! other interpolation
          WRITE (numout,*) ' this kind of interpolation do not exist at the moment : we stop'
          STOP 'dtadyn'         
       END IF
-      
-      ! In any case, we need rhop
-      CALL eos( tn, sn, rhd, rhop ) 
-      
-#if ! defined key_off_degrad && defined key_traldf_c2d
-      ! In case of 2D varying coefficients, we need aeiv and aeiu
-      IF( lk_traldf_eiv )   CALL ldf_eiv( kt )      ! eddy induced velocity coefficient
-#endif
-   
+      !
+      CALL eos( tsn, rhd, rhop )       ! In any case, we need rhop
+      !
+#if ! defined key_degrad && defined key_traldf_c2d
+      !                                ! In case of 2D varying coefficients, we need aeiv and aeiu
+      IF( lk_traldf_eiv )   CALL dta_eiv( kt )      ! eddy induced velocity coefficient
+#endif
+      !
+      IF( .NOT. l_offbbl ) THEN       ! Compute bbl coefficients if needed
+         tsb(:,:,:,:) = tsn(:,:,:,:)
+         CALL bbl( kt, 'TRC')
+      END IF
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
+         CALL prt_ctl(tab3d_1=un               , clinfo1=' un      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
+         CALL prt_ctl(tab3d_1=vn               , clinfo1=' vn      - : ', mask1=tmask, ovlap=1, kdim=jpk   )
+         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=emps             , clinfo1=' emps    - : ', 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
+      !
    END SUBROUTINE dta_dyn
+
 
    SUBROUTINE dynrea( kt, kenr )
@@ -354,48 +303,24 @@
       !! ** Purpose : READ dynamics fiels from OPA9 netcdf output
       !! 
-      !! ** Method : READ the kenr records of DATA and store in
-      !!             in udta(...,2), ....  
-      !! 
-      !! ** History : additions : M. Levy et M. Benjelloul jan 2001 
-      !!              (netcdf FORMAT) 
-      !!              05-03 (O. Aumont and A. El Moussaoui) F90
-      !!              06-07 : (C. Ethe) use of iom module
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE iom
-
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt, kenr       ! time index
-      !! * Local declarations
+      !! ** Method : READ the kenr records of DATA and store in udta(...,2), ....  
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt, kenr   ! time index
+      !!
       INTEGER ::  jkenr
-
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-        zu, zv, zw, zt, zs, zavt ,   &     ! 3-D dynamical fields
-        zhdiv                              ! horizontal divergence
-
-      REAL(wp), DIMENSION(jpi,jpj) :: &
-         zemp, zqsr, zmld, zice, zwspd, &
-         ztaux, ztauy
-#if defined key_trcbbl_dif   ||   defined key_trcbbl_adv
-      REAL(wp), DIMENSION(jpi,jpj) :: zbblx, zbbly
-#endif
-
-#if ! defined key_off_degrad && defined key_traldf_c2d
-      REAL(wp), DIMENSION(jpi,jpj) :: zahtw 
-#   if defined key_trcldf_eiv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zu, zv, zw, zt, zs, zavt , zhdiv              ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::  zemp, zqsr, zmld, zice, zwspd, ztaux, ztauy   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::  zbblx, zbbly
+
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
       REAL(wp), DIMENSION(jpi,jpj) :: zaeiw 
-#  endif
-#endif
-
-#if defined key_off_degrad
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-      zahtu, zahtv, zahtw  !  Lateral diffusivity
-# if defined key_trcldf_eiv
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
-      zaeiu, zaeiv, zaeiw  ! G&M coefficient
+#endif
+#if defined key_degrad
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zahtu, zahtv, zahtw  !  Lateral diffusivity
+# if defined key_traldf_eiv
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zaeiu, zaeiv, zaeiw  ! G&M coefficient
 # endif
 #endif
-
-      !---------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
       ! 0. Initialization
       
@@ -407,7 +332,7 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'Dynrea : reading dynamical fields, kenr = ', jkenr
-         WRITE(numout,*) ' ~~~~~~~'
-#if defined key_off_degrad
+         WRITE(numout,*) 'Dynrea : read dynamical fields, kenr = ', jkenr
+         WRITE(numout,*) '~~~~~~~'
+#if defined key_degrad
          WRITE(numout,*) ' Degraded fields'
 #endif
@@ -442,14 +367,12 @@
       CALL iom_get( numfl_u, jpdom_data, 'vozocrtx', zu   (:,:,:), jkenr )
       CALL iom_get( numfl_v, jpdom_data, 'vomecrty', zv   (:,:,:), jkenr )
-
-#if defined key_trcbbl_dif || defined key_trcbbl_adv
-      IF( iom_varid( numfl_u, 'sobblcox', ldstop = .FALSE. ) > 0  .AND. &
-      &   iom_varid( numfl_v, 'sobblcoy', ldstop = .FALSE. ) > 0 ) THEN
-         CALL iom_get( numfl_u, jpdom_data, 'sobblcox', zbblx(:,:), jkenr )
-         CALL iom_get( numfl_v, jpdom_data, 'sobblcoy', zbbly(:,:), jkenr )
-      ELSE
-         CALL bbl_sign( zt, zs, zbblx, zbbly )    
-      ENDIF
-#endif
+      IF( lk_trabbl .AND. .NOT. lk_c1d .AND. nn_bbl_ldf == 1 ) THEN
+         IF( iom_varid( numfl_u, 'sobblcox', ldstop = .FALSE. ) > 0  .AND. &
+         &   iom_varid( numfl_v, 'sobblcoy', ldstop = .FALSE. ) > 0 ) THEN
+             CALL iom_get( numfl_u, jpdom_data, 'sobblcox', zbblx(:,:), jkenr )
+             CALL iom_get( numfl_v, jpdom_data, 'sobblcoy', zbbly(:,:), jkenr )
+             l_offbbl = .TRUE.
+         ENDIF
+      ENDIF
 
       ! file grid-W
@@ -458,22 +381,19 @@
       CALL wzv( zu, zv, zw, zhdiv )
 
-# if defined key_zdfddm
-      CALL iom_get( numfl_w, jpdom_data, 'voddmavs', zavt (:,:,:), jkenr )
-#else
-      CALL iom_get( numfl_w, jpdom_data, 'votkeavt', zavt (:,:,:), jkenr )
-#endif 
-
-#if ! defined key_off_degrad && defined key_traldf_c2d
-      CALL iom_get( numfl_w, jpdom_data, 'soleahtw', zahtw (:,: ), jkenr )
-#  if   defined key_trcldf_eiv 
+      IF( iom_varid( numfl_w, 'voddmavs', ldstop = .FALSE. ) > 0 ) THEN          ! avs exist: it is used
+         CALL iom_get( numfl_w, jpdom_data, 'voddmavs', zavt (:,:,:), jkenr )
+      ELSE                                                                       ! no avs: use avt
+         CALL iom_get( numfl_w, jpdom_data, 'votkeavt', zavt (:,:,:), jkenr )
+      ENDIF
+
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv 
       CALL iom_get( numfl_w, jpdom_data, 'soleaeiw', zaeiw (:,: ), jkenr )
-#  endif
-#endif
-
-#if defined key_off_degrad
+#endif
+
+#if defined key_degrad
       CALL iom_get( numfl_u, jpdom_data, 'vozoahtu', zahtu(:,:,:), jkenr )
       CALL iom_get( numfl_v, jpdom_data, 'vomeahtv', zahtv(:,:,:), jkenr )
       CALL iom_get( numfl_w, jpdom_data, 'voveahtw', zahtw(:,:,:), jkenr )
-#  if defined key_trcldf_eiv
+#  if defined key_traldf_eiv
       CALL iom_get( numfl_u, jpdom_data, 'vozoaeiu', zaeiu(:,:,:), jkenr )
       CALL iom_get( numfl_v, jpdom_data, 'vomeaeiv', zaeiv(:,:,:), jkenr )
@@ -486,24 +406,17 @@
       wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:)
 
-#if defined key_trc_diatrd
-      hdivdta(:,:,:,2) = zhdiv(:,:,:) * tmask(:,:,:)
-#endif
-
       tdta(:,:,:,2)   = zt  (:,:,:) * tmask(:,:,:)
       sdta(:,:,:,2)   = zs  (:,:,:) * tmask(:,:,:)
       avtdta(:,:,:,2) = zavt(:,:,:) * tmask(:,:,:)
 
-#if ! defined key_off_degrad && defined key_traldf_c2d
-      ahtwdta(:,:,2)  = zahtw(:,:) * tmask(:,:,1)
-#if defined key_trcldf_eiv
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
       aeiwdta(:,:,2)  = zaeiw(:,:) * tmask(:,:,1)
 #endif
-#endif
-
-#if defined key_off_degrad
+
+#if defined key_degrad
         ahtudta(:,:,:,2) = zahtu(:,:,:) * umask(:,:,:)
         ahtvdta(:,:,:,2) = zahtv(:,:,:) * vmask(:,:,:)
         ahtwdta(:,:,:,2) = zahtw(:,:,:) * tmask(:,:,:)
-#  if defined key_trcldf_eiv
+#  if defined key_traldf_eiv
         aeiudta(:,:,:,2) = zaeiu(:,:,:) * umask(:,:,:)
         aeivdta(:,:,:,2) = zaeiv(:,:,:) * vmask(:,:,:)
@@ -519,13 +432,12 @@
       qsrdta (:,:,2)  = zqsr(:,:) * tmask(:,:,1)
       hmlddta(:,:,2)  = zmld(:,:) * tmask(:,:,1)
-      
-#if defined key_trcbbl_dif   ||   defined key_trcbbl_adv
-      bblxdta(:,:,2) = MAX( 0., zbblx(:,:) )
-      bblydta(:,:,2) = MAX( 0., zbbly(:,:) )
-
-      WHERE( bblxdta(:,:,2) > 2. ) bblxdta(:,:,2) = 0.
-      WHERE( bblydta(:,:,2) > 2. ) bblydta(:,:,2) = 0.
-#endif
-
+
+      IF( l_offbbl ) THEN 
+         bblxdta(:,:,2) = MAX( 0., zbblx(:,:) )
+         bblydta(:,:,2) = MAX( 0., zbbly(:,:) )
+         WHERE( bblxdta(:,:,2) > 2. ) bblxdta(:,:,2) = 0.
+         WHERE( bblydta(:,:,2) > 2. ) bblydta(:,:,2) = 0.
+      ENDIF
+      
       IF( kt == nitend ) THEN
          CALL iom_close ( numfl_t )
@@ -534,7 +446,8 @@
          CALL iom_close ( numfl_w )
       ENDIF
-      
+      !      
    END SUBROUTINE dynrea
 
+
    SUBROUTINE dta_dyn_init
       !!----------------------------------------------------------------------
@@ -544,18 +457,8 @@
       !!
       !! ** Method :
-      !!
-      !! History :
-      !!    ! original  : 92-01 (M. Imbard: sub domain)
-      !!    ! 98-04 (L.Bopp MA Foujols: slopes for isopyc.)
-      !!    ! 98-05 (L. Bopp read output of coupled run)
-      !!    ! 05-03 (O. Aumont and A. El Moussaoui) F90
-      !!----------------------------------------------------------------------
-      !! * Modules used
-
-      !! * Local declarations
-
+      !!----------------------------------------------------------------------
       REAL(wp) ::   znspyr   !: number of time step per year
-
-      NAMELIST/namdyn/ ndtadyn, ndtatot, nsptint, nficdyn, lperdyn,  &
+      !!
+      NAMELIST/namdyn/ ndtadyn, ndtatot, nsptint, lperdyn,  &
       &                cfile_grid_T, cfile_grid_U, cfile_grid_V, cfile_grid_W
       !!----------------------------------------------------------------------
@@ -564,9 +467,8 @@
       ! ======================================
 
-      ! Read Namelist namdyn : Lateral physics on tracers
-      REWIND( numnam )
+      REWIND( numnam )              ! Read Namelist namdyn : Lateral physics on tracers
       READ  ( numnam, namdyn )
 
-      IF(lwp) THEN
+      IF(lwp) THEN                  ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'namdyn : offline dynamical selection'
@@ -577,5 +479,4 @@
          WRITE(numout,*) ' total number of elements in the FILE       ndtatot = ' , ndtatot
          WRITE(numout,*) ' type of interpolation                      nsptint = ' , nsptint
-         WRITE(numout,*) ' number of dynamics FILE                    nficdyn = ' , nficdyn
          WRITE(numout,*) ' loop on the same FILE                      lperdyn = ' , lperdyn
          WRITE(numout,*) '  '
@@ -586,11 +487,14 @@
          WRITE(numout,*) ' '
       ENDIF
-
+      !
       znspyr   = nyear_len(1) * rday / rdt  
       rnspdta  = znspyr / FLOAT( ndtadyn )
       rnspdta2 = rnspdta * 0.5 
-
+      !
+      CALL dta_dyn( nit000 )
+      !
    END SUBROUTINE dta_dyn_init
 
+
    SUBROUTINE wzv( pu, pv, pw, phdiv )
       !!----------------------------------------------------------------------
@@ -599,29 +503,20 @@
       !! ** Purpose :   Compute the now vertical velocity after the array swap
       !!
-      !! ** Method  :
-      !! ** Method  : - Divergence:
-      !!      - compute the now divergence given by :
-      !!         * z-coordinate
+      !! ** 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) and,
-      !!      in regid-lid case, w=0 at the sea surface.
-      !!
-      !!
-      !! History :
-      !!   9.0  !  02-07  (G. Madec)  Vector optimization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in  )   :: pu, pv    !:  horizontal velocities
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out )   :: pw        !:  verticla velocity
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: phdiv     !:  horizontal divergence
-
-      !! * Local declarations
+      !!        The boundary conditions are w=0 at the bottom (no flux).
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pu, pv    !:  horizontal velocities
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: pw        !:  verticla velocity
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv     !:  horizontal divergence
+      !!
       INTEGER  ::  ji, jj, jk
       REAL(wp) ::  zu, zu1, zv, zv1, zet
-
-
+      !!----------------------------------------------------------------------
+      !
       ! Computation of vertical velocity using horizontal divergence
       phdiv(:,:,:) = 0.
@@ -629,11 +524,4 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-#if defined key_zco
-               zu  = pu(ji  ,jj  ,jk) * umask(ji  ,jj  ,jk) * e2u(ji  ,jj  )
-               zu1 = pu(ji-1,jj  ,jk) * umask(ji-1,jj  ,jk) * e2u(ji-1,jj  )
-               zv  = pv(ji  ,jj  ,jk) * vmask(ji  ,jj  ,jk) * e1v(ji  ,jj  )
-               zv1 = pv(ji  ,jj-1,jk) * vmask(ji  ,jj-1,jk) * e1v(ji  ,jj-1)
-               zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) )
-#else
                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)
@@ -641,21 +529,51 @@
                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) )
-#endif
                phdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 
             END DO
          END DO
-      ENDDO
-
-      ! Lateral boundary conditions on phdiv
-      CALL lbc_lnk( phdiv, 'T', 1. )
-
-
+      END DO
+      CALL lbc_lnk( phdiv, 'T', 1. )      ! Lateral boundary conditions on phdiv
+      !
       ! computation of vertical velocity from the bottom
-      pw(:,:,jpk) = 0.
+      pw(:,:,jpk) = 0._wp
       DO jk = jpkm1, 1, -1
          pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * phdiv(:,:,jk)
       END DO
-
+      !
    END SUBROUTINE wzv
+
+
+   SUBROUTINE dta_eiv( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dta_eiv  ***
+      !!
+      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
+      !!      growth rate of baroclinic instability.
+      !!
+      !! ** Method : Specific to the offline model. Computes the horizontal
+      !!             values from the vertical value
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
+      !!
+      INTEGER ::   ji, jj           ! dummy loop indices
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dta_eiv : eddy induced velocity coefficients'
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+      !
+      ! Average the diffusive coefficient at u- v- points
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
+            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
+         END DO
+      END DO
+      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )    ! lateral boundary condition
+      !
+   END SUBROUTINE dta_eiv
+
 
    SUBROUTINE tau2wnd( ptaux, ptauy, pwspd )
@@ -667,13 +585,11 @@
       !! ** Method  : |tau|=rhoa*Cd*|U|^2
       !!---------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( in  ) ::  &
-         ptaux, ptauy                              !: wind stress in i-j direction resp.
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::  &
-         pwspd                                     !: wind speed 
-      REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3
-      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient
-      REAL(wp) ::   ztx, zty, ztau, zcoef ! temporary variables
-      INTEGER  ::   ji, jj                ! dummy indices
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptaux, ptauy   ! wind stress in i-j direction resp.
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pwspd          ! wind speed
+      !! 
+      REAL(wp) ::   zrhoa  = 1.22_wp       ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3_wp     ! drag coefficient
+      REAL(wp) ::   ztx, zty, ztau, zcoef  ! temporary variables
+      INTEGER  ::   ji, jj                 ! dummy indices
       !!---------------------------------------------------------------------
       zcoef = 1. / ( zrhoa * zcdrag )
@@ -689,185 +605,7 @@
       END DO
       CALL lbc_lnk( pwspd(:,:), 'T', 1. )
-
+      !
    END SUBROUTINE tau2wnd
 
-#if defined key_trcbbl_dif   ||   defined key_trcbbl_adv
-
-   SUBROUTINE bbl_sign( ptn, psn, pbblx, pbbly )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE bbl_sign  ***
-      !!
-      !! ** Purpose :   Compute the sign of local gradient of density multiplied by the slope
-      !!                along the bottom slope gradient : grad( rho) * grad(h)
-      !!                Need to compute the diffusive bottom boundary layer
-      !!
-      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad
-      !!      is an along bottom slope gradient) an additional lateral diffu-
-      !!      sive trend along the bottom slope is added to the general tracer
-      !!      trend, otherwise nothing is done. See trcbbl.F90
-      !!
-      !!
-      !! History :
-      !!   9.0  !  02-07  (G. Madec)  Vector optimization
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in  ) ::  &
-         ptn             ,  &                           !: temperature 
-         psn                                            !: salinity 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::  &
-         pbblx , pbbly                                  !: sign of bbl in i-j direction resp. 
-      
-      !! * Local declarations
-      INTEGER  ::   ji, jj                   ! dummy loop indices
-      INTEGER  ::   ik
-      REAL(wp) ::   &
-         ztx, zsx, zhx, zalbetx, zgdrhox,     &  ! temporary scalars
-         zty, zsy, zhy, zalbety, zgdrhoy 
-      REAL(wp), DIMENSION(jpi,jpj) ::    &
-        ztnb, zsnb, zdep
-      REAL(wp) ::    fsalbt, pft, pfs, pfh   ! statement function
-      !!----------------------------------------------------------------------
-      ! ratio alpha/beta
-      ! ================
-      !  fsalbt: ratio of thermal over saline expension coefficients
-      !       pft :  potential temperature in degrees celcius
-      !       pfs :  salinity anomaly (s-35) in psu
-      !       pfh :  depth in meters
-
-      fsalbt( pft, pfs, pfh ) =                                              &
-         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
-                                   - 0.203814e-03 ) * pft                    &
-                                   + 0.170907e-01 ) * pft                    &
-                                   + 0.665157e-01                            &
-         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
-         +  ( ( - 0.302285e-13 * pfh                                         &
-                - 0.251520e-11 * pfs                                         &
-                + 0.512857e-12 * pft * pft          ) * pfh                  &
-                                     - 0.164759e-06   * pfs                  &
-             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
-                                     + 0.380374e-04 ) * pfh
-
-      ! 0. 2D fields of bottom temperature and salinity, and bottom slope
-      ! -----------------------------------------------------------------
-      ! mbathy= number of w-level, minimum value=1 (cf domrea.F90)
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-#  endif
-            ik          =  MAX( mbathy(ji,jj) - 1, 1 )    ! vertical index of the bottom ocean T-level
-            ztnb(ji,jj) = ptn(ji,jj,ik) * tmask(ji,jj,1)  ! masked T and S at ocean bottom
-            zsnb(ji,jj) = psn(ji,jj,ik) * tmask(ji,jj,1)
-            zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
-#  if ! defined key_vectopt_loop
-         END DO
-#  endif
-      END DO
-
-      !!----------------------------------------------------------------------
-      ! 1. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
-      ! --------------------------------------------
-      ! Sign of the local density gradient along the i- and j-slopes
-      ! multiplied by the slope of the ocean bottom
-
-      SELECT CASE ( neos )
-
-      CASE ( 0 )                 ! Jackett and McDougall (1994) formulation
-
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-#  endif
-            ! temperature, salinity anomalie and depth
-            ztx = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
-            zsx = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
-            zhx = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
-            !
-            zty = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
-            zsy = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
-            zhy = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
-            ! masked ratio alpha/beta
-            zalbetx = fsalbt( ztx, zsx, zhx ) * umask(ji,jj,1)
-            zalbety = fsalbt( zty, zsy, zhy ) * vmask(ji,jj,1)
-            ! local density gradient along i-bathymetric slope
-            zgdrhox = zalbetx * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
-                   -            ( zsnb(ji+1,jj) - zsnb(ji,jj) )
-            ! local density gradient along j-bathymetric slope
-            zgdrhoy = zalbety * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
-                   -            ( zsnb(ji,jj+1) - zsnb(ji,jj) )
-            ! sign of local i-gradient of density multiplied by the i-slope
-            pbblx(ji,jj) = 0.5 - SIGN( 0.5, -zgdrhox * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
-            ! sign of local j-gradient of density multiplied by the j-slope
-            pbbly(ji,jj) = 0.5 - SIGN( 0.5, -zgdrhoy * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
-#  if ! defined key_vectopt_loop
-         END DO
-#  endif
-      END DO
-
-      CASE ( 1 )               ! Linear formulation function of temperature only
-                               !
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-#  endif
-            ! local 'density/temperature' gradient along i-bathymetric slope
-            zgdrhox =  ztnb(ji+1,jj) - ztnb(ji,jj)
-            ! local density gradient along j-bathymetric slope
-            zgdrhoy =  ztnb(ji,jj+1) - ztnb(ji,jj)
-            ! sign of local i-gradient of density multiplied by the i-slope
-            pbblx(ji,jj) = 0.5 - SIGN( 0.5, -zgdrhox * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
-            ! sign of local j-gradient of density multiplied by the j-slope
-            pbbly(ji,jj) = 0.5 - SIGN( 0.5, -zgdrhoy * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
-#  if ! defined key_vectopt_loop
-         END DO
-#  endif
-      END DO
-
-      CASE ( 2 )               ! Linear formulation function of temperature and salinity
-
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-#  endif     
-            ! local density gradient along i-bathymetric slope
-            zgdrhox = - ( rbeta*( zsnb(ji+1,jj) - zsnb(ji,jj) )   &
-                      -  ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
-            ! local density gradient along j-bathymetric slope
-            zgdrhoy = - ( rbeta*( zsnb(ji,jj+1) - zsnb(ji,jj) )   &
-                      -  ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )
-            ! sign of local i-gradient of density multiplied by the i-slope
-            pbblx(ji,jj) = 0.5 - SIGN( 0.5, - zgdrhox * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
-            ! sign of local j-gradient of density multiplied by the j-slope
-            pbbly(ji,jj) = 0.5 - SIGN( 0.5, -zgdrhoy * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
-#  if ! defined key_vectopt_loop
-         END DO
-#  endif
-      END DO
-
-      CASE DEFAULT
-
-         WRITE(ctmp1,*) '          bad flag value for neos = ', neos
-         CALL ctl_stop(ctmp1)
-
-      END SELECT
-   
-      ! Lateral boundary conditions
-      CALL lbc_lnk( pbblx, 'U', 1. )
-      CALL lbc_lnk( pbbly, 'V', 1. )
-
-   END SUBROUTINE bbl_sign
-
-#endif
 
    SUBROUTINE swap_dyn_data
@@ -876,10 +614,6 @@
       !!
       !! ** Purpose :   swap array data
-      !!
-      !! History :
-      !!   9.0  !  07-09  (C. Ethe)
-      !!----------------------------------------------------------------------
-
-
+      !!----------------------------------------------------------------------
+      !
       ! swap from record 2 to 1
       tdta   (:,:,:,1) = tdta   (:,:,:,2)
@@ -889,9 +623,5 @@
       vdta   (:,:,:,1) = vdta   (:,:,:,2)
       wdta   (:,:,:,1) = wdta   (:,:,:,2)
-#if defined key_trc_diatrd
-      hdivdta(:,:,:,1) = hdivdta(:,:,:,2)
-#endif
-
-#if defined key_ldfslp
+#if defined key_ldfslp && ! defined key_c1d
       uslpdta (:,:,:,1) = uslpdta (:,:,:,2)
       vslpdta (:,:,:,1) = vslpdta (:,:,:,2)
@@ -904,17 +634,18 @@
       empdta (:,:,1) = empdta (:,:,2) 
       qsrdta (:,:,1) = qsrdta (:,:,2) 
-
-#if ! defined key_off_degrad && defined key_traldf_c2d
-      ahtwdta(:,:,1) = ahtwdta(:,:,2)
-#  if defined key_trcldf_eiv
+      IF( l_offbbl ) THEN
+         bblxdta(:,:,1) = bblxdta(:,:,2)
+         bblydta(:,:,1) = bblydta(:,:,2) 
+      ENDIF
+
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
       aeiwdta(:,:,1) = aeiwdta(:,:,2)
-#  endif
-#endif
-
-#if defined key_off_degrad
+#endif
+
+#if defined key_degrad
       ahtudta(:,:,:,1) = ahtudta(:,:,:,2)
       ahtvdta(:,:,:,1) = ahtvdta(:,:,:,2)
       ahtwdta(:,:,:,1) = ahtwdta(:,:,:,2)
-#  if defined key_trcldf_eiv
+#  if defined key_traldf_eiv
       aeiudta(:,:,:,1) = aeiudta(:,:,:,2)
       aeivdta(:,:,:,1) = aeivdta(:,:,:,2)
@@ -922,11 +653,7 @@
 #  endif
 #endif
-
-#if defined key_trcbbl_dif || defined key_trcbbl_adv
-      bblxdta(:,:,1) = bblxdta(:,:,2)
-      bblydta(:,:,1) = bblydta(:,:,2)
-#endif
-
+      !
    END SUBROUTINE swap_dyn_data
+
 
    SUBROUTINE assign_dyn_data
@@ -939,22 +666,13 @@
       !!----------------------------------------------------------------------
       
-      tn (:,:,:) = tdta  (:,:,:,2)
-      sn (:,:,:) = sdta  (:,:,:,2)
-      avt(:,:,:) = avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = sdta  (:,:,:,2)
+      avt(:,:,:)        = avtdta(:,:,:,2)
       
       un (:,:,:) = udta  (:,:,:,2) 
       vn (:,:,:) = vdta  (:,:,:,2)
       wn (:,:,:) = wdta  (:,:,:,2)
-
-#if defined key_trc_diatrd
-      hdivn(:,:,:) = hdivdta(:,:,:,2)
-#endif
-
-#if defined key_zdfddm
-      avs(:,:,:)   = avtdta (:,:,:,2)
-#endif
-
-      
-#if defined key_ldfslp
+      
+#if defined key_ldfslp && ! defined key_c1d
       uslp (:,:,:) = uslpdta (:,:,:,2) 
       vslp (:,:,:) = vslpdta (:,:,:,2) 
@@ -969,42 +687,34 @@
       emps(:,:) = emp(:,:) 
       qsr (:,:) = qsrdta (:,:,2) 
-
-#if ! defined key_off_degrad && defined key_traldf_c2d    
-      ahtw(:,:) = ahtwdta(:,:,2)
-#  if defined key_trcldf_eiv
+      IF( l_offbbl ) THEN
+         ahu_bbl(:,:) = bblxdta(:,:,2)
+         ahv_bbl(:,:) = bblydta(:,:,2) 
+      ENDIF
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
       aeiw(:,:) = aeiwdta(:,:,2)
-#  endif
-#endif
-      
-#if defined key_off_degrad
+#endif
+      
+#if defined key_degrad
       ahtu(:,:,:) = ahtudta(:,:,:,2)
       ahtv(:,:,:) = ahtvdta(:,:,:,2)
       ahtw(:,:,:) = ahtwdta(:,:,:,2)
-#  if defined key_trcldf_eiv
+#  if defined key_traldf_eiv
       aeiu(:,:,:) = aeiudta(:,:,:,2)
       aeiv(:,:,:) = aeivdta(:,:,:,2)
       aeiw(:,:,:) = aeiwdta(:,:,:,2)
 #  endif
-      
-#endif
-      
-#if defined key_trcbbl_dif ||  defined key_trcbbl_adv
-      bblx(:,:) = bblxdta(:,:,2)
-      bbly(:,:) = bblydta(:,:,2)
-#endif
-
+#endif
+      !
    END SUBROUTINE assign_dyn_data
 
+
    SUBROUTINE linear_interp_dyn_data( pweigh )
       !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE linear_interp_dyn_data  ***
+      !!               ***  ROUTINE linear_interp_dyn_data  ***
       !!
       !! ** Purpose :   linear interpolation of data
-      !!
-      !!----------------------------------------------------------------------
-      !! * Argument
-      REAL(wp), INTENT( in ) ::   pweigh       ! weigh
-
-      !! * Local declarations
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   pweigh   ! weigh
+      !!
       REAL(wp) :: zweighm1
       !!----------------------------------------------------------------------
@@ -1012,22 +722,13 @@
       zweighm1 = 1. - pweigh
       
-      tn (:,:,:) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
-      sn (:,:,:) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
-      avt(:,:,:) = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
+      tsn(:,:,:,jp_tem) = zweighm1 * tdta  (:,:,:,1) + pweigh * tdta  (:,:,:,2)
+      tsn(:,:,:,jp_sal) = zweighm1 * sdta  (:,:,:,1) + pweigh * sdta  (:,:,:,2)
+      avt(:,:,:)        = zweighm1 * avtdta(:,:,:,1) + pweigh * avtdta(:,:,:,2)
       
       un (:,:,:) = zweighm1 * udta  (:,:,:,1) + pweigh * udta  (:,:,:,2) 
       vn (:,:,:) = zweighm1 * vdta  (:,:,:,1) + pweigh * vdta  (:,:,:,2)
       wn (:,:,:) = zweighm1 * wdta  (:,:,:,1) + pweigh * wdta  (:,:,:,2)
-
-#if defined key_trc_diatrd
-      hdivn(:,:,:) = zweighm1 * hdivdta(:,:,:,1) + pweigh * hdivdta(:,:,:,2)
-#endif
-
-#if defined key_zdfddm
-      avs(:,:,:)   = zweighm1 * avtdta (:,:,:,1) + pweigh * avtdta (:,:,:,2)
-#endif
-
-      
-#if defined key_ldfslp
+      
+#if defined key_ldfslp && ! defined key_c1d
       uslp (:,:,:) = zweighm1 * uslpdta (:,:,:,1) + pweigh * uslpdta (:,:,:,2) 
       vslp (:,:,:) = zweighm1 * vslpdta (:,:,:,1) + pweigh * vslpdta (:,:,:,2) 
@@ -1042,17 +743,18 @@
       emps(:,:) = emp(:,:) 
       qsr (:,:) = zweighm1 * qsrdta (:,:,1) + pweigh  * qsrdta (:,:,2) 
-
-#if ! defined key_off_degrad && defined key_traldf_c2d    
-      ahtw(:,:) = zweighm1 * ahtwdta(:,:,1) + pweigh * ahtwdta(:,:,2)
-#  if defined key_trcldf_eiv
+      IF( l_offbbl ) THEN
+         ahu_bbl(:,:) = zweighm1 * bblxdta(:,:,1) +  pweigh  * bblxdta(:,:,2)
+         ahv_bbl(:,:) = zweighm1 * bblydta(:,:,1) +  pweigh  * bblydta(:,:,2)
+      ENDIF
+
+#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv 
       aeiw(:,:) = zweighm1 * aeiwdta(:,:,1) + pweigh * aeiwdta(:,:,2)
-#  endif
-#endif
-      
-#if defined key_off_degrad
+#endif
+      
+#if defined key_degrad
       ahtu(:,:,:) = zweighm1 * ahtudta(:,:,:,1) + pweigh * ahtudta(:,:,:,2)
       ahtv(:,:,:) = zweighm1 * ahtvdta(:,:,:,1) + pweigh * ahtvdta(:,:,:,2)
       ahtw(:,:,:) = zweighm1 * ahtwdta(:,:,:,1) + pweigh * ahtwdta(:,:,:,2)
-#  if defined key_trcldf_eiv
+#  if defined key_traldf_eiv
       aeiu(:,:,:) = zweighm1 * aeiudta(:,:,:,1) + pweigh * aeiudta(:,:,:,2)
       aeiv(:,:,:) = zweighm1 * aeivdta(:,:,:,1) + pweigh * aeivdta(:,:,:,2)
@@ -1060,11 +762,7 @@
 #  endif
 #endif
-      
-#if defined key_trcbbl_dif   ||   defined key_trcbbl_adv
-      bblx(:,:) = zweighm1 * bblxdta(:,:,1) + pweigh * bblxdta(:,:,2)
-      bbly(:,:) = zweighm1 * bblydta(:,:,1) + pweigh * bblydta(:,:,2)
-#endif
-
+      !      
    END SUBROUTINE linear_interp_dyn_data
 
+   !!======================================================================
 END MODULE dtadyn
