Index: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90
===================================================================
--- branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5845)
+++ branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90	(revision 5883)
@@ -1,39 +1,42 @@
-MODULE traldf_lap
+MODULE traldf_lap_blp
    !!==============================================================================
-   !!                       ***  MODULE  traldf_lap  ***
+   !!                       ***  MODULE  traldf_lap_blp  ***
    !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
    !!==============================================================================
-   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
-   !!                 ! 1991-11  (G. Madec)
-   !!                 ! 1995-11  (G. Madec)  suppress volumetric scale factors
-   !!                 ! 1996-01  (G. Madec)  statement function for e3
-   !!            NEMO ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            1.0  ! 2004-08  (C. Talandier) New trends organization
-   !!                 ! 2005-11  (G. Madec)  add zps case
-   !!            3.0  ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!            3.7  ! 2014-01  (G. Madec, S. Masson) re-entrant laplacian 
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   tra_ldf_lap : update the tracer trend with the lateral diffusion : iso-level laplacian operator
-   !!   tra_ldf_blp : update the tracer trend with the lateral diffusion : iso-level bilaplacian operator
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE ldftra          ! lateral physics: eddy diffusivity
-   USE diaptr          ! poleward transport diagnostics
-   USE trc_oce         ! share passive tracers/Ocean variables
-   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
+   !! History :  3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_ldf_lap   : tracer trend update with iso-level laplacian diffusive operator
+   !!   tra_ldf_blp   : tracer trend update with iso-level or iso-neutral bilaplacian operator
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE ldftra         ! lateral physics: eddy diffusivity
+   USE traldf_iso     ! iso-neutral lateral diffusion (standard operator)     (tra_ldf_iso   routine)
+   USE traldf_triad   ! iso-neutral lateral diffusion (triad    operator)     (tra_ldf_triad routine)
+   USE diaptr         ! poleward transport diagnostics
+   USE trc_oce        ! share passive tracers/Ocean variables
+   USE zpshde         ! partial step: hor. derivative     (zps_hde routine)
    !
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp         ! distribued memory computing library
-   USE timing          ! Timing
-   USE wrk_nemo        ! Memory allocation
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! distribued memory computing library
+   USE timing         ! Timing
+   USE wrk_nemo       ! Memory allocation
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_ldf_lap   ! routine called by traldf.F90
+   PUBLIC   tra_ldf_lap   ! called by traldf.F90
+   PUBLIC   tra_ldf_blp   ! called by traldf.F90
+
+   !                      ! Flag to control the type of lateral diffusive operator
+   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10   ! error in specification of lateral diffusion
+   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00   ! without operator (i.e. no lateral diffusive trend)
+   !                          !!      laplacian     !    bilaplacian    !
+   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  ! iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11   ,   np_blp_i  = 21  ! standard iso-neutral or geopotential operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_it = 12   ,   np_blp_it = 22  ! triad    iso-neutral or geopotential operator
 
    !! * Substitutions
@@ -162,4 +165,89 @@
    END SUBROUTINE tra_ldf_lap
    
+
+   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                    pgui, pgvi,   &
+      &                                                    ptb , pta , kjpt, kldf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tra_ldf_blp  ***
+      !!                    
+      !! ** Purpose :   Compute the before lateral tracer diffusive 
+      !!      trend and add it to the general trend of tracer equation.
+      !!
+      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
+      !!      operator applied to before field (forward in time).
+      !!      It is computed by two successive calls to laplacian routine
+      !!
+      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !!----------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000     ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
+      INTEGER                              , INTENT(in   ) ::   kldf       ! type of operator used
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      !
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zlap         ! laplacian at t-point
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_blp')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,kjpt,   zlap ) 
+      CALL wrk_alloc( jpi,jpj,    kjpt,   zglu, zglv, zgui, zgvi ) 
+      !
+      IF( kt == kit000 .AND. lwp )  THEN
+         WRITE(numout,*)
+         SELECT CASE ( kldf )
+         CASE ( np_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
+         CASE ( np_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
+         CASE ( np_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
+         END SELECT
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+
+      zlap(:,:,:,:) = 0._wp
+      !
+      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
+      !
+      CASE ( np_blp    )               ! iso-level bilaplacian
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
+      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+      END SELECT
+      !
+      DO jn = 1, kjpt
+         CALL lbc_lnk( zlap(:,:,:,jn) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
+      END DO
+      !                                               ! Partial top/bottom cell: GRADh( zlap )  
+      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
+      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, kjpt, zlap, zglu, zglv )              ! only bottom 
+      ENDIF
+      !
+      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
+      !
+      CASE ( np_blp    )               ! iso-level bilaplacian
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
+      CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+      CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,kjpt,   zlap ) 
+      CALL wrk_dealloc( jpi,jpj    ,kjpt,   zglu, zglv, zgui, zgvi ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_blp')
+      !
+   END SUBROUTINE tra_ldf_blp
+
    !!==============================================================================
-END MODULE traldf_lap
+END MODULE traldf_lap_blp
