Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 2392)
@@ -5,4 +5,10 @@
    !!                          assimilation
    !!======================================================================
+   !! History :       ! 2007-03  (M. Martin)  Met Office version
+   !!                 ! 2007-04  (A. Weaver)  calc_date original code
+   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
+   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
+   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init 
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -15,99 +21,31 @@
    !!   ssh_asm_inc  : Apply the SSH increment
    !!----------------------------------------------------------------------
-   !! * Modules used   
-   USE par_kind, ONLY : &       ! Precision variables
-      & wp
-   USE in_out_manager, ONLY : & ! I/O manager 
-      & lwp,      &
-      & numnam,   &
-      & numout,   &
-      & ctl_warn, &
-      & ctl_stop, &
-      & nit000,   &
-      & nstop,    &
-      & ln_rstart
-   USE par_oce, ONLY : & ! Ocean space and time domain variables
-      & jpi,  & 
-      & jpj,  &
-      & jpk,  &
-      & jpkm1
-   USE dom_oce, ONLY : & ! Ocean space and time domain
-      & rdt,    &
-      & n_cla,  &
-      & neuler, &
-      & ln_zps, &
-      & tmask,  &
-      & umask,  &
-      & vmask
+   USE in_out_manager   ! I/O manager
+   USE par_oce          ! Ocean space and time domain variables
+   USE dom_oce          ! Ocean space and time domain
+   USE oce              ! Dynamics and active tracers defined in memory
+   USE divcur           ! Horizontal divergence and relative vorticity
+   USE eosbn2           ! Equation of state - in situ and potential density
+   USE zpshde           ! Partial step : Horizontal Derivative
+   USE iom              ! Library to read input files
+   USE asmpar           ! Parameters for the assmilation interface
 #if defined key_c1d
-   USE c1d, ONLY : &    ! 1D initialization
-      & lk_c1d
+   USE c1d, ONLY :   lk_c1d    ! 1D initialization
 #endif
-   USE oce, ONLY : &      ! Dynamics and active tracers defined in memory
-      & ub, un, ua,    &
-      & vb, vn, va,    &
-      & tsb, tsn, tsa, &
-      & sshb, sshn,    &
-      & rhd, rhop,     &
-      & rotb, rotn,    &
-      & hdivb, hdivn,  &
-      & gtsu, gru,     &
-      & gtsv, grv 
-   USE divcur, ONLY : &   ! Horizontal divergence and relative vorticity
-      & div_cur
-   USE cla_div, ONLY : &  ! Specific update of the horizontal divergence
-      & div_cla           ! (specific to ORCA_R2)
-   USE eosbn2, ONLY : &   ! Equation of state - in situ and potential density
-      & eos
-   USE zpshde, ONLY : &   ! Partial step : Horizontal Derivative
-      & zps_hde
-!   USE phycst, ONLY : &   ! Calendar parameters
-!      & rjjss  
-   USE iom                ! Library to read input files
-   USE asmpar             ! Parameters for the assmilation interface
-   USE dom_oce, ONLY : &
-      & ndastp
-!   USE daymod, ONLY : &
-!      & nmonth_len         ! length of month in days
-!      & nbiss, &
-!      & nobis, &
- 
 
    IMPLICIT NONE
-
-   !! * Routine accessibility
    PRIVATE
-   PUBLIC asm_inc_init, & !: Initialize the increment arrays and IAU weights
-      &   calc_date,    & !: Compute the calendar date YYYYMMDD on a given step
-      &   tra_asm_inc,  & !: Apply the tracer (T and S) increments
-      &   dyn_asm_inc,  & !: Apply the dynamic (u and v) increments
-      &   ssh_asm_inc     !: Apply the SSH increment
-
-   !! * Private Module variables
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: &
-      & t_bkg,    &       !: Background temperature 
-      & s_bkg,    &       !: Background salinity
-      & u_bkg,    &       !: Background u-component velocity
-      & v_bkg,    &       !: Background v-component velocity
-      & t_bkginc, &       !: Increment to the background temperature
-      & s_bkginc, &       !: Increment to the background salinity
-      & u_bkginc, &       !: Increment to the u-component velocity 
-      & v_bkginc          !: Increment to the v-component velocity 
-         
-   REAL(wp), PRIVATE, DIMENSION(:,:), ALLOCATABLE :: &
-      & ssh_bkg,   &      !: Background sea surface height
-      & ssh_bkginc        !: Increment to the background sea surface height
-
-   REAL(wp), PUBLIC, DIMENSION(:), ALLOCATABLE :: &
-      & wgtiau            !: IAU weights for each time step
-
-   !! * Shared Module variables
-   LOGICAL, PUBLIC, PARAMETER :: &
+   
+   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
+   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
+   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
+   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
+   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
+
 #if defined key_asminc
-      & lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
+    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
 #else
-      & lk_asminc = .FALSE.  !: No assimilation increments
+    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
 #endif
-
    LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields
    LOGICAL, PUBLIC :: ln_trjwri = .FALSE. !: No output of the state trajectory fields
@@ -119,29 +57,29 @@
    LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check
 
-   REAL, PUBLIC :: salfixmin   !: Ensure that the salinity is larger than 
-   !                           !: this  value if (ln_salfix)
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components 
+   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
 #if defined key_asminc
-   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: &
-      & ssh_iau           !: IAU-weighted sea surface height increment
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment
 #endif
-
-   INTEGER, PUBLIC :: nitbkg     !: Time step of the background state used in the Jb term
-                                 !: (relative to the cycle interval [0,nitend-nit000-1])
-   INTEGER, PUBLIC :: nitdin     !: Time step of the background state for direct initialization
-                                 !: (relative to the cycle interval [0,nitend-nit000-1])
-   INTEGER, PUBLIC :: nitiaustr  !: Time step of the start of the IAU interval 
-                                 !: (relative to the cycle interval [0,nitend-nit000-1])
-   INTEGER, PUBLIC :: nitiaufin  !: Time step of the end of the IAU interval 
-                                 !: (relative to the cycle interval [0,nitend-nit000-1])
-   INTEGER, PUBLIC :: niaufn     !: Type of IAU weighing function
-                                 !:   0 = Constant weighting
-                                 !:   1 = Linear hat-like, centred in middle of IAU interval 
+   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
+   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
+   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
+   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval 
+   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
+   ! 
+   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
+   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval 
+   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
+
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
 
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -155,22 +93,6 @@
       !!
       !! ** Action  : 
-      !!
-      !! History :
-      !!        !  07-03  (M. Martin) Met Office version
-      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
-      !!----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-      !! * Modules used
-      NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri, &
-         &                 ln_trainc, ln_dyninc, ln_sshinc, &
-         &                 ln_asmdin, ln_asmiau, &
-         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, &
-         &                 nittrjfrq, ln_salfix, salfixmin
-
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER :: jt
-
       INTEGER :: imid
       INTEGER :: inum
@@ -190,4 +112,11 @@
       REAL(wp) :: zdate_bkg    ! Date in background state file for DI
       REAL(wp) :: zdate_inc    ! Time axis in increments file
+      !!
+      NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri,                           &
+         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
+         &                 ln_asmdin, ln_asmiau,                           &
+         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
+         &                 nittrjfrq, ln_salfix, salfixmin
+      !!----------------------------------------------------------------------
 
       !-----------------------------------------------------------------------
@@ -220,35 +149,20 @@
          WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '          Namelist namasm : set assimilation increment parameters'
-         WRITE(numout,*) '             Logical switch for writing out background state         ', &
-            &            ' ln_bkgwri = ', ln_bkgwri
-         WRITE(numout,*) '             Logical switch for writing out state trajectory         ', &
-            &            ' ln_trjwri = ', ln_trjwri
-         WRITE(numout,*) '             Logical switch for applying tracer increments           ', &
-            &            ' ln_trainc = ', ln_trainc
-         WRITE(numout,*) '             Logical switch for applying velocity increments         ', &
-            &            ' ln_dyninc = ', ln_dyninc
-         WRITE(numout,*) '             Logical switch for applying SSH increments              ', &
-            &            ' ln_sshinc = ', ln_sshinc
-         WRITE(numout,*) '             Logical switch for Direct Initialization (DI)           ', &
-            &            ' ln_asmdin = ', ln_asmdin
-         WRITE(numout,*) '             Logical switch for Incremental Analysis Updating (IAU)  ', &
-            &            ' ln_asmiau = ', ln_asmiau
-         WRITE(numout,*) '             Timestep of background in [0,nitend-nit000-1]           ', &
-            &            '    nitbkg = ', nitbkg
-         WRITE(numout,*) '             Timestep of background for DI in [0,nitend-nit000-1]    ', &
-            &            '    nitdin = ', nitdin
-         WRITE(numout,*) '             Timestep of start of IAU interval in [0,nitend-nit000-1]', &
-            &            ' nitiaustr = ', nitiaustr
-         WRITE(numout,*) '             Timestep of end of IAU interval in [0,nitend-nit000-1]  ', &
-            &            ' nitiaufin = ', nitiaufin
-         WRITE(numout,*) '             Type of IAU weighting function                          ', &
-            &            '    niaufn = ', niaufn
-         WRITE(numout,*) '             Frequency of trajectory output for 4D-VAR               ', &
-            &            ' nittrjfrq = ', nittrjfrq
-         WRITE(numout,*) '             Logical switch for ensuring that the sa > salfixmin     ', &
-            &            ' ln_salfix = ', ln_salfix
-         WRITE(numout,*) '             Minimum salinity after applying the increments          ', &
-            &            ' salfixmin = ', salfixmin
+         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
+         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
+         WRITE(numout,*) '      Logical switch for writing out state trajectory          ln_trjwri = ', ln_trjwri
+         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
+         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
+         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
+         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
+         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
+         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
+         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
+         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
+         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
+         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
+         WRITE(numout,*) '      Frequency of trajectory output for 4D-VAR                nittrjfrq = ', nittrjfrq
+         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
+         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
       ENDIF
 
@@ -565,6 +479,7 @@
 
       ENDIF
-
+      !
    END SUBROUTINE asm_inc_init
+
 
    SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
@@ -577,21 +492,10 @@
       !!
       !! ** Action  : 
-      !!
-      !! History :
-      !!        !  07-04  (A. Weaver) 
-      !!        !  10-05  (D. Lea)        Update to work with NEMO vn3.2
-      !!----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-      !! * Arguments
-
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kit000  ! Initial time step
       INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
       INTEGER, INTENT(IN) :: kdate0  ! Initial date
       INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
-
-      !! * Local declarations
-
+      !
       INTEGER :: iyea0    ! Initial year
       INTEGER :: imon0    ! Initial month
@@ -648,8 +552,9 @@
          idaycnt = idaycnt + 1
       END DO
-
+      !
       kdate = iyea * 10000 + imon * 100 + iday
-
+      !
    END SUBROUTINE
+
 
    SUBROUTINE calc_month_len( iyear, imonth_len )
@@ -660,14 +565,9 @@
       !!
       !! ** Method  : 
-      !!
-      !! ** Action  : 
-      !!
-      !! History :
-      !!        !  10-05  (D. Lea)   New routine based on day_init 
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
       INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
       INTEGER :: iyear         !: year
-
+      !!----------------------------------------------------------------------
+      !
       ! length of the month of the current year (from nleapy, read in namelist)
       IF ( nleapy < 2 ) THEN 
@@ -681,6 +581,7 @@
          imonth_len(:) = nleapy   ! all months with nleapy days per year
       ENDIF
-
+      !
    END SUBROUTINE
+
 
    SUBROUTINE tra_asm_inc( kt )
@@ -693,19 +594,11 @@
       !!
       !! ** Action  : 
-      !!
-      !! History :
-      !!        !  07-03  (M. Martin) Met Office version
-      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
-      !!----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt               ! Current time step
-
-      !! * Local declarations
+      !
       INTEGER :: ji,jj,jk
       INTEGER :: it
       REAL(wp) :: zincwgt  ! IAU weight for current time step
+      !!----------------------------------------------------------------------
 
       IF ( ln_asmiau ) THEN
@@ -793,8 +686,9 @@
             DEALLOCATE( s_bkg    )
          ENDIF
-         
-      ENDIF
-
+         !  
+      ENDIF
+      !
    END SUBROUTINE tra_asm_inc
+
 
    SUBROUTINE dyn_asm_inc( kt )
@@ -807,19 +701,11 @@
       !!
       !! ** Action  : 
-      !!
-      !! History :
-      !!        !  07-03  (M. Martin) Met Office version
-      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
-      !!----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt   ! Current time step
-
-      !! * Local declarations
+      !
       INTEGER :: jk
       INTEGER :: it
       REAL(wp) :: zincwgt  ! IAU weight for current time step
+      !!----------------------------------------------------------------------
 
       IF ( ln_asmiau ) THEN
@@ -871,6 +757,5 @@
             vb(:,:,:) = vn(:,:,:)
  
-            CALL div_cur( kt )                  ! Compute divergence and curl for now fields
-            IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence)
+            CALL div_cur( kt )            ! Compute divergence and curl for now fields
 
             rotb (:,:,:) = rotn (:,:,:)   ! Update before fields
@@ -881,10 +766,10 @@
             DEALLOCATE( u_bkginc )
             DEALLOCATE( v_bkginc )
-
-         ENDIF
-
-      ENDIF
-
+         ENDIF
+         !
+      ENDIF
+      !
    END SUBROUTINE dyn_asm_inc
+
 
    SUBROUTINE ssh_asm_inc( kt )
@@ -897,18 +782,10 @@
       !!
       !! ** Action  : 
-      !!
-      !! History :
-      !!        !  07-03  (M. Martin) Met Office version
-      !!        !  07-04  (A. Weaver) Merge with OPAVAR/NEMOVAR
-      !!----------------------------------------------------------------------
-
-      IMPLICIT NONE
-
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT(IN) :: kt   ! Current time step
-
-      !! * Local declarations
+      !
       INTEGER :: it
       REAL(wp) :: zincwgt  ! IAU weight for current time step
+      !!----------------------------------------------------------------------
 
       IF ( ln_asmiau ) THEN
@@ -960,8 +837,9 @@
 
          ENDIF
-
-      ENDIF
-
+         !
+      ENDIF
+      !
    END SUBROUTINE ssh_asm_inc
 
+   !!======================================================================
 END MODULE asminc
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2392)
@@ -55,7 +55,4 @@
    !                                         !!* Namelist namcla : cross land advection
    INTEGER, PUBLIC ::   nn_cla = 0            !: =1 cross land advection for exchanges through some straits (ORCA2)
-
-   !                                          ! old non-DOCTOR names still used in the model
-   INTEGER, PUBLIC ::   n_cla = 0             !: =1 cross land advection for exchanges through some straits (ORCA2)
 
    !!----------------------------------------------------------------------
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 2392)
@@ -243,9 +243,4 @@
       ENDIF
 
-      n_cla = nn_cla                ! conversion DOCTOR names into model names (this should disappear soon)
-
-      IF( lk_mpp_rep .AND. n_cla /= 0 )   CALL ctl_stop( ' Reproductibility tests (nbit_cmp=1) require n_cla = 0' )
-      !
-
 #if defined key_netcdf4
       !                             ! NetCDF 4 case   ("key_netcdf4" defined)
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90	(revision 2392)
@@ -4,43 +4,37 @@
    !! Ocean initialization : domain initialization
    !!==============================================================================
-   !! History :       !  88-03  (G. Madec)
-   !!                 !  91-11  (G. Madec)
-   !!                 !  92-06  (M. Imbard)
-   !!                 !  96-01  (G. Madec)  terrain following coordinates
-   !!                 !  97-02  (G. Madec)  print mesh informations
-   !!                 !  99-11  (M. Imbard) NetCDF format with IO-IPSL
-   !!                 !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
-   !!                 !  01-09  (M. Levy)  eel config: grid in km, beta-plane
-   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module, namelist
-   !!            9.0  !  04-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
-   !!                           use of parameters in par_CONFIG-Rxx.h90, not in namelist
-   !!                 !  04-05  (A. Koch-Larrouy) Add Gyre configuration 
+   !! History :  OPA  ! 1988-03  (G. Madec) Original code
+   !!            7.0  ! 1996-01  (G. Madec)  terrain following coordinates
+   !!            8.0  ! 1997-02  (G. Madec)  print mesh informations
+   !!            8.1  ! 1999-11  (M. Imbard) NetCDF format with IO-IPSL
+   !!            8.2  ! 2000-08  (D. Ludicone) Reduced section at Bab el Mandeb
+   !!             -   ! 2001-09  (M. Levy)  eel config: grid in km, beta-plane
+   !!  NEMO      1.0  ! 2002-08  (G. Madec)  F90: Free form and module, namelist
+   !!             -   ! 2004-01  (A.M. Treguier, J.M. Molines) Case 4 (Mercator mesh)
+   !!                            use of parameters in par_CONFIG-Rxx.h90, not in namelist
+   !!             -   ! 2004-05  (A. Koch-Larrouy) Add Gyre configuration 
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   dom_hgr        : initialize the horizontal mesh 
-   !!   hgr_read       : read "coordinate" NetCDF file 
+   !!   dom_hgr       : initialize the horizontal mesh 
+   !!   hgr_read      : read "coordinate" NetCDF file 
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Module variables
-   REAL(wp) ::   glam0, gphi0           ! variables corresponding to parameters
-      !                                 ! ppglam0 ppgphi0 set in par_oce
-
-   !! * Routine accessibility
-   PUBLIC dom_hgr        ! called by domain.F90
+   REAL(wp) ::   glam0, gphi0   ! variables corresponding to parameters ppglam0 ppgphi0 set in par_oce
+
+   PUBLIC   dom_hgr   ! called by domain.F90
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -100,15 +94,13 @@
       !!                Madec, Imbard, 1996, Clim. Dyn.
       !!----------------------------------------------------------------------
-      INTEGER  ::   ji, jj              ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
-      INTEGER  ::   ijeq                ! index of equator T point (used in case 4)
-      REAL(wp) ::   &
-         zti, zui, zvi, zfi,         &  ! temporary scalars
-         ztj, zuj, zvj, zfj,         &  !
-         zphi0, zbeta, znorme,       &  !
-         zarg, zf0, zminff, zmaxff
-      REAL(wp) ::   &
-         zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg,   &
-         zphi1, zsin_alpha, zim05, zjm05
+      INTEGER  ::   ji, jj               ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
+      INTEGER  ::   ijeq                 ! index of equator T point (used in case 4)
+      REAL(wp) ::   zti, zui, zvi, zfi   ! local scalars
+      REAL(wp) ::   ztj, zuj, zvj, zfj   !   -      -
+      REAL(wp) ::   zphi0, zbeta, znorme !
+      REAL(wp) ::   zarg, zf0, zminff, zmaxff
+      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg
+      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05
       !!----------------------------------------------------------------------
 
@@ -138,5 +130,5 @@
          IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
             !                                             ! =====================
-            IF( n_cla == 0 ) THEN
+            IF( nn_cla == 0 ) THEN
                !
                ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u = 20 km)
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 2392)
@@ -5,12 +5,12 @@
    !!======================================================================
    !! History :  OPA  ! 1987-07  (G. Madec)  Original code
-   !!             -   ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
-   !!             -   ! 1996-01  (G. Madec)  suppression of common work arrays
+   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
+   !!            7.0  ! 1996-01  (G. Madec)  suppression of common work arrays
    !!             -   ! 1996-05  (G. Madec)  mask computed from tmask and sup-
    !!                 !                      pression of the double computation of bmask
-   !!             -   ! 1997-02  (G. Madec)  mesh information put in domhgr.F
-   !!             -   ! 1997-07  (G. Madec)  modification of mbathy and fmask
+   !!            8.0  ! 1997-02  (G. Madec)  mesh information put in domhgr.F
+   !!            8.1  ! 1997-07  (G. Madec)  modification of mbathy and fmask
    !!             -   ! 1998-05  (G. Roullet)  free surface
-   !!             -   ! 2000-03  (G. Madec)  no slip accurate
+   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
    !!             -   ! 2001-09  (J.-M. Molines)  Open boundaries
    !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
@@ -44,7 +44,6 @@
    !! NEMO/OPA 3.2 , LODYC-IPSL  (2009)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
    
@@ -132,5 +131,5 @@
       ENDIF
 
-      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral free-slip '
+      IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  free-slip '
       ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  no-slip '
       ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ocean lateral  partial-slip '
@@ -308,5 +307,5 @@
       IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA_R2 configuration
          !                                                 ! Increased lateral friction near of some straits
-         IF( n_cla == 0 ) THEN
+         IF( nn_cla == 0 ) THEN
             !                                ! Gibraltar strait  : partial slip (fmask=0.5)
             ij0 = 101   ;   ij1 = 101
@@ -322,5 +321,4 @@
             !
          ENDIF
-
          !                                ! Danish straits  : strong slip (fmask > 2)
 ! We keep this as an example but it is instable in this case 
@@ -331,51 +329,42 @@
          !
       ENDIF
-            !                                               ! =====================
-         IF( cp_cfg == "orca" .AND. jp_cfg .eq. 1 ) THEN    ! ORCA R1 configuration
-            !                                               ! =====================
-
-            ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait 
-            ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : '
-
-            ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait 
-            ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : '
-
-            ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top) 
-            ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: Makassar (Top) : '
-
-            ii0 =  44   ;   ii1 =  44        ! Lombok Strait 
-            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : '
-
-            ii0 =  53   ;   ii1 =  53        ! Ombai Strait 
-            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : '
-
-            ii0 =  56   ;   ii1 =  56        ! Timor Passage 
-            ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: '
-
-            ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait 
-            ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: West Halmahera : '
-
-            ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait 
-            ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0  
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) '             orca_r1: East Halmahera : '
-
-            !
-            !
-         ENDIF
+      !
+      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration
+         !                                                 ! Increased lateral friction near of some straits
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : '
+         IF(lwp) WRITE(numout,*) '      Gibraltar '
+         ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait 
+         ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
+
+         IF(lwp) WRITE(numout,*) '      Bhosporus '
+         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait 
+         ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
+
+         IF(lwp) WRITE(numout,*) '      Makassar (Top) '
+         ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top) 
+         ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  3.0  
+
+         IF(lwp) WRITE(numout,*) '      Lombok '
+         ii0 =  44   ;   ii1 =  44        ! Lombok Strait 
+         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  2.0  
+
+         IF(lwp) WRITE(numout,*) '      Ombai '
+         ii0 =  53   ;   ii1 =  53        ! Ombai Strait 
+         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0  
+
+         IF(lwp) WRITE(numout,*) '      Timor Passage '
+         ii0 =  56   ;   ii1 =  56        ! Timor Passage 
+         ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2.0  
+
+         IF(lwp) WRITE(numout,*) '      West Halmahera '
+         ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait 
+         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0  
+
+         IF(lwp) WRITE(numout,*) '      East Halmahera '
+         ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait 
+         ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3.0  
+         !
+      ENDIF
       !
       CALL lbc_lnk( fmask, 'F', 1. )      ! Lateral boundary conditions on fmask
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 2392)
@@ -412,6 +412,5 @@
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
                !                                             ! =====================
-               IF( n_cla == 0 ) THEN
-                  !
+               IF( nn_cla == 0 ) THEN
                   ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open 
                   ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
@@ -422,5 +421,5 @@
                   END DO
                   IF(lwp) WRITE(numout,*)
-                  IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
+                  IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
                   !
                   ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
@@ -432,6 +431,5 @@
                   END DO
                   IF(lwp) WRITE(numout,*)
-                  IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
-                  !
+                  IF(lwp) WRITE(numout,*) '      orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
                ENDIF
                !
@@ -443,9 +441,9 @@
             CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
             CALL iom_close (inum)
-!                                                            ! =====================
+            !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
-               ii0 = 142   ;   ii1 = 142                     ! Close Halmera Strait  
-               ij0 =  51   ;   ij1 =  53                     ! =====================
-               DO ji = mi0(ii0), mi1(ii1)
+               ii0 = 142   ;   ii1 = 142                     ! =====================
+               ij0 =  51   ;   ij1 =  53                     
+               DO ji = mi0(ii0), mi1(ii1)                    ! Close Halmera Strait
                   DO jj = mj0(ij0), mj1(ij1)
                      bathy(ji,jj) = 0.0 
@@ -453,13 +451,12 @@
                END DO
                IF(lwp) WRITE(numout,*)
-               IF(lwp) WRITE(numout,*) '             orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
+               IF(lwp) WRITE(numout,*) '      orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
             ENDIF
             !                                                ! =====================
             IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
                !                                             ! =====================
-              IF( n_cla == 0 ) THEN
-                 !
-                 ii0 = 140   ;   ii1 = 140                  ! Gibraltar Strait open 
-                 ij0 = 102   ;   ij1 = 102                  ! (Thomson, Ocean Modelling, 1995)
+              IF( nn_cla == 0 ) THEN
+                 ii0 = 140   ;   ii1 = 140                   ! Gibraltar Strait open 
+                 ij0 = 102   ;   ij1 = 102                   ! (Thomson, Ocean Modelling, 1995)
                  DO ji = mi0(ii0), mi1(ii1)
                     DO jj = mj0(ij0), mj1(ij1)
@@ -468,8 +465,8 @@
                  END DO
                  IF(lwp) WRITE(numout,*)
-                 IF(lwp) WRITE(numout,*) '             orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
+                 IF(lwp) WRITE(numout,*) '      orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
                  !
-                 ii0 = 160   ;   ii1 = 160                  ! Bab el mandeb Strait open
-                 ij0 = 88    ;   ij1 = 88                   ! (Thomson, Ocean Modelling, 1995)
+                 ii0 = 160   ;   ii1 = 160                   ! Bab el mandeb Strait open
+                 ij0 = 88    ;   ij1 = 88                    ! (Thomson, Ocean Modelling, 1995)
                  DO ji = mi0(ii0), mi1(ii1)
                     DO jj = mj0(ij0), mj1(ij1)
@@ -479,5 +476,4 @@
                  IF(lwp) WRITE(numout,*)
                  IF(lwp) WRITE(numout,*) '             orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
-                 !
               ENDIF
               !
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtasal.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtasal.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtasal.F90	(revision 2392)
@@ -4,4 +4,10 @@
    !! Ocean data  :  read ocean salinity data from monthly atlas data
    !!=====================================================================
+   !! History :  OPA  ! 1991-03  ()  Original code
+   !!             -   ! 1992-07  (M. Imbard)
+   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT 
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
+   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
+   !!----------------------------------------------------------------------
 #if defined key_dtasal   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -10,5 +16,4 @@
    !!   dta_sal      : read ocean salinity data
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -23,12 +28,9 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dta_sal   ! called by step.F90 and inidta.F90
+   PUBLIC   dta_sal   ! called by step.F90 and inidta.F90
    
-   !! * Shared module variables
    LOGICAL , PUBLIC, PARAMETER ::   lk_dtasal = .TRUE.    !: salinity data flag
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   s_dta    !: salinity data at given time-step
 
-   !! * Module variables
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sal       ! structure of input SST (file informations, fields read)
 
@@ -38,12 +40,7 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   !!----------------------------------------------------------------------
-   !!   Default option:                                         NetCDF file
-   !!----------------------------------------------------------------------
 
    SUBROUTINE dta_sal( kt )
@@ -54,15 +51,10 @@
       !!             
       !! ** Method  : - Read on unit numsdt the monthly salinity data interpo-
-      !!     lated onto the model grid.
+      !!                lated onto the model grid.
       !!              - At each time step, a linear interpolation is applied
-      !!     between two monthly values.
-      !!
-      !! History :
-      !!        !  91-03  ()  Original code
-      !!        !  92-07  (M. Imbard)
-      !!   9.0  !  02-06  (G. Madec)  F90: Free form and module 
+      !!                between two monthly values.
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt             ! ocean time step
-      
+      !
       INTEGER ::   ji, jj, jk, jl, jkk            ! dummy loop indicies
       INTEGER ::   ik, ierror                     ! temporary integers
@@ -71,5 +63,4 @@
 #endif
       REAL(wp)::   zl
-      
 #if defined key_orca_lev10
       INTEGER ::   ikr, ikw, ikt, jjk 
@@ -80,6 +71,7 @@
       TYPE(FLD_N)              :: sn_sal
       LOGICAL , SAVE           :: linit_sal = .FALSE.
+      !!
+      NAMELIST/namdta_sal/   cn_dir, sn_sal
       !!----------------------------------------------------------------------
-      NAMELIST/namdta_sal/cn_dir,sn_sal
      
       ! 1. Initialization
@@ -91,9 +83,9 @@
          cn_dir = './'             ! directory in which the model is executed
          ! ... default values (NB: frequency positive => hours, negative => months)
-         !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
-         !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
-         sn_sal = FLD_N( 'salinity',  -1.  ,  'vosaline',  .false.   , .true.  ,  'monthly'  , ''       , ''         )
-
-         REWIND ( numnam )         ! ... read in namlist namdta_sal 
+         !            !   file    ! frequency ! variable ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
+         !            !   name    !  (hours)  !  name    !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
+         sn_sal = FLD_N( 'salinity',  -1.     ,'vosaline',  .false.   , .true.  ,  'monthly'  , ''       , ''         )
+
+         REWIND ( numnam )         ! read in namlist namdta_sal 
          READ( numnam, namdta_sal ) 
 
@@ -115,5 +107,5 @@
          IF( sn_sal%ln_tint )   ALLOCATE( sf_sal(1)%fdta(jpi,jpj,jpk,2) )
 #endif
-         ! fill sf_sal with sn_sal and control print
+         !                         ! fill sf_sal with sn_sal and control print
          CALL fld_fill( sf_sal, (/ sn_sal /), cn_dir, 'dta_sal', 'Salinity data', 'namdta_sal' )
          linit_sal = .TRUE.        
@@ -132,9 +124,6 @@
 
 #if defined key_tradmp
-      IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-   
-         !                                        ! =======================
-         !                                        !  ORCA_R2 configuration
-         !                                        ! =======================
+      IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN     !  ORCA_R2 configuration
+         !
          ij0 = 101   ;   ij1 = 109
          ii0 = 141   ;   ii1 = 155   
@@ -147,6 +136,6 @@
             END DO
          END DO
-
-         IF( n_cla == 1 ) THEN 
+         !
+         IF( nn_cla == 1 ) THEN 
             !                                         ! New salinity profile at Gibraltar
             il0 = 138   ;   il1 = 138   
@@ -230,5 +219,5 @@
          !                                  ! Mask
          s_dta(:,:,:) = s_dta(:,:,:) * tmask(:,:,:)
-         s_dta(:,:,jpk) = 0. 
+         s_dta(:,:,jpk) = 0.e0
          IF( ln_zps ) THEN               ! z-coord. partial steps
             DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
@@ -254,5 +243,5 @@
          CALL prihre(s_dta(:,:,jpkm1),jpi,jpj,1,jpi,20,1,jpj,20,1.,numout)
       ENDIF
-
+      !
    END SUBROUTINE dta_sal
 
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtatem.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtatem.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DTA/dtatem.F90	(revision 2392)
@@ -4,4 +4,10 @@
    !! Ocean data  :  read ocean temperature data from monthly atlas data
    !!=====================================================================
+   !! History :  OPA  ! 1991-03  ()  Original code
+   !!             -   ! 1992-07  (M. Imbard)
+   !!            8.0  ! 1999-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT 
+   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
+   !!            3.3  ! 2010-10  (C. Bricaud, S. Masson)  use of fldread
+   !!----------------------------------------------------------------------
 #if defined key_dtatem   ||   defined key_esopa
    !!----------------------------------------------------------------------
@@ -10,5 +16,4 @@
    !!   dta_tem      : read ocean temperature data
    !!---l-------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -22,12 +27,9 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC dta_tem   ! called by step.F90 and inidta.F90
-
-   !! * Shared module variables
+   PUBLIC   dta_tem    ! called by step.F90 and inidta.F90
+
    LOGICAL , PUBLIC, PARAMETER ::   lk_dtatem = .TRUE.   !: temperature data flag
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::  t_dta    !: temperature data at given time-step
 
-   !! * Module variables
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tem      ! structure of input SST (file informations, fields read)
 
@@ -37,12 +39,7 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   !!----------------------------------------------------------------------
-   !!   Default case                                            NetCDF file
-   !!----------------------------------------------------------------------
 
    SUBROUTINE dta_tem( kt )
@@ -62,13 +59,7 @@
       !!
       !! ** Action  :   define t_dta array at time-step kt
-      !!
-      !! History :
-      !!        !  91-03  ()  Original code
-      !!        !  92-07  (M. Imbard)
-      !!        !  99-10  (M.A. Foujols, M. Imbard)  NetCDF FORMAT 
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form and module
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-
+      !
       INTEGER ::   ji, jj, jk, jl, jkk            ! dummy loop indicies
       INTEGER ::   ik, ierror                     ! temporary integers
@@ -85,6 +76,7 @@
       TYPE(FLD_N)              ::   sn_tem
       LOGICAL , SAVE           ::   linit_tem = .FALSE.
+      !!
+      NAMELIST/namdta_tem/   cn_dir, sn_tem
       !!----------------------------------------------------------------------
-      NAMELIST/namdta_tem/cn_dir,sn_tem
  
       ! 1. Initialization 
@@ -96,9 +88,9 @@
          cn_dir = './'       ! directory in which the model is executed
          ! ... default values (NB: frequency positive => hours, negative => months)
-         !            !   file    ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
-         !            !   name    !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
-         sn_tem = FLD_N( 'temperature',  -1.  ,  'votemper',  .false.   , .true.  ,  'yearly'  , ''       , ''         )
-
-         REWIND( numnam )         ! ... read in namlist namdta_tem 
+         !            !   file    ! frequency ! variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !   name    !  (hours)  !  name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_tem = FLD_N( 'temperature',  -1.  , 'votemper',  .false.   , .true.  ,  'yearly'   , ''       , ''       )
+
+         REWIND( numnam )          ! read in namlist namdta_tem 
          READ( numnam, namdta_tem ) 
 
@@ -120,5 +112,5 @@
          IF( sn_tem%ln_tint )   ALLOCATE( sf_tem(1)%fdta(jpi,jpj,jpk,2) )
 #endif
-         ! fill sf_tem with sn_tem and control print
+         !                         ! fill sf_tem with sn_tem and control print
          CALL fld_fill( sf_tem, (/ sn_tem /), cn_dir, 'dta_tem', 'Temperature data', 'namdta_tem' )
          linit_tem = .TRUE.
@@ -138,8 +130,6 @@
          
 #if defined key_tradmp
-      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN
-         !                                        ! =======================
-         !                                        !  ORCA_R2 configuration
-         !                                        ! ======================= 
+      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN      !  ORCA_R2 configuration
+         !
          ij0 = 101   ;   ij1 = 109
          ii0 = 141   ;   ii1 = 155
@@ -151,6 +141,6 @@
             END DO
          END DO
-            
-         IF( n_cla == 1 ) THEN 
+         !
+         IF( nn_cla == 1 ) THEN 
             !                                         ! New temperature profile at Gibraltar
             il0 = 138   ;   il1 = 138
@@ -175,5 +165,4 @@
                END DO
             END DO
-            !
          ELSE
             !                                         ! Reduced temperature at Red Sea
@@ -251,21 +240,21 @@
                      t_dta(ji,jj,ik) = (1.-zl) * t_dta(ji,jj,ik) + zl * t_dta(ji,jj,ik-1)
                   ENDIF
-            END DO
-         END DO
-      ENDIF
-
-   ENDIF
-         
-   IF( lwp .AND. kt == nit000 ) THEN
-      WRITE(numout,*) ' temperature Levitus '
-      WRITE(numout,*)
-      WRITE(numout,*)'  level = 1'
-      CALL prihre( t_dta(:,:,1    ), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-      WRITE(numout,*)'  level = ', jpk/2
-      CALL prihre( t_dta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-      WRITE(numout,*)'  level = ', jpkm1
-      CALL prihre( t_dta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-   ENDIF
-
+               END DO
+            END DO
+         ENDIF
+         !
+      ENDIF
+         
+      IF( lwp .AND. kt == nit000 ) THEN
+         WRITE(numout,*) ' temperature Levitus '
+         WRITE(numout,*)
+         WRITE(numout,*)'  level = 1'
+         CALL prihre( t_dta(:,:,1    ), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+         WRITE(numout,*)'  level = ', jpk/2
+         CALL prihre( t_dta(:,:,jpk/2), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+         WRITE(numout,*)'  level = ', jpkm1
+         CALL prihre( t_dta(:,:,jpkm1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
+      ENDIF
+      !
    END SUBROUTINE dta_tem
 
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 2392)
@@ -4,4 +4,18 @@
    !! Ocean diagnostic variable : horizontal divergence and relative vorticity
    !!==============================================================================
+   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
+   !!            4.0  ! 1991-11  (G. Madec)
+   !!            6.0  ! 1993-03  (M. Guyon)  symetrical conditions
+   !!            7.0  ! 1996-01  (G. Madec)  s-coordinates
+   !!            8.0  ! 1997-06  (G. Madec)  lateral boundary cond., lbc
+   !!            8.1  ! 1997-08  (J.M. Molines)  Open boundaries
+   !!            8.2  ! 2000-03  (G. Madec)  no slip accurate
+   !!  NEMO      1.0  ! 2002-09  (G. Madec, E. Durand)  Free form, F90
+   !!             -   ! 2005-01  (J. Chanut) Unstructured open boundaries
+   !!             -   ! 2003-08  (G. Madec)  merged of cur and div, free form, F90
+   !!             -   ! 2005-01  (J. Chanut, A. Sellar) unstructured open boundaries
+   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
+   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -9,18 +23,17 @@
    !!                vorticity fields
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
+   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoff 
+   USE obc_oce         ! ocean lateral open boundary condition
+   USE cla             ! cross land advection             (cla_div routine)
    USE in_out_manager  ! I/O manager
-   USE obc_oce         ! ocean lateral open boundary condition
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE sbcrnf         ! river runoff 
-   USE sbc_oce, ONLY : ln_rnf   ! surface boundary condition: ocean
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC div_cur    ! routine called by step.F90 and istate.F90
+   PUBLIC   div_cur    ! routine called by step.F90 and istate.F90
 
    !! * Substitutions
@@ -30,7 +43,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -48,13 +60,11 @@
       !!      vorticity at before and now time-step
       !!
-      !! ** Method  : 
-      !!      I.  divergence :
+      !! ** Method  : I.  divergence :
       !!         - save the divergence computed at the previous time-step
       !!      (note that the Asselin filter has not been applied on hdivb)
       !!         - compute the now divergence given by :
       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      above expression
-      !!         - apply lateral boundary conditions on hdivn 
-      !!      II. vorticity :
+      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 
+      !!              II. vorticity :
       !!         - save the curl computed at the previous time-step
       !!            rotb = rotn
@@ -62,6 +72,4 @@
       !!         - compute the now curl in tensorial formalism:
       !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
-      !!         - apply lateral boundary conditions on rotn through a call
-      !!      of lbc_lnk routine.
       !!         - Coastal boundary condition: 'key_noslip_accurate' defined,
       !!      the no-slip boundary condition is computed using Schchepetkin
@@ -69,20 +77,11 @@
       !!      For example, along east coast, the one-sided finite difference
       !!      approximation used for di[v] is:
-      !!         di[e2v vn] =  1/(e1f*e2f)
-      !!                    * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
+      !!         di[e2v vn] =  1/(e1f*e2f) * ( (e2v vn)(i) + (e2v vn)(i-1) + (e2v vn)(i-2) )
       !!
       !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
       !!              - update rotb , rotn , the before & now rel. vorticity
-      !!
-      !! History :
-      !!   8.2  !  00-03  (G. Madec)  no slip accurate
-      !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
-      !!        !  05-01  (J. Chanut, A. Sellar) unstructured open boundaries
-      !! NEMO 3.3  !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      
-      !! * Local declarations
+      !
       INTEGER ::   ji, jj, jk     ! dummy loop indices
       INTEGER ::   ii, ij, jl     ! temporary integer
@@ -102,8 +101,8 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-
+         !
          hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
          rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
-
+         !
          !                                             ! --------
          ! Horizontal divergence                       !   div
@@ -198,6 +197,5 @@
             DO ji = 1, fs_jpim1   ! vector opt.
                rotn(ji,jj,jk) = (  zwv(ji+1,jj  ) - zwv(ji,jj)      &
-                                 - zwu(ji  ,jj+1) + zwu(ji,jj)  )   &
-                              * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
+                  &              - zwu(ji  ,jj+1) + zwu(ji,jj)  ) * fmask(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
             END DO
          END DO
@@ -228,16 +226,15 @@
                * ( -4. * zwu(ii,ij) + zwu(ii,ij-1) - 0.2 * zwu(ii,ij-2) )
          END DO
-
          !                                             ! ===============
       END DO                                           !   End of slab
       !                                                ! ===============
 
-      IF( ln_rnf )  CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
       
       ! 4. Lateral boundary conditions on hdivn and rotn
       ! ---------------------------------=======---======
-      CALL lbc_lnk( hdivn, 'T', 1. )     ! T-point, no sign change
-      CALL lbc_lnk( rotn , 'F', 1. )     ! F-point, no sign change
-
+      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )    ! lateral boundary cond. (no sign change)
+      !
    END SUBROUTINE div_cur
    
@@ -259,6 +256,5 @@
       !!      - compute the now divergence given by :
       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      above expression
-      !!      - apply lateral boundary conditions on hdivn 
+      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 
       !!              - Relavtive Vorticity :
       !!      - save the curl computed at the previous time-step (rotb = rotn)
@@ -266,6 +262,4 @@
       !!      - compute the now curl in tensorial formalism:
       !!            rotn = 1/(e1f*e2f) ( di[e2v vn] - dj[e1u un] )
-      !!      - apply lateral boundary conditions on rotn through a call to
-      !!      routine lbc_lnk routine.
       !!      Note: Coastal boundary condition: lateral friction set through
       !!      the value of fmask along the coast (see dommsk.F90) and shlat
@@ -274,19 +268,7 @@
       !! ** Action  : - update hdivb, hdivn, the before & now hor. divergence
       !!              - update rotb , rotn , the before & now rel. vorticity
-      !!
-      !! History :
-      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
-      !!   4.0  !  91-11  (G. Madec)
-      !!   6.0  !  93-03  (M. Guyon)  symetrical conditions
-      !!   7.0  !  96-01  (G. Madec)  s-coordinates
-      !!   8.0  !  97-06  (G. Madec)  lateral boundary cond., lbc
-      !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
-      !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
-      !!        !  05-01  (J. Chanut) Unstructured open boundaries
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
-      
-      !! * Local declarations
+      !
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
       REAL(wp) ::  zraur, zdep
@@ -302,8 +284,8 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-
+         !
          hdivb(:,:,jk) = hdivn(:,:,jk)    ! time swap of div arrays
          rotb (:,:,jk) = rotn (:,:,jk)    ! time swap of rot arrays
-
+         !
          !                                             ! --------
          ! Horizontal divergence                       !   div 
@@ -312,6 +294,6 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                hdivn(ji,jj,jk) =   &
-                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj  )*fse3u(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)       &
-                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji  ,jj-1)*fse3v(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  )    &
+                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk)       &
+                   + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk)  )    &
                   / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
             END DO  
@@ -349,11 +331,11 @@
       !                                                ! ===============
 
-      IF( ln_rnf )  CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
 
       ! 4. Lateral boundary conditions on hdivn and rotn
       ! ---------------------------------=======---======
-      CALL lbc_lnk( hdivn, 'T', 1. )       ! T-point, no sign change
-      CALL lbc_lnk( rotn , 'F', 1. )       ! F-point, no sign change
-
+      CALL lbc_lnk( hdivn, 'T', 1. )   ;   CALL lbc_lnk( rotn , 'F', 1. )     ! lateral boundary cond. (no sign change)
+      !
    END SUBROUTINE div_cur
    
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 2392)
@@ -164,12 +164,10 @@
       !                        ! Control of timestep choice
       IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN
-         IF( n_cla == 1 )   &
-           &   CALL ctl_stop( ' Crossland advection not implemented for this free surface formulation ' )
+         IF( nn_cla == 1 )   CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )
       ENDIF
 
-      !                        ! Control of momentum forulation
+      !                        ! Control of momentum formulation
       IF( lk_dynspg_ts .AND. lk_vvl ) THEN
-         IF( .NOT. ln_dynadv_vec )   &
-           &   CALL ctl_stop( ' Flux formulae not implemented for this free surface formulation ' )
+         IF( .NOT.ln_dynadv_vec )   CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )
       ENDIF
 
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 2392)
@@ -39,5 +39,5 @@
    USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine) 
    USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine)
-   USE cla_dynspg      ! cross land advection
+   USE cla             ! cross land advection
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing library
@@ -199,7 +199,5 @@
       CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
 #endif
-#if defined key_orca_r2
-      IF( n_cla == 1 )   CALL dyn_spg_cla( kt )      ! Cross Land Advection (update (ua,va))
-#endif
+      IF( nn_cla == 1 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
 
       ! compute the next vertically averaged velocity (effect of the additional force not included)
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 2392)
@@ -19,5 +19,4 @@
    USE domvvl          ! Variable volume
    USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE cla_div         ! cross land: hor. divergence      (div_cla routine)
    USE iom             ! I/O library
    USE restart         ! only for lrst_oce
@@ -147,8 +146,6 @@
       ENDIF
       !
-
-                         CALL div_cur( kt )           ! Horizontal divergence & Relative vorticity
-      IF( n_cla == 1 )   CALL div_cla( kt )           ! Cross Land Advection (Update Hor. divergence)
-
+      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      !
       z2dt = 2. * rdt                                 ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt =rdt
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 2392)
@@ -21,4 +21,5 @@
    USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
    USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine)
+   USE cla             ! cross land advection      (cla_traadv     routine)
    USE ldftra_oce      ! lateral diffusion coefficient on tracers
    USE in_out_manager  ! I/O manager
@@ -29,19 +30,18 @@
    PRIVATE
 
-   PUBLIC   tra_adv         ! routine called by step module
-   PUBLIC   tra_adv_init    ! routine called by opa module
+   PUBLIC   tra_adv        ! routine called by step module
+   PUBLIC   tra_adv_init   ! routine called by opa module
  
-   !                                                   !!* Namelist namtra_adv *
-   LOGICAL ::   ln_traadv_cen2   = .TRUE.       ! 2nd order centered scheme flag
-   LOGICAL ::   ln_traadv_tvd    = .FALSE.      ! TVD scheme flag
-   LOGICAL ::   ln_traadv_muscl  = .FALSE.      ! MUSCL scheme flag
-   LOGICAL ::   ln_traadv_muscl2 = .FALSE.      ! MUSCL2 scheme flag
-   LOGICAL ::   ln_traadv_ubs    = .FALSE.      ! UBS scheme flag
-   LOGICAL ::   ln_traadv_qck    = .FALSE.      ! QUICKEST scheme flag
+   !                                        !!* Namelist namtra_adv *
+   LOGICAL ::   ln_traadv_cen2   = .TRUE.    ! 2nd order centered scheme flag
+   LOGICAL ::   ln_traadv_tvd    = .FALSE.   ! TVD scheme flag
+   LOGICAL ::   ln_traadv_muscl  = .FALSE.   ! MUSCL scheme flag
+   LOGICAL ::   ln_traadv_muscl2 = .FALSE.   ! MUSCL2 scheme flag
+   LOGICAL ::   ln_traadv_ubs    = .FALSE.   ! UBS scheme flag
+   LOGICAL ::   ln_traadv_qck    = .FALSE.   ! QUICKEST scheme flag
 
    INTEGER ::   nadv   ! choice of the type of advection scheme
 
-   REAL(wp), DIMENSION(jpk) ::   r2dt  ! vertical profile time-step, = 2 rdttra
-      !                                ! except at nit000 (=rdttra) if neuler=0
+   REAL(wp), DIMENSION(jpk) ::   r2dt   ! vertical profile time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
 
    !! * Substitutions
@@ -51,7 +51,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -65,5 +64,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
+      !
       INTEGER ::   jk   ! dummy loop index
       REAL(wp), DIMENSION(jpi,jpj,jpk)   ::  zun, zvn, zwn   ! 3D workspace: effective transport
@@ -75,4 +74,7 @@
          r2dt(:) = 2. * rdttra(:)                      ! = 2 rdttra (leapfrog)
       ENDIF
+      !
+      IF( nn_cla == 1 )   CALL cla_traadv( kt )       !==  Cross Land Advection  ==! (hor. advection)
+      !
       !                                               !==  effective transport  ==!
       DO jk = 1, jpkm1
@@ -84,5 +86,5 @@
       !
       IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &
-         &                  CALL tra_adv_eiv( kt, zun, zvn, zwn, 'TRA' )      ! add the eiv transport (if necessary)
+         &              CALL tra_adv_eiv( kt, zun, zvn, zwn, 'TRA' )          ! add the eiv transport (if necessary)
       !
       CALL iom_put( "uoce_eff", zun )                                         ! output effective transport      
@@ -97,5 +99,4 @@
       CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS 
       CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, 'TRA', r2dt, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST 
-
       !
       CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
@@ -118,7 +119,6 @@
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask,               &
             &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         !
       END SELECT
-
+      !
       !                                              ! print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv  - Ta: ', mask1=tmask,               &
@@ -137,6 +137,6 @@
       INTEGER ::   ioptio
       !!
-      NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd,    &
-         &                 ln_traadv_muscl, ln_traadv_muscl2, &
+      NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd,     &
+         &                 ln_traadv_muscl, ln_traadv_muscl2,  &
          &                 ln_traadv_ubs  , ln_traadv_qck
       !!----------------------------------------------------------------------
@@ -156,5 +156,5 @@
          WRITE(numout,*) '      UBS    advection scheme        ln_traadv_ubs    = ', ln_traadv_ubs
          WRITE(numout,*) '      QUICKEST advection scheme      ln_traadv_qck    = ', ln_traadv_qck
-    ENDIF
+      ENDIF
 
       ioptio = 0                      ! Parameter control
@@ -168,7 +168,4 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' )
-
-      IF( n_cla == 1 .AND. .NOT. ln_traadv_cen2 )   &
-         &                CALL ctl_stop( 'cross-land advection only with 2nd order advection scheme' )
 
       !                              ! Set nadv
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRD/trdmld.F90	(revision 2392)
@@ -757,5 +757,5 @@
       END IF
 
-      IF( ( lk_trdmld ) .AND. ( n_cla == 1 ) ) THEN
+      IF( ( lk_trdmld ) .AND. ( nn_cla == 1 ) ) THEN
          WRITE(numout,cform_war)
          WRITE(numout,*) '                You set n_cla = 1. Note that the Mixed-Layer diagnostics  '
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla.F90	(revision 2392)
@@ -1,59 +1,107 @@
 MODULE cla
-   !!==============================================================================
-   !!                         ***  MODULE  cla  ***
-   !! Cross Land Advection : parameterize ocean exchanges through straits by a
-   !!                        specified advection across land.
-   !!==============================================================================
+   !!======================================================================
+   !!                    ***  MODULE  cla  ***
+   !! Cross Land Advection : specific update of the horizontal divergence,
+   !!                        tracer trends and after velocity
+   !!
+   !!                 ---   Specific to ORCA_R2   ---
+   !!
+   !!======================================================================
+   !! History :  1.0  ! 2002-11 (A. Bozec)  Original code
+   !!            3.2  ! 2009-07 (G. Madec)  merge cla, cla_div, tra_cla, cla_dynspg
+   !!                 !                     and correct a mpp bug reported by A.R. Porter
+   !!----------------------------------------------------------------------
 #if defined key_orca_r2
    !!----------------------------------------------------------------------
-   !!   'key_orca_r2'   :                             ORCA R2 configuration
+   !!   'key_orca_r2'                                 global ocean model R2
    !!----------------------------------------------------------------------
-   !!   tra_cla           : update the tracer trend with the horizontal 
-   !!                       and vertical advection trends at straits
-   !!   tra_bab_el_mandeb : 
-   !!   tra_gibraltar     :
-   !!   tra_hormuz        :
-   !!   tra_cla_init      :
+   !!   cla_div           : update of horizontal divergence at cla straits
+   !!   tra_cla           : update of tracers at cla straits
+   !!   cla_dynspg        : update of after horizontal velocities at cla straits
+   !!   cla_init          : initialisation - control check
+   !!   cla_bab_el_mandeb : cross land advection for Bab-el-mandeb strait
+   !!   cla_gibraltar     : cross land advection for Gibraltar strait
+   !!   cla_hormuz        : cross land advection for Hormuz strait
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! distributed memory computing
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE dynspg_oce     ! ocean dynamics: surface pressure gradient variables
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! distributed memory computing library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
    PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC tra_cla        ! routine called by step.F90
-   PUBLIC tra_cla_init   ! routine called by opa.F90
-
-   !! * Modules variables   
-   REAL(wp) :: zempmed, zempred
-
-   REAL(wp) :: zisw_rs, zurw_rs, zbrw_rs          ! Imposed transport Red Sea 
-   REAL(wp) :: zisw_ms, zmrw_ms, zurw_ms, zbrw_ms ! Imposed transport Med Sea 
-   REAL(wp) :: zisw_pg, zbrw_pg                   ! Imposed transport Persic Gulf 
-
-   REAL(wp), DIMENSION(jpk) ::   &
-      zu1_rs_i, zu2_rs_i, zu3_rs_i,              &  ! Red Sea velocities
-      zu1_ms_i, zu2_ms_i, zu3_ms_i,              &  ! Mediterranean Sea velocities
-      zu_pg                                         ! Persic Gulf velocities
-   REAL(wp), DIMENSION (jpk) :: zthor, zshor        ! Temperature, salinity Hormuz 
+   
+   PUBLIC   cla_init     ! routine called by opa.F90
+   PUBLIC   cla_div      ! routine called by divcur.F90
+   PUBLIC   cla_traadv   ! routine called by traadv.F90
+   PUBLIC   cla_dynspg   ! routine called by dynspg_flt.F90
+
+   INTEGER  ::   nbab, ngib, nhor   ! presence or not of required grid-point on local domain
+   !                                ! for Bab-el-Mandeb, Gibraltar, and Hormuz straits
+   
+   !                                                              !!! profile of hdiv for some straits
+   REAL(wp), DIMENSION (jpk) ::   hdiv_139_101, hdiv_139_101_kt    ! Gibraltar     strait, fixed & time evolving part (i,j)=(172,101)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_139_102                     ! Gibraltar     strait, fixed part only            (i,j)=(139,102)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_141_102, hdiv_141_102_kt    ! Gibraltar     strait, fixed & time evolving part (i,j)=(141,102)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_161_88 , hdiv_161_88_kt     ! Bab-el-Mandeb strait, fixed & time evolving part (i,j)=(161,88)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_161_87                      ! Bab-el-Mandeb strait, fixed part only            (i,j)=(161,87)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_160_89 , hdiv_160_89_kt     ! Bab-el-Mandeb strait, fixed & time evolving part (i,j)=(160,89)
+   REAL(wp), DIMENSION (jpk) ::   hdiv_172_94                      ! Hormuz        strait, fixed part only            (i,j)=(172, 94)
+
+   REAL(wp), DIMENSION (jpk) ::   t_171_94_hor, s_171_94_hor       ! Temperature, salinity in the Hormuz strait
    
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
-   SUBROUTINE tra_cla( kt )
+   SUBROUTINE cla_div( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE div_cla  ***
+      !!
+      !! ** Purpose :   update the horizontal divergence of the velocity field
+      !!              at some straits ( Gibraltar, Bab el Mandeb and Hormuz ).
+      !!
+      !! ** Method  : - first time-step: initialisation of cla
+      !!              - all   time-step: using imposed transport at each strait, 
+      !!              the now horizontal divergence is updated
+      !!
+      !! ** Action  :   phdivn   updted now horizontal divergence at cla straits
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !!----------------------------------------------------------------------
+      !     
+      IF( kt == nit000 ) THEN
+         !
+         CALL cla_init                                        ! control check 
+         !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'div_cla : cross land advection on hdiv '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~'
+         !
+         IF( nbab == 1 )   CALL cla_bab_el_mandeb('ini')    ! Bab el Mandeb ( Red Sea - Indian ocean )
+         IF( ngib == 1 )   CALL cla_gibraltar    ('ini')    ! Gibraltar strait (Med Sea - Atlantic ocean)
+         IF( nhor == 1 )   CALL cla_hormuz       ('ini')    ! Hormuz Strait ( Persian Gulf - Indian ocean )
+         !
+      ENDIF                           
+      !
+      IF( nbab == 1    )   CALL cla_bab_el_mandeb('div')    ! Bab el Mandeb ( Red Sea - Indian ocean )
+      IF( ngib == 1    )   CALL cla_gibraltar    ('div')    ! Gibraltar strait (Med Sea - Atlantic ocean)
+      IF( nhor == 1    )   CALL cla_hormuz       ('div')    ! Hormuz Strait ( Persian Gulf - Indian ocean )
+      !
+!!gm  lbc useless here, no?
+!!gm      CALL lbc_lnk( hdivn, 'T', 1. )                    ! Lateral boundary conditions on hdivn
+      !
+   END SUBROUTINE cla_div
+   
+   
+   SUBROUTINE cla_traadv( kt )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_cla  ***
@@ -63,716 +111,632 @@
       !!      at some straits ( Bab el Mandeb, Gibraltar, Hormuz ).
       !!
-      !! ** Method  :   ...
-      !!         Add this trend now to the general trend of tracer (ta,sa):
-      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
-      !!
-      !! ** Action  :   update (ta,sa) with the now advective tracer trends
-      !!
-      !! History :
-      !!        !         (A. Bozec)  original code  
-      !!   8.5  !  02-11  (A. Bozec)  F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Arguments
+      !! ** Method  :   using both imposed transport at each strait and T & S
+      !!              budget, the now tracer trends is updated
+      !!
+      !! ** Action  :   (ta,sa)   updated now tracer trends at cla straits
+      !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
       !!----------------------------------------------------------------------
- 
-      ! Bab el Mandeb strait horizontal advection
-
-      CALL tra_bab_el_mandeb 
-
-      ! Gibraltar strait horizontal advection
-
-      CALL tra_gibraltar
-
-      ! Hormuz Strait ( persian Gulf) horizontal advection
-
-      CALL tra_hormuz 
-
-   END SUBROUTINE tra_cla
-
-
-   SUBROUTINE tra_bab_el_mandeb
+      !
+      IF( kt == nit000 ) THEN 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_cla : cross land advection on tracers '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~'
+      ENDIF
+      !
+      IF( nbab == 1    )   CALL cla_bab_el_mandeb('tra')      ! Bab el Mandeb strait
+      IF( ngib == 1    )   CALL cla_gibraltar    ('tra')      ! Gibraltar strait
+      IF( nhor == 1    )   CALL cla_hormuz       ('tra')      ! Hormuz Strait ( Persian Gulf)
+      !
+   END SUBROUTINE cla_traadv
+
+   
+   SUBROUTINE cla_dynspg( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE cla_dynspg  ***
+      !!                   
+      !! ** Purpose :   Update the after velocity at some straits 
+      !!              (Bab el Mandeb, Gibraltar, Hormuz).
+      !!
+      !! ** Method  :   required to compute the filtered surface pressure gradient 
+      !!
+      !! ** Action  :   (ua,va)   after velocity at the cla straits
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 ) THEN 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'cla_dynspg : cross land advection on (ua,va) '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+      ENDIF
+      !
+      IF( nbab == 1    )   CALL cla_bab_el_mandeb('spg')      ! Bab el Mandeb strait
+      IF( ngib == 1    )   CALL cla_gibraltar    ('spg')      ! Gibraltar strait
+      IF( nhor == 1    )   CALL cla_hormuz       ('spg')      ! Hormuz Strait ( Persian Gulf)
+      !
+!!gm lbc is needed here, not?
+!!gm      CALL lbc_lnk( hdivn, 'U', -1. )   ;   CALL lbc_lnk( hdivn, 'V', -1. )      ! Lateral boundary conditions 
+      !
+   END SUBROUTINE cla_dynspg
+
+
+   SUBROUTINE cla_init
+      !! -------------------------------------------------------------------
+      !!                   ***  ROUTINE cla_init  ***
+      !!           
+      !! ** Purpose :   control check for mpp computation  
+      !!
+      !! ** Method  : - All the strait grid-points must be inside one of the 
+      !!              local domain interior for the cla advection to work
+      !!              properly in mpp (i.e. inside (2:jpim1,2:jpjm1) ).
+      !!              Define the corresponding indicators (nbab, ngib, nhor)
+      !!              - The profiles of cross-land fluxes are currently hard
+      !!              coded for L31 levels. Stop if jpk/=31
+      !!
+      !! ** Action  :   nbab, ngib, nhor   strait inside the local domain or not
       !!---------------------------------------------------------------------
-      !!             ***  ROUTINE tra_bab_el_mandeb  ***
-      !!
-      !! ** Purpose :   Update the horizontal advective trend of tracers
-      !!      correction in Bab el Mandeb strait and
-      !!      add it to the general trend of tracer equations.
+      REAL(wp) ::   ztemp
+      !!---------------------------------------------------------------------
+      !
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'cla_init : cross land advection initialisation '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      !
+      IF( .NOT.lk_dynspg_flt )   CALL ctl_stop( 'cla_init: Cross Land Advection works only with lk_dynspg_flt=T ' )
+      !
+      IF( lk_vvl    )   CALL ctl_stop( 'cla_init: Cross Land Advection does not work with lk_vvl=T option' )
+      !
+      IF( jpk /= 31 )   CALL ctl_stop( 'cla_init: Cross Land Advection hard coded for ORCA_R2_L31' )
+      !
+      !                                        _|_______|_______|_
+      !                                     89  |       |///////|  
+      !                                        _|_______|_______|_
+      ! ------------------------ !          88  |///////|       | 
+      !   Bab el Mandeb strait   !             _|_______|_______|_
+      ! ------------------------ !          87  |///////|       | 
+      !                                        _|_______|_______|_
+      !                                         |  160  |  161  |  
+      !
+      ! The 6 Bab el Mandeb grid-points must be inside one of the interior of the
+      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
+      nbab = 0
+      IF(  ( 1 <= mj0( 88) .AND. mj1( 89) <= jpj ) .AND.    &  !* (161,89), (161,88) and (161,88) on the local pocessor
+         & ( 1 <= mi0(160) .AND. mi1(161) <= jpi )       )    nbab = 1 
+      !
+      ! test if there is no local domain that includes all required grid-points
+      ztemp = REAL( nbab )
+      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
+      IF( ztemp == 0 ) THEN                     ! Only 2 points in each direction, this should never be a problem
+         CALL ctl_stop( ' cross land advection at Bab-el_Mandeb does not work with your processor cutting: change it' )
+      ENDIF
+      !                                        ___________________________
+      ! ------------------------ !         102  |       |///////|       |
+      !     Gibraltar strait     !             _|_______|_______|_______|_
+      ! ------------------------ !         101  |       |///////|       |
+      !                                        _|_______|_______|_______|_ 
+      !                                         |  139  |  140  |  141  |
+      !
+      ! The 6 Gibraltar grid-points must be inside one of the interior of the
+      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
+      ngib = 0
+      IF(  ( 2 <= mj0(101) .AND. mj1(102) <= jpjm1 ) .AND.    &  !* (139:141,101:102) on the local pocessor
+         & ( 2 <= mi0(139) .AND. mi1(141) <= jpim1 )       )    ngib = 1 
+      !
+      ! test if there is no local domain that includes all required grid-points
+      ztemp = REAL( ngib )
+      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
+      IF( ztemp == 0 ) THEN                     ! 3 points in i-direction, this may be a problem with some cutting
+           CALL ctl_stop( ' cross land advection at Gibraltar does not work with your processor cutting: change it' )
+      ENDIF
+      !                                        _______________
+      ! ------------------------ !          94  |/////|     | 
+      !       Hormuz strait      !             _|_____|_____|_
+      ! ------------------------ !                171   172     
+      !           
+      ! The 2 Hormuz grid-points must be inside one of the interior of the
+      ! local domain for the cla advection to work properly (i.e. (2:jpim1,2:jpjm1)
+      nhor = 0
+      IF(    2 <= mj0( 94) .AND. mj1( 94) <= jpjm1  .AND.  & 
+         &   2 <= mi0(171) .AND. mi1(172) <= jpim1         )   nhor = 1 
+      !
+      ! test if there is no local domain that includes all required grid-points
+      ztemp = REAL( nhor )
+      IF( lk_mpp )   CALL mpp_sum( ztemp )      ! sum with other processors value
+      IF( ztemp == 0 ) THEN                     ! 3 points in i-direction, this may be a problem with some cutting
+           CALL ctl_stop( ' cross land advection at Hormuz does not work with your processor cutting: change it' )
+      ENDIF
+      !
+   END SUBROUTINE cla_init
+
+
+   SUBROUTINE cla_bab_el_mandeb( cd_td )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE cla_bab_el_mandeb  ***
+      !!       
+      !! ** Purpose :   update the now horizontal divergence, the tracer tendancy
+      !!              and the after velocity in vicinity of Bab el Mandeb ( Red Sea - Indian ocean).
+      !!
+      !! ** Method  :   compute the exchanges at each side of the strait :
+      !!
+      !!       surf. zio_flow
+      !! (+ balance of emp) /\  |\\\\\\\\\\\|
+      !!                    ||  |\\\\\\\\\\\|  
+      !!    deep zio_flow   ||  |\\\\\\\\\\\|  
+      !!            |  ||   ||  |\\\\\\\\\\\|  
+      !!        89  |  ||   ||  |\\\\\\\\\\\|  
+      !!            |__\/_v_||__|____________ 
+      !!            !\\\\\\\\\\\|          surf. zio_flow
+      !!            |\\\\\\\\\\\|<===    (+ balance of emp)
+      !!            |\\\\\\\\\\\u
+      !!        88  |\\\\\\\\\\\|<---      deep  zrecirc (upper+deep at 2 different levels)
+      !!            |___________|__________   
+      !!            !\\\\\\\\\\\|         
+      !!            |\\\\\\\\\\\| ---\     deep  zrecirc (upper+deep) 
+      !!        87  !\\\\\\\\\\\u ===/   + deep  zio_flow   (all at the same level)
+      !!            !\\\\\\\\\\\|  
+      !!            !___________|__________ 
+      !!                160         161
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='div' update the divergence
+      !                                         ! ='tra' update the tracers
+      !                                         ! ='spg' update after velocity
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zemp_red     ! temporary scalar
+      REAL(wp) ::   zio_flow, zrecirc_upp, zrecirc_mid, zrecirc_bot
+      !!---------------------------------------------------------------------
+      !
+      SELECT CASE( cd_td ) 
+      !                     ! ---------------- !
+      CASE( 'ini' )         !  initialisation  ! 
+         !                  ! ---------------- ! 
+         !                                   
+         zio_flow    = 0.4e6                       ! imposed in/out flow
+         zrecirc_upp = 0.2e6                       ! imposed upper recirculation water
+         zrecirc_bot = 0.5e6                       ! imposed bottom  recirculation water
+
+         hdiv_161_88(:) = 0.e0                     ! (161,88) Gulf of Aden side, north point
+         hdiv_161_87(:) = 0.e0                     ! (161,87) Gulf of Aden side, south point
+         hdiv_160_89(:) = 0.e0                     ! (160,89) Red sea side
+
+         DO jj = mj0(88), mj1(88)              !** profile of hdiv at (161,88)   (Gulf of Aden side, north point)
+            DO ji = mi0(161), mi1(161)         !------------------------------
+               DO jk = 1, 8                        ! surface in/out flow   (Ind -> Red)   (div >0)
+                  hdiv_161_88(jk) = + zio_flow / ( 8. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               !                                   ! recirculation water   (Ind -> Red)   (div >0)
+               hdiv_161_88(20) =                 + zrecirc_upp   / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,20) )
+               hdiv_161_88(21) = + ( zrecirc_bot - zrecirc_upp ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,21) )
+            END DO
+         END DO
+         !
+         DO jj = mj0(87), mj1(87)              !** profile of hdiv at (161,88)   (Gulf of Aden side, south point)
+            DO ji = mi0(161), mi1(161)         !------------------------------
+               !                                   ! deep out flow + recirculation   (Red -> Ind)   (div <0)
+               hdiv_161_87(21) = - ( zio_flow + zrecirc_bot ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,21) )
+            END DO
+         END DO
+         !
+         DO jj = mj0(89), mj1(89)              !** profile of hdiv at (161,88)   (Red sea side)
+            DO ji = mi0(160), mi1(160)         !------------------------------
+               DO jk = 1, 8                        ! surface inflow    (Ind -> Red)   (div <0)
+                  hdiv_160_89(jk) = - zio_flow / ( 8. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               !                                   ! deep    outflow   (Red -> Ind)   (div >0)
+               hdiv_160_89(16)    = + zio_flow / (      e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,16) )
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
+         !                  ! ---------=====-- ! 
+         !                                     !** emp on the Red Sea   (div >0) 
+         zemp_red = 0.e0                       !---------------------
+         DO jj = mj0(87), mj1(96)                  ! sum over the Red sea
+            DO ji = mi0(148), mi1(160) 
+               zemp_red = zemp_red + emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( zemp_red )   ! sum with other processors value
+         zemp_red = zemp_red * 1.e-3               ! convert in m3
+         !
+         !                                     !** Correct hdivn (including emp adjustment)
+         !                                     !-------------------------------------------
+         DO jj = mj0(88), mj1(88)                  !* profile of hdiv at (161,88)   (Gulf of Aden side, north point)
+            DO ji = mi0(161), mi1(161) 
+               hdiv_161_88_kt(:) = hdiv_161_88(:)
+               DO jk = 1, 8                              ! increase the inflow from the Indian   (div >0) 
+                  hdiv_161_88_kt(jk) = hdiv_161_88(jk) + zemp_red / (8. * e2u(ji,jj) * fse3u(ji,jj,jk) )
+               END DO
+               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_161_88_kt(:)
+            END DO
+         END DO
+         DO jj = mj0(87), mj1(87)                  !* profile of divergence at (161,87)   (Gulf of Aden side, south point)
+            DO ji = mi0(161), mi1(161) 
+               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_161_87(:)
+            END DO
+         END DO
+         DO jj = mj0(89), mj1(89)                  !* profile of divergence at (160,89)   (Red sea side)
+            DO ji = mi0(160), mi1(160) 
+               hdiv_160_89_kt(:) = hdiv_160_89(:)
+               DO jk = 1, 18                              ! increase the inflow from the Indian   (div <0) 
+                  hdiv_160_89_kt(jk) = hdiv_160_89(jk) - zemp_red / (10. * e1v(ji,jj) * fse3v(ji,jj,jk) )
+               END DO
+               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_160_89_kt(:)
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
+         !                  ! --------=======- !
+         !
+         DO jj = mj0(88), mj1(88)              !** (161,88)   (Gulf of Aden side, north point)
+            DO ji = mi0(161), mi1(161) 
+               DO jk = 1, jpkm1                         ! surf inflow + reciculation (from Gulf of Aden)
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_161_88_kt(jk) * tn(ji,jj,jk)
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_161_88_kt(jk) * sn(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         DO jj = mj0(87), mj1(87)              !** (161,87)   (Gulf of Aden side, south point)
+            DO ji = mi0(161), mi1(161) 
+               jk =  21                                 ! deep outflow + recirulation (combined flux)
+               ta(ji,jj,jk) = ta(ji,jj,jk) + hdiv_161_88(20) * tn(ji  ,jj+1,20)   &  ! upper recirculation from Gulf of Aden
+                  &                        + hdiv_161_88(21) * tn(ji  ,jj+1,21)   &  ! deep  recirculation from Gulf of Aden
+                  &                        + hdiv_160_89(16) * tn(ji-1,jj+2,16)      ! deep inflow from Red sea
+               sa(ji,jj,jk) = sa(ji,jj,jk) + hdiv_161_88(20) * sn(ji  ,jj+1,20)   &
+                  &                        + hdiv_161_88(21) * sn(ji  ,jj+1,21)   &
+                  &                        + hdiv_160_89(16) * sn(ji-1,jj+2,16)   
+            END DO
+         END DO
+         DO jj = mj0(89), mj1(89)              !** (161,88)   (Red sea side)
+            DO ji = mi0(160), mi1(160)
+               DO jk = 1, 14                            ! surface inflow (from Gulf of Aden)
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_160_89_kt(jk) * tn(ji+1,jj-1,jk)
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_160_89_kt(jk) * sn(ji+1,jj-1,jk)
+               END DO
+               !                                        ! deep    outflow (from Red sea)
+               ta(ji,jj,16) = ta(ji,jj,16) - hdiv_160_89(jk) * tn(ji,jj,jk)
+               sa(ji,jj,16) = sa(ji,jj,16) - hdiv_160_89(jk) * sn(ji,jj,jk)
+            END DO
+         END DO
+         !
+         !                  ! ---------------- !
+      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
+         !                  ! --------=======- !
+         ! at this stage, (ua,va) are the after velocity, not the tendancy
+         ! compute the velocity from the divergence at T-point
+         !
+         DO jj = mj0(88), mj1(88)              !** (160,88)   (Gulf of Aden side, north point)
+            DO ji = mi0(160), mi1(160)                   ! 160, not 161 as it is a U-point) 
+               ua(ji,jj,:) = - hdiv_161_88_kt(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
+                  &                              * e2u(ji,jj) * fse3u(ji,jj,:)
+            END DO
+         END DO
+         DO jj = mj0(87), mj1(87)              !** (160,87)   (Gulf of Aden side, south point)
+            DO ji = mi0(160), mi1(160)                   ! 160, not 161 as it is a U-point) 
+               ua(ji,jj,:) = - hdiv_161_87(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
+                  &                           * e2u(ji,jj) * fse3u(ji,jj,:)
+            END DO
+         END DO
+         DO jj = mj0(88), mj1(88)              !** profile of divergence at (160,89)   (Red sea side)
+            DO ji = mi0(160), mi1(160)                   ! 88, not 89 as it is a V-point)
+               va(ji,jj,:) = - hdiv_160_89_kt(:) / ( e1t(ji,jj+1) * e2t(ji,jj+1) * fse3t(ji,jj+1,:) )   &
+                  &                              * e1v(ji,jj) * fse3v(ji,jj,:)
+            END DO
+         END DO
+      END SELECT
+      !
+   END SUBROUTINE cla_bab_el_mandeb
+   
+
+   SUBROUTINE cla_gibraltar( cd_td )
+      !! -------------------------------------------------------------------
+      !!                 ***  ROUTINE cla_gibraltar  ***
+      !!        
+      !! ** Purpose :   update the now horizontal divergence, the tracer 
+      !!              tendancyand the after velocity in vicinity of Gibraltar 
+      !!              strait ( Persian Gulf - Indian ocean ).
       !!
       !! ** Method :
-      !!     We impose transport at Bab el Mandeb and knowing T and S in
-      !!     surface and depth at each side of the  strait, we deduce T and S
-      !!     of the deep outflow of the Red Sea in the Indian ocean . 
-      !!                                          |
-      !!            |/ \|            N          |\ /|
-      !!            |_|_|______      |          |___|______
-      !!        88  |   |<-       W - - E    88 |   |<-
-      !!        87  |___|______      |       87 |___|->____
-      !!             160 161         S           160 161 
-      !!       horizontal view                horizontal view
-      !!          surface                        depth
-      !!    
-      !!     The horizontal advection is evaluated by a second order cen-
-      !!     tered scheme using now fields (leap-frog scheme). In specific
-      !!     areas (vicinity of major river mouths, some straits, or tn
-      !!     approaching the freezing point) it is mixed with an upstream 
-      !!     scheme for stability reasons. 
-      !!
-      !!         C A U T I O N : the trend saved is the centered trend only.
-      !!      It doesn't take into account the upstream part of the scheme.
-      !!
-      !! ** history :
-      !!           !  02-11  (A. Bozec) Original code 
-      !!      8.5  !  02-11  (A. Bozec) F90: Free form and module
+      !!                     _______________________            
+      !!      deep  zio_flow /====|///////|====> surf. zio_flow
+      !!    + deep  zrecirc  \----|///////|     (+balance of emp)
+      !! 102                      u///////u
+      !!      mid.  recicul    <--|///////|<==== deep  zio_flow
+      !!                     _____|_______|_____  
+      !!      surf. zio_flow ====>|///////|       
+      !!    (+balance of emp)     |///////|
+      !! 101                      u///////|             
+      !!      mid.  recicul    -->|///////|               Caution: zrecirc split into 
+      !!      deep  zrecirc  ---->|///////|                  upper & bottom recirculation
+      !!                   _______|_______|_______ 
+      !!                     139     140     141  
+      !!
       !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::  ji, jj, jk               ! dummy loop indices
-      REAL(wp) :: zsu, zvt
-      REAL(wp) :: zsumt, zsumt1, zsumt2, zsumt3, zsumt4
-      REAL(wp) :: zsums, zsums1, zsums2, zsums3, zsums4
-      REAL(wp) :: zt, zs
-      REAL(wp) :: zwei
-      REAL(wp), DIMENSION (jpk) ::  zu1_rs, zu2_rs, zu3_rs
+      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='div' update the divergence
+      !                                         ! ='tra' update the tracers
+      !                                         ! ='spg' update after velocity
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zemp_med     ! temporary scalar
+      REAL(wp) ::   zio_flow, zrecirc_upp, zrecirc_mid, zrecirc_bot
       !!---------------------------------------------------------------------
-
-      ! Initialization of vertical sum for T and S transport
-      ! ----------------------------------------------------
-
-      zsumt  = 0.e0       ! East  Bab el Mandeb surface north point (T)
-      zsums  = 0.e0       ! East  Bab el Mandeb surface north point (S)
-      zsumt1 = 0.e0       ! East  Bab el Mandeb depth   south point (T)
-      zsums1 = 0.e0       ! East  Bab el Mandeb depth   south point (S)
-      zsumt2 = 0.e0       ! West  Bab el Mandeb surface             (T)
-      zsums2 = 0.e0       ! West  Bab el Mandeb surface             (S)
-      zsumt3 = 0.e0       ! West  Bab el Mandeb depth               (T)
-      zsums3 = 0.e0       ! West  Bab el Mandeb depth               (S)
-      zsumt4 = 0.e0       ! East  Bab el Mandeb depth   north point (T) 
-      zsums4 = 0.e0       ! East  Bab el Mandeb depth   north point (S) 
-      
-      ! EMP of the Red Sea 
-      ! ------------------
-
-      zempred = 0.e0
-      zwei = 0.e0
-      DO jj = mj0(87), mj1(96)
-         DO ji = mi0(148), mi1(160)
-            zwei    = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempred = zempred + ( emp(ji,jj) - rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempred )      ! sum with other processors value
-
-      ! convert in m3
-      zempred = zempred * 1.e-3
-
-      ! Velocity profile at each point
-      ! ------------------------------
-
-      zu1_rs(:) = zu1_rs_i(:)
-      zu2_rs(:) = zu2_rs_i(:)
-      zu3_rs(:) = zu3_rs_i(:)
-
-      ! velocity profile at 161,88 East Bab el Mandeb North point 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1, 8                                      
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu1_rs(jk) = zu1_rs(jk) - ( zempred / 8. ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )     
-            END DO
-         END DO
-      END DO
-
-      ! velocity profile at 161, 88 West Bab el Mandeb 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  10                                     
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu3_rs(jk) = zu3_rs(jk) + ( zempred / 10. ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-      
-      ! Balance of temperature and salinity
-      ! -----------------------------------
-
-      ! east Bab el Mandeb surface vertical sum of transport* S,T
-      DO jk =  1, 19
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(161), mi1(161) 
-               zsumt  = zsumt  + tn(ji,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zu1_rs(jk)  
-               zsums  = zsums  + sn(ji,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zu1_rs(jk) 
-            END DO
-         END DO
-      END DO
-
-      ! west Bab el Mandeb surface vertical sum of transport* S,T
-      DO jk =  1, 10
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(161), mi1(161) 
-               zsumt2 = zsumt2 + tn(ji,jj,jk) * e1v(ji-1,jj) * fse3v(ji-1,jj,jk) * zu3_rs(jk)
-               zsums2 = zsums2 + sn(ji,jj,jk) * e1v(ji-1,jj) * fse3v(ji-1,jj,jk) * zu3_rs(jk)
-            END DO
-         END DO
-      END DO
-
-      ! west Bab el Mandeb deeper
-      DO jj = mj0(89), mj1(89) 
-         DO ji = mi0(160), mi1(160) 
-            zsumt3 = tn(ji,jj,16) * e1v(ji,jj-1) * fse3v(ji,jj-1,16) * zu3_rs(16)
-            zsums3 = sn(ji,jj,16) * e1v(ji,jj-1) * fse3v(ji,jj-1,16) * zu3_rs(16)
-         END DO
-      END DO
-
-      ! east  Bab el Mandeb deeper  
-      DO jk =  20, 21
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(161), mi1(161) 
-               zsumt4 =  zsumt4 + tn(ji,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zu1_rs(jk)
-               zsums4 =  zsums4 + sn(ji,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * zu1_rs(jk) 
-            END DO
-         END DO
-      END DO
-
-      ! Total transport  
-      zsumt1 = -( zsumt3 + zsumt2 + zsumt + zsumt4 )
-      zsums1 = -( zsums3 + zsums2 + zsums + zsums4 )
-
-      ! Temperature and Salinity at East Bab el Mandeb, Level 21
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zt = zsumt1 / ( zu2_rs(21) * e2u(ji,jj-1) * fse3u(ji,jj-1,21) )
-            zs = zsums1 / ( zu2_rs(21) * e2u(ji,jj-1) * fse3u(ji,jj-1,21) )
-         END DO
-      END DO
-      
-      ! New Temperature and Salinity at East Bab el Mandeb
-      ! --------------------------------------------------
-
-      ! north point 	
-      DO jk = 1, jpk
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(161), mi1(161) 
-               zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e2u(ji-1,jj) * fse3u(ji-1,jj,jk)
-               ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu1_rs(jk) * tn(ji,jj,jk)
-               sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu1_rs(jk) * sn(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! south point
-      jk =  21
-      DO jj = mj0(87), mj1(87) 
-         DO ji = mi0(161), mi1(161) 
-            zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-            zsu = e2u(ji-1,jj) * fse3u(ji-1,jj,jk)
-            ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu2_rs(jk) * zt
-            sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu2_rs(jk) * zs
-         END DO
-      END DO
-
-
-      ! New Temperature and Salinity at West Bab el Mandeb 
-      ! --------------------------------------------------
-
-      ! surface	
-      DO jk = 1, 10
-         DO jj = mj0(89), mj1(89) 
-            DO ji = mi0(160), mi1(160) 
-               zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e1v(ji,jj-1) * fse3v(ji,jj-1,jk)
-               ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu3_rs(jk) * tn(ji+1,jj-1,jk)
-               sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu3_rs(jk) * sn(ji+1,jj-1,jk)
-            END DO
-         END DO
-      END DO
-      ! deeper
-      jk =  16
-      DO jj = mj0(89), mj1(89) 
-         DO ji = mi0(160), mi1(160) 
-            zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-            zsu = e1v(ji,jj-1) * fse3v(ji,jj-1,jk)
-            ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu3_rs(jk) * tn(ji,jj,jk)
-            sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu3_rs(jk) * sn(ji,jj,jk)
-         END DO
-      END DO
-
-   END SUBROUTINE tra_bab_el_mandeb
-
-
-   SUBROUTINE tra_gibraltar
+      !
+      SELECT CASE( cd_td ) 
+      !                     ! ---------------- !
+      CASE( 'ini' )         !  initialisation  ! 
+         !                  ! ---------------- ! 
+         !                                     !** initialization of the velocity
+         hdiv_139_101(:) = 0.e0                     !  139,101 (Atlantic side, south point)
+         hdiv_139_102(:) = 0.e0                     !  139,102 (Atlantic side, north point)
+         hdiv_141_102(:) = 0.e0                     !  141,102 (Med sea  side)
+            
+         !                                     !** imposed transport
+         zio_flow    = 0.8e6                        ! inflow surface  water
+         zrecirc_mid = 0.7e6                        ! middle recirculation water
+         zrecirc_upp = 2.5e6                        ! upper  recirculation water
+         zrecirc_bot = 3.5e6                        ! bottom recirculation water
+         !
+         DO jj = mj0(101), mj1(101)            !** profile of hdiv at 139,101 (Atlantic side, south point)
+            DO ji = mi0(139), mi1(139)         !-----------------------------
+               DO jk = 1, 14                        ! surface in/out flow (Atl -> Med)   (div >0)
+                  hdiv_139_101(jk) = + zio_flow / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               DO jk = 15, 20                       ! middle  reciculation (Atl 101 -> Atl 102)   (div >0)   
+                  hdiv_139_101(jk) = + zrecirc_mid / ( 6. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               !                                    ! upper reciculation (Atl 101 -> Atl 101)   (div >0)
+               hdiv_139_101(21) =               + zrecirc_upp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               !
+               !                                    ! upper & bottom reciculation (Atl 101 -> Atl 101 & 102)   (div >0)
+               hdiv_139_101(22) = ( zrecirc_bot - zrecirc_upp ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)            !** profile of hdiv at 139,102 (Atlantic side, north point)
+            DO ji = mi0(139), mi1(139)         !-----------------------------
+               DO jk = 15, 20                       ! middle reciculation (Atl 101 -> Atl 102)   (div <0)                
+                  hdiv_139_102(jk) = - zrecirc_mid / ( 6. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               !                                    ! outflow of Mediterranean sea + deep recirculation   (div <0) 
+               hdiv_139_102(22) = - ( zio_flow + zrecirc_bot ) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)            !** velocity profile at 141,102  (Med sea side)
+            DO ji = mi0(141), mi1(141)         !------------------------------
+               DO  jk = 1, 14                       ! surface inflow in the Med     (div <0)
+                  hdiv_141_102(jk) = - zio_flow / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               !                                    ! deep    outflow toward the Atlantic    (div >0) 
+               hdiv_141_102(21)    = + zio_flow / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
+         !                  ! ---------=====-- ! 
+         !                                     !** emp on the Mediterranean Sea  (div >0) 
+         zemp_med = 0.e0                       !-------------------------------
+         DO jj = mj0(96), mj1(110)                  ! sum over the Med sea
+            DO ji = mi0(141),mi1(181)
+               zemp_med = zemp_med + emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj) 
+            END DO
+         END DO
+         DO jj = mj0(96), mj1(96)                   ! minus 2 points in Red Sea 
+            DO ji = mi0(148),mi1(148)
+               zemp_med = zemp_med - emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
+            END DO
+            DO ji = mi0(149),mi1(149)
+               zemp_med = zemp_med - emp(ji,jj) * e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( zemp_med )    ! sum with other processors value
+         zemp_med = zemp_med * 1.e-3                ! convert in m3
+         !
+         !                                     !** Correct hdivn (including emp adjustment)
+         !                                     !-------------------------------------------
+         DO jj = mj0(101), mj1(101)                 !* 139,101 (Atlantic side, south point)
+            DO ji = mi0(139), mi1(139) 
+               hdiv_139_101_kt(:) = hdiv_139_101(:)      
+               DO jk = 1, 14                              ! increase the inflow from the Atlantic   (div >0) 
+                  hdiv_139_101_kt(jk) = hdiv_139_101(jk) + zemp_med / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_139_101_kt(:)
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)                 !* 139,102 (Atlantic side, north point)
+            DO ji = mi0(139), mi1(139) 
+               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_139_102(:)
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)                 !* 141,102 (Med side)
+            DO ji = mi0(141), mi1(141) 
+               hdiv_141_102(:) = hdiv_141_102(:)
+               DO jk = 1, 14                              ! increase the inflow from the Atlantic   (div <0)
+                  hdiv_141_102_kt(jk) = hdiv_141_102(jk) - zemp_med / ( 14. * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               hdivn(ji, jj,:) = hdivn(ji, jj,:) + hdiv_141_102_kt(:)
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
+         !                  ! --------=======- !
+         !
+         DO jj = mj0(101), mj1(101)            !** 139,101 (Atlantic side, south point)   (div >0)
+            DO ji = mi0(139), mi1(139) 
+               DO jk = 1, jpkm1                         ! surf inflow + mid. & bottom reciculation (from Atlantic)	
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_139_101_kt(jk) * tn(ji,jj,jk)
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_139_101_kt(jk) * sn(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         DO jj = mj0(102), mj1(102)            !** 139,102 (Atlantic side, north point)   (div <0)
+            DO ji = mi0(139), mi1(139) 
+               DO jk = 15, 20                            ! middle  reciculation (Atl 101 -> Atl 102)   (div <0)
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_139_102(jk) * tn(ji,jj-1,jk)  ! middle Atlantic recirculation
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_139_102(jk) * sn(ji,jj-1,jk)
+               END DO
+               !                                         ! upper & bottom Atl. reciculation (Atl 101 -> Atl 102) - (div <0)
+               !                                         ! deep Med flow                    (Med 102 -> Atl 102) - (div <0)
+               ta(ji,jj,22) = ta(ji,jj,22) + hdiv_141_102(21) * tn(ji+2,jj  ,21)   &  ! deep Med flow  
+                  &                        + hdiv_139_101(21) * tn(ji  ,jj-1,21)   &  ! upper  Atlantic recirculation  
+                  &                        + hdiv_139_101(22) * tn(ji  ,jj-1,22)      ! bottom Atlantic recirculation  
+               sa(ji,jj,22) = sa(ji,jj,22) + hdiv_141_102(21) * sn(ji+2,jj  ,21)   &
+                  &                        + hdiv_139_101(21) * sn(ji  ,jj-1,21)   &
+                  &                        + hdiv_139_101(22) * sn(ji  ,jj-1,22) 
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)                 !* 141,102 (Med side)   (div <0)
+            DO ji = mi0(141), mi1(141) 
+               DO jk = 1, 14                             ! surface flow from Atlantic to Med sea
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_141_102_kt(jk) * tn(ji-2,jj-1,jk)
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_141_102_kt(jk) * sn(ji-2,jj-1,jk)
+               END DO
+               !                                         ! deeper flow from Med sea to Atlantic
+               ta(ji,jj,21) = ta(ji,jj,21) - hdiv_141_102(21) * tn(ji,jj,21)
+               sa(ji,jj,21) = sa(ji,jj,21) - hdiv_141_102(21) * sn(ji,jj,21)
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
+         !                  ! --------=======- !
+         ! at this stage, (ua,va) are the after velocity, not the tendancy
+         ! compute the velocity from the divergence at T-point
+         !
+         DO jj = mj0(101), mj1(101)            !** 139,101 (Atlantic side, south point)
+            DO ji = mi0(139), mi1(139)                    ! div >0 => ua >0, same sign
+               ua(ji,jj,:) = hdiv_139_101_kt(:) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,:) )   &
+                  &                             * e2u(ji,jj) * fse3u(ji,jj,:)
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)            !** 139,102 (Atlantic side, north point)
+            DO ji = mi0(139), mi1(139)                    ! div <0 => ua <0, same sign
+               ua(ji,jj,:) = hdiv_139_102(:) / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,:) )   &
+                  &                          * e2u(ji,jj) * fse3u(ji,jj,:)   
+            END DO
+         END DO
+         DO jj = mj0(102), mj1(102)            !** 140,102 (Med side) (140 not 141 as it is a U-point)
+            DO ji = mi0(140), mi1(140)                    ! div >0 => ua <0, opposite sign
+               ua(ji,jj,:) = - hdiv_141_102(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
+                  &                            * e2u(ji,jj) * fse3u(ji,jj,:)
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+   END SUBROUTINE cla_gibraltar
+
+
+   SUBROUTINE cla_hormuz( cd_td )
+      !! -------------------------------------------------------------------
+      !!                   ***  ROUTINE div_hormuz  ***
+      !!              
+      !! ** Purpose :   update the now horizontal divergence, the tracer 
+      !!              tendancyand the after velocity in vicinity of Hormuz 
+      !!              strait ( Persian Gulf - Indian ocean ).
+      !!
+      !! ** Method  :   Hormuz strait
+      !!            ______________   
+      !!            |/////|<==      surface inflow
+      !!        94  |/////|     
+      !!            |/////|==>      deep    outflow
+      !!            |_____|_______
+      !!              171    172     
       !!---------------------------------------------------------------------
-      !!               ***  ROUTINE tra_gibraltar  ***
-      !!
-      !! ** Purpose :
-      !!        Update the horizontal advective trend of tracers (t & s)
-      !!        correction in Gibraltar  and
-      !!        add it to the general trend of tracer equations.
-      !!
-      !! ** Method :
-      !!      We impose transport at Gibraltar and knowing T and S in
-      !!      surface and deeper at each side of the strait, we deduce T and S
-      !!      of the outflow of the Mediterranean Sea in the Atlantic ocean .
-      !!                                
-      !!          ________________      N        ________________
-      !! 102           |    |->         |           <-|    |<-
-      !! 101      ___->|____|_____   W - - E     ___->|____|_____
-      !!           139   140  141       |         139   140  141
-      !!          horizontal view       S        horizontal view
-      !!            surface                          depth
-      !!         C A U T I O N : the trend saved is the centered trend only.
-      !!      It doesn't take into account the upstream part of the scheme.
-      !!
-      !! ** history :
-      !!           !  02-06  (A. Bozec) Original code 
-      !!      8.5  !  02-11  (A. Bozec) F90: Free form and module
+      CHARACTER(len=1), INTENT(in) ::   cd_td   ! ='ini' initialisation
+      !!                                        ! ='div' update the divergence
+      !!                                        ! ='tra' update the tracers
+      !!                                        ! ='spg' update after velocity
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zio_flow     ! temporary scalar
       !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::  ji, jj, jk               ! dummy loop indices
-      REAL(wp) :: zsu, zvt
-      REAL(wp) :: zsumt, zsumt1, zsumt2, zsumt3, zsumt4
-      REAL(wp) :: zsums, zsums1, zsums2, zsums3, zsums4
-      REAL(wp) :: zt, zs
-      REAL(wp) :: zwei
-      REAL(wp), DIMENSION (jpk) ::  zu1_ms, zu2_ms, zu3_ms
-      !!---------------------------------------------------------------------
-      
-      ! Initialization of vertical sum for T and S transport
-      ! ----------------------------------------------------
-
-      zsumt  = 0.e0        ! West Gib. surface south point ( T )
-      zsums  = 0.e0        ! West Gib. surface south point ( S )
-      zsumt1 = 0.e0        ! East Gib. surface north point ( T )
-      zsums1 = 0.e0        ! East Gib. surface north point ( S )
-      zsumt2 = 0.e0        ! East Gib. depth   north point ( T )
-      zsums2 = 0.e0        ! East Gib. depth   north point ( S )
-      zsumt3 = 0.e0        ! West Gib. depth   south point ( T ) 
-      zsums3 = 0.e0        ! West Gib. depth   south point ( S ) 
-      zsumt4 = 0.e0        ! West Gib. depth   north point ( T ) 
-      zsums4 = 0.e0        ! West Gib. depth   north point ( S ) 
-      
-      ! EMP of Mediterranean Sea 
-      ! ------------------------
- 
-      zempmed = 0.e0
-      zwei = 0.e0
-      DO jj = mj0(96),mj1(110)
-         DO ji = mi0(141),mi1(181)
-            zwei    = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempmed = zempmed + ( emp(ji,jj) - rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempmed )      ! sum with other processors value
-
-
-      ! minus 2 points in Red Sea and 3 in Atlantic ocean
-      DO jj = mj0(96),mj1(96)
-         DO ji = mi0(148),mi1(148)
-            zempmed = zempmed  -  ( emp(ji  ,jj)-rnf(ji  ,jj) ) * tmask(ji  ,jj,1) * e1t(ji  ,jj) * e2t(ji  ,jj)   & 
-                               -  ( emp(ji+1,jj)-rnf(ji+1,jj) ) * tmask(ji+1,jj,1) * e1t(ji+1,jj) * e2t(ji+1,jj)    
-         END DO
-      END DO
-
-      ! convert in m3
-      zempmed = zempmed * 1.e-3
-
-      ! Velocity profile at each point
-      ! ------------------------------
-
-      zu1_ms(:) = zu1_ms_i(:)
-      zu2_ms(:) = zu2_ms_i(:)
-      zu3_ms(:) = zu3_ms_i(:)
-
-      ! velocity profile at 139,101  South point + (emp-rnf) on surface 
-      DO jk = 1, 14                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms(jk) = zu1_ms(jk) + ( zempmed / 14. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      ! profile at East Gibraltar    
-      ! velocity profile at 141,102  + (emp-rnf) on surface 
-      DO  jk = 1, 14                     
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu3_ms(jk) = zu3_ms(jk) +  ( zempmed / 14. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) ) 
-            END DO
-         END DO
-      END DO
-     
-      ! Balance of temperature and salinity
-      ! -----------------------------------
-
-      ! west gibraltar surface vertical sum of transport* S,T
-      DO  jk =  1, 14 
-         DO jj = mj0(101), mj1(101) 
-            DO ji = mi0(139), mi1(139) 
-               zsumt  = zsumt + tn(ji, jj,jk) * e2u(ji, jj) * fse3u(ji, jj,jk) * zu1_ms(jk)  
-               zsums  = zsums + sn(ji, jj,jk) * e2u(ji, jj) * fse3u(ji, jj,jk) * zu1_ms(jk) 
-            END DO
-         END DO
-      END DO
-
-      ! east Gibraltar surface  vertical sum of transport* S,T
-      DO  jk =  1, 14 
-         DO jj = mj0(101), mj1(101) 
-            DO ji = mi0(139), mi1(139) 
-               zsumt1 = zsumt1 + tn(ji, jj,jk) * e2u(ji+1, jj+1) * fse3u(ji+1, jj+1,jk) * zu3_ms(jk)
-               zsums1 = zsums1 + sn(ji, jj,jk) * e2u(ji+1, jj+1) * fse3u(ji+1, jj+1,jk) * zu3_ms(jk)
-            END DO
-         END DO
-      END DO
-
-      ! east Gibraltar deeper  vertical sum of transport* S,T
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(141), mi1(141) 
-            zsumt2 = tn(ji, jj,21) * e2u(ji-1, jj) * fse3u(ji-1, jj,21) * zu3_ms(21)
-            zsums2 = sn(ji, jj,21) * e2u(ji-1, jj) * fse3u(ji-1, jj,21) * zu3_ms(21)
-         END DO
-      END DO
-      
-      ! west Gibraltar deeper vertical sum of transport* S,T
-      DO  jk =  21, 22 
-         DO jj = mj0(101), mj1(101) 
-            DO ji = mi0(139), mi1(139) 
-               zsumt3 = zsumt3 + tn(ji, jj,jk) * e2u(ji, jj) * fse3u(ji, jj,jk) * zu1_ms(jk)
-               zsums3 = zsums3 + sn(ji, jj,jk) * e2u(ji, jj) * fse3u(ji, jj,jk) * zu1_ms(jk)
-            END DO
-         END DO
-      END DO
-
-      ! Total transport = 0.
-      zsumt4 = zsumt2 + zsumt1 - zsumt - zsumt3
-      zsums4 = zsums2 + zsums1 - zsums - zsums3
-
-      ! Temperature and Salinity at West gibraltar , Level 22
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zt = zsumt4 / ( zu2_ms(22) * e2u(ji-1, jj) * fse3u(ji-1, jj, 22) )
-            zs = zsums4 / ( zu2_ms(22) * e2u(ji-1, jj) * fse3u(ji-1, jj, 22) )
-         END DO
-      END DO
-      
-      ! New Temperature and Salinity trend at West Gibraltar
-      ! ----------------------------------------------------
-
-      ! south point 	
-      DO jk = 1, 22
-         DO jj = mj0(101), mj1(101) 
-            DO ji = mi0(139), mi1(139) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji, jj) * fse3u(ji, jj,jk)
-               ta(ji, jj,jk) = ta(ji, jj,jk) - ( 1. / zvt ) * zsu * zu1_ms(jk) * tn(ji, jj,jk)
-               sa(ji, jj,jk) = sa(ji, jj,jk) - ( 1. / zvt ) * zsu * zu1_ms(jk) * sn(ji, jj,jk)
-            END DO
-         END DO
-      END DO
-
-      ! north point 
-      DO jk = 15, 20
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(139), mi1(139) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji, jj) * fse3u(ji, jj,jk)
-               ta(ji, jj,jk) = ta(ji, jj,jk) - ( 1. / zvt ) * zsu * zu2_ms(jk) * tn(ji, jj-1,jk)
-               sa(ji, jj,jk) = sa(ji, jj,jk) - ( 1. / zvt ) * zsu * zu2_ms(jk) * sn(ji, jj-1,jk)
-            END DO
-         END DO
-      END DO
-
-      ! Gibraltar outflow, north point deeper
-      jk =  22
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(139), mi1(139) 
-            zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-            zsu = e2u(ji, jj) * fse3u(ji, jj,jk)
-            ta(ji, jj,jk) = ta(ji, jj,jk) - ( 1. / zvt ) * zsu * zu2_ms(jk) * zt
-            sa(ji, jj,jk) = sa(ji, jj,jk) - ( 1. / zvt ) * zsu * zu2_ms(jk) * zs
-         END DO
-      END DO
-
-
-      ! New Temperature and Salinity at East Gibraltar 
-      ! ----------------------------------------------
-
-      ! surface	
-      DO jk = 1, 14
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(141), mi1(141) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji-1, jj) * fse3u(ji-1, jj,jk)
-               ta(ji, jj,jk) = ta(ji, jj,jk) + ( 1. / zvt ) * zsu * zu3_ms(jk) * tn(ji-2, jj-1,jk)
-               sa(ji, jj,jk) = sa(ji, jj,jk) + ( 1. / zvt ) * zsu * zu3_ms(jk) * sn(ji-2, jj-1,jk)
-            END DO
-         END DO
-      END DO
-      ! deeper
-      jk =  21
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(141), mi1(141) 
-            zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-            zsu = e2u(ji-1, jj) * fse3u(ji-1, jj,jk)
-            ta(ji, jj,jk) = ta(ji, jj,jk) + ( 1. / zvt ) * zsu * zu3_ms(jk) * tn(ji, jj,jk)
-            sa(ji, jj,jk) = sa(ji, jj,jk) + ( 1. / zvt ) * zsu * zu3_ms(jk) * sn(ji, jj,jk)
-         END DO
-      END DO
-
-   END SUBROUTINE tra_gibraltar
-
-
-   SUBROUTINE tra_hormuz
-      !!---------------------------------------------------------------------
-      !!               ***  ROUTINE tra_hormuz  ***
-      !!
-      !! ** Purpose :   Update the horizontal advective trend of tracers
-      !!        correction in Hormuz.
-      !!
-      !! ** Method :   We impose transport at Hormuz .
-      !!                                
-      !! ** history :
-      !!           !  02-11  (A. Bozec) Original code 
-      !!      8.5  !  02-11  (A. Bozec) F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::  ji, jj, jk              ! dummy loop indices
-      REAL(wp) :: zsu, zvt
-      !!---------------------------------------------------------------------
-      
-      ! New trend at Hormuz strait
-      ! --------------------------
-      DO jk = 1, 8   
-         DO jj = mj0(94), mj1(94) 
+      !
+      SELECT CASE( cd_td ) 
+      !                     ! ---------------- !
+      CASE( 'ini' )         !  initialisation  ! 
+         !                  ! ---------------- ! 
+         !                                     !** profile of horizontal divergence due to cross-land advection
+         zio_flow  = 1.e6                          ! imposed in/out flow
+         !
+         hdiv_172_94(:) = 0.e0         
+         !
+         DO jj = mj0(94), mj1(94)                  ! in/out flow at (i,j) = (172,94)
             DO ji = mi0(172), mi1(172) 
-               zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e2u(ji-1,jj) * fse3u(ji-1,jj,jk)
-               ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu_pg(jk) * tn(ji,jj,jk) 
-               sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu_pg(jk) * sn(ji,jj,jk) 
-            END DO
-         END DO
-      END DO
-      DO jk = 16, 18
-         DO jj = mj0(94), mj1(94) 
+               DO jk = 1, 8                            ! surface inflow  (Indian ocean to Persian Gulf) (div<0)
+                  hdiv_172_94(jk) = - ( zio_flow / 8.e0 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+               DO jk = 16, 18                          ! deep    outflow (Persian Gulf to Indian ocean) (div>0)
+                  hdiv_172_94(jk) = + ( zio_flow / 3.e0 * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         !                                     !** T & S profile in the Hormuz strait (use in deep outflow)
+         !      Temperature       and         Salinity
+         t_171_94_hor(:)  = 0.e0   ;   s_171_94_hor(:)  = 0.e0
+         t_171_94_hor(16) = 18.4   ;   s_171_94_hor(16) = 36.27
+         t_171_94_hor(17) = 17.8   ;   s_171_94_hor(17) = 36.4
+         t_171_94_hor(18) = 16.    ;   s_171_94_hor(18) = 36.27
+         !
+         !                  ! ---------------- !
+      CASE( 'div' )         !   update hdivn   ! (call by divcur module)
+         !                  ! ---------=====-- ! 
+         !                                   
+         DO jj = mj0(94), mj1(94)              !** 172,94 (Indian ocean side)
             DO ji = mi0(172), mi1(172) 
-               zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e2u(ji-1,jj) * fse3u(ji-1,jj,jk)
-               ta(ji,jj,jk) = ta(ji,jj,jk) + ( 1. / zvt ) * zsu * zu_pg(jk) * zthor(jk)
-               sa(ji,jj,jk) = sa(ji,jj,jk) + ( 1. / zvt ) * zsu * zu_pg(jk) * zshor(jk)
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE tra_hormuz
-
-
-   SUBROUTINE tra_cla_init
-      !!---------------------------------------------------------------------
-      !!               ***  ROUTINE tra_cla_init  ***
-      !!
-      !! ** Purpose :   Initialization of variables
-      !!
-      !! ** history :
-      !!      9.0  !  02-11  (A. Bozec) Original code
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER ::  ji, jj, jk              ! dummy loop indices
-      !!---------------------------------------------------------------------
-
-      ! Control print
-      ! -------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'tra_cla_init : cross land advection on tracer '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
-
-      ! Initialization at Bab el Mandeb
-      ! -------------------------------
-
-      ! imposed transport
-      zisw_rs = 0.4e6        ! inflow surface water
-      zurw_rs = 0.2e6        ! upper recirculation water
-!!Alex      zbrw_rs = 1.2e6        ! bottom  recirculation water
-      zbrw_rs = 0.5e6        ! bottom  recirculation water
-
-      ! initialization of the velocity at Bab el Mandeb
-      zu1_rs_i(:) = 0.e0      ! velocity profile at 161,88 South point
-      zu2_rs_i(:) = 0.e0      ! velocity profile at 161,87 North point
-      zu3_rs_i(:) = 0.e0      ! velocity profile at 160,88 East  point
-
-      ! velocity profile at 161,88 East Bab el Mandeb North point 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1, 8                                      
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu1_rs_i(jk) = -( zisw_rs / 8. ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )     
-            END DO
-         END DO
-      END DO
-
-      ! recirculation water 
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zu1_rs_i(20) = -(           zurw_rs ) / ( e2u(ji,jj) * fse3u(ji,jj,20) )
-            zu1_rs_i(21) = -( zbrw_rs - zurw_rs ) / ( e2u(ji,jj) * fse3u(ji,jj,21) )
-         END DO
-      END DO
-      
-      ! velocity profile at 161,87 East Bab el Mandeb South point
-      DO jj = mj0(87), mj1(87) 
-         DO ji = mi0(160), mi1(160) 
-            zu2_rs_i(21) =  ( zbrw_rs + zisw_rs ) / ( e2u(ji,jj) * fse3u(ji,jj,21) )
-         END DO
-      END DO
-
-      ! velocity profile at 161, 88 West Bab el Mandeb 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  10                                     
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu3_rs_i(jk) =  ( zisw_rs / 10. ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )
-            END DO
-         END DO
-      END DO
-
-      ! deeper
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zu3_rs_i(16)  = - zisw_rs /( e1v(ji,jj) * fse3v(ji,jj,16) )
-         END DO
-      END DO
-
-
-      ! Initialization at Gibraltar
-      ! ---------------------------
-
-      ! imposed transport
-      zisw_ms = 0.8e6         ! atlantic-mediterranean  water
-      zmrw_ms = 0.7e6         ! middle recirculation water
-      zurw_ms = 2.5e6         ! upper  recirculation water 
-      zbrw_ms = 3.5e6         ! bottom recirculation water 
-
-      ! initialization of the velocity
-      zu1_ms_i(:) = 0.e0       ! velocity profile at 139,101 South point
-      zu2_ms_i(:) = 0.e0       ! velocity profile at 139,102 North point
-      zu3_ms_i(:) = 0.e0       ! velocity profile at 141,102 East  point
-
-      ! velocity profile at 139,101  South point
-      DO jk = 1, 14                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms_i(jk) = ( zisw_ms / 14. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk)) 
-            END DO
-         END DO
-      END DO
-
-      ! middle recirculation ( uncounting in the balance )
-      DO jk = 15, 20                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms_i(jk) = ( zmrw_ms / 6. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu1_ms_i(21) =  (           zurw_ms ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,21) )
-            zu1_ms_i(22) =  ( zbrw_ms - zurw_ms ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,22) )
-         END DO
-      END DO
-
-      ! velocity profile at 139,102  North point
-      ! middle recirculation ( uncounting in the balance )
-      DO jk = 15, 20                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu2_ms_i(jk) = -( zmrw_ms / 6. ) / ( e2u(ji-1, jj) * fse3u(ji-1, jj,jk) ) 
-            END DO
-         END DO
-      END DO 
-
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu2_ms_i(22) = -( zisw_ms + zbrw_ms ) / ( e2u(ji-1, jj) * fse3u(ji-1, jj,22) )
-         END DO
-      END DO 
-
-      ! profile at East Gibraltar   
-      ! velocity profile at 141,102 
-      DO  jk = 1, 14                     
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu3_ms_i(jk) =  ( zisw_ms / 14. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      ! deeper 
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu3_ms_i(21) = -zisw_ms / ( e2u(ji, jj) * fse3u(ji, jj,21) )
-         END DO
-      END DO
-
-
-      ! Initialization at Hormuz
-      ! ------------------------
-
-      ! imposed transport
-      zisw_pg = 4. * 0.25e6      ! surface and bottom  water
-
-      ! initialization of the velocity
-      zu_pg(:) = 0.e0       ! velocity profile at 139,101 South point
-
-      ! Velocity profile 
-      DO jk = 1, 8 
-         DO jj = mj0(94), mj1(94) 
+               hdivn(ji,jj,:) = hdivn(ji,jj,:) + hdiv_172_94(:)
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'tra' )         !  update (ta,sa)  ! (call by traadv module)
+         !                  ! --------=======- !
+         !                          
+         DO jj = mj0(94), mj1(94)              !** 172,94 (Indian ocean side)
             DO ji = mi0(172), mi1(172) 
-               zu_pg(jk) = -( zisw_pg / 8. ) /  ( e2u(ji-1,jj) * fse3u(ji-1,jj,jk) )
-            END DO
-         END DO
-      END DO
-     DO jk = 16, 18
-         DO jj = mj0(94), mj1(94) 
-            DO ji = mi0(172), mi1(172) 
-               zu_pg(jk) =  ( zisw_pg / 3. )  / ( e2u(ji-1,jj) * fse3u(ji-1,jj,jk) )
-            END DO
-         END DO
-      END DO
-
-      ! Temperature and Salinity at Hormuz
-      zthor(:) = 0.e0
-      zshor(:) = 0.e0
-
-      zthor(16) = 18.4
-      zshor(16) = 36.27
-      !
-      zthor(17) = 17.8
-      zshor(17) = 36.4
-      !
-      zthor(18) = 16.
-      zshor(18) = 36.27
- 
-   END SUBROUTINE tra_cla_init
-
+               DO jk = 1, 8                          ! surface inflow   (Indian ocean to Persian Gulf) (div<0)
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_172_94(jk) * tn(ji,jj,jk) 
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_172_94(jk) * sn(ji,jj,jk) 
+               END DO
+               DO jk = 16, 18                        ! deep outflow     (Persian Gulf to Indian ocean) (div>0)
+                  ta(ji,jj,jk) = ta(ji,jj,jk) - hdiv_172_94(jk) * t_171_94_hor(jk)
+                  sa(ji,jj,jk) = sa(ji,jj,jk) - hdiv_172_94(jk) * s_171_94_hor(jk)
+               END DO
+            END DO
+         END DO
+         !                  ! ---------------- !
+      CASE( 'spg' )         !  update (ua,va)  ! (call by dynspg module)
+         !                  ! --------=======- !
+         ! No barotropic flow through Hormuz strait
+         ! at this stage, (ua,va) are the after velocity, not the tendancy
+         ! compute the velocity from the divergence at T-point
+         DO jj = mj0(94), mj1(94)              !** 171,94 (Indian ocean side) (171 not 172 as it is the western U-point)
+            DO ji = mi0(171), mi1(171)                ! div >0 => ua >0, opposite sign
+               ua(ji,jj,:) = - hdiv_172_94(:) / ( e1t(ji+1,jj) * e2t(ji+1,jj) * fse3t(ji+1,jj,:) )   &
+                  &                           * e2u(ji,jj) * fse3u(ji,jj,:)
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+   END SUBROUTINE cla_hormuz
+   
 #else
    !!----------------------------------------------------------------------
-   !!   Default option                              NO cross land advection
+   !!   Default key                                            Dummy module
    !!----------------------------------------------------------------------
-   USE in_out_manager  ! I/O manager
+   USE in_out_manager ! I/O manager
 CONTAINS
-   SUBROUTINE tra_cla_init 
-   END SUBROUTINE tra_cla_init
-   SUBROUTINE tra_cla( kt ) 
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step indice
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra_cla : No use of cross land advection'
-         WRITE(numout,*) '~~~~~~~'
-      ENDIF
-   END SUBROUTINE tra_cla
+   SUBROUTINE cla_init
+      CALL ctl_stop( 'cla_init: Cross Land Advection hard coded for ORCA_R2 with 31 levels' )
+   END SUBROUTINE cla_init
+   SUBROUTINE cla_div( kt )
+      WRITE(*,*) 'cla_div: You should have not see this print! error?', kt
+   END SUBROUTINE cla_div
+   SUBROUTINE cla_traadv( kt ) 
+      WRITE(*,*) 'cla_traadv: You should have not see this print! error?', kt
+   END SUBROUTINE cla_traadv
+   SUBROUTINE cla_dynspg( kt ) 
+      WRITE(*,*) 'dyn_spg_cla: You should have not see this print! error?', kt
+   END SUBROUTINE cla_dynspg
 #endif
-
+   
    !!======================================================================
 END MODULE cla
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla_div.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla_div.F90	(revision 2391)
+++ 	(revision )
@@ -1,541 +1,0 @@
-MODULE cla_div
-   !!==============================================================================
-   !!                    ***  MODULE  cla_div  ***
-   !! Ocean diagnostic variable : specific update of the horizontal divergence
-   !!                             CAUTION: Specific to ORCA_R2
-   !!==============================================================================
-#if defined key_orca_r2
-   !!----------------------------------------------------------------------
-   !!   'key_orca_r2'                                 global ocean model R2
-   !!----------------------------------------------------------------------
-   !!   div_cla      : 
-   !!   div_bab_el_mandeb
-   !!   div_gibraltar
-   !!   div_hormuz
-   !!   div_cla_init :
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distributed memory computing library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Module variables
-   REAL(wp) :: zempmed, zempred       ! EMP of Mediterranean and Red Sea 
-
-   REAL(wp) :: zisw_rs, zbrw_rs, zurw_rs          ! imposed transport at Red Sea
-   REAL(wp) :: zisw_ms, zbrw_ms, zurw_ms, zmrw_ms ! imposed transport at Mediterranean Sea
-   REAL(wp) :: zisw_pg, zbrw_pg                   ! imposed transport at Persic Gulf
-
-   REAL(wp), DIMENSION (jpk) ::   &
-      zu1_rs_i, zu2_rs_i, zu3_rs_i,                   & ! Red Sea velocities
-      zu1_ms_i, zu2_ms_i, zu3_ms_i,                   & ! Mediterranean Sea velocities
-      zu_pg                                             ! Persic Gulf velocities
-   
-   !! * Routine accessibility
-   PUBLIC div_cla     ! routine called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE div_cla ( kt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE div_cla  ***
-      !!
-      !! ** Purpose :   update the horizontal divergence of the velocity field
-      !!      for at some straits ( Gibraltar, Bab el Mandeb and Hormuz ).
-      !!
-      !! ** Method  :   With imposed transport at each strait, we compute 
-      !!      corresponding velocities and update horizontal divergence.
-      !!        Apply lateral boundary conditions on hdivn through a call
-      !!      to routine lbc_lnk.
-      !!
-      !! ** Action  :   update hdivn array : the now horizontal divergence
-      !!
-      !! History :
-      !!   8.5  !  02-11 (A. Bozec)  Free form, F90
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !!----------------------------------------------------------------------
-
-      ! Correction of the Divergence at some straits
-            
-      IF( kt == nit000 )     CALL div_cla_init    ! Initialization
-
-      CALL div_bab_el_mandeb                        ! New divergence at Bab el Mandeb 
-
-      CALL div_gibraltar                            ! New divergence at Gibraltar
-
-      CALL div_hormuz                               ! Hormuz Strait ( persian Gulf)
-
-      ! Lateral boundary conditions on hdivn
-      CALL lbc_lnk( hdivn, 'T', 1. )
-
-   END SUBROUTINE div_cla
-
-   SUBROUTINE div_bab_el_mandeb
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE div_bab_el_mandeb  ***
-      !!       
-      !! ** Purpose :   Update  the now horizontal divergence of the velocity 
-      !!     field in Bab el Mandeb ( Red Sea strait ).
-      !!
-      !! ** Method  :   Set the velocity field at each side of the strait :
-      !!                                          |
-      !!            |/ \|            N          |\ /|
-      !!            |_|_|______      |          |___|______
-      !!        88  |   |<-       W - - E    88 |   |<-
-      !!        87  |___|______      |       87 |___|->____
-      !!             160 161         S           160 161 
-      !!       horizontal view                horizontal view
-      !!          surface                        depth
-      !!      The now divergence is given by :
-      !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
-      !!
-      !! ** History :
-      !!           !         (A. Bozec) Original code 
-      !!      8.5  !  02-11  (A. Bozec) F90: Free form and module
-      !!----------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  :: ji, jj, jk   ! dummy loop indices
-      REAL(wp) :: zsu, zvt, zwei   ! temporary scalar
-      REAL(wp), DIMENSION (jpk) ::  zu1_rs, zu2_rs, zu3_rs
-      !!---------------------------------------------------------------------
-      
-      ! EMP on the Red Sea 
-      ! ------------------
-
-      zempred = 0.e0
-      zwei = 0.e0
-      DO jj = mj0(87), mj1(96)
-         DO ji = mi0(148), mi1(160) 
-            zwei    = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempred = zempred + ( emp(ji,jj) - rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempred )      ! sum with other processors value
-
-
-      ! convert in m3
-      zempred = zempred * 1.e-3         
-
-      ! Velocity profile at each point
-      ! ------------------------------
-
-      zu1_rs(:) = zu1_rs_i(:)
-      zu2_rs(:) = zu2_rs_i(:)
-      zu3_rs(:) = zu3_rs_i(:)
-
-      ! velocity profile at 161,88 North point 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  8                                      
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu1_rs(jk) = zu1_rs(jk) - ( zempred / 8. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) )
-            END DO
-         END DO
-      END DO
-
-      ! velocity profile at 160,88 North  point
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  10                                     
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu3_rs(jk) = zu3_rs(jk) + ( zempred / 10. ) / ( e1v(ji, jj) * fse3v(ji, jj,jk) )
-            END DO
-         END DO
-      END DO
-       
-      ! Divergence at each point of the straits
-      ! ---------------------------------------
-
-      ! compute the new divergence at 161,88 
-      DO jk = 1, 21
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(161), mi1(161) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji-1, jj) * fse3u(ji-1, jj,jk)
-               hdivn(ji, jj  ,jk) = hdivn(ji, jj  ,jk) - ( 1. / zvt ) * zsu * zu1_rs(jk)
-            END DO
-         END DO
-      END DO
-
-      ! compute the new divergence at 161,87 
-      DO jk = 1, 21
-         DO jj = mj0(87), mj1(87) 
-            DO ji = mi0(161), mi1(161) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji-1, jj) * fse3u(ji-1, jj,jk)
-               hdivn(ji, jj,jk) = hdivn(ji, jj,jk) - ( 1. / zvt ) * zsu * zu2_rs(jk)
-            END DO
-         END DO
-      END DO
-
-      ! compute the divergence at 160,89 
-      DO jk = 1, 18
-         DO jj = mj0(89), mj1(89) 
-            DO ji = mi0(160), mi1(160) 
-               zvt = e1t(ji, jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e1v(ji, jj-1) * fse3v(ji, jj-1,jk)
-               hdivn(ji, jj,jk) = hdivn(ji, jj,jk) - ( 1. / zvt ) * zsu * zu3_rs(jk)
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE div_bab_el_mandeb
-
-   SUBROUTINE div_gibraltar
-      !! -------------------------------------------------------------------
-      !!                 ***  ROUTINE div_gibraltar  ***
-      !!        
-      !! ** Purpose :   update the now horizontal divergence of the velocity 
-      !!     field in Gibraltar.
-      !!
-      !! ** Method :
-      !!          ________________      N        ________________
-      !! 102           |    |->         |           <-|    |<-
-      !! 101      ___->|____|_____   W - - E     ___->|____|_____
-      !!           139   140  141       |         139   140  141
-      !!          horizontal view       S        horizontal view
-      !!            surface                          depth
-      !!      The now divergence is given by :
-      !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
-      !!
-      !! ** History :
-      !!           !         (A. Bozec) Original code 
-      !!      8.5  !  02-10  (A. Bozec) F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  :: ji, jj, jk   ! dummy loop indices
-      REAL(wp) :: zsu, zvt
-      REAL(wp) :: zwei
-      REAL(wp), DIMENSION (jpk) ::  zu1_ms, zu2_ms, zu3_ms
-      !!---------------------------------------------------------------------
-      
-      ! EMP on the Mediterranean Sea 
-      ! ----------------------------
-
-      zempmed = 0.e0
-      zwei = 0.e0
-      DO jj = mj0(96), mj1(110)
-         DO ji = mi0(141),mi1(181)
-            zwei    = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempmed = zempmed + ( emp(ji,jj) - rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempmed )      ! sum with other processors value
-
-      ! minus 2 points in Red Sea and 3 in Atlantic 
-      DO jj = mj0(96), mj1(96)
-         DO ji = mi0(148),mi1(148)
-            zempmed = zempmed -  ( emp(ji  , jj) - rnf(ji  ,jj) ) * tmask(ji  , jj,1) * e1t(ji  , jj) * e2t(ji  , jj)   & 
-                              -  ( emp(ji+1, jj) - rnf(ji+1,jj) ) * tmask(ji+1, jj,1) * e1t(ji+1, jj) * e2t(ji+1, jj)    
-         END DO
-      END DO
-
-      ! convert in m3
-      zempmed = zempmed * 1.e-3
-
-      ! Velocity profile at each point
-      ! ------------------------------
-
-      zu1_ms(:) = zu1_ms_i(:)
-      zu2_ms(:) = zu2_ms_i(:)
-      zu3_ms(:) = zu3_ms_i(:)
-
-      ! velocity profile at 139,101 South point 
-      ! we imposed zisw + EMP above the Mediterranean Sea 
-      DO jk = 1, 14                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms(jk) =  zu1_ms(jk) + ( zempmed / 14. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk) ) 
-            END DO
-         END DO
-      END DO
-      
-      ! velocity profile at 141,102  East point
-      ! flux in surface inflow of the Atlantic ocean + EMP   
-      DO  jk = 1, 14                     
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu3_ms(jk) = zu3_ms(jk) +  ( zempmed / 14. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) ) 
-            END DO
-         END DO
-      END DO
-     
-      ! Divergence at each point of the straits
-      ! ---------------------------------------
-
-      ! compute the new divergence at 139,101 South point  
-      DO jk = 1, jpk
-         DO jj = mj0(101), mj1(101) 
-            DO ji = mi0(139), mi1(139) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji, jj) * fse3u(ji, jj,jk)
-               hdivn(ji, jj,jk) = hdivn(ji, jj,jk) + ( 1. / zvt ) * zsu * zu1_ms(jk) 
-            END DO
-         END DO
-      END DO
-
-      ! compute the new divergence at 139,102 deep North point
-      DO jk = 1, jpk
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(139), mi1(139) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji, jj) * fse3u(ji, jj,jk)
-               hdivn(ji, jj,jk) = hdivn(ji, jj,jk) + ( 1. / zvt ) * zsu * zu2_ms(jk) 
-            END DO
-         END DO
-      END DO
-
-      ! compute the divergence at 141,102 East point
-      DO jk = 1, jpk
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(141), mi1(141) 
-               zvt = e1t(ji, jj) * e2t(ji, jj) * fse3t(ji, jj,jk)
-               zsu = e2u(ji-1, jj) * fse3u(ji-1, jj,jk)
-               hdivn(ji, jj,jk) = hdivn(ji, jj,jk) - ( 1. / zvt ) * zsu * zu3_ms(jk)
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE div_gibraltar
-
-   SUBROUTINE div_hormuz
-      !! -------------------------------------------------------------------
-      !!                   ***  ROUTINE div_hormuz  ***
-      !!              
-      !! ** Purpose :   update the now horizontal divergence of the velocity 
-      !!     field in Hormuz ( Persic Gulf strait ) .
-      !!
-      !! ** Method :
-      !!      The now divergence is given by :
-      !!         hdivn = 1/(e1t*e2t) [ di(e2u  un) + dj(e1v  vn) ]
-      !!
-      !! ** History :
-      !!           !         (A. Bozec) Original code 
-      !!      8.5  !  02-10  (A. Bozec) F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  :: ji, jj, jk   ! dummy loop indices
-      REAL(wp) :: zsu, zvt     ! temporary scalars
-      !!---------------------------------------------------------------------
-
-      ! New divergence at Hormuz 
-      ! ------------------------
-      DO jk = 1, jpk
-         DO jj = mj0(94), mj1(94) 
-            DO ji = mi0(172), mi1(172) 
-               zvt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
-               zsu = e2u(ji-1,jj) * fse3u(ji-1,jj,jk)
-               hdivn(ji,jj,jk) = hdivn(ji,jj,jk) - ( 1. / zvt ) * zsu * zu_pg(jk)
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE div_hormuz
-
-
-   SUBROUTINE div_cla_init
-      !! -------------------------------------------------------------------
-      !!                   ***  ROUTINE div_cla_init  ***
-      !!           
-      !! ** Purpose :   Initialization of variables at all straits  
-      !!
-      !! ** History :
-      !!           !         (A. Bozec) Original code 
-      !!      8.5  !  02-10  (A. Bozec) F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Local declarations
-      INTEGER  :: ji, jj, jk   ! dummy loop indices
-      !!---------------------------------------------------------------------
-
-      ! Control print
-      ! -------------
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'divmod_cross_land : cross land advection on divergence '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
-      IF(lwp) WRITE(numout,*) ' '
-
-      ! Initialization at Bab el Mandeb
-      ! -------------------------------
-
-      ! imposed transport
-      zisw_rs = 0.4e6             ! inflow surface water
-      zurw_rs = 0.2e6             ! upper recirculation water
-!Alex      zbrw = 1.2e6        ! bottom  recirculation water
-      zbrw_rs = 0.5e6             ! bottom  recirculation water
-
-      ! initialization of the velocity
-      zu1_rs_i(:) = 0.e0          ! velocity profile at 161,88 South point
-      zu2_rs_i(:) = 0.e0          ! velocity profile at 161,87 North point
-      zu3_rs_i(:) = 0.e0          ! velocity profile at 160,88 East  point
-
-      ! velocity profile at 161,88 North point 
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  8                                      
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu1_rs_i(jk) = -( zisw_rs / 8. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) )
-            END DO
-         END DO
-      END DO
-
-      ! recirculation water 
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zu1_rs_i(20) =            - zurw_rs   / ( e2u(ji, jj) * fse3u(ji, jj,20) )
-            zu1_rs_i(21) = -( zbrw_rs - zurw_rs ) / ( e2u(ji, jj) * fse3u(ji, jj,21) )
-         END DO
-      END DO
-      
-      ! velocity profile at 161,87 South point
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zu2_rs_i(21) =  ( zbrw_rs + zisw_rs ) / ( e2u(ji, jj-1 ) * fse3u(ji, jj-1,21) )
-         END DO
-      END DO
-
-      ! velocity profile at 160,88 North  point
-      ! we imposed zisw_rs + EMP above the Red Sea 
-      DO jk = 1,  10                                     
-         DO jj = mj0(88), mj1(88) 
-            DO ji = mi0(160), mi1(160) 
-               zu3_rs_i(jk) = ( zisw_rs / 10. ) / ( e1v(ji, jj) * fse3v(ji, jj,jk) )
-            END DO
-         END DO
-      END DO
-
-      ! deeper
-      DO jj = mj0(88), mj1(88) 
-         DO ji = mi0(160), mi1(160) 
-            zu3_rs_i(16)  = - zisw_rs / ( e1v(ji, jj) * fse3v(ji, jj,16) )
-         END DO
-      END DO
-      
-
-      ! Initialization at Gibraltar
-      ! ---------------------------
-
-      ! imposed transport
-      zisw_ms = 0.8e6             ! inflow surface  water
-      zmrw_ms = 0.7e6             ! middle recirculation water
-      zurw_ms = 2.5e6             ! upper  recirculation water
-      zbrw_ms = 3.5e6             ! bottom recirculation water
-
-      ! initialization of the velocity
-      zu1_ms_i(:) = 0.e0          ! velocity profile at 139,101 South point
-      zu2_ms_i(:) = 0.e0          ! velocity profile at 139,102 North point
-      zu3_ms_i(:) = 0.e0          ! velocity profile at 141,102 East  point
-
-      ! velocity profile at 139,101 South point 
-      DO jk = 1, 14                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms_i(jk) =  ( zisw_ms / 14. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      ! recirculation water 
-      DO jk = 15, 20                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu1_ms_i(jk) =  ( zmrw_ms / 6. ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu1_ms_i(21) =  (           zurw_ms ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,21) )
-            zu1_ms_i(22) =  ( zbrw_ms - zurw_ms ) / ( e2u(ji-1, jj-1) * fse3u(ji-1, jj-1,22) )
-         END DO
-      END DO
-      
-      ! velocity profile at 139,102 North point
-      DO jk = 15, 20                      
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu2_ms_i(jk) = -( zmrw_ms / 6. ) / ( e2u(ji-1, jj) * fse3u(ji-1, jj,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      ! outflow of Mediterranean sea + recirculation
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu2_ms_i(22) = -( zisw_ms + zbrw_ms ) / ( e2u(ji-1, jj) * fse3u(ji-1, jj,22) )
-         END DO
-      END DO
-      
-      ! velocity profile at 141,102  East point
-      ! flux in surface inflow of the Atlantic ocean    
-      DO  jk = 1, 14                     
-         DO jj = mj0(102), mj1(102) 
-            DO ji = mi0(140), mi1(140) 
-               zu3_ms_i(jk) =  ( zisw_ms / 14. ) / ( e2u(ji, jj) * fse3u(ji, jj,jk) ) 
-            END DO
-         END DO
-      END DO
-
-      ! deeper 
-      DO jj = mj0(102), mj1(102) 
-         DO ji = mi0(140), mi1(140) 
-            zu3_ms_i(21) = - zisw_ms / ( e2u(ji, jj) * fse3u(ji, jj,21) )
-         END DO
-      END DO
-
-      ! Initialization at Hormuz
-      ! ------------------------
-
-      ! imposed transport
-      zisw_pg = 4. * 0.25e6       ! inflow surface  water
-      zbrw_pg = 4. * 0.25e6       ! bottom recirculation water
-
-      ! initialization of the velocity
-      zu_pg(:) = 0.e0           ! velocity profile at 172,94 
-
-      ! velocity profile 
-      DO jk = 1, 8 
-         DO jj = mj0(94), mj1(94) 
-            DO ji = mi0(172), mi1(172) 
-               zu_pg(jk) = -( zisw_pg / 8. ) / ( e2u(ji-1,jj) * fse3u(ji-1,jj, jk) )
-            END DO
-         END DO
-      END DO
-
-      DO jk = 16, 18
-         DO jj = mj0(94), mj1(94) 
-            DO ji = mi0(172), mi1(172) 
-               zu_pg(jk) =  ( zbrw_pg / 3. ) / ( e2u(ji-1,jj) * fse3u(ji-1,jj, jk) )
-            END DO
-         END DO
-      END DO
-
-   END SUBROUTINE div_cla_init
-#else
-   !!----------------------------------------------------------------------
-   !!   Default key                                            Dummy module
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE div_cla( kt )
-      WRITE(*,*) 'div_cla: You should have not see this print! error?', kt
-   END SUBROUTINE div_cla
-#endif
-   
-   !!======================================================================
-END MODULE cla_div
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla_dynspg.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/cla_dynspg.F90	(revision 2391)
+++ 	(revision )
@@ -1,241 +1,0 @@
-MODULE cla_dynspg
-   !!======================================================================
-   !!                       ***  cla_dynspg  ***
-   !!======================================================================
-   !!   dyn_spg      : update the momentum trend with the surface pressure
-   !!                  gradient in the free surface constant volume case
-   !!                  with vector optimization
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain 
-   USE zdf_oce         ! ocean vertical physics
-   USE obc_oce         ! Lateral open boundary condition
-   USE sol_oce         ! solver variables
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE phycst          ! physical constants
-   USE solpcg          ! preconditionned conjugate gradient solver
-   USE solsor          ! Successive Over-relaxation solver
-   USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
-   USE obcvol          ! ocean open boundary condition (obc_vol routines)
-   USE in_out_manager  ! I/O manager
-   USE lib_mpp         ! distribued memory computing
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Accessibility
-   PUBLIC dyn_spg_cla   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE dyn_spg_cla( kt ) 
-      !!----------------------------------------------------------------------
-      !!              ***  routine dyn_spg_cross_land  ***
-      !!
-      !! ** Purpose :
-      !!
-      !! ** Method :
-      !!
-      !! ** Action :
-      !!
-      !! History :
-      !!        !         (A. Bozec)  Original code
-      !!   8.5  !  02-11  (A. Bozec)  F90: Free form and module
-      !!---------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt           ! ocean time-step
-      !! * Local declarations
-      INTEGER  ::   ji, jj, jk                ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1        ! temporary integer
-      REAL(wp) ::    &    
-         zempmed, zempred,   &                ! EMP on Med Sea ans Red Sea
-         zwei,   &                            !              
-         zisw_rs, zurw_rs, zbrw_rs,      &    ! imposed transport Red sea
-         zisw_ms, zurw_ms, zbrw_ms, zmrw_ms   ! imposed transport Med Sea
-      !!----------------------------------------------------------------------
-
-      ! Different velocities for straits ( Gibraltar, Bab el Mandeb...)
-         
-      ! Control print
-      ! -------------
-      IF( kt == nit000 ) THEN 
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dynspg_cross_land : cross land advection on surface '
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~   pressure '
-         IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      ! EMP on Mediterranean Sea and Red Sea 
-      ! ------------------------------------
-      ! compute the (emp-rnf) in Mediterranean Sea
-      zempmed = 0.e0
-      zwei = 0.e0
-      ij0 =  96   ;   ij1 = 110
-      ii0 = 141   ;   ii1 = 181
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            zwei    = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempmed = zempmed + ( emp(ji,jj)-rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempmed )      ! sum with other processors value
-
-      ! minus 2 points in Red Sea and 3 in Atlantic 
-      ij0 =  96   ;   ij1 =  96
-      ii0 = 148   ;   ii1 = 148
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            zempmed = zempmed - ( emp(ji  ,jj) - rnf(ji  ,jj) ) * tmask(ji  ,jj,1) * e1t(ji  ,jj) * e2t(ji  ,jj)   & 
-               &              - ( emp(ji+1,jj) - rnf(ji+1,jj) ) * tmask(ji+1,jj,1) * e1t(ji+1,jj) * e2t(ji+1,jj)    
-         END DO
-      END DO
-      ! we convert in m3
-      zempmed = zempmed * 1.e-3
-
-      ! compute the (emp-rnf) in Red Sea   
-      zempred = 0.e0
-      zwei = 0.e0
-      ij0 =  87   ;   ij1 =  96
-      ii0 = 148   ;   ii1 = 160
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            zwei      = tmask(ji,jj,1) * e1t(ji,jj) * e2t(ji,jj)
-            zempred   = zempred + ( emp(ji,jj) - rnf(ji,jj) ) * zwei 
-         END DO
-      END DO
-      IF( lk_mpp )   CALL mpp_sum( zempred )      ! sum with other processors value
-
-      ! we convert in m3
-      zempred = zempred * 1.e-3
-
-      ! New Transport at Bab el Mandeb and Gibraltar
-      ! --------------------------------------------
-
-      ! imposed transport at Bab el Mandeb
-      zisw_rs = 0.4e6        ! inflow surface water
-      zurw_rs = 0.2e6        ! upper recirculation water
-!!Alex      zbrw_rs = 1.2e6        ! bottom  recirculation water
-      zbrw_rs = 0.5e6        ! bottom  recirculation water
-
-      ! imposed transport at Gibraltar
-      zisw_ms  = 0.8e6          ! atlantic-mediterranean  water
-      zmrw_ms  = 0.7e6          ! middle recirculation water
-      zurw_ms  = 2.5e6          ! upper  recirculation water 
-      zbrw_ms  = 3.5e6          ! bottom recirculation water 
-
-      ! Different velocities for straits ( Gibraltar, Bab el Mandeb )
-      ! -------------------------------------------------------------
-
-      ! Bab el Mandeb
-      ! -------------
-      ! 160,88 north point Bab el Mandeb
-      ij0 =  88   ;   ij1 =  88
-      ii0 = 160   ;   ii1 = 160
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj  ,: ) = 0.e0  !  North East Bab el Mandeb 
-         END DO
-      END DO
-      !                              ! surface
-      DO jk = 1,  8                                      
-         DO jj = mj0(ij0), mj1(ij1)
-            DO ji = mi0(ii0),mi1(ii1)
-               ua(ji, jj,jk) = -( ( zisw_rs + zempred ) / 8. ) / ( e2u(ji, jj) * fse3t(ji, jj,jk) )     
-            END DO
-         END DO
-      END DO
-      !                              ! deeper
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji, jj,21) = - zbrw_rs / ( e2u(ji, jj) * fse3t(ji, jj,21) )
-         END DO
-      END DO
-
-      ! 160,87 south point Bab el Mandeb
-      ij0 =  87   ;   ij1 =  87
-      ii0 = 160   ;   ii1 = 160
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj  ,: ) = 0.e0  !  South East Bab el Mandeb 
-         END DO
-      END DO
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji, jj,21) =  ( zisw_rs + zbrw_rs ) / ( e2u(ji,jj )*fse3t(ji, jj,21) )      
-         END DO
-      END DO
-
-      ! Gibraltar
-      ! ---------
-
-      ! initialisation of velocity at concerned points 
-      ! 139, 101 south point in Gibraltar 
-      ij0 = 101   ;   ij1 = 101
-      ii0 = 139   ;   ii1 = 139
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj  ,: ) = 0.e0  !  South West Gibraltar
-            ua(ji,jj+1,: ) = 0.e0  !  North West Gibraltar
-         END DO
-      END DO
-      !                            ! surface
-      DO jk = 1, 14                      
-         DO jj = mj0(ij0), mj1(ij1)
-            DO ji = mi0(ii0),mi1(ii1)
-               ua(ji,jj,jk) =  ( ( zisw_ms + zempmed ) / 14. ) / ( e2u(ji,jj) * fse3t(ji,jj,jk) ) 
-            END DO
-         END DO
-      END DO
-      !                            ! middle circulation
-      DO jk = 15, 20                      
-         DO jj = mj0(ij0), mj1(ij1)
-            DO ji = mi0(ii0),mi1(ii1)
-               ua(ji,jj,jk) =  ( zmrw_ms / 6. ) / ( e2u(ji,jj) * fse3t(ji,jj,jk) ) 
-            END DO
-         END DO
-      END DO
-      !                            ! deeper 
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj,21) =             zurw_ms   / ( e2u(ji,jj) * fse3t(ji,jj,21) )
-            ua(ji,jj,22) = ( zbrw_ms - zurw_ms ) / ( e2u(ji,jj) * fse3t(ji,jj,22) )
-         END DO
-      END DO
-
-      ! 139,102 north point in Gibraltar
-      ij0 = 102   ;   ij1 = 102
-      ii0 = 139   ;   ii1 = 139
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj  ,: ) = 0.e0  !  North West Gibraltar
-         END DO
-      END DO
-      DO jk = 15, 20                      
-         DO jj = mj0(ij0), mj1(ij1)
-            DO ji = mi0(ii0),mi1(ii1)
-               ua(ji,jj,jk) = -( zmrw_ms / 6. ) / ( e2u(ji,jj) * fse3t(ji,jj,jk) ) 
-            END DO
-         END DO
-      END DO
-      !                            ! deeper
-      DO jj = mj0(ij0), mj1(ij1)
-         DO ji = mi0(ii0),mi1(ii1)
-            ua(ji,jj,22) = -( zisw_ms + zbrw_ms ) / ( e2u(ji,jj) * fse3t(ji,jj,22) )
-         END DO
-      END DO
-
-   END SUBROUTINE dyn_spg_cla
-
-   !!======================================================================
-END MODULE cla_dynspg
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/opa.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/opa.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/opa.F90	(revision 2392)
@@ -38,4 +38,5 @@
    USE step_oce        ! module used in the ocean time stepping module
    USE sbc_oce         ! surface boundary condition: ocean
+   USE cla             ! cross land advection               (tra_cla routine)
    USE domcfg          ! domain configuration               (dom_cfg routine)
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
@@ -46,5 +47,5 @@
    USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
    USE ldftra          ! lateral diffusivity setting    (ldftra_init routine)
-   USE zdfini
+   USE zdfini          ! vertical physics setting          (zdf_init routine)
    USE phycst          ! physical constant                  (par_cst routine)
    USE trdmod          ! momentum/tracers trends       (trd_mod_init routine)
@@ -67,5 +68,4 @@
    USE trcini          ! passive tracer initialisation
 #endif
-
    USE lib_mpp         ! distributed memory computing
 #if defined key_iomput
@@ -82,5 +82,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id $
+   !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -277,5 +277,5 @@
                             CALL ldf_tra_init      ! Lateral ocean tracer physics
                             CALL ldf_dyn_init      ! Lateral ocean momentum physics
-      IF( lk_ldfslp )       CALL ldf_slp_init      ! slope of lateral mixing
+      IF( lk_ldfslp     )   CALL ldf_slp_init      ! slope of lateral mixing
 
       !                                     ! Active tracers
@@ -285,5 +285,4 @@
       IF( lk_tradmp     )   CALL tra_dmp_init   ! internal damping trends
                             CALL tra_adv_init   ! horizontal & vertical advection
-      IF( n_cla == 1    )   CALL tra_cla_init   ! Cross Land Advection (Update Hor. advection)
                             CALL tra_ldf_init   ! lateral mixing
                             CALL tra_zdf_init   ! vertical mixing and after tracer fields
@@ -296,4 +295,8 @@
                             CALL dyn_zdf_init   ! vertical diffusion
                             CALL dyn_spg_init   ! surface pressure gradient
+                            
+      !                                     ! Misc. options
+      IF( nn_cla == 1   )   CALL cla_init       ! Cross Land Advection
+      
 #if defined key_top
       !                                     ! Passive tracers
@@ -302,15 +305,15 @@
       !                                     ! Diagnostics
                             CALL     iom_init   ! iom_put initialization
-      IF( lk_floats    )    CALL     flo_init   ! drifting Floats
-      IF( lk_diaar5    )    CALL dia_ar5_init   ! ar5 diag
+      IF( lk_floats     )   CALL     flo_init   ! drifting Floats
+      IF( lk_diaar5     )   CALL dia_ar5_init   ! ar5 diag
                             CALL dia_ptr_init   ! Poleward TRansports initialization
                             CALL dia_hsb_init   ! heat content, salt content and volume budgets
                             CALL trd_mod_init   ! Mixed-layer/Vorticity/Integral constraints trends
-      IF( lk_diaobs    ) THEN                   ! Observation & model comparison
+      IF( lk_diaobs     ) THEN                  ! Observation & model comparison
                             CALL dia_obs_init            ! Initialize observational data
                             CALL dia_obs( nit000 - 1 )   ! Observation operator for restart
       ENDIF      
       !                                     ! Assimilation increments
-      IF( lk_asminc    )    CALL asm_inc_init   ! Initialize assimilation increments
+      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
       !
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2392)
@@ -186,5 +186,4 @@
       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
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2391)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 2392)
@@ -31,5 +31,4 @@
    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)
