Index: NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90	(revision 14789)
+++ NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90	(revision 14852)
@@ -266,11 +266,13 @@
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::  zhpi, zhpj
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::  zhpi, zhpj
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
+         ENDIF
       ENDIF
       !
@@ -318,13 +320,15 @@
       INTEGER  ::   iku, ikv                         ! temporary integers
       REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj
-      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zgtsu, zgtsv
-      REAL(wp), DIMENSION(jpi,jpj)     :: zgru, zgrv
-      !!----------------------------------------------------------------------
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj
+      REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zgtsu, zgtsv
+      REAL(wp), DIMENSION(A2D(nn_hls)     ) :: zgru, zgrv
+      !!----------------------------------------------------------------------
+      !
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
+         ENDIF
       ENDIF
 
@@ -410,14 +414,16 @@
       REAL(wp) ::   zcoef0, zuap, zvap, ztmp       ! local scalars
       LOGICAL  ::   ll_tmp1, ll_tmp2               ! local logical variables
-      REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk)  ::   zhpi, zhpj
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
       !
-      IF( ln_wd_il ) ALLOCATE(zcpx(jpi,jpj), zcpy(jpi,jpj))
-      !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OCE original scheme used'
+      IF( ln_wd_il ) ALLOCATE(zcpx(A2D(nn_hls)), zcpy(A2D(nn_hls)))
+      !
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OCE original scheme used'
+         ENDIF
       ENDIF
       !
@@ -462,5 +468,4 @@
           END IF
         END_2D
-        CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
       !
@@ -548,7 +553,7 @@
       REAL(wp) ::   ze3w, ze3wi1, ze3wj1   ! local scalars
       REAL(wp) ::   zcoef0, zuap, zvap     !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
-      REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
-      REAL(wp), DIMENSION(jpi,jpj)      ::  zrhdtop_oce
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk ) ::  zhpi, zhpj
+      REAL(wp), DIMENSION(A2D(nn_hls),jpts) ::  zts_top
+      REAL(wp), DIMENSION(A2D(nn_hls))      ::  zrhdtop_oce
       !!----------------------------------------------------------------------
       !
@@ -560,11 +565,9 @@
       ! compute rhd at the ice/oce interface (ocean side)
       ! usefull to reduce residual current in the test case ISOMIP with no melting
-      DO ji = 1, jpi
-        DO jj = 1, jpj
-          ikt = mikt(ji,jj)
-          zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
-          zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
-        END DO
-      END DO
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         ikt = mikt(ji,jj)
+         zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
+         zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
+      END_2D
       CALL eos( zts_top, risfdep, zrhdtop_oce )
 
@@ -636,21 +639,21 @@
       INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points 
       REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
-      REAL(wp) ::   z_grav_10, z1_12
+      REAL(wp) ::   z_grav_10, z1_12, z1_cff
       REAL(wp) ::   cffu, cffx          !    "         "
       REAL(wp) ::   cffv, cffy          !    "         "
       LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
- 
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz')
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz')
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho')
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho')
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays 
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zhpi, zhpj
+
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz')
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz')
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho')
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho')
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
       !
       IF( ln_wd_il ) THEN
-         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
+         ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) )
         DO_2D( 0, 0, 0, 0 )
           ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
@@ -689,11 +692,12 @@
           END IF
         END_2D
-        CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
 
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
+         ENDIF
       ENDIF
 
@@ -723,9 +727,8 @@
       zdz_k  (:,:,:) = 0._wp
 
-      DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 
-         cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1)
-         IF( cffw > zep) THEN
-            zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) )
-         ENDIF
+      DO_3D( 1, 1, 1, 1, 2, jpk-2 )
+         cffw = MAX( 2._wp * zdrhoz(ji,jj,jk) * zdrhoz(ji,jj,jk+1), 0._wp )
+         z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1)
+         zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
          zdz_k(ji,jj,jk) = 2._wp *   zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1)   &
             &                  / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) )
@@ -737,6 +740,8 @@
 
 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition
-      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd    (:,:,2) - rhd    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)
-      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2)
+      DO_2D( 1, 1, 1, 1 )
+         zdrho_k(ji,jj,1) = aco_bc_vrt * ( rhd  (ji,jj,2) - rhd  (ji,jj,1) ) - bco_bc_vrt * zdrho_k(ji,jj,2)
+         zdz_k  (ji,jj,1) = aco_bc_vrt * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_vrt * zdz_k  (ji,jj,2)
+      END_2D
 
       DO_2D( 1, 1, 1, 1 )
@@ -785,13 +790,17 @@
       !  5. compute and store elementary horizontal differences in provisional arrays 
       !----------------------------------------------------------------------------------------
-
-      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
-         zdrhox(ji,jj,jk) =   rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
-         zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  )
-         zdrhoy(ji,jj,jk) =   rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
-         zdzy  (ji,jj,jk) = - gde3w(ji  ,jj+1,jk) + gde3w(ji,jj,jk  )
-      END_3D
-
-      CALL lbc_lnk( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 
+      zdrhox(:,:,:) = 0._wp
+      zdzx  (:,:,:) = 0._wp
+      zdrhoy(:,:,:) = 0._wp
+      zdzy  (:,:,:) = 0._wp
+
+      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
+         zdrhox(ji,jj,jk) = rhd  (ji+1,jj  ,jk) - rhd  (ji  ,jj  ,jk)
+         zdzx  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji+1,jj  ,jk)
+         zdrhoy(ji,jj,jk) = rhd  (ji  ,jj+1,jk) - rhd  (ji  ,jj  ,jk)
+         zdzy  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji  ,jj+1,jk)
+      END_3D
+
+      IF( nn_hls == 1 ) CALL lbc_lnk( 'dynhpg', zdrhox, 'U', -1._wp, zdzx, 'U', -1._wp, zdrhoy, 'V', -1._wp, zdzy, 'V', -1._wp )
 
       !-------------------------------------------------------------------------
@@ -800,31 +809,19 @@
 
       DO_3D( 0, 1, 0, 1, 1, jpkm1 )
-         cffu = 2._wp * zdrhox(ji-1,jj  ,jk) * zdrhox(ji,jj,jk  )
-         IF( cffu > zep ) THEN
-            zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) )
-         ELSE
-            zdrho_i(ji,jj,jk ) = 0._wp
-         ENDIF
-
-         cffx = 2._wp * zdzx  (ji-1,jj  ,jk) * zdzx  (ji,jj,jk  )
-         IF( cffx > zep ) THEN
-            zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) )
-         ELSE
-            zdz_i(ji,jj,jk) = 0._wp
-         ENDIF
-
-         cffv = 2._wp * zdrhoy(ji  ,jj-1,jk) * zdrhoy(ji,jj,jk  )
-         IF( cffv > zep ) THEN
-            zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) )
-         ELSE
-            zdrho_j(ji,jj,jk) = 0._wp
-         ENDIF
-
-         cffy = 2._wp * zdzy  (ji  ,jj-1,jk) * zdzy  (ji,jj,jk  )
-         IF( cffy > zep ) THEN
-            zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) )
-         ELSE
-            zdz_j(ji,jj,jk) = 0._wp
-         ENDIF
+         cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp )
+         z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk)
+         zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
+
+         cffx = MAX( 2._wp * zdzx(ji-1,jj,jk)   * zdzx(ji,jj,jk), 0._wp )
+         z1_cff = zdzx(ji-1,jj,jk)   + zdzx(ji,jj,jk)
+         zdz_i(ji,jj,jk)   = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
+
+         cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp )
+         z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk)
+         zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
+
+         cffy = MAX( 2._wp * zdzy(ji,jj-1,jk)   * zdzy(ji,jj,jk), 0._wp )
+         z1_cff = zdzy(ji,jj-1,jk)   + zdzy(ji,jj,jk)
+         zdz_j(ji,jj,jk)   = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
       END_3D
       
@@ -840,32 +837,30 @@
          zz_drho_j(:,:) = zdrho_j(:,:,jk)
          zz_dz_j  (:,:) = zdz_j  (:,:,jk)
-         DO_2D( 0, 1, 0, 1)
-            ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
-            IF (ji < jpi) THEN
-               IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN  
-                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 
-                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk)
-               END IF
+         ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
+         DO_2D( 0, 0, 0, 1 )
+            IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN
+               zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)
+               zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk)
             END IF
-            ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
-            IF (ji > 2) THEN
-               IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN
-                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)  
-                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk)
-               END IF
+         END_2D
+         ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
+         DO_2D( -1, 1, 0, 1 )
+            IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN
+               zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)
+               zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk)
             END IF
-            ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
-            IF (jj < jpj) THEN
-               IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN
-                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
-                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk)
-               END IF
-            END IF 
-            ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
-            IF (jj > 2) THEN
-               IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 
-                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 
-                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk)
-               END IF
+         END_2D
+         ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
+         DO_2D( 0, 1, 0, 0 )
+            IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN
+               zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
+               zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk)
+            END IF
+         END_2D
+         ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
+         DO_2D( 0, 1, -1, 1 )
+            IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN
+               zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)
+               zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk)
             END IF
          END_2D
@@ -974,15 +969,17 @@
       REAL(wp) :: zrhdt1
       REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zsshu_n, zsshv_n
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zpgu, zpgv   ! 2D workspace
+      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zsshu_n, zsshv_n
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zdept, zrhh
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
       !
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
+      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
+         ENDIF
       ENDIF
 
@@ -1001,64 +998,61 @@
       !
       IF( ln_wd_il ) THEN
-         ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
+         ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) )
          DO_2D( 0, 0, 0, 0 )
-          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
-               &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
-               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
-               &                                                      > rn_wdmin1 + rn_wdmin2
-          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (         &
-               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                &
-               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
-
-          IF(ll_tmp1) THEN
-            zcpx(ji,jj) = 1.0_wp
-          ELSE IF(ll_tmp2) THEN
-            ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
-            zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
-                        &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
-           
-             zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
-          ELSE
-            zcpx(ji,jj) = 0._wp
-          END IF
-   
-          ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
-               &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-               &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
-               &                                                      > rn_wdmin1 + rn_wdmin2
-          ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (      &
-               &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
-               &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
-
-          IF(ll_tmp1) THEN
-            zcpy(ji,jj) = 1.0_wp
-          ELSE IF(ll_tmp2) THEN
-            ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
-            zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
-                        &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
-             zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
-
+            ll_tmp1 = MIN(   ssh(ji,jj,Kmm)              ,   ssh(ji+1,jj,Kmm)                 ) >       &
+               &      MAX( -ht_0(ji,jj)                  , -ht_0(ji+1,jj)                     ) .AND.   &
+               &      MAX(   ssh(ji,jj,Kmm) + ht_0(ji,jj),   ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) >       &
+               &      rn_wdmin1 + rn_wdmin2
+            ll_tmp2 = ( ABS(   ssh(ji,jj,Kmm) -   ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND.                   &
+               &      ( MAX(   ssh(ji,jj,Kmm) ,   ssh(ji+1,jj,Kmm) ) >                                  &
+               &        MAX( -ht_0(ji,jj)     , -ht_0(ji+1,jj)     ) + rn_wdmin1 + rn_wdmin2 )
+
+            IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+            ELSE IF(ll_tmp2) THEN
+               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
+               zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp)
+            ELSE
+               zcpx(ji,jj) = 0._wp
+            END IF
+
+            ll_tmp1 = MIN(   ssh(ji,jj,Kmm)              ,   ssh(ji,jj+1,Kmm)                 ) >       &
+               &      MAX( -ht_0(ji,jj)                  , -ht_0(ji,jj+1)                     ) .AND.   &
+               &      MAX(   ssh(ji,jj,Kmm) + ht_0(ji,jj),   ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) >       &
+               &      rn_wdmin1 + rn_wdmin2
+            ll_tmp2 = ( ABS(   ssh(ji,jj,Kmm) -   ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND.                   &
+               &      ( MAX(   ssh(ji,jj,Kmm) ,   ssh(ji,jj+1,Kmm) ) >                                  &
+               &        MAX( -ht_0(ji,jj)     , -ht_0(ji,jj+1)     ) + rn_wdmin1 + rn_wdmin2 )
+
+            IF(ll_tmp1) THEN
+               zcpy(ji,jj) = 1.0_wp
+            ELSE IF(ll_tmp2) THEN
+               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm)) )
+               zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp)
             ELSE
                zcpy(ji,jj) = 0._wp
             ENDIF
          END_2D
-         CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       ENDIF
 
       ! Clean 3-D work arrays
       zhpi(:,:,:) = 0._wp
-      zrhh(:,:,:) = rhd(:,:,:)
+      zrhh(:,:,:) = rhd(A2D(nn_hls),:)
 
       ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
       DO_2D( 1, 1, 1, 1 )
-       jk = mbkt(ji,jj)
-       IF(     jk <=  1   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
-       ELSEIF( jk ==  2   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
-       ELSEIF( jk < jpkm1 ) THEN
-          DO jkk = jk+1, jpk
-             zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   &
-                &                      gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2))
-          END DO
-       ENDIF
+         jk = mbkt(ji,jj)
+         IF(     jk <=  1   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
+         ELSEIF( jk ==  2   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
+         ELSEIF( jk < jpkm1 ) THEN
+            DO jkk = jk+1, jpk
+               zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   &
+                  &                      gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2))
+            END DO
+         ENDIF
       END_2D
 
@@ -1082,17 +1076,17 @@
       ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
       DO_2D( 0, 1, 0, 1 )
-       zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
-          &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
-
-       ! assuming linear profile across the top half surface layer
-       zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1
+         zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
+            &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
+
+         ! assuming linear profile across the top half surface layer
+         zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1
       END_2D
 
       ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
       DO_3D( 0, 1, 0, 1, 2, jpkm1 )
-      zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  &
-         &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   &
-         &                           asp  (ji,jj,jk-1), bsp  (ji,jj,jk-1), &
-         &                           csp  (ji,jj,jk-1), dsp  (ji,jj,jk-1)  )
+         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  &
+            &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   &
+            &                           asp  (ji,jj,jk-1), bsp  (ji,jj,jk-1), &
+            &                           csp  (ji,jj,jk-1), dsp  (ji,jj,jk-1)  )
       END_3D
 
@@ -1107,158 +1101,153 @@
 !                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
 !!gm not this:
-       zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
-                      & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-       zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
-                      & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
-      END_2D
-
-      CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
+         zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
+                        & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
+         zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
+                        & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
+      END_2D
 
       DO_2D( 0, 0, 0, 0 )
-       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 
-       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )
+         zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) )
+         zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )
       END_2D
 
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-      zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
-      zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
+         zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
+         zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
       END_3D
 
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-      zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
-      zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
+         zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
+         zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
       END_3D
 
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-      zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
-      zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
-      zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
-      zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
+         zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+         zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
+         zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
+         zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
       END_3D
 
 
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-      zpwes = 0._wp; zpwed = 0._wp
-      zpnss = 0._wp; zpnsd = 0._wp
-      zuijk = zu(ji,jj,jk)
-      zvijk = zv(ji,jj,jk)
-
-      !!!!!     for u equation
-      IF( jk <= mbku(ji,jj) ) THEN
-         IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
-           jis = ji + 1; jid = ji
-         ELSE
-           jis = ji;     jid = ji +1
+         zpwes = 0._wp; zpwed = 0._wp
+         zpnss = 0._wp; zpnsd = 0._wp
+         zuijk = zu(ji,jj,jk)
+         zvijk = zv(ji,jj,jk)
+
+         !!!!!     for u equation
+         IF( jk <= mbku(ji,jj) ) THEN
+            IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
+              jis = ji + 1; jid = ji
+            ELSE
+              jis = ji;     jid = ji +1
+            ENDIF
+
+            ! integrate the pressure on the shallow side
+            jk1 = jk
+            DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
+               IF( jk1 == mbku(ji,jj) ) THEN
+                  zuijk = -zdept(jis,jj,jk1)
+                  EXIT
+               ENDIF
+               zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
+               zpwes = zpwes +                                      &
+                  integ_spline(zdept(jis,jj,jk1), zdeps,            &
+                                 asp(jis,jj,jk1), bsp(jis,jj,jk1),  &
+                                 csp(jis,jj,jk1), dsp(jis,jj,jk1))
+               jk1 = jk1 + 1
+            END DO
+
+            ! integrate the pressure on the deep side
+            jk1 = jk
+            DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
+               IF( jk1 == 1 ) THEN
+                  zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad)
+                  zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
+                                                    bsp(jid,jj,1)  , csp(jid,jj,1), &
+                                                    dsp(jid,jj,1)) * zdeps
+                  zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
+                  EXIT
+               ENDIF
+               zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
+               zpwed = zpwed +                                        &
+                  integ_spline(zdeps,             zdept(jid,jj,jk1),  &
+                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
+                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
+               jk1 = jk1 - 1
+            END DO
+
+            ! update the momentum trends in u direction
+            zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
+            IF( .NOT.ln_linssh ) THEN
+               zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
+                  &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) )
+            ELSE
+               zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
+            ENDIF
+            IF( ln_wd_il ) THEN
+               zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)
+               zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
+            ENDIF
+            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk)
          ENDIF
 
-         ! integrate the pressure on the shallow side
-         jk1 = jk
-         DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
-           IF( jk1 == mbku(ji,jj) ) THEN
-             zuijk = -zdept(jis,jj,jk1)
-             EXIT
-           ENDIF
-           zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
-           zpwes = zpwes +                                    &
-                integ_spline(zdept(jis,jj,jk1), zdeps,            &
-                       asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
-                       csp(jis,jj,jk1),    dsp(jis,jj,jk1))
-           jk1 = jk1 + 1
-         END DO
-
-         ! integrate the pressure on the deep side
-         jk1 = jk
-         DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
-           IF( jk1 == 1 ) THEN
-             zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad)
-             zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
-                                               bsp(jid,jj,1),   csp(jid,jj,1), &
-                                               dsp(jid,jj,1)) * zdeps
-             zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
-             EXIT
-           ENDIF
-           zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
-           zpwed = zpwed +                                        &
-                  integ_spline(zdeps,              zdept(jid,jj,jk1), &
-                         asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
-                         csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
-           jk1 = jk1 - 1
-         END DO
-
-         ! update the momentum trends in u direction
-
-         zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
-         IF( .NOT.ln_linssh ) THEN
-           zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
-              &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) )
-          ELSE
-           zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
+         !!!!!     for v equation
+         IF( jk <= mbkv(ji,jj) ) THEN
+            IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
+               jjs = jj + 1; jjd = jj
+            ELSE
+               jjs = jj    ; jjd = jj + 1
+            ENDIF
+
+            ! integrate the pressure on the shallow side
+            jk1 = jk
+            DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
+               IF( jk1 == mbkv(ji,jj) ) THEN
+                  zvijk = -zdept(ji,jjs,jk1)
+                  EXIT
+               ENDIF
+               zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
+               zpnss = zpnss +                                       &
+                  integ_spline(zdept(ji,jjs,jk1), zdeps,             &
+                               asp(ji,jjs,jk1),   bsp(ji,jjs,jk1),   &
+                               csp(ji,jjs,jk1),   dsp(ji,jjs,jk1) )
+              jk1 = jk1 + 1
+            END DO
+
+            ! integrate the pressure on the deep side
+            jk1 = jk
+            DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
+               IF( jk1 == 1 ) THEN
+                  zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)
+                  zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
+                                                    bsp(ji,jjd,1)  , csp(ji,jjd,1), &
+                                                    dsp(ji,jjd,1) ) * zdeps
+                  zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
+                  EXIT
+               ENDIF
+               zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
+               zpnsd = zpnsd +                                        &
+                  integ_spline(zdeps,             zdept(ji,jjd,jk1),  &
+                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1),  &
+                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
+               jk1 = jk1 - 1
+            END DO
+
+            ! update the momentum trends in v direction
+            zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
+            IF( .NOT.ln_linssh ) THEN
+               zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
+                       ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )
+            ELSE
+               zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
+            ENDIF
+            IF( ln_wd_il ) THEN
+               zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)
+               zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)
+            ENDIF
+
+            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)
          ENDIF
-         IF( ln_wd_il ) THEN
-            zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)
-            zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
-         ENDIF
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 
-      ENDIF
-
-      !!!!!     for v equation
-      IF( jk <= mbkv(ji,jj) ) THEN
-         IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
-           jjs = jj + 1; jjd = jj
-         ELSE
-           jjs = jj    ; jjd = jj + 1
-         ENDIF
-
-         ! integrate the pressure on the shallow side
-         jk1 = jk
-         DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
-           IF( jk1 == mbkv(ji,jj) ) THEN
-             zvijk = -zdept(ji,jjs,jk1)
-             EXIT
-           ENDIF
-           zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
-           zpnss = zpnss +                                      &
-                  integ_spline(zdept(ji,jjs,jk1), zdeps,            &
-                         asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
-                         csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
-           jk1 = jk1 + 1
-         END DO
-
-         ! integrate the pressure on the deep side
-         jk1 = jk
-         DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
-           IF( jk1 == 1 ) THEN
-             zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)
-             zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
-                                               bsp(ji,jjd,1),   csp(ji,jjd,1), &
-                                               dsp(ji,jjd,1) ) * zdeps
-             zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
-             EXIT
-           ENDIF
-           zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
-           zpnsd = zpnsd +                                        &
-                  integ_spline(zdeps,              zdept(ji,jjd,jk1), &
-                         asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
-                         csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
-           jk1 = jk1 - 1
-         END DO
-
-
-         ! update the momentum trends in v direction
-
-         zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
-         IF( .NOT.ln_linssh ) THEN
-            zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
-                    ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )
-         ELSE
-            zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
-         ENDIF
-         IF( ln_wd_il ) THEN
-            zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 
-            zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 
-         ENDIF
-
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)
-      ENDIF
          !
       END_3D
@@ -1279,91 +1268,79 @@
       !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   fsp, xsp           ! value and coordinate
-      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   asp, bsp, csp, dsp ! coefficients of the interpoated function
-      INTEGER                   , INTENT(in   ) ::   polynomial_type    ! 1: cubic spline   ;   2: Linear
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   fsp, xsp           ! value and coordinate
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(  out) ::   asp, bsp, csp, dsp ! coefficients of the interpoated function
+      INTEGER                             , INTENT(in   ) ::   polynomial_type    ! 1: cubic spline   ;   2: Linear
       !
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      INTEGER  ::   jpi, jpj, jpkm1
       REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
       REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha
-      REAL(wp) ::   zdf(size(fsp,3))
-      !!----------------------------------------------------------------------
-      !
-!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!!  
-      jpi   = size(fsp,1)
-      jpj   = size(fsp,2)
-      jpkm1 = MAX( 1, size(fsp,3) - 1 )
+      REAL(wp) ::   zdf(jpk)
+      !!----------------------------------------------------------------------
       !
       IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-           !!Fritsch&Butland's method, 1984 (preferred, but more computation)
-           !    DO jk = 2, jpkm1-1
-           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)
-           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
-           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
-           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
-           !
-           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
-           !
-           !       IF(zdf1 * zdf2 <= 0._wp) THEN
-           !           zdf(jk) = 0._wp
-           !       ELSE
-           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
-           !       ENDIF
-           !    END DO
-
-           !!Simply geometric average
-               DO jk = 2, jpkm1-1
-                  zdf1 = (fsp(ji,jj,jk  ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk  ) - xsp(ji,jj,jk-1))
-                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk  )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk  ))
-
-                  IF(zdf1 * zdf2 <= 0._wp) THEN
-                     zdf(jk) = 0._wp
-                  ELSE
-                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
-                  ENDIF
-               END DO
-
-               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
-                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) )           -  0.5_wp * zdf(2)
-               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
-                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1)
-
-               DO jk = 1, jpkm1 - 1
-                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
-                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
-                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
-                 zddf1  = -2._wp * ztmp1 + ztmp2
-                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
-                 zddf2  =  2._wp * ztmp1 - ztmp2
-
-                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
-                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
-                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
-                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
-                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
-                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
-                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
-                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
-                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
-               END DO
+         DO_2D( 1, 1, 1, 1 )
+            !!Fritsch&Butland's method, 1984 (preferred, but more computation)
+            !    DO jk = 2, jpkm1-1
+            !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)
+            !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
+            !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
+            !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
+            !
+            !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
+            !
+            !       IF(zdf1 * zdf2 <= 0._wp) THEN
+            !           zdf(jk) = 0._wp
+            !       ELSE
+            !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
+            !       ENDIF
+            !    END DO
+
+            !!Simply geometric average
+            DO jk = 2, jpk-2
+               zdf1 = (fsp(ji,jj,jk  ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk  ) - xsp(ji,jj,jk-1))
+               zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk  )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk  ))
+
+               IF(zdf1 * zdf2 <= 0._wp) THEN
+                  zdf(jk) = 0._wp
+               ELSE
+                  zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
+               ENDIF
             END DO
-         END DO
+
+            zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
+                       &          ( xsp(ji,jj,2) - xsp(ji,jj,1) )           -  0.5_wp * zdf(2)
+            zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
+                       &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2)
+
+            DO jk = 1, jpk-2
+               zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
+               ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
+               ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
+               zddf1  = -2._wp * ztmp1 + ztmp2
+               ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
+               zddf2  =  2._wp * ztmp1 - ztmp2
+
+               dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
+               csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
+               bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
+                             & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
+                             & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
+                             &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
+               asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
+                             &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
+                             &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
+            END DO
+         END_2D
 
       ELSEIF ( polynomial_type == 2 ) THEN     ! Linear
-         DO ji = 1, jpi
-            DO jj = 1, jpj
-               DO jk = 1, jpkm1-1
-                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
-                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
-
-                  dsp(ji,jj,jk) = 0._wp
-                  csp(ji,jj,jk) = 0._wp
-                  bsp(ji,jj,jk) = ztmp1 / zdxtmp
-                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+         DO_3D( 1, 1, 1, 1, 1, jpk-2 )
+            zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
+            ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
+
+            dsp(ji,jj,jk) = 0._wp
+            csp(ji,jj,jk) = 0._wp
+            bsp(ji,jj,jk) = ztmp1 / zdxtmp
+            asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
+         END_3D
          !
       ELSE
