Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10740)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10806)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
+   SUBROUTINE dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs, kpass )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_lap  ***
@@ -45,10 +45,11 @@
       !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
       !!
-      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
+      !! ** Action : - pu_rhs, pva_rhs increased by the harmonic operator applied on pu, pv.
       !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   kt              ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
       INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! velocity trend   [m/s2]
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -76,12 +77,12 @@
 !!gm open question here : e3f  at before or now ?    probably now...
 !!gm note that ahmf has already been multiplied by fmask
-               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
-                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
-                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
+               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
+                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
                !                                      ! ahm * div        (computed from 2 to jpi/jpj)
 !!gm note that ahmt has already been multiplied by tmask
-               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
-                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
-                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
+               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev1)                                         &
+                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,ktlev1) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,ktlev1) * pu(ji-1,jj,jk)  &
+                  &        + e1v(ji,jj)*e3v(ji,jj,jk,ktlev1) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,ktlev1) * pv(ji,jj-1,jk)  )
             END DO  
          END DO  
@@ -89,10 +90,10 @@
          DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
-                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)   &
                   &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
                   !
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
-                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
+               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,ktlev2)   &
                   &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
             END DO
@@ -105,5 +106,5 @@
 
 
-   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
+   SUBROUTINE dyn_ldf_blp( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dyn_ldf_blp  ***
@@ -116,9 +117,10 @@
       !!      It is computed by two successive calls to dyn_ldf_lap routine
       !!
-      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
+      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! momentum trend
       !
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
@@ -134,9 +136,9 @@
       zvlap(:,:,:) = 0._wp
       !
-      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
+      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap)
       !
       CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
       !
-      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
+      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, zulap, zvlap, pu_rhs, pva_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt_rhs)
       !
    END SUBROUTINE dyn_ldf_blp
