Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90	(revision 5038)
@@ -25,6 +25,7 @@
    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 sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean
    USE sbcrnf          ! river runoff 
+   USE sbcisf          ! ice shelf 
    USE cla             ! cross land advection             (cla_div routine)
    USE in_out_manager  ! I/O manager
@@ -66,5 +67,6 @@
       !!         - compute the now divergence given by :
       !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
-      !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 
+      !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf)
+      !!      and cross land flow (div_cla) 
       !!              II. vorticity :
       !!         - save the curl computed at the previous time-step
@@ -225,6 +227,7 @@
       !                                                ! ===============
 
-      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
-      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence)
+      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field)
+      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (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
@@ -323,5 +326,6 @@
       !                                                ! ===============
 
-      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( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field)
       IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
       !
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5038)
@@ -30,4 +30,5 @@
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
+   LOGICAL, PUBLIC ::   ln_dynzad_zts   !: vertical advection with sub-timestepping (requires vector form)
    
    INTEGER ::   nadv   ! choice of the formulation and scheme for the advection
@@ -64,7 +65,10 @@
                       CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
                       CALL dyn_zad     ( kt )    ! vector form : vertical advection
-      CASE ( 1 ) 
+      CASE ( 1 )     
+                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping
+      CASE ( 2 ) 
                       CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme
-      CASE ( 2 )   
+      CASE ( 3 )   
                       CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme
       !
@@ -91,5 +95,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs
+      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
       !!----------------------------------------------------------------------
 
@@ -111,4 +115,5 @@
          WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
          WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -120,9 +125,12 @@
 
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
+      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
+          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
 
       !                               ! Set nadv
       IF( ln_dynadv_vec  )   nadv =  0 
-      IF( ln_dynadv_cen2 )   nadv =  1
-      IF( ln_dynadv_ubs  )   nadv =  2
+      IF( ln_dynzad_zts  )   nadv =  1
+      IF( ln_dynadv_cen2 )   nadv =  2
+      IF( ln_dynadv_ubs  )   nadv =  3
       IF( lk_esopa       )   nadv = -1
 
@@ -130,6 +138,7 @@
          WRITE(numout,*)
          IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
-         IF( nadv ==  1 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
-         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
+         IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
          IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection formulation'
       ENDIF
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90	(revision 5038)
@@ -15,11 +15,12 @@
    USE oce            ! ocean dynamics and tracers
    USE dom_oce        ! ocean space and time domain
-   USE trdmod_oce     ! ocean variables trends
-   USE trdmod         ! ocean dynamics trends
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   !
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -103,5 +104,5 @@
          zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
          zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
-         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
+         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
          zfu_t(:,:,:) = ua(:,:,:)
          zfv_t(:,:,:) = va(:,:,:)
@@ -153,5 +154,5 @@
          zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
          zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
-         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
+         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
       ENDIF
       !                                            ! Control print
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90	(revision 5038)
@@ -16,12 +16,13 @@
    USE oce            ! ocean dynamics and tracers
    USE dom_oce        ! ocean space and time domain
-   USE trdmod         ! ocean dynamics trends
-   USE trdmod_oce     ! ocean variables trends
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   !
    USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -196,5 +197,5 @@
          zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:)
          zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:)
-         CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )
+         CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt )
          zfu_t(:,:,:) = ua(:,:,:)
          zfv_t(:,:,:) = va(:,:,:)
@@ -245,5 +246,5 @@
          zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:)
          zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:)
-         CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )
+         CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt )
       ENDIF
       !                                            ! Control print
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90	(revision 5038)
@@ -10,21 +10,21 @@
 
    !!----------------------------------------------------------------------
-   !!   dyn_bfr      : Update the momentum trend with the bottom friction contribution
+   !!   dyn_bfr       : Update the momentum trend with the bottom friction contribution
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE zdf_oce         ! ocean vertical physics variables
-   USE zdfbfr          ! ocean bottom friction variables
-   USE trdmod          ! ocean active dynamics and tracers trends 
-   USE trdmod_oce      ! ocean variables trends
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE timing          ! Timing
-   USE wrk_nemo        ! Memory Allocation
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables 
+   USE zdf_oce        ! ocean vertical physics variables
+   USE zdfbfr         ! ocean bottom friction variables
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
+   USE timing         ! Timing
+   USE wrk_nemo       ! Memory Allocation
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   dyn_bfr    !  routine called by step.F90
+   PUBLIC   dyn_bfr   !  routine called by step.F90
 
    !! * Substitutions
@@ -57,7 +57,10 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_bfr')
       !
+!!gm issue: better to put the logical in step to control the call of zdf_bfr
+!!          ==> change the logical from ln_bfrimp to ln_bfr_exp !!
       IF( .NOT.ln_bfrimp) THEN     ! only for explicit bottom friction form
                                     ! implicit bfr is implemented in dynzdf_imp
 
+!!gm bug : time step is only rdt (not 2 rdt if euler start !)
         zm1_2dt = - 1._wp / ( 2._wp * rdt )
 
@@ -69,11 +72,6 @@
 
 
-# if defined key_vectopt_loop
-        DO jj = 1, 1
-           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
         DO jj = 2, jpjm1
            DO ji = 2, jpim1
-# endif
               ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels
               ikbv = mbkv(ji,jj)
@@ -82,4 +80,16 @@
               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu)
               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv)
+              
+              ! (ISF) stability criteria for top friction
+              ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels
+              ikbv = mikv(ji,jj)
+              !
+              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)
+              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) &
+                 &             * (1.-umask(ji,jj,1))
+              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) &
+                 &             * (1.-vmask(ji,jj,1))
+              ! (ISF)
+              
            END DO
         END DO
@@ -89,5 +99,5 @@
            ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
            ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-           CALL trd_mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt )
+           CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt )
            CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
         ENDIF
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 5038)
@@ -29,12 +29,15 @@
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
+   USE sbc_oce         ! surface variable (only for the flag with ice shelf)
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
-   USE trdmod          ! ocean dynamics trends
-   USE trdmod_oce      ! ocean variables trends
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   !
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lbclnk          ! lateral boundary condition
+   USE lbclnk          ! lateral boundary condition 
    USE lib_mpp         ! MPP library
+   USE eosbn2          ! compute density
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
@@ -74,5 +77,5 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!             - Save the trend (l_trddyn=T)
+      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
@@ -99,5 +102,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
          CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
       ENDIF
@@ -175,4 +178,6 @@
       IF( ln_hpg_prj )   ioptio = ioptio + 1
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
+      IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   &
+          &  CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' )
       !
    END SUBROUTINE dyn_hpg_init
@@ -315,11 +320,6 @@
 
       ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
-# if defined key_vectopt_loop
-         jj = 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-# endif
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
@@ -338,7 +338,5 @@
                va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
             ENDIF
-# if ! defined key_vectopt_loop
-         END DO
-# endif
+         END DO
       END DO
       !
@@ -368,12 +366,16 @@
       INTEGER, INTENT(in) ::   kt    ! ocean time-step index
       !!
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
-      !
-      IF( kt == nit000 ) THEN
+      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop
+      REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj, 2, ztstop) 
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
+      CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
+      !
+     IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
@@ -384,17 +386,74 @@
       zcoef0 = - grav * 0.5_wp
       ! To use density and not density anomaly
-      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
-      ELSE                 ;     znad = 0._wp         ! Fixed volume
-      ENDIF
-
-      ! Surface value
+!      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
+!      ELSE                 ;     znad = 0._wp         ! Fixed volume
+!      ENDIF
+      znad=1._wp
+      ! iniitialised to 0. zhpi zhpi 
+      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
+
+!==================================================================================     
+!=====Compute iceload and contribution of the half first wet layer ================= 
+!===================================================================================
+
+      ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
+      ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
+
+      ! compute density of the water displaced by the ice shelf 
+      zrhd = rhd ! save rhd
+      DO jk = 1, jpk
+           zdept(:,:)=gdept_1d(jk)
+           CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))
+      END DO
+      WHERE ( tmask(:,:,:) == 1._wp)
+        rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
+      END WHERE
+      
+      ! compute rhd at the ice/oce interface (ice shelf side)
+      CALL eos(ztstop,risfdep,zrhdtop_isf)
+
+      ! compute rhd at the ice/oce interface (ocean side)
+      DO ji=1,jpi
+        DO jj=1,jpj
+          ikt=mikt(ji,jj)
+          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
+          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
+        END DO
+      END DO
+      CALL eos(ztstop,risfdep,zrhdtop_oce)
+      !
+      ! Surface value + ice shelf gradient
+      ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
+      ziceload = 0._wp
+      DO jj = 1, jpj
+         DO ji = 1, jpi   ! vector opt.
+            ikt=mikt(ji,jj)
+            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+            DO jk=2,ikt-1
+               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
+                  &                              * (1._wp - tmask(ji,jj,jk))
+            END DO
+            IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
+                               &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
+         END DO
+      END DO
+      riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
+      ! compute zp from z=0 to first T wet point (correction due to zps not yet applied)
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ! hydrostatic pressure gradient along s-surfaces
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-            ! s-coordinate pressure gradient correction
+            ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1)
+            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
+            ! we assume ISF is in isostatic equilibrium
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    &
+               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   &
+               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &
+               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              ) 
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    &
+               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   &
+               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
+               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   &
+               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             ) 
+            ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
                &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
@@ -402,35 +461,155 @@
                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
-         END DO
-      END DO
-
-      ! interior value (2=<jk=<jpkm1)
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
+            va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
+         END DO
+      END DO
+!==================================================================================     
+!===== Compute partial cell contribution for the top cell ========================= 
+!==================================================================================
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            iku = miku(ji,jj) ; 
+            zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp
+            ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
+            ! u direction
+            IF ( iku .GT. 1 ) THEN
+               ! case iku
+               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            &
+                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       &
+                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
+               ! corrective term ( = 0 if z coordinate )
+               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
+               ! zhpi will be added in interior loop
+               ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap
+               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
+               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku)
+
+               ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
+               zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &    
+                  &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          &
+                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   &
+                  &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          &
+                  &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   )
+               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 
+            END IF
+               
+            ! v direction
+            ikv = mikv(ji,jj)
+            ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
+            IF ( ikv .GT. 1 ) THEN
+               ! case ikv
+               zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            &
+                     &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           &
+                     &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
+               ! corrective term ( = 0 if z coordinate )
+               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
+               ! zhpi will be added in interior loop
+               va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap
+               ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure  
+               IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv) 
+               
+               ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
+               zhpjint        =  zcoef0 / e2v(ji,jj)                                                              &
+                  &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         &
+                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    &
+                  &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         &
+                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  )
+               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
+            END IF
+         END DO
+      END DO
+
+!==================================================================================     
+!===== Compute interior value ===================================================== 
+!==================================================================================
+
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            iku=miku(ji,jj); ikv=mikv(ji,jj)
+            DO jk = 2, jpkm1
                ! hydrostatic pressure gradient along s-surfaces
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
+               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
+               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              &
+                  &                            + zcoef0 / e1u(ji,jj)                                                           &
+                  &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 &
+                  &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   &
+                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 &
+                  &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   ) 
                ! s-coordinate pressure gradient correction
-               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
-               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
+               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
+               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    &
+                  &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
+
+               ! hydrostatic pressure gradient along s-surfaces
+               ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc)
+               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              &
+                  &                            + zcoef0 / e2v(ji,jj)                                                           &
+                  &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 &
+                  &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   &
+                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 &
+                  &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   )
+               ! s-coordinate pressure gradient correction
+               ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
+               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     &
+                  &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+
+!==================================================================================     
+!===== Compute bottom cell contribution (partial cell) ============================ 
+!==================================================================================
+
+# if defined key_vectopt_loop
+         jj = 1
+         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
+# else
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+# endif
+            iku = mbku(ji,jj)
+            ikv = mbkv(ji,jj)
+
+            IF (iku .GT. 1) THEN
+               ! remove old value (interior case)
+               zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   &
+                     &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj)
+               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
+               ! put new value
+               ! -zpshpi to avoid double contribution of the partial step in the top layer 
+               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj)
+               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 
+               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
+            END IF
+            ! v direction
+            IF (ikv .GT. 1) THEN
+               ! remove old value (interior case)
+               zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   &
+                     &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj)
+               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
+               ! put new value
+               ! -zpshpj to avoid double contribution of the partial step in the top layer 
+               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj)
+               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)   
+               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
+            END IF
+# if ! defined key_vectopt_loop
+         END DO
+# endif
+      END DO
+     
+      ! set back to original density value into the ice shelf cell (maybe useless because it is masked)
+      rhd = zrhd
+      !
+      CALL wrk_dealloc( jpi,jpj,2, ztstop)
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
+      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
       !
    END SUBROUTINE hpg_sco
+
 
    SUBROUTINE hpg_djc( kt )
@@ -671,6 +850,4 @@
       !!
       !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!             - Save the trend (l_trddyn=T)
-      !!
       !!----------------------------------------------------------------------
       INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
@@ -724,13 +901,19 @@
 
       ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
-      DO jj = 1, jpj;   DO ji = 1, jpi
-          zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad
-      END DO        ;   END DO
-
-      DO jk = 2, jpk;   DO jj = 1, jpj;   DO ji = 1, jpi
-          zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)
-      END DO        ;   END DO        ;   END DO
-
-      fsp(:,:,:) = zrhh(:,:,:)
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad
+         END DO
+      END DO
+
+      DO jk = 2, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      fsp(:,:,:) = zrhh (:,:,:)
       xsp(:,:,:) = zdept(:,:,:)
 
@@ -933,4 +1116,5 @@
    END SUBROUTINE hpg_prj
 
+
    SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
       !!----------------------------------------------------------------------
@@ -940,6 +1124,6 @@
       !!
       !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3
+      !!
       !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
-      !!
       !!----------------------------------------------------------------------
       IMPLICIT NONE
@@ -949,6 +1133,5 @@
       INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline
                                                                     ! 2: Linear
-
-      ! Local Variables
+      !
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
       INTEGER  ::   jpi, jpj, jpkm1
@@ -1040,5 +1223,4 @@
       ENDIF
 
-
    END SUBROUTINE cspline
 
@@ -1050,8 +1232,6 @@
       !! ** Purpose :   1-d linear interpolation
       !!
-      !! ** Method  :
-      !!                interpolation is straight forward
+      !! ** Method  :   interpolation is straight forward
       !!                extrapolation is also permitted (no value limit)
-      !!
       !!----------------------------------------------------------------------
       IMPLICIT NONE
@@ -1070,4 +1250,5 @@
    END FUNCTION interp1
 
+
    FUNCTION interp2(x, a, b, c, d)  RESULT(f)
       !!----------------------------------------------------------------------
@@ -1133,5 +1314,4 @@
    END FUNCTION integ_spline
 
-
    !!======================================================================
 END MODULE dynhpg
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5038)
@@ -14,6 +14,7 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -52,8 +53,8 @@
       !!
       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
-      !!             - save this trends (l_trddyn=T) for post-processing
+      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
+      !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zu, zv       ! temporary scalars
@@ -131,5 +132,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
          CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
       ENDIF
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 5038)
@@ -15,4 +15,5 @@
    USE phycst         ! physical constants
    USE ldfdyn_oce     ! ocean dynamics lateral physics
+   USE ldftra_oce     ! ocean tracers  lateral physics
    USE ldfslp         ! lateral mixing: slopes of mixing orientation
    USE dynldf_bilapg  ! lateral mixing            (dyn_ldf_bilapg routine)
@@ -20,7 +21,7 @@
    USE dynldf_iso     ! lateral mixing            (dyn_ldf_iso    routine)
    USE dynldf_lap     ! lateral mixing            (dyn_ldf_lap    routine)
-   USE ldftra_oce, ONLY: ln_traldf_hor     ! ocean tracers lateral physics
-   USE trdmod         ! ocean dynamics and tracer trends
-   USE trdmod_oce     ! ocean variables trends
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics   (trd_dyn        routine)
+   !
    USE prtctl         ! Print control
    USE in_out_manager ! I/O manager
@@ -30,5 +31,4 @@
    USE timing          ! Timing
 
-
    IMPLICIT NONE
    PRIVATE
@@ -55,5 +55,4 @@
       !! ** Purpose :   compute the lateral ocean dynamics physics.
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
@@ -107,5 +106,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv )
       ENDIF
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90	(revision 5038)
@@ -19,7 +19,6 @@
    USE dom_oce         ! ocean space and time domain
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
+   !
    USE in_out_manager  ! I/O manager
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE wrk_nemo        ! Memory Allocation
@@ -70,11 +69,8 @@
       !!      Add this before trend to the general trend (ua,va):
       !!            (ua,va) = (ua,va) + (diffu,diffv)
-      !!      'key_trddyn' defined: the two components of the horizontal
-      !!                               diffusion trend are saved.
       !!
       !! ** Action : - Update (ua,va) with the before iso-level biharmonic
       !!               mixing trend.
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90	(revision 5038)
@@ -19,9 +19,8 @@
    USE dom_oce         ! ocean space and time domain
    USE ldfdyn_oce      ! ocean dynamics lateral physics
+   USE zdf_oce         ! ocean vertical physics
+   USE ldfslp          ! iso-neutral slopes available
    USE ldftra_oce, ONLY: ln_traldf_iso
-   USE zdf_oce         ! ocean vertical physics
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
-   USE ldfslp          ! iso-neutral slopes available
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -81,9 +80,7 @@
       !!         -3- Add this trend to the general trend (ta,sa):
       !!            (ua,va) = (ua,va) + (zwk3,zwk4)
-      !!      'key_trddyn' defined: the trend is saved for diagnostics.
       !!
       !! ** Action  : - Update (ua,va) arrays with the before geopotential
       !!                biharmonic mixing trend.
-      !!              - save the trend in (zwk3,zwk4) ('key_trddyn')
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
@@ -201,7 +198,5 @@
       !!                          pu and pv (all the components except
       !!                          second order vertical derivative term)
-      !!      'key_trddyn' defined: the trend is saved for diagnostics.
-      !!----------------------------------------------------------------------
-      !!
+      !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity 
       !                                                               ! 2nd call: ahm x these fields
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 5038)
@@ -22,7 +22,6 @@
    USE ldftra_oce      ! ocean tracer   lateral physics
    USE zdf_oce         ! ocean vertical physics
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
    USE ldfslp          ! iso-neutral slopes 
+   !
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90	(revision 5038)
@@ -19,8 +19,6 @@
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
    USE zdf_oce         ! ocean vertical physics
+   !
    USE in_out_manager  ! I/O manager
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
-   USE ldfslp          ! iso-neutral slopes 
    USE timing          ! Timing
 
@@ -57,9 +55,6 @@
       !!      Add this before trend to the general trend (ua,va):
       !!            (ua,va) = (ua,va) + (diffu,diffv)
-      !!      'key_trddyn' activated: the two components of the horizontal
-      !!                                 diffusion trend are saved.
       !!
-      !! ** Action : - Update (ua,va) with the before iso-level harmonic 
-      !!               mixing trend.
+      !! ** Action : - Update (ua,va) with the iso-level harmonic mixing trend
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 5038)
@@ -18,4 +18,5 @@
    !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
    !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
+   !!            3.7  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
    !!-------------------------------------------------------------------------
   
@@ -34,14 +35,18 @@
    USE bdydyn          ! ocean open boundary conditions
    USE bdyvol          ! ocean open boundary condition (bdy_vol routines)
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   USE trdken          ! trend manager: kinetic energy
+   !
    USE in_out_manager  ! I/O manager
+   USE iom             ! I/O manager library
    USE lbclnk          ! lateral boundary condition (or mpp link)
    USE lib_mpp         ! MPP library
    USE wrk_nemo        ! Memory Allocation
    USE prtctl          ! Print control
-
+   USE timing          ! Timing
 #if defined key_agrif
    USE agrif_opa_interp
 #endif
-   USE timing          ! Timing
 
    IMPLICIT NONE
@@ -79,5 +84,5 @@
       !!             at the local domain boundaries through lbc_lnk call,
       !!             at the one-way open boundaries (lk_bdy=T),
-      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
+      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
       !!
       !!              * Apply the time filter applied and swap of the dynamics
@@ -99,14 +104,14 @@
       REAL(wp) ::   z2dt         ! temporary scalar
 #endif
-      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
-      REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
+      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec      ! local scalars
+      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zue, zve
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f, zua, zva 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_nxt')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva )
+      IF( nn_timing == 1 )   CALL timing_start('dyn_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
+      IF( lk_dynspg_ts )   CALL wrk_alloc( jpi,jpj, zue, zve )
       !
       IF( kt == nit000 ) THEN
@@ -152,15 +157,16 @@
 
 # if defined key_dynspg_ts
+!!gm IF ( lk_dynspg_ts ) THEN ....
       ! Ensure below that barotropic velocities match time splitting estimate
       ! Compute actual transport and replace it with ts estimate at "after" time step
-      zua(:,:) = 0._wp
-      zva(:,:) = 0._wp
-      DO jk = 1, jpkm1
-         zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+      zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
+      zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
+      DO jk = 2, jpkm1
+         zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
+         zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
       END DO
       DO jk = 1, jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
+         ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
       END DO
 
@@ -175,4 +181,5 @@
          END DO  
       ENDIF
+!!gm ENDIF
 # endif
 
@@ -195,4 +202,23 @@
 # endif
 #endif
+
+      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
+         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step 
+         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
+         !
+         !                                  ! Kinetic energy and Conversion
+         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
+         !
+         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
+            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
+            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
+            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
+            CALL iom_put( "vtrd_tot", zva )
+         ENDIF
+         !
+         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
+         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
+         !                                  !  computation of the asselin filter trends)
+      ENDIF
 
       ! Time filter and swap of dynamics arrays
@@ -217,6 +243,6 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi    
-                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )
-                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )
+                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
+                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                      !
                      ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
@@ -301,13 +327,13 @@
             ! Revert "before" velocities to time split estimate
             ! Doing it here also means that asselin filter contribution is removed  
-            zua(:,:) = 0._wp
-            zva(:,:) = 0._wp
-            DO jk = 1, jpkm1
-               zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
-               zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
+            zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
+            zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)    
+            DO jk = 2, jpkm1
+               zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
+               zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)    
             END DO
             DO jk = 1, jpkm1
-               ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
-               vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
+               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)
+               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)
             END DO
          ENDIF
@@ -327,6 +353,6 @@
             hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
          END DO
-         hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
-         hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
+         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) )
+         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) )
       ENDIF
       !
@@ -335,11 +361,6 @@
       !
       DO jk = 1, jpkm1
-#if defined key_vectopt_loop
-         DO jj = 1, 1         !Vector opt. => forced unrolling
-            DO ji = 1, jpij
-#else 
          DO jj = 1, jpj
             DO ji = 1, jpi
-#endif                  
                un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
                vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
@@ -358,9 +379,16 @@
       !
       !
+
+      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
+         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
+         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
+         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
+      ENDIF
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
       ! 
-      CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f )
-      IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva )
+      CALL wrk_dealloc( jpi,jpj,jpk,  ze3u_f, ze3v_f, zua, zva )
+      IF( lk_dynspg_ts )   CALL wrk_dealloc( jpi,jpj, zue, zve )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_nxt')
@@ -370,3 +398,2 @@
    !!=========================================================================
 END MODULE dynnxt
-
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 5038)
@@ -26,12 +26,13 @@
    USE sbctide
    USE updtide
-   USE trdmod         ! ocean dynamics trends
-   USE trdmod_oce     ! ocean variables trends
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   !
    USE prtctl         ! Print control                     (prt_ctl routine)
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
-   USE solver          ! solver initialization
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE solver         ! solver initialization
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
 
@@ -163,5 +164,8 @@
                END DO
             END DO
-         END DO         
+         END DO    
+         
+!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
+              
       ENDIF
 
@@ -191,9 +195,9 @@
          CASE( 2 )
             z2dt = 2. * rdt
-            IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
+            IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
             ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:)
             ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:)
          END SELECT
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt )
          !
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
@@ -246,4 +250,6 @@
       IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   &
            &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' )
+      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0 )   &
+           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' )
       !
       IF( lk_esopa     )   nspg = -1
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 5038)
@@ -19,4 +19,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE phycst          ! physical constants
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing library
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 5038)
@@ -13,4 +13,5 @@
    !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines. 
    !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
+   !!            3.7  !  2014-04  (F. Roquet, G. Madec)  add some trends diag
    !!----------------------------------------------------------------------
 #if defined key_dynspg_flt   ||   defined key_esopa  
@@ -36,4 +37,7 @@
    USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
    USE cla             ! cross land advection
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! distributed memory computing library
@@ -43,8 +47,8 @@
    USE iom
    USE lib_fortran
+   USE timing          ! Timing
 #if defined key_agrif
    USE agrif_opa_interp
 #endif
-   USE timing          ! Timing
 
    IMPLICIT NONE
@@ -99,15 +103,16 @@
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
       !!
-      !! References : Roullet and Madec 1999, JGR.
+      !! References : Roullet and Madec, JGR, 2000.
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
       INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge)
-      !!                                   
+      !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpw
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
-      !
       !
       IF( kt == nit000 ) THEN
@@ -179,4 +184,17 @@
          END DO
          !
+         IF( l_trddyn )   THEN                      ! temporary save of spg trends
+            CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
+            DO jk = 1, jpkm1              ! unweighted time stepping 
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk)
+                     ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt )
+         ENDIF
+         !
       ENDIF
 
@@ -194,31 +212,18 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            spgu(ji,jj) = 0._wp
-            spgv(ji,jj) = 0._wp
-         END DO
-      END DO
-
-      ! vertical sum
-!CDIR NOLOOPCHG
-      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
-         DO jk = 1, jpkm1
-            DO ji = 1, jpij
-               spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk)
-               spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk)
-            END DO
-         END DO
-      ELSE                        ! No  vector opt.
-         DO jk = 1, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
-                  spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-
-      ! transport: multiplied by the horizontal scale factor
-      DO jj = 2, jpjm1
+            spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1)
+            spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1)
+         END DO
+      END DO
+      DO jk = 2, jpkm1                     ! vertical sum
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
+               spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                     ! transport: multiplied by the horizontal scale factor
          DO ji = fs_2, fs_jpim1   ! vector opt.
             spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
@@ -322,4 +327,26 @@
       ENDIF
 #endif      
+
+      IF( l_trddyn )   THEN                     
+         ztrdu(:,:,:) = ua(:,:,:)                 ! save the after velocity before the filtered SPG
+         ztrdv(:,:,:) = va(:,:,:)
+         !
+         CALL wrk_alloc( jpi, jpj, zpw )
+         !
+         zpw(:,:) = - z2dt * gcx(:,:)
+         CALL iom_put( "ssh_flt" , zpw )          ! output equivalent ssh modification due to implicit filter
+         !
+         !                                        ! save surface pressure flux: -pw at z=0
+         zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1)
+         CALL iom_put( "pw0_exp" , zpw )
+         zpw(:,:) = wn(:,:,1)
+         CALL iom_put( "w0" , zpw )
+         zpw(:,:) =  rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1)
+         CALL iom_put( "pw0_flt" , zpw )
+         !
+         CALL wrk_dealloc( jpi, jpj, zpw ) 
+         !                                   
+      ENDIF
+      
       ! Add the trends multiplied by z2dt to the after velocity
       ! -------------------------------------------------------
@@ -336,10 +363,15 @@
       END DO
 
-      ! write filtered free surface arrays in restart file
-      ! --------------------------------------------------
-      IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' )
-      !
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_flt')
+      IF( l_trddyn )   THEN                      ! save the explicit SPG trends for further diagnostics
+         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt
+         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt )
+         !
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+      ENDIF
+
+      IF( lrst_oce )   CALL flt_rst( kt, 'WRITE' )      ! write filtered free surface arrays in restart file
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_spg_flt')
       !
    END SUBROUTINE dyn_spg_flt
@@ -352,6 +384,6 @@
       !! ** Purpose : Read or write filtered free surface arrays in restart file
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
+      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !!----------------------------------------------------------------------
       !
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 5038)
@@ -97,6 +97,6 @@
       ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) )
 
-      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
-                             &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
+      IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 
+                                                    &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
 
       dyn_spg_ts_alloc = MAXVAL(ierr(:))
@@ -218,5 +218,25 @@
       !
       IF ( kt == nit000 .OR. lk_vvl ) THEN
-         IF ( ln_dynvor_een ) THEN
+         IF ( ln_dynvor_een_old ) THEN
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    &
+                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp  
+                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( zwz, 'F', 1._wp )
+            zwz(:,:) = ff(:,:) * zwz(:,:)
+
+            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp
+            DO jj = 2, jpj
+               DO ji = fs_2, jpi   ! vector opt.
+                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
+                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
+                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
+                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
+               END DO
+            END DO
+         ELSE IF ( ln_dynvor_een ) THEN
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
@@ -339,5 +359,5 @@
          END DO
          !
-      ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme
+      ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN  ! enstrophy and energy conserving scheme
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -687,5 +707,5 @@
             END DO
             !
-         ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==!
+         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 5038)
@@ -15,4 +15,5 @@
    !!            3.2  ! 2009-04  (R. Benshila)  vvl: correction of een scheme
    !!            3.3  ! 2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.7  ! 2014-04  (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 
    !!----------------------------------------------------------------------
 
@@ -29,6 +30,6 @@
    USE dommsk         ! ocean mask
    USE dynadv         ! momentum advection (use ln_dynadv_vec value)
-   USE trdmod         ! ocean dynamics trends 
-   USE trdmod_oce     ! ocean variables trends
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE prtctl         ! Print control
@@ -50,4 +51,5 @@
    LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme
    LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme
+   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation)
 
    INTEGER ::   nvor = 0   ! type of vorticity trend used
@@ -73,5 +75,6 @@
       !! ** Action : - Update (ua,va) with the now vorticity term trend
       !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) ('key_trddyn')
+      !!               and planetary vorticity trends) and send them to trd_dyn 
+      !!               for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
@@ -108,5 +111,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -114,6 +117,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
             CALL vor_ene( kt, ntot, ua, va )                ! total vorticity
@@ -127,5 +129,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -133,6 +135,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
             CALL vor_ens( kt, ntot, ua, va )                ! total vorticity
@@ -146,5 +147,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -152,6 +153,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
             CALL vor_mix( kt )                               ! total vorticity (mix=ens-ene)
@@ -165,5 +165,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt )
             ztrdu(:,:,:) = ua(:,:,:)
             ztrdv(:,:,:) = va(:,:,:)
@@ -171,6 +171,5 @@
             ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
             ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt )
-            CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt )
+            CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt )
          ELSE
             CALL vor_een( kt, ntot, ua, va )                ! total vorticity
@@ -211,6 +210,4 @@
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) ('key_trddyn')
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
@@ -328,6 +325,4 @@
       !!
       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) ('key_trddyn')
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
@@ -444,6 +439,4 @@
       !!
       !! ** Action : - Update (ua,va) arrays with the now vorticity term trend
-      !!             - Save the trends in (ztrdu,ztrdv) in 2 parts (relative 
-      !!               and planetary vorticity trends) ('key_trddyn')
       !!
       !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.
@@ -557,6 +550,4 @@
       !!
       !! ** Action : - Update (ua,va) with the now vorticity term trend
-      !!             - save the trends in (ztrdu,ztrdv) in 2 parts (relative
-      !!               and planetary vorticity trends) ('key_trddyn')
       !!
       !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
@@ -601,20 +592,34 @@
             IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )
          ENDIF
-         ze3f(:,:,:) = 0.d0
+         ze3f(:,:,:) = 0._wp
 #endif
       ENDIF
 
       IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
-         DO jk = 1, jpk
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
-                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
-                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
-                  IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
-               END DO
-            END DO
-         END DO
+
+         IF( ln_dynvor_een_old ) THEN ! original formulation
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3
+                  END DO
+               END DO
+            END DO
+         ELSE ! new formulation from NEMO 3.6
+            DO jk = 1, jpk
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1
+                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
+                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
+                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
+                  END DO
+               END DO
+            END DO
+         ENDIF
+
          CALL lbc_lnk( ze3f, 'F', 1. )
       ENDIF
@@ -715,5 +720,5 @@
       INTEGER ::   ios             ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een
+      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old
       !!----------------------------------------------------------------------
 
@@ -736,4 +741,5 @@
          WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix
          WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een
+         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old
       ENDIF
 
@@ -759,4 +765,5 @@
       IF( ln_dynvor_mix )   ioptio = ioptio + 1
       IF( ln_dynvor_een )   ioptio = ioptio + 1
+      IF( ln_dynvor_een_old )   ioptio = ioptio + 1
       IF( lk_esopa      )   ioptio =          1
 
@@ -767,5 +774,5 @@
       IF( ln_dynvor_ens )   nvor =  1
       IF( ln_dynvor_mix )   nvor =  2
-      IF( ln_dynvor_een )   nvor =  3
+      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3
       IF( lk_esopa      )   nvor = -1
       
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90	(revision 5038)
@@ -16,16 +16,18 @@
    USE dom_oce        ! ocean space and time domain
    USE sbc_oce        ! surface boundary condition: ocean
-   USE trdmod_oce     ! ocean variables trends
-   USE trdmod         ! ocean dynamics trends 
+   USE trd_oce        ! trends: ocean variables
+   USE trddyn         ! trend manager: dynamics
+   !
    USE in_out_manager ! I/O manager
-   USE lib_mpp         ! MPP library
+   USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
    PRIVATE
    
-   PUBLIC   dyn_zad   ! routine called by step.F90
+   PUBLIC   dyn_zad       ! routine called by dynadv.F90
+   PUBLIC   dyn_zad_zts   ! routine called by dynadv.F90
 
    !! * Substitutions
@@ -53,5 +55,5 @@
       !!
       !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
-      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+      !!              - Send the trends to trddyn for diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
@@ -83,5 +85,5 @@
          DO jj = 2, jpj                   ! vertical fluxes 
             DO ji = fs_2, jpi             ! vector opt.
-               zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+               zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
             END DO
          END DO
@@ -95,8 +97,8 @@
       DO jj = 2, jpjm1              ! Surface and bottom values set to zero
          DO ji = fs_2, fs_jpim1           ! vector opt.
-            zwuw(ji,jj, 1 ) = 0.e0
-            zwvw(ji,jj, 1 ) = 0.e0
-            zwuw(ji,jj,jpk) = 0.e0
-            zwvw(ji,jj,jpk) = 0.e0
+            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
+            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
          END DO  
       END DO
@@ -118,5 +120,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
       ENDIF
@@ -132,4 +134,146 @@
    END SUBROUTINE dyn_zad
 
+   SUBROUTINE dyn_zad_zts ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dynzad_zts  ***
+      !! 
+      !! ** Purpose :   Compute the now vertical momentum advection trend and 
+      !!      add it to the general trend of momentum equation. This version
+      !!      uses sub-timesteps for improved numerical stability with small
+      !!      vertical grid sizes. This is especially relevant when using 
+      !!      embedded ice with thin surface boxes.
+      !!
+      !! ** Method  :   The now vertical advection of momentum is given by:
+      !!         w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]
+      !!         w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]
+      !!      Add this trend to the general trend (ua,va):
+      !!         (ua,va) = (ua,va) + w dz(u,v)
+      !!
+      !! ** Action  : - Update (ua,va) with the vert. momentum adv. trends
+      !!              - Save the trends in (ztrdu,ztrdv) ('key_trddyn')
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
+      !
+      INTEGER  ::   ji, jj, jk, jl  ! dummy loop indices
+      INTEGER  ::   jnzts = 5       ! number of sub-timesteps for vertical advection
+      INTEGER  ::   jtb, jtn, jta   ! sub timestep pointers for leap-frog/euler forward steps
+      REAL(wp) ::   zua, zva        ! temporary scalars
+      REAL(wp) ::   zr_rdt          ! temporary scalar
+      REAL(wp) ::   z2dtzts         ! length of Euler forward sub-timestep for vertical advection
+      REAL(wp) ::   zts             ! length of sub-timestep for vertical advection
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zwuw , zwvw, zww
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  zus , zvs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_zad_zts')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps'
+      ENDIF
+
+      IF( l_trddyn )   THEN         ! Save ua and va trends
+         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
+      ENDIF
+      
+      IF( neuler == 0 .AND. kt == nit000 ) THEN
+          z2dtzts =         rdt / REAL( jnzts, wp )   ! = rdt (restart with Euler time stepping)
+      ELSE
+          z2dtzts = 2._wp * rdt / REAL( jnzts, wp )   ! = 2 rdt (leapfrog)
+      ENDIF
+      
+      DO jk = 2, jpkm1                    ! Calculate and store vertical fluxes
+         DO jj = 2, jpj                   
+            DO ji = fs_2, jpi             ! vector opt.
+               zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero
+         DO ji = fs_2, fs_jpim1           ! vector opt.
+            zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp
+            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp
+            zwuw(ji,jj,jpk) = 0._wp
+            zwvw(ji,jj,jpk) = 0._wp
+         END DO  
+      END DO
+
+! Start with before values and use sub timestepping to reach after values
+
+      zus(:,:,:,1) = ub(:,:,:)
+      zvs(:,:,:,1) = vb(:,:,:)
+
+      DO jl = 1, jnzts                   ! Start of sub timestepping loop
+
+         IF( jl == 1 ) THEN              ! Euler forward to kick things off
+           jtb = 1   ;   jtn = 1   ;   jta = 2
+           zts = z2dtzts
+         ELSEIF( jl == 2 ) THEN          ! First leapfrog step
+           jtb = 1   ;   jtn = 2   ;   jta = 3
+           zts = 2._wp * z2dtzts
+         ELSE                            ! Shuffle pointers for subsequent leapfrog steps
+           jtb = MOD(jtb,3) + 1
+           jtn = MOD(jtn,3) + 1
+           jta = MOD(jta,3) + 1
+         ENDIF
+
+         DO jk = 2, jpkm1           ! Vertical momentum advection at level w and u- and v- vertical
+            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point
+               DO ji = fs_2, fs_jpim1        ! vector opt.
+                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) )
+                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) )
+               END DO  
+            END DO   
+         END DO
+         DO jk = 1, jpkm1           ! Vertical momentum advection at u- and v-points
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1       ! vector opt.
+                  !                         ! vertical momentum advective trends
+                  zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
+                  zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
+                  zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts
+                  zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts
+               END DO  
+            END DO  
+         END DO
+
+      END DO      ! End of sub timestepping loop
+
+      zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts )
+      DO jk = 1, jpkm1              ! Recover trends over the outer timestep
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1       ! vector opt.
+               !                         ! vertical momentum advective trends
+               !                         ! add the trends to the general momentum trends
+               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt
+               va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt
+            END DO  
+         END DO  
+      END DO
+
+      IF( l_trddyn ) THEN           ! save the vertical advection trends for diagnostic
+         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt )
+         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
+      ENDIF
+      !                             ! Control print
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dyn_zad_zts')
+      !
+   END SUBROUTINE dyn_zad_zts
+
    !!======================================================================
 END MODULE dynzad
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 5038)
@@ -20,6 +20,6 @@
 
    USE ldfdyn_oce      ! ocean dynamics: lateral physics
-   USE trdmod          ! ocean active dynamics and tracers trends 
-   USE trdmod_oce      ! ocean variables trends
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -91,6 +91,5 @@
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt )
-         !
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
       ENDIF
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 5038)
@@ -70,6 +70,4 @@
       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars
       REAL(wp) ::   ze3ua, ze3va
-      !!----------------------------------------------------------------------
-
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwd, zws
       !!----------------------------------------------------------------------
@@ -101,15 +99,14 @@
 
       IF( ln_bfrimp ) THEN
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-# endif
                ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
                ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
                avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
                avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
+               ikbu = miku(ji,jj)       ! ocean top level at u- and v-points 
+               ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
+               IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
+               IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
             END DO
          END DO
@@ -138,6 +135,6 @@
             ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
             va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
-         ENDDO
-         ! Add bottom stress due to barotropic component only:
+         END DO
+         ! Add bottom/top stress due to barotropic component only:
          DO jj = 2, jpjm1        
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -148,4 +145,10 @@
                ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
                va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
+               ikbu = miku(ji,jj)         ! top ocean level at u- and v-points 
+               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu)
+               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv)
+               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -166,5 +169,5 @@
                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  )
                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk)
-               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
+               zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1)
                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws
@@ -194,7 +197,8 @@
       !-----------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      DO jj = 2, jpjm1   
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = miku(ji,jj)+1, jpkm1
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
@@ -204,17 +208,14 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1) 
+            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj)) 
 #if defined key_dynspg_ts
-            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+            ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
                &                                      / ( ze3ua * rau0 ) 
 #else
-            ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
-               &                                                     / ( fse3u(ji,jj,1) * rau0     ) ) 
-#endif
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) &
+               &                   + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+               &                                  / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) ) 
+#endif
+            DO jk = miku(ji,jj)+1, jpkm1
 #if defined key_dynspg_ts
                zrhs = ua(ji,jj,jk)   ! zrhs=right hand side
@@ -230,9 +231,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = jpk-2, miku(ji,jj), -1
                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
@@ -292,7 +289,8 @@
       !-----------------------------------------------------------------------
       !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+      DO jj = 2, jpjm1   
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = mikv(ji,jj)+1, jpkm1        
                zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
             END DO
@@ -302,17 +300,14 @@
       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1) 
+            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj)) 
 #if defined key_dynspg_ts            
-            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+            va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
                &                                      / ( ze3va * rau0 ) 
 #else
-            va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
-               &                                                       / ( fse3v(ji,jj,1) * rau0     )  )
-#endif
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) &
+               &                   + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+               &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  )
+#endif
+            DO jk = mikv(ji,jj)+1, jpkm1
 #if defined key_dynspg_ts
                zrhs = va(ji,jj,jk)   ! zrhs=right hand side
@@ -328,9 +323,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
+            DO jk = jpk-2, mikv(ji,jj), -1
                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
             END DO
@@ -352,17 +343,16 @@
       !! restore bottom layer avmu(v) 
       IF( ln_bfrimp ) THEN
-# if defined key_vectopt_loop
-      DO jj = 1, 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-# endif
-            ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-            ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-            avmu(ji,jj,ikbu+1) = 0.e0
-            avmv(ji,jj,ikbv+1) = 0.e0
-         END DO
-      END DO
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+              ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+              avmu(ji,jj,ikbu+1) = 0.e0
+              avmv(ji,jj,ikbv+1) = 0.e0
+              ikbu = miku(ji,jj)         ! ocean top level at u- and v-points 
+              ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+              IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
+              IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
+           END DO
+        END DO
       ENDIF
       !
Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4792)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 5038)
@@ -31,5 +31,4 @@
    USE bdy_par         
    USE bdydyn2d        ! bdy_ssh routine
-   USE diaar5, ONLY:   lk_diaar5
    USE iom
 #if defined key_agrif
@@ -111,5 +110,5 @@
       ! 
       z1_rau0 = 0.5_wp * r1_rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
 
 #if ! defined key_dynspg_ts
@@ -138,6 +137,6 @@
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
-      CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
+      CALL iom_put( "ssh" , sshn )   ! sea surface height
+      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
       !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask, ovlap=1 )
@@ -233,5 +232,5 @@
       !                                           !------------------------------!
       CALL iom_put( "woce", wn )   ! vertical velocity
-      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
+      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          CALL wrk_alloc( jpi, jpj, z2d ) 
          CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
@@ -241,6 +240,6 @@
             z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
          END DO
-         CALL iom_put( "w_masstr" , z3d                     )  
-         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
+         CALL iom_put( "w_masstr" , z3d )  
+         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
          CALL wrk_dealloc( jpi, jpj, z2d  ) 
          CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
@@ -291,5 +290,5 @@
       ELSE                                         !** Leap-Frog time-stepping: Asselin filter + swap
          sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )     ! before <-- now filtered
-         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)
+         IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * ssmask(:,:)
          sshn(:,:) = ssha(:,:)                           ! now <-- after
       ENDIF
