Index: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5598)
+++ branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5955)
@@ -2,40 +2,44 @@
    !!==============================================================================
    !!                       ***  MODULE  traldf_lap  ***
-   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
+   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
    !!==============================================================================
-   !! History :  OPA  !  87-06  (P. Andrich, D. L Hostis)  Original code
-   !!                 !  91-11  (G. Madec)
-   !!                 !  95-11  (G. Madec)  suppress volumetric scale factors
-   !!                 !  96-01  (G. Madec)  statement function for e3
-   !!            NEMO !  02-06  (G. Madec)  F90: Free form and module
-   !!            1.0  !  04-08  (C. Talandier) New trends organization
-   !!                 !  05-11  (G. Madec)  add zps case
-   !!            3.0  !  10-06  (C. Ethe, G. Madec) Merge TRA-TRC
+   !! 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 horizontal diffusion
-   !!                 using a iso-level harmonic (laplacien) operator.
+   !!   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_oce      ! ocean active tracers: lateral physics
-   USE in_out_manager  ! I/O manager
+   USE ldftra          ! lateral physics: eddy diffusivity
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-   USE lib_mpp         ! MPP library
+   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
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_ldf_lap   ! routine called by step.F90
+   PUBLIC   tra_ldf_lap   ! routine called by traldf.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -43,7 +47,7 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu , pgv ,    &
-      &                                        pgui, pgvi,    &
-      &                                ptb, pta, kjpt ) 
+   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                    pgui, pgvi,   &
+      &                                        ptb , pta , kjpt, kpass ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
@@ -55,6 +59,6 @@
       !!      fields (forward time scheme). The horizontal diffusive trends of 
       !!      the tracer is given by:
-      !!          difft = 1/(e1t*e2t*e3t) {  di-1[ aht e2u*e3u/e1u di(tb) ]
-      !!                                   + dj-1[ aht e1v*e3v/e2v dj(tb) ] }
+      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
+      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
       !!      Add this trend to the general tracer trend pta :
       !!          pta = pta + difft
@@ -63,106 +67,100 @@
       !!                harmonic mixing trend.
       !!----------------------------------------------------------------------
-      USE oce, ONLY:   ztu => ua , ztv => va  ! (ua,va) used as workspace
-      !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000          ! first 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   ) ::   kpass      ! =1/2 first or second passage
+      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,    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
-      INTEGER  ::   iku, ikv, ierr       ! local integers
-      REAL(wp) ::   zabe1, zabe2, zbtr   ! local scalars
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zsign            ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztu, ztv, zaheeu, zaheev
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('tra_ldf_lap')
+      IF( nn_timing == 1 )   CALL timing_start('tra_ldf_lap')
       !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+      IF( kt == nit000 .AND. lwp )  THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
+         WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
-
-      !                                                          ! =========== !
-      DO jn = 1, kjpt                                            ! tracer loop !
-         !                                                       ! =========== !    
-         DO jk = 1, jpkm1                                            ! slab loop
-            !                                           
-            ! 1. First derivative (gradient)
-            ! -------------------
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
+      !
+      !                                !==  Initialization of metric arrays used for all tracers  ==!
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
+      ELSE                    ;   zsign = -1._wp
+      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked!
+               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      !                             ! =========== !
+      DO jn = 1, kjpt               ! tracer loop !
+         !                          ! =========== !    
+         !                               
+         DO jk = 1, jpkm1              !== First derivative (gradient)  ==!
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk)
-                  ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+               DO ji = 1, fs_jpim1
+                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
+                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
                END DO
             END DO
-            IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell
+         END DO  
+         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
+            DO jj = 1, jpjm1                    ! bottom
+               DO ji = 1, fs_jpim1
+                  ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
+                  ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
+               END DO
+            END DO  
+            IF( ln_isfcav ) THEN                ! top in ocean cavities only
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     ! last level
-                     iku = mbku(ji,jj)
-                     ikv = mbkv(ji,jj)
-                     IF( iku == jk ) THEN
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku)
-                        ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn)
-                     ENDIF
-                     IF( ikv == jk ) THEN
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv)
-                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn)
-                     ENDIF
+                     IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
+                     IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
                   END DO
                END DO
             ENDIF
-            ! (ISH)
-            IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell
-                                                   ! into a cavity
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     ! ice shelf level level MAX(2,jk) => only where ice shelf
-                     iku = miku(ji,jj) 
-                     ikv = mikv(ji,jj) 
-                     IF( iku == MAX(2,jk) ) THEN 
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku) 
-                        ztu(ji,jj,jk) = zabe1 * pgui(ji,jj,jn) 
-                     ENDIF 
-                     IF( ikv == MAX(2,jk) ) THEN 
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 
-                        ztv(ji,jj,jk) = zabe2 * pgvi(ji,jj,jn) 
-                     END IF 
-                  END DO
-               END DO
-            ENDIF
-         
-         
-            ! 2. Second derivative (divergence) added to the general tracer trends
-            ! ---------------------------------------------------------------------
+         ENDIF
+         !
+         DO jk = 1, jpkm1              !== Second derivative (divergence) added to the general tracer trends  ==!
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ! horizontal diffusive trends added to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
-                     &                                          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
+               DO ji = fs_2, fs_jpim1
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
+                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
+                     &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                END DO
             END DO
-            !
-         END DO                                             !  End of slab  
+         END DO  
          !
-         ! "Poleward" diffusive heat or salt transports
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-            IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( ztv(:,:,:) )
-            IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( ztv(:,:,:) )
+         !                             !== "Poleward" diffusive heat or salt transports  ==!
+         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
+             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
+            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+               IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( -ztv(:,:,:) )
+               IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( -ztv(:,:,:) )
+            ENDIF
          ENDIF
-         !                                                  ! ==================
-      END DO                                                ! end of tracer loop
-      !                                                     ! ==================
-      IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_lap')
+         !                          ! ==================
+      END DO                        ! end of tracer loop
+      !                             ! ==================
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_lap')
       !
    END SUBROUTINE tra_ldf_lap
-
+   
    !!==============================================================================
 END MODULE traldf_lap
