Index: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 4705)
+++ branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 6225)
@@ -12,34 +12,32 @@
 
    !!----------------------------------------------------------------------
-   !!   ssh_nxt        : after ssh
-   !!   ssh_swp        : filter ans swap the ssh arrays
-   !!   wzv            : compute now vertical velocity
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain variables 
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE domvvl          ! Variable volume
-   USE divcur          ! hor. divergence and curl      (div & cur routines)
-   USE iom             ! I/O library
-   USE restart         ! only for lrst_oce
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
-   USE phycst
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE lib_mpp         ! MPP library
-   USE bdy_oce
-   USE bdy_par         
-   USE bdydyn2d        ! bdy_ssh routine
-   USE diaar5, ONLY:   lk_diaar5
-   USE iom
+   !!   ssh_nxt       : after ssh
+   !!   ssh_swp       : filter ans swap the ssh arrays
+   !!   wzv           : compute now vertical velocity
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables 
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE domvvl         ! Variable volume
+   USE divhor         ! horizontal divergence
+   USE phycst         ! physical constants
+   USE bdy_oce        ! 
+   USE bdy_par        !
+   USE bdydyn2d       ! bdy_ssh routine
 #if defined key_agrif
-   USE agrif_opa_update
    USE agrif_opa_interp
 #endif
 #if defined key_asminc   
-   USE asminc          ! Assimilation increment
-#endif
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
+   USE   asminc       ! Assimilation increment
+#endif
+   !
+   USE in_out_manager ! I/O manager
+   USE restart        ! only for lrst_oce
+   USE prtctl         ! Print control
+   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! Memory Allocation
+   USE timing         ! Timing
+   USE wet_dry         ! Wetting/Drying flux limting
 
    IMPLICIT NONE
@@ -51,5 +49,4 @@
 
    !! * Substitutions
-#  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -70,29 +67,26 @@
       !!      by the time step.
       !!
-      !! ** action  :   ssha    : after sea surface height
+      !! ** action  :   ssha, after sea surface height
       !!
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zhdiv
-      INTEGER, INTENT(in) ::   kt                      ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
       ! 
-      INTEGER             ::   jk                      ! dummy loop indice
-      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ssh_nxt')
-      !
-      CALL wrk_alloc( jpi, jpj, zhdiv ) 
+      INTEGER  ::   jk            ! dummy loop indice
+      REAL(wp) ::   z2dt, zcoef   ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhdiv   ! 2D workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('ssh_nxt')
+      !
+      CALL wrk_alloc( jpi,jpj,   zhdiv ) 
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
-         !
-      ENDIF
-      !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      ENDIF
+      !
+      CALL div_hor( kt )                              ! Horizontal divergence
       !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
@@ -104,5 +98,5 @@
       zhdiv(:,:) = 0._wp
       DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
-        zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
+        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
       END DO
       !                                                ! Sea surface elevation time stepping
@@ -110,34 +104,33 @@
       ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
       ! 
-      z1_rau0 = 0.5_wp * r1_rau0
-      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
-
-#if ! defined key_dynspg_ts
-      ! These lines are not necessary with time splitting since
-      ! boundary condition on sea level is set during ts loop
-#if defined key_agrif
-      CALL agrif_ssh( kt )
-#endif
-#if defined key_bdy
-      IF (lk_bdy) THEN
-         CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
-         CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
-      ENDIF
-#endif
-#endif
+      zcoef = 0.5_wp * r1_rau0
+
+      IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
+
+      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
+
+      IF ( .NOT.ln_dynspg_ts ) THEN
+         ! These lines are not necessary with time splitting since
+         ! boundary condition on sea level is set during ts loop
+# if defined key_agrif
+         CALL agrif_ssh( kt )
+# endif
+# if defined key_bdy
+         IF( lk_bdy ) THEN
+            CALL lbc_lnk( ssha, 'T', 1. )    ! Not sure that's necessary
+            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
+         ENDIF
+# endif
+      ENDIF
 
 #if defined key_asminc
-      !                                                ! Include the IAU weighted SSH increment
-      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
+      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN     ! Include the IAU weighted SSH increment
          CALL ssh_asm_inc( kt )
          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
       ENDIF
 #endif
-
       !                                           !------------------------------!
       !                                           !           outputs            !
       !                                           !------------------------------!
-      CALL iom_put( "ssh" , sshn                  )   ! sea surface height
-      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 )
@@ -165,17 +158,15 @@
       !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
       !!----------------------------------------------------------------------
-      !
-      INTEGER, INTENT(in) ::   kt           ! time step
+      INTEGER, INTENT(in) ::   kt   ! time step
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   z1_2dt       ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  z3d, zhdiv
-      !
-      INTEGER             ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)            ::   z1_2dt       ! local scalars
-      !!----------------------------------------------------------------------
-      
-      IF( nn_timing == 1 )  CALL timing_start('wzv')
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('wzv')
       !
       IF( kt == nit000 ) THEN
-         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
@@ -183,5 +174,4 @@
          !
          wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
-         !
       ENDIF
       !                                           !------------------------------!
@@ -199,5 +189,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
+                  zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
                END DO
             END DO
@@ -208,6 +198,6 @@
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)                    &
-               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
+               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
          END DO
          !          IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
@@ -216,11 +206,11 @@
          DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
             ! computation of w
-            wn(:,:,jk) = wn(:,:,jk+1) - (   fse3t_n(:,:,jk) * hdivn(:,:,jk)                                   &
-               &                          + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
+            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
+               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
          END DO
       ENDIF
 
 #if defined key_bdy
-      IF (lk_bdy) THEN
+      IF( lk_bdy ) THEN
          DO jk = 1, jpkm1
             wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
@@ -229,26 +219,8 @@
 #endif
       !
-      !                                           !------------------------------!
-      !                                           !           outputs            !
-      !                                           !------------------------------!
-      CALL iom_put( "woce", wn )   ! vertical velocity
-      IF( lk_diaar5 ) THEN                            ! vertical mass transport & its square value
-         CALL wrk_alloc( jpi, jpj, z2d ) 
-         CALL wrk_alloc( jpi, jpj, jpk, z3d ) 
-         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
-         z2d(:,:) = rau0 * e12t(:,:)
-         DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
-         END DO
-         CALL iom_put( "w_masstr" , z3d                     )  
-         CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
-         CALL wrk_dealloc( jpi, jpj, z2d  ) 
-         CALL wrk_dealloc( jpi, jpj, jpk, z3d ) 
-      ENDIF
-      !
       IF( nn_timing == 1 )  CALL timing_stop('wzv')
-
-
+      !
    END SUBROUTINE wzv
+
 
    SUBROUTINE ssh_swp( kt )
@@ -272,4 +244,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      REAL(wp) ::   zcoef   ! local scalar
       !!----------------------------------------------------------------------
       !
@@ -281,26 +255,25 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
       ENDIF
-
-# if defined key_dynspg_ts
-      IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN    !** Euler time-stepping: no filter
-# else
-      IF ( neuler == 0 .AND. kt == nit000 ) THEN   !** Euler time-stepping at first time-step : no filter
-#endif
-         sshb(:,:) = sshn(:,:)                           ! before <-- now
-         sshn(:,:) = ssha(:,:)                           ! now    <-- after  (before already = now)
-      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)
-         sshn(:,:) = ssha(:,:)                           ! now <-- after
-      ENDIF
-      !
-      ! Update velocity at AGRIF zoom boundaries
-#if defined key_agrif
-      IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
-#endif
+      !              !==  Euler time-stepping: no filter, just swap  ==!
+      IF(  ( neuler == 0 .AND. kt == nit000 ) .OR.    &
+         & ( ln_bt_fw    .AND. ln_dynspg_ts )      ) THEN 
+         sshb(:,:) = sshn(:,:)                              ! before <-- now
+         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
+         !
+      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
+         !                                                  ! before <-- now filtered
+         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
+         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
+            zcoef = atfp * rdt * r1_rau0
+            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
+               &                             -    rnf_b(:,:) + rnf   (:,:)   &
+               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
+         ENDIF
+         sshn(:,:) = ssha(:,:)                              ! now <-- after
+      ENDIF
       !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask, ovlap=1 )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('ssh_swp')
+      IF( nn_timing == 1 )   CALL timing_stop('ssh_swp')
       !
    END SUBROUTINE ssh_swp
