Index: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90
===================================================================
--- branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5017)
+++ branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/wadlmt.F90	(revision 5017)
@@ -0,0 +1,338 @@
+
+MODULE wadlmt
+   !!==============================================================================
+   !!                       ***  MODULE  wadlmt  ***
+   !! compute and apply flux limiters and preserve water depth positivity
+   !! only effects if wetting/drying is on (ln_wd == .true.
+   !!==============================================================================
+   !! History :   
+   !!  NEMO      3.6  ! 2014-09  ((H.Liu)  Original code
+   !!                 ! will add the runoff and periodic BC case later
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   wad_lmt    : Compute the horizontal flux limiter and the limited velocity
+   !!                when wetting and drying happens 
+   !!----------------------------------------------------------------------
+   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 sbcrnf          ! river runoff 
+   USE cla             ! cross land advection             (cla_div routine)
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! MPP library
+   USE wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+
+   PUBLIC   wad_lmt, wad_lmt_bt    ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+CONTAINS
+
+   SUBROUTINE wad_lmt( sshb1, sshemp, z2dt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   generate flux limiters for wetting/drying
+      !!
+      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif) 
+      !!
+      !! ** Action  : - calculate flux limiter and W/D flag
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   sshb1
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   sshemp
+      REAL(wp), INTENT(in) :: z2dt
+      !
+      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
+      INTEGER  ::   zflag               ! local scalar
+      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
+      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
+      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
+      REAL(wp) ::   ztmp                ! local scalars
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu,  zflxv            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
+
+      !!----------------------------------------------------------------------
+      !
+
+      IF( nn_timing == 1 )  CALL timing_start('wad_lmt')
+
+      IF(ln_wd) THEN
+
+        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
+        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        !
+       
+        IF(lwp) WRITE(numout,*)
+        IF(lwp) WRITE(numout,*) 'wad_lmt : wetting/drying limiters and velocity limiting'
+       
+        zflag  = 0
+        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
+
+       
+        zflxp(:,:)   = 0._wp
+        zflxn(:,:)   = 0._wp
+        zflxu(:,:)   = 0._wp
+        zflxv(:,:)   = 0._wp
+
+        zwdlmtu(:,:)  = 1._wp
+        zwdlmtv(:,:)  = 1._wp
+       
+        ! Horizontal Flux in u and v direction
+        DO jk = 1, jpkm1  
+           DO jj = 1, jpjm1
+              DO ji = 1, jpim1
+                 zflxu(ji,jj) = zflxu(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+                 zflxv(ji,jj) = zflxv(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+              END DO  
+           END DO  
+        END DO
+       
+        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+       
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+
+             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
+             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
+
+              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
+                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
+              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
+                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
+
+              zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1
+              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
+                !zdep2 = 0._wp
+                sshb1(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+              END IF
+           ENDDO
+        END DO
+
+      
+        !! start limiter iterations 
+        DO jk1 = 1, nn_wdit + 1
+       
+          
+           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+          
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+        
+                 wdmask(ji,jj) = 0
+                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
+                 IF(bathy(ji,jj) > zdepwd) CYCLE
+        
+                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
+                 ztmp = e12t(ji,jj)
+
+                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
+                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
+                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
+                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
+          
+                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+                 zdep2 = bathy(ji,jj) + sshb1(ji,jj) - rn_wdmin1 - z2dt * sshemp(ji,jj)  ! this one can be moved out of the loop
+          
+                 IF(zdep1 > zdep2) THEN
+                   zflag = 1
+                   wdmask(ji, jj) = 1
+                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                   zcoef = max(zcoef, 0._wp)
+                   IF(jk1 > nn_wdit) zcoef = 0._wp
+                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
+                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
+                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
+                 END IF
+              END DO ! ji loop
+           END DO  ! jj loop
+
+           CALL lbc_lnk( zwdlmtu, 'U', 1. )
+           CALL lbc_lnk( zwdlmtv, 'V', 1. )
+
+           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
+
+           IF(zflag == 0) EXIT
+          
+           zflag = 0     ! flag indicating if any further iteration is needed?
+        END DO  ! jk1 loop
+       
+        un(:,:,:) = un(:,:,:) * zwdlmtu(ji, jj) 
+        vn(:,:,:) = vn(:,:,:) * zwdlmtv(ji, jj) 
+
+        CALL lbc_lnk( un, 'U', -1. )
+        CALL lbc_lnk( vn, 'V', -1. )
+       
+        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt!!!'
+       
+        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+        !
+        !
+        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu, zflxv, zflxu1, zflxv1 )
+        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
+      !
+      END IF
+
+      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+   END SUBROUTINE wad_lmt
+
+   SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rdtbt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE wad_lmt  ***
+      !!                    
+      !! ** Purpose :   limiting flux in the barotropic stepping (dynspg_ts)
+      !!
+      !! ** Method  : - Prevent negative depth occurring (Not ready for Agrif) 
+      !!
+      !! ** Action  : - calculate flux limiter and W/D flag
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in)    ::   rdtbt    ! ocean time-step index
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zflxu,  zflxv, sshn_e, zssh_frc  
+      !
+      INTEGER  ::   ji, jj, jk, jk1     ! dummy loop indices
+      INTEGER  ::   zflag         ! local scalar
+      REAL(wp) ::   z2dt
+      REAL(wp) ::   zcoef, zdep1, zdep2 ! local scalars
+      REAL(wp) ::   zzflxp, zzflxn      ! local scalars
+      REAL(wp) ::   zdepwd              ! local scalar, always wet cell depth
+      REAL(wp) ::   ztmp                ! local scalars
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zwdlmtu, zwdlmtv         !: W/D flux limiters
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxp,  zflxn            ! local 2D workspace
+      REAL(wp), POINTER,  DIMENSION(:,:) ::   zflxu1, zflxv1           ! local 2D workspace
+
+      !!----------------------------------------------------------------------
+      !
+
+      IF( nn_timing == 1 )  CALL timing_start('wad_lmt_bt')
+
+      IF(ln_wd) THEN
+
+        CALL wrk_alloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
+        CALL wrk_alloc( jpi, jpj, zwdlmtu, zwdlmtv)
+        !
+       
+        IF(lwp) WRITE(numout,*)
+        IF(lwp) WRITE(numout,*) 'wad_lmt_bt : wetting/drying limiters and velocity limiting'
+       
+        zflag  = 0
+        zdepwd = 50._wp   !maximum depth that ocean cells can have W/D processes
+
+        z2dt = rdtbt   
+       
+        zflxp(:,:)   = 0._wp
+        zflxn(:,:)   = 0._wp
+        zflxu(:,:)   = 0._wp
+        zflxv(:,:)   = 0._wp
+
+        zwdlmtu(:,:)  = 1._wp
+        zwdlmtv(:,:)  = 1._wp
+       
+        ! Horizontal Flux in u and v direction
+       
+        zflxu(:,:) = zflxu(:,:) * e2u(:,:)
+        zflxv(:,:) = zflxv(:,:) * e1v(:,:)
+       
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+
+             IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE   ! we don't care about land cells
+             IF(bathy(ji,jj) > zdepwd) CYCLE       ! and cells which will unlikely go dried out
+
+              zflxp(ji,jj) = max(zflxu(ji,jj), 0._wp) - min(zflxu(ji-1,jj),   0._wp) + &
+                           & max(zflxv(ji,jj), 0._wp) - min(zflxv(ji,  jj-1), 0._wp) 
+              zflxn(ji,jj) = min(zflxu(ji,jj), 0._wp) - max(zflxu(ji-1,jj),   0._wp) + &
+                           & min(zflxv(ji,jj), 0._wp) - max(zflxv(ji,  jj-1), 0._wp) 
+
+              zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1
+              IF(zdep2 < 0._wp) THEN  !add more safty, but not necessary
+                !zdep2 = 0._wp
+               sshn_e(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+              END IF
+           ENDDO
+        END DO
+
+      
+        !! start limiter iterations 
+        DO jk1 = 1, nn_wdit + 1
+       
+          
+           zflxu1(:,:) = zflxu(:,:) * zwdlmtu(:,:)
+           zflxv1(:,:) = zflxv(:,:) * zwdlmtv(:,:)
+          
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+        
+                 wdmask(ji,jj) = 0
+                 IF(tmask(ji, jj, 1) < 0.5_wp) CYCLE 
+                 IF(bathy(ji,jj) > zdepwd) CYCLE
+        
+                 !ztmp = e1t(ji,jj) * e2t(ji,jj)     !there must be an array ready for this
+                 ztmp = e12t(ji,jj)
+
+                 zzflxp = max(zflxu1(ji,jj), 0._wp) - min(zflxu1(ji-1,jj),   0._wp) + &
+                        & max(zflxv1(ji,jj), 0._wp) - min(zflxv1(ji,  jj-1), 0._wp) 
+                 zzflxn = min(zflxu1(ji,jj), 0._wp) - max(zflxu1(ji-1,jj),   0._wp) + &
+                        & min(zflxv1(ji,jj), 0._wp) - max(zflxv1(ji,  jj-1), 0._wp) 
+          
+                 zdep1 = (zzflxp + zzflxn) * z2dt / ztmp
+                 zdep2 = bathy(ji,jj) + sshn_e(ji,jj) - rn_wdmin1   ! this one can be moved out of the loop
+                 zdep2 = zdep2 - z2dt * zssh_frc(ji,jj)
+          
+                 IF(zdep1 > zdep2) THEN
+                   zflag = 1
+                   !wdmask(ji, jj) = 1
+                   zcoef = ( ( zdep2 - rn_wdmin2 ) * ztmp - zzflxn * z2dt ) / ( zflxp(ji,jj) * z2dt )
+                   zcoef = max(zcoef, 0._wp)
+                   IF(jk1 > nn_wdit) zcoef = 0._wp
+                   IF(zflxu1(ji,  jj) > 0._wp) zwdlmtu(ji  ,jj) = zcoef
+                   IF(zflxu1(ji-1,jj) < 0._wp) zwdlmtu(ji-1,jj) = zcoef
+                   IF(zflxv1(ji,  jj) > 0._wp) zwdlmtv(ji  ,jj) = zcoef
+                   IF(zflxv1(ji,jj-1) < 0._wp) zwdlmtv(ji-1,jj) = zcoef
+                 END IF
+              END DO ! ji loop
+           END DO  ! jj loop
+
+           CALL lbc_lnk( zwdlmtu, 'U', 1. )
+           CALL lbc_lnk( zwdlmtv, 'V', 1. )
+
+           IF(lk_mpp) CALL mpp_max(zflag)   !max over the global domain
+
+           IF(zflag == 0) EXIT
+          
+           zflag = 0     ! flag indicating if any further iteration is needed?
+        END DO  ! jk1 loop
+       
+        zflxu(:,:) = zflxu(:,:) * zwdlmtu(ji, jj) 
+        zflxv(:,:) = zflxv(:,:) * zwdlmtv(ji, jj) 
+
+        CALL lbc_lnk( zflxu, 'U', -1. )
+        CALL lbc_lnk( zflxv, 'V', -1. )
+       
+        IF(zflag == 1 .AND. lwp) WRITE(numout,*) 'Need more iterations in wad_lmt_bt!!!'
+       
+        !IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field)
+        !IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field)
+        !
+        !
+        CALL wrk_dealloc( jpi, jpj, zflxp, zflxn, zflxu1, zflxv1 )
+        CALL wrk_dealloc( jpi, jpj, zwdlmtu, zwdlmtv)
+      !
+      END IF
+
+      IF( nn_timing == 1 )  CALL timing_stop('wad_lmt')
+   END SUBROUTINE wad_lmt_bt
+   !!======================================================================
+END MODULE wadlmt
