Changeset 5014

Timestamp:

2015-01-07T19:03:53+01:00 (9 years ago)

Author:

hliu

Message:

upload the modifications for W/D based on r:4826

Location:

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC

Files:

• DOM/dom_oce.F90 (modified) (5 diffs)
• DOM/domain.F90 (modified) (1 diff)
• DOM/domvvl.F90 (modified) (3 diffs)
• DOM/domzgr.F90 (modified) (21 diffs)
• DYN/dynhpg.F90 (modified) (17 diffs)
• DYN/dynspg_ts.F90 (modified) (62 diffs)
• DYN/sshwzv.F90 (modified) (4 diffs)
• nemogcm.F90 (modified) (3 diffs)
• step_oce.F90 (modified) (1 diff)

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

@@ -10 +10 @@
    !!            3.5  ! 2012     (S. Mocavero, I. Epicoco) Add arrays associated
    !!                             to the optimization of BDY communications
+   !!            3.6.?! 2014     (H. Liu) Add arrays associated with Wetting and Drying
    !!----------------------------------------------------------------------
 
@@ -257 +258 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol           !: north fold mask (jperio= 3 or 4)
+
 
 #if defined key_noslip_accurate
@@ -263 +265 @@
    INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: nicoa, njcoa       !: ???
 #endif
+
+   !!----------------------------------------------------------------------
+   !! critical depths,limiters,and masks for  Wetting and Drying
+   !! ---------------------------------------------------------------------
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wdmask             !: u- and v- limiter 
+
+   LOGICAL,  PUBLIC, SAVE ::   ln_wd       !: key to turn on/off wetting/drying (T: on, F: off)
+   REAL(wp), PUBLIC, SAVE ::   rn_wdmin1   !: minimum water depth on dried cells
+   REAL(wp), PUBLIC, SAVE ::   rn_wdmin2   !: tolerrance of minimum water depth on dried cells
+   REAL(wp), PUBLIC, SAVE ::   rn_wdld     !: land elevation below which wetting/drying will be considered
+   INTEGER , PUBLIC, SAVE ::   nn_wdit     !: maximum number of iteration for W/D limiter
+
+   !LOGICAL,  PUBLIC, SAVE ::   ln_wd     =  .FALSE.  !: turn on wetting/drying (T: on, F: off)
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin1 = 0.10_wp   !: minimum water depth on dried cells
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdmin2 = 0.01_wp   !: tolerrance of minimum water depth on dried cells
+   !REAL(wp), PUBLIC, SAVE ::   rn_wdld   = 20.0_wp   !: land elevation below which wetting/drying will be considered
+   !INTEGER , PUBLIC, SAVE ::   nn_wdit   =   10      !: maximum number of iteration for W/D limiter
 
    !!----------------------------------------------------------------------
@@ -326 +346 @@
    INTEGER FUNCTION dom_oce_alloc()
       !!----------------------------------------------------------------------
-      INTEGER, DIMENSION(11) :: ierr
+      INTEGER, DIMENSION(12) :: ierr
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -391 +411 @@
       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) )
 #endif
+
+      IF(ln_wd) &
+      ALLOCATE( wdmask(jpi,jpj), STAT=ierr(12) )
       !
       dom_oce_alloc = MAXVAL(ierr)

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

@@ -86 +86 @@
                              CALL dom_zgr      ! Vertical mesh and bathymetry
                              CALL dom_msk      ! Masks
-      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
+      IF( ln_sco.AND.(.NOT.ln_wd) ) &
+                         &   CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency
       !
       ht_0(:,:) = 0.0_wp                       ! Reference ocean depth at T-points

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

@@ -145 +145 @@
       fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
 
+      IF(ln_wd) THEN
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            IF(e3t_0(ji,jj,1) < 0._wp) fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+          ENDDO
+        ENDDO
+      END IF
+
       ! Reconstruction of all vertical scale factors at now and before time steps
       ! =============================================================================
@@ -793 +801 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !! * Local declarations
-      INTEGER ::   jk
+      INTEGER ::   ji, jj, jk
       INTEGER ::   id1, id2, id3, id4, id5     ! local integers
       !!----------------------------------------------------------------------
@@ -867 +875 @@
             fse3t_n(:,:,:) = e3t_0(:,:,:)
             sshn(:,:) = 0.0_wp
+
+            IF(ln_wd) THEN
+              DO jj = 1, jpj
+                DO ji = 1, jpi
+                  IF(e3t_0(ji,jj,1) < 0._wp) THEN
+                    fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+                    fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1 
+                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)
+                  ENDIF
+                ENDDO
+              ENDDO
+            END IF
+
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
                tilde_e3t_b(:,:,:) = 0.0_wp

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

@@ -17 +17 @@
    !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case  
+   !!            3.6  ! 2014-09  (H. Liu) Modifications for Wetting/Drying
    !!----------------------------------------------------------------------
 
@@ -451 +452 @@
          END DO
          !
-         !                                            ! ================ !
       ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain)
          !                                            ! ================ !
@@ -1121 +1121 @@
       REAL(wp) ::   zrmax, ztaper   ! temporary scalars
       REAL(wp) ::   zrfact
+      REAL(wp) ::   zsmth
       !
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
@@ -1176 +1177 @@
       bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
 
+      IF(.NOT.ln_wd) THEN                      
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -1181 +1183 @@
          END DO
       END DO
+      END IF
       !                                        ! =============================
       !                                        ! Define the envelop bathymetry   (hbatt)
@@ -1187 +1190 @@
       zenv(:,:) = bathy(:,:)
       !
+      IF(.NOT.ln_wd) THEN    
       ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
       DO jj = 1, jpj
@@ -1202 +1206 @@
          END DO
       END DO
+      END IF
+
       ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
       CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' )
@@ -1241 +1247 @@
          zri(:,:) = 0._wp
          zrj(:,:) = 0._wp
+         
+         !IF(ln_wd) THEN                     !extend the smoothed region to cover the W/D zones
+         !  zsmth = -rn_wdld    
+         !ELSE
+           zsmth = 0._wp                     ! The original form (only smooth ocean points)
+         !ENDIF
+
          DO jj = 1, nlcj
             DO ji = 1, nlci
                iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
                ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
+               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(iip1,jj) > zsmth)) THEN
                   zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
                END IF
-               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
+               IF( (zenv(ji,jj) > zsmth) .AND. (zenv(ji,ijp1) > zsmth)) THEN
                   zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
                END IF
@@ -1257 +1270 @@
             END DO
          END DO
+
          IF( lk_mpp )   CALL mpp_max( zrmax )   ! max over the global domain
          !
@@ -1271 +1285 @@
       END DO                                                !     End loop     !
       !                                                     ! ================ !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
-         END DO
-      END DO
+      IF(ln_wd) THEN
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              zenv(ji,jj) = MAX( zenv(ji,jj), -rn_wdld )    ! filt out land bathy data 
+           END DO
+        END DO
+      ELSE
+        DO jj = 1, jpj
+           DO ji = 1, jpi
+              zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
+           END DO
+        END DO
+      ENDIF
       !
       ! Envelope bathymetry saved in hbatt
@@ -1305 +1327 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         IF(.NOT.ln_wd) THEN
+           WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
+         ELSE
+           WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min
+           WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld
+         ENDIF
       ENDIF
       hbatu(:,:) = rn_sbot_min
@@ -1318 +1345 @@
         END DO
       END DO
+
+      IF(ln_wd) THEN               !avoid the zero depth on T- (U-,V-,F-) points
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            IF(ABS(hbatt(ji,jj)) < rn_wdmin1) &
+              & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatu(ji,jj)) < rn_wdmin1) &
+              & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatv(ji,jj)) < rn_wdmin1) &
+              & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1
+            IF(ABS(hbatf(ji,jj)) < rn_wdmin1) &
+              & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1
+          END DO
+        END DO
+      END IF
       ! 
       ! Apply lateral boundary condition
 !!gm  ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
+
       zhbat(:,:) = hbatu(:,:)   ;   CALL lbc_lnk( hbatu, 'U', 1._wp )
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF( hbatu(ji,jj) == 0._wp ) THEN
+               !No worries about the following line when ln_wd == .true.
                IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
                IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
@@ -1330 +1374 @@
          END DO
       END DO
+
       zhbat(:,:) = hbatv(:,:)   ;   CALL lbc_lnk( hbatv, 'V', 1._wp )
       DO jj = 1, jpj
@@ -1339 +1384 @@
          END DO
       END DO
+
       zhbat(:,:) = hbatf(:,:)   ;   CALL lbc_lnk( hbatf, 'F', 1._wp )
       DO jj = 1, jpj
@@ -1350 +1396 @@
 
 !!bug:  key_helsinki a verifer
-      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
-      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
-      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
-      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      IF(.NOT.ln_wd) THEN
+       hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
+       hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
+       hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
+       hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
+      END IF
 
       IF( nprint == 1 .AND. lwp )   THEN
@@ -1395 +1443 @@
       CALL lbc_lnk( e3vw_0, 'V', 1._wp )
 
-      fsdepw(:,:,:) = gdepw_0 (:,:,:)
-      fsde3w(:,:,:) = gdep3w_0(:,:,:)
-      !
-      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
-      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
-      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
-      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
-      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
-      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
-      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      !fsdepw(:,:,:) = gdepw_0 (:,:,:)
+      !fsde3w(:,:,:) = gdep3w_0(:,:,:)
+      !
+      IF(.NOT.ln_wd) THEN
+        where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
+        where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
+        where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
+        where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
+        where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
+        where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
+        where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
+      END IF
 
 #if defined key_agrif
@@ -1446 +1496 @@
                IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
             END DO
-            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
-         END DO
-      END DO
+
+            IF(ln_wd) THEN
+              IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN
+                mbathy(ji,jj) = 0
+              ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN
+                mbathy(ji,jj) = 2
+              ENDIF
+            ELSE
+              IF( scobot(ji,jj) == 0._wp )   mbathy(ji,jj) = 0
+            ENDIF
+         END DO
+      END DO
+
       IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )
@@ -1568 +1628 @@
       INTEGER  ::   ji, jj, jk           ! dummy loop argument
       REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
+      REAL(wp) ::   ztmpu, ztmpv, ztmpf
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -1621 +1682 @@
       DO ji = 1, jpim1
          DO jj = 1, jpjm1
+            ! extended for Wetting/Drying case
+            ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
+            ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
+            ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
             DO jk = 1, jpk
-               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
-                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
-                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               IF(ABS(ztmpu) < rn_wdmin1.AND.ln_wd) THEN
+                 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+               ELSE
+                 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+                 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                        &              / ztmpu
+               END IF
+
+               IF(ABS(ztmpv) < rn_wdmin1.AND.ln_wd) THEN
+                 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+               ELSE
+                 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+                 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                        &              / ztmpv
+               END IF
+
+               IF(ABS(ztmpf) < rn_wdmin1.AND.ln_wd) THEN
+                 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
+                    &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+               ELSE
+                 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  &
+                    &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
+                    &              / ztmpf
+               END IF
+
+               !
+               !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+               !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
+               !   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+               !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+               !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
+               !   &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
                !
                e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
@@ -1672 +1767 @@
       REAL(wp) ::   zsmth               ! smoothing around critical depth
       REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
+      REAL(wp) ::   ztmpu, ztmpv, ztmpf
       !
       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
@@ -1750 +1846 @@
         DO jj=1,jpj-1
 
-          DO jk = 1, jpk
-                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
-                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
-                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
-                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
-                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
-                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+           ! extend to suit for Wetting/Drying case
+           ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
+           ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
+           ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
+           DO jk = 1, jpk
+              IF(ABS(ztmpu) < 1.e-5.AND.ln_wd) THEN
+                z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) )
+                z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) )
+              ELSE
+                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
+                       &              / ztmpu
+              END IF
+
+              IF(ABS(ztmpv) < 1.e-5.AND.ln_wd) THEN
+                z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) )
+                z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) )
+              ELSE
+                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
+                       &              / ztmpv
+              END IF
+
+              IF(ABS(ztmpf) < 1.e-5.AND.ln_wd) THEN
+                z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk)   + z_esigt3(ji+1,jj,jk)  &
+                   &                            + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) )
+              ELSE
+                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)  &
+                   &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
+                   &              / ztmpf
+              END IF
+
+             !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
+             !z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
+             !                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
+             !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
+             !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
+             !                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
 
              e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

@@ -16 +16 @@
    !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
+   !!            3.6? !  2014-09  (H. Liu) add Wetting/Drying pressure filter
    !!----------------------------------------------------------------------
 
@@ -369 +370 @@
       !!
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
+      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp   ! temporary scalars
+      LOGICAL  ::   ll_tmp1, ll_tmp2           ! local logical variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -386 +390 @@
       IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume
       ELSE                 ;     znad = 0._wp         ! Fixed volume
+      ENDIF
+
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+             ELSE
+               zcpx(ji,jj) = 0._wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpy(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+             ELSE
+               zcpy(ji,jj) = 0._wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
       ENDIF
 
@@ -401 +440 @@
             zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
                &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
+
+            IF(ln_wd) THEN
+              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
+              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+              zuap = zuap * zcpx(ji,jj)
+              zvap = zvap * zcpy(ji,jj)
+            ENDIF
+
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
@@ -423 +470 @@
                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)
+
+               IF(ln_wd) THEN
+                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
+                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+                 zuap = zuap * zcpx(ji,jj)
+                 zvap = zvap * zcpy(ji,jj)
+               ENDIF
+
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
@@ -431 +486 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_sco
@@ -448 +504 @@
       REAL(wp) ::   z1_10, cffu, cffx   !    "         "
       REAL(wp) ::   z1_12, cffv, cffy   !    "         "
+      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  dzx, dzy, dzz, dzu, dzv, dzw
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  drhox, drhoy, drhoz, drhou, drhov, drhow
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  rho_i, rho_j, rho_k
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
       !!----------------------------------------------------------------------
       !
@@ -457 +515 @@
       CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
       CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
-      !
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
+      !
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+             ELSE
+               zcpx(ji,jj) = 0._wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpy(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+             ELSE
+               zcpy(ji,jj) = 0._wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
+
 
       IF( kt == nit000 ) THEN
@@ -628 +725 @@
             zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
             zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
+            IF(ln_wd) THEN
+              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
+              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
+            ENDIF
             ! add to the general momentum trend
             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
@@ -647 +748 @@
                   &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
                   &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) / e2v(ji,jj)
+               IF(ln_wd) THEN
+                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
+                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
+               ENDIF
                ! add to the general momentum trend
                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
@@ -657 +762 @@
       CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
       CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k,  zhpi,  zhpj        )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_djc
@@ -682 +788 @@
       !! The local variables for the correction term
       INTEGER  :: jk1, jis, jid, jjs, jjd
+      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables
       REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
       REAL(wp) :: zrhdt1
@@ -687 +794 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
-      !!----------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zcpx, zcpy    !W/D pressure filter
+      !!----------------------------------------------------------------------
+      !
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
+      IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
       !
       IF( kt == nit000 ) THEN
@@ -703 +813 @@
       znad = 0.0_wp
       IF( lk_vvl ) znad = 1._wp
+
+      IF(ln_wd) THEN
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1 
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1) THEN
+               zcpx(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji+1,jj) - sshn(ji,jj)))
+             ELSE
+               zcpx(ji,jj) = 0._wp
+             END IF
+      
+             ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+             ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
+                                                       & rn_wdmin1 + rn_wdmin2
+
+             IF(ll_tmp1.OR.ll_tmp2) THEN
+               zcpy(ji,jj) = 1.0_wp
+             ELSE IF(ll_tmp2) THEN
+               ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
+                           &     (sshn(ji,jj+1) - sshn(ji,jj)))
+             ELSE
+               zcpy(ji,jj) = 0._wp
+             END IF
+           END DO
+        END DO
+        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+      ENDIF
 
       ! Clean 3-D work arrays
@@ -862 +1007 @@
                ENDIF
 
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * &
-               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
+               IF(ln_wd) THEN
+                  zdpdx1 = zdpdx1 * zcpx(ji,jj)
+                  zdpdx2 = zdpdx2 * zcpx(ji,jj)
+                ENDIF
+                ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
             ENDIF
 
@@ -919 +1067 @@
                ENDIF
 
-               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*&
-               &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
+               IF(ln_wd) THEN
+                  zdpdy1 = zdpdy1 * zcpy(ji,jj)
+                  zdpdy2 = zdpdy2 * zcpy(ji,jj)
+                ENDIF
+
+               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
             ENDIF
 
@@ -930 +1082 @@
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
       CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
+      IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
       !
    END SUBROUTINE hpg_prj

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -11 +11 @@
    !!             3.5  ! 2013-07  (J. Chanut) Switch to Forward-backward time stepping
    !!             3.6  ! 2013-11  (A. Coward) Update for z-tilde compatibility
+   !!             3.?  ! 2014-09  (H. Liu) Add Wetting/Drying components
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -17 +18 @@
    !!----------------------------------------------------------------------
    !!   dyn_spg_ts  : compute surface pressure gradient trend using a time-
-   !!                 splitting scheme and add to the general trend 
+   !!                 splitting scheme and add to the general trend
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -26 +27 @@
    USE dynvor          ! vorticity term
    USE bdy_par         ! for lk_bdy
-   USE bdytides        ! open boundary condition data     
+   USE bdytides        ! open boundary condition data
    USE bdydyn2d        ! open boundary conditions on barotropic variables
    USE sbctide         ! tides
@@ -38 +39 @@
    USE zdf_oce         ! Bottom friction coefts
    USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing    
+   USE timing          ! Timing
    USE sbcapr          ! surface boundary condition: atmospheric pressure
+   USE wadlmt          ! wetting/drying flux limter
    USE dynadv, ONLY: ln_dynadv_vec
 #if defined key_agrif
    USE agrif_opa_interp ! agrif
 #endif
-#if defined key_asminc   
+#if defined key_asminc
    USE asminc          ! Assimilation increment
 #endif
@@ -51 +53 @@
    PRIVATE
 
-   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90 
+   PUBLIC dyn_spg_ts        ! routine called in dynspg.F90
    PUBLIC dyn_spg_ts_alloc  !    "      "     "    "
    PUBLIC dyn_spg_ts_init   !    "      "     "    "
@@ -97 +99 @@
       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) , & 
+      IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &
                              &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )
 
@@ -110 +112 @@
       !!----------------------------------------------------------------------
       !!
-      !! ** Purpose :   
+      !! ** Purpose :
       !!      -Compute the now trend due to the explicit time stepping
-      !!      of the quasi-linear barotropic system. 
+      !!      of the quasi-linear barotropic system.
       !!
-      !! ** Method  :  
+      !! ** Method  :
       !!      Barotropic variables are advanced from internal time steps
       !!      "n"   to "n+1" if ln_bt_fw=T
-      !!      or from 
+      !!      or from
       !!      "n-1" to "n+1" if ln_bt_fw=F
       !!      thanks to a generalized forward-backward time stepping (see ref. below).
@@ -126 +128 @@
       !!      -Compute barotropic advective velocities at step "n" : un_adv, vn_adv
       !!      These are used to advect tracers and are compliant with discrete
-      !!      continuity equation taken at the baroclinic time steps. This 
+      !!      continuity equation taken at the baroclinic time steps. This
       !!      ensures tracers conservation.
       !!      -Update 3d trend (ua, va) with barotropic component.
       !!
-      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 
-      !!              The regional oceanic modeling system (ROMS): 
+      !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005:
+      !!              The regional oceanic modeling system (ROMS):
       !!              a split-explicit, free-surface,
-      !!              topography-following-coordinate oceanic model. 
-      !!              Ocean Modelling, 9, 347-404. 
+      !!              topography-following-coordinate oceanic model.
+      !!              Ocean Modelling, 9, 347-404.
       !!---------------------------------------------------------------------
       !
       INTEGER, INTENT(in)  ::   kt   ! ocean time-step index
       !
-      LOGICAL  ::   ll_fw_start        ! if true, forward integration 
-      LOGICAL  ::   ll_init             ! if true, special startup of 2d equations
-      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices
-      INTEGER  ::   ikbu, ikbv, noffset      ! local integers
-      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
+      LOGICAL  ::   ll_fw_start                 ! if true, forward integration
+      LOGICAL  ::   ll_init                         ! if true, special startup of 2d equations
+      LOGICAL  ::   ll_tmp1, ll_tmp2            ! local logical variables used in W/D
+      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
+      INTEGER  ::   ikbu, ikbv, noffset         ! local integers
+      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars
       REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      -
-      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2    !   -      -
+      REAL(wp) ::   z1_12, z1_8, z1_4, z1_2       !   -      -
       REAL(wp) ::   zu_spg, zv_spg             !   -      -
-      REAL(wp) ::   zhura, zhvra               !   -      -
-      REAL(wp) ::   za0, za1, za2, za3           !   -      -
+      REAL(wp) ::   zhura, zhvra                     !   -      -
+      REAL(wp) ::   za0, za1, za2, za3                 !   -      -
+
+      REAL(wp) ::   ztmp                       ! temporary vaialbe
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e
@@ -156 +161 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a
       REAL(wp), POINTER, DIMENSION(:,:) :: zhf
+      REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
       !!----------------------------------------------------------------------
       !
@@ -167 +173 @@
       CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a                                   )
       CALL wrk_alloc( jpi, jpj, zhf )
+      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
       !
       !                                         !* Local constant initialization
-      z1_12 = 1._wp / 12._wp 
-      z1_8  = 0.125_wp                                   
+      z1_12 = 1._wp / 12._wp
+      z1_8  = 0.125_wp
       z1_4  = 0.25_wp
-      z1_2  = 0.5_wp     
+      z1_2  = 0.5_wp
       zraur = 1._wp / rau0
       !
-      IF( kt == nit000 .AND. neuler == 0 ) THEN    ! reciprocal of baroclinic time step 
+      IF( kt == nit000 .AND. neuler == 0 ) THEN         ! reciprocal of baroclinic time step
         z2dt_bf = rdt
       ELSE
         z2dt_bf = 2.0_wp * rdt
       ENDIF
-      z1_2dt_b = 1.0_wp / z2dt_bf 
-      !
-      ll_init = ln_bt_av                           ! if no time averaging, then no specific restart 
+      z1_2dt_b = 1.0_wp / z2dt_bf
+      !
+      ll_init = ln_bt_av                                ! if no time averaging, then no specific restart
       ll_fw_start = .FALSE.
       !
-                                                       ! time offset in steps for bdy data update
+                                                        ! time offset in steps for bdy data update
       IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ;  noffset = 0 ; ENDIF
       !
-      IF( kt == nit000 ) THEN                !* initialisation
+      IF( kt == nit000 ) THEN                   !* initialisation
          !
          IF(lwp) WRITE(numout,*)
@@ -245 +252 @@
             IF ( .not. ln_sco ) THEN
 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case
-!     Set it to zero for the time being 
+!     Set it to zero for the time being
 !              IF( rn_hmin < 0._wp ) THEN    ;   jk = - INT( rn_hmin )                                      ! from a nb of level
 !              ELSE                          ;   jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 )  ! from a depth
@@ -264 +271 @@
             END DO
             CALL lbc_lnk( zhf, 'F', 1._wp )
-            ! JC: TBC. hf should be greater than 0 
+            ! JC: TBC. hf should be greater than 0
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -274 +281 @@
       ENDIF
       !
-      ! If forward start at previous time step, and centered integration, 
+      ! If forward start at previous time step, and centered integration,
       ! then update averaging weights:
       IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN
@@ -280 +287 @@
          CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2)
       ENDIF
-                          
+
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
       ! -----------------------------------------------------------------------------
-      !      
+      !
       !
       !                                   !* e3*d/dt(Ua) (Vertically integrated)
@@ -293 +300 @@
       DO jk = 1, jpkm1
          zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)         
+         zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
       END DO
       !
@@ -311 +318 @@
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
-      zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes 
+      zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes
       zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
       !
@@ -353 +360 @@
          END DO
          !
-      ENDIF 
+      ENDIF
       !
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
       IF( lk_vvl ) THEN                         ! Variable volume : remove surface pressure gradient
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
-               zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
-            END DO
-         END DO
+        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+          DO jj = 2, jpjm1
+             DO ji = 2, jpim1
+                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
+                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
+                                                          & rn_wdmin1 + rn_wdmin2
+                IF(ll_tmp1) THEN
+                  zcpx(ji,jj) = 1.0_wp
+                ELSE IF(ll_tmp2) THEN
+                   ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
+                  zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
+                              &     (sshn(ji+1,jj) - sshn(ji,jj)))
+                ELSE
+                   zcpx(ji,jj) = 0._wp
+                END IF
+
+                ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+                ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
+                                                          & rn_wdmin1 + rn_wdmin2
+                IF(ll_tmp1) THEN
+                   zcpy(ji,jj) = 1.0_wp
+                ELSE IF(ll_tmp2) THEN
+                   ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
+                  zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /&
+                              &     (sshn(ji,jj+1) - sshn(ji,jj)))
+                ELSE
+                  zcpy(ji,jj) = 0._wp
+                END IF
+             END DO
+           END DO
+           CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+        ENDIF
+
+        IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1
+                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj  ) - sshn(ji  ,jj ) ) / &
+                               &                 e1u(ji,jj) * zcpx(ji,jj)
+                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji  ,jj+1) - sshn(ji  ,jj ) ) / &
+                               &                 e2v(ji,jj) * zcpy(ji,jj)
+              END DO
+           END DO
+         ELSE
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) / e1u(ji,jj)
+                 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) / e2v(ji,jj)
+              END DO
+           END DO
+        END IF
+
       ENDIF
 
@@ -371 +423 @@
              zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)
           END DO
-      END DO 
-      !
-      !                 ! Add bottom stress contribution from baroclinic velocities:      
+      END DO
+      !
+      !                                         ! Add bottom stress contribution from baroclinic velocities:
       IF (ln_bt_fw) THEN
-         DO jj = 2, jpjm1                          
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikbu = mbku(ji,jj)       
-               ikbv = mbkv(ji,jj)    
+               ikbu = mbku(ji,jj)
+               ikbv = mbkv(ji,jj)
                zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities
                zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj)
@@ -386 +438 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikbu = mbku(ji,jj)       
-               ikbv = mbkv(ji,jj)    
+               ikbu = mbku(ji,jj)
+               ikbv = mbkv(ji,jj)
                zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities
                zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj)
@@ -397 +449 @@
       zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)
       zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)
-      !       
+      !
       IF (ln_bt_fw) THEN                        ! Add wind forcing
          zu_frc(:,:) =  zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)
@@ -404 +456 @@
          zu_frc(:,:) =  zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)
          zv_frc(:,:) =  zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)
-      ENDIF  
+      ENDIF
       !
       IF ( ln_apr_dyn ) THEN                    ! Add atm pressure forcing
          IF (ln_bt_fw) THEN
-            DO jj = 2, jpjm1              
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zu_spg =  grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) /e1u(ji,jj)
@@ -417 +469 @@
             END DO
          ELSE
-            DO jj = 2, jpjm1              
+            DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   zu_spg =  grav * z1_2 * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)    &
@@ -427 +479 @@
                END DO
             END DO
-         ENDIF 
+         ENDIF
       ENDIF
       !                                   !* Right-Hand-Side of the barotropic ssh equation
@@ -450 +502 @@
       !
       ! -----------------------------------------------------------------------
-      !  Phase 2 : Integration of the barotropic equations 
+      !  Phase 2 : Integration of the barotropic equations
       ! -----------------------------------------------------------------------
       !
       !                                             ! ==================== !
       !                                             !    Initialisations   !
-      !                                             ! ==================== !  
-      ! Initialize barotropic variables:      
+      !                                             ! ==================== !
+      ! Initialize barotropic variables:
       IF( ll_init )THEN
          sshbb_e(:,:) = 0._wp
@@ -465 +517 @@
          vb_e   (:,:) = 0._wp
       ENDIF
-      !
-      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                    
-         sshn_e(:,:) = sshn (:,:)            
-         zun_e (:,:) = un_b (:,:)            
+
+      IF(ln_wd) THEN      !preserve the positivity of water depth
+                          !ssh[b,n,a] should have already been processed for this
+         sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:))
+         sshb_e(:,:)  = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:))
+      ENDIF
+      !
+      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields
+         sshn_e(:,:) = sshn (:,:)
+         zun_e (:,:) = un_b (:,:)
          zvn_e (:,:) = vn_b (:,:)
          !
-         hu_e  (:,:) = hu   (:,:)       
-         hv_e  (:,:) = hv   (:,:) 
-         hur_e (:,:) = hur  (:,:)    
+         hu_e  (:,:) = hu   (:,:)
+         hv_e  (:,:) = hv   (:,:)
+         hur_e (:,:) = hur  (:,:)
          hvr_e (:,:) = hvr  (:,:)
       ELSE                                ! CENTRED integration: start from BEFORE fields
          sshn_e(:,:) = sshb (:,:)
-         zun_e (:,:) = ub_b (:,:)         
+         zun_e (:,:) = ub_b (:,:)
          zvn_e (:,:) = vb_b (:,:)
          !
-         hu_e  (:,:) = hu_b (:,:)       
-         hv_e  (:,:) = hv_b (:,:) 
-         hur_e (:,:) = hur_b(:,:)    
+         hu_e  (:,:) = hu_b (:,:)
+         hv_e  (:,:) = hv_b (:,:)
+         hur_e (:,:) = hur_b(:,:)
          hvr_e (:,:) = hvr_b(:,:)
       ENDIF
@@ -489 +547 @@
       !
       ! Initialize sums:
-      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
+      ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)
       va_b  (:,:) = 0._wp
       ssha  (:,:) = 0._wp       ! Sum for after averaged sea level
@@ -506 +564 @@
          !
          ! Set extrapolation coefficients for predictor step:
-         IF ((jn<3).AND.ll_init) THEN      ! Forward           
-           za1 = 1._wp                                          
-           za2 = 0._wp                        
-           za3 = 0._wp                        
-         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105 
+         IF ((jn<3).AND.ll_init) THEN      ! Forward
+           za1 = 1._wp
+           za2 = 0._wp
+           za3 = 0._wp
+         ELSE                              ! AB3-AM4 Coefficients: bet=0.281105
            za1 =  1.781105_wp              ! za1 =   3/2 +   bet
            za2 = -1.06221_wp               ! za2 = -(1/2 + 2*bet)
@@ -524 +582 @@
             ! Extrapolate Sea Level at step jit+0.5:
             zsshp2_e(:,:) = za1 * sshn_e(:,:)  + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)
+
+
             !
             DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points
@@ -535 +595 @@
                END DO
             END DO
+
             CALL lbc_lnk( zwx, 'U', 1._wp )    ;   CALL lbc_lnk( zwy, 'V', 1._wp )
             !
             zhup2_e (:,:) = hu_0(:,:) + zwx(:,:)                ! Ocean depth at U- and V-points
             zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:)
+            IF(ln_wd) THEN
+              zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1)
+              zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1)
+            END IF
          ELSE
             zhup2_e (:,:) = hu(:,:)
@@ -552 +617 @@
          !
 #if defined key_agrif
-         ! Set fluxes during predictor step to ensure 
+         ! Set fluxes during predictor step to ensure
          ! volume conservation
          IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN
@@ -577 +642 @@
          ENDIF
 #endif
+         IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
          !
          ! Sum over sub-time-steps to compute advective velocities
@@ -584 +650 @@
          !
          ! Set next sea level:
-         DO jj = 2, jpjm1                                 
+         DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zhdiv(ji,jj) = (   zwx(ji,jj) - zwx(ji-1,jj)   &
@@ -590 +656 @@
             END DO
          END DO
+         !IF(ln_wd) THEN
+
          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+         IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
@@ -600 +669 @@
          IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )
 #endif
-         !  
+         !
          ! Sea Surface Height at u-,v-points (vvl case only)
-         IF ( lk_vvl ) THEN                                
+         IF ( lk_vvl ) THEN
             DO jj = 2, jpjm1
                DO ji = 2, jpim1      ! NO Vector Opt.
@@ -613 +682 @@
                END DO
             END DO
+
             CALL lbc_lnk( zsshu_a, 'U', 1._wp )   ;   CALL lbc_lnk( zsshv_a, 'V', 1._wp )
-         ENDIF   
-         !                                 
+         ENDIF
+         !
          ! Half-step back interpolation of SSH for surface pressure computation:
          !----------------------------------------------------------------------
          IF ((jn==1).AND.ll_init) THEN
            za0=1._wp                        ! Forward-backward
-           za1=0._wp                           
+           za1=0._wp
            za2=0._wp
            za3=0._wp
@@ -627 +697 @@
            za1=-0.1666666666666_wp          ! za1 = gam
            za2= 0.0833333333333_wp          ! za2 = eps
-           za3= 0._wp              
-         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
-           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps    
-           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps 
+           za3= 0._wp
+         ELSE                               ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880
+           za0=0.614_wp                     ! za0 = 1/2 +   gam + 2*eps
+           za1=0.285_wp                     ! za1 = 1/2 - 2*gam - 3*eps
            za2=0.088_wp                     ! za2 = gam
            za3=0.013_wp                     ! za3 = eps
@@ -640 +710 @@
          !
          ! Compute associated depths at U and V points:
-         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN      
-            !                                        
-            DO jj = 2, jpjm1                            
+         IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN
+            !
+            DO jj = 2, jpjm1
                DO ji = 2, jpim1
                   zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e12u(ji  ,jj)    &
@@ -650 +720 @@
                      &       * ( e12t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  &
                      &       +   e12t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) )
-                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1 
+                  zhust_e(ji,jj) = hu_0(ji,jj) + zx1
                   zhvst_e(ji,jj) = hv_0(ji,jj) + zy1
                END DO
             END DO
+
+            IF(ln_wd) THEN
+              zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 )
+              zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 )
+            END IF
+            !shall we call lbc_lnk for zhu(v)st_e() here?
+
          ENDIF
          !
@@ -692 +769 @@
                   zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) &
                      &                                    + ftnw(ji+1,jj) * zwy(ji+1,jj  ) &
-                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) & 
+                     &                                    + ftse(ji,jj  ) * zwy(ji  ,jj-1) &
                      &                                    + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
-                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 
+                  zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) &
                      &                                    + ftse(ji,jj+1) * zwx(ji  ,jj+1) &
-                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) & 
+                     &                                    + ftnw(ji,jj  ) * zwx(ji-1,jj  ) &
                      &                                    + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
                END DO
             END DO
-            ! 
+            !
          ENDIF
          !
@@ -720 +797 @@
          !
          ! Surface pressure trend:
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! Add surface pressure gradient
-               zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
-               zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
-               zwx(ji,jj) = zu_spg
-               zwy(ji,jj) = zv_spg
-            END DO
-         END DO
+         IF(ln_wd) THEN                    ! Calculating and applying W/D gravity filters
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1
+                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
+                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
+                                                           & rn_wdmin1 + rn_wdmin2
+                 IF(ll_tmp1) THEN
+                   zcpx(ji,jj) = 1.0_wp
+                 ELSE IF(ll_tmp2) THEN
+                    ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here
+                   zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
+                               & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)))
+                 ELSE
+                    zcpx(ji,jj) = 0._wp
+                 END IF
+
+                 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
+                 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
+                                                           & rn_wdmin1 + rn_wdmin2
+                 IF(ll_tmp1) THEN
+                    zcpy(ji,jj) = 1.0_wp
+                 ELSE IF(ll_tmp2) THEN
+                    ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here
+                   zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /&
+                               & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)))
+                 ELSE
+                   zcpy(ji,jj) = 0._wp
+                 END IF
+              END DO
+            END DO
+            CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
+         ENDIF
+
+         IF(ln_wd) THEN
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1 
+                 ! Add surface pressure gradient
+                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+                 zwx(ji,jj) = zu_spg * zcpx(ji,jj)
+                 zwy(ji,jj) = zv_spg * zcpy(ji,jj)
+              END DO
+           END DO
+         ELSE
+           DO jj = 2, jpjm1
+              DO ji = fs_2, fs_jpim1   ! vector opt.
+                 ! Add surface pressure gradient
+                 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj)
+                 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj)
+                 zwx(ji,jj) = zu_spg
+                 zwy(ji,jj) = zv_spg
+              END DO
+           END DO
+         END IF
+
          !
          ! Set next velocities:
-         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN    ! Vector form
+         IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN   ! Vector form
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua_e(ji,jj) = (                                zun_e(ji,jj)   & 
+                  ua_e(ji,jj) = (                                zun_e(ji,jj)   &
                             &     + rdtbt * (                      zwx(ji,jj)   &
                             &                                 + zu_trd(ji,jj)   &
-                            &                                 + zu_frc(ji,jj) ) & 
+                            &                                 + zu_frc(ji,jj) ) &
                             &   ) * umask(ji,jj,1)
 
@@ -748 +871 @@
             END DO
 
-         ELSE                 ! Flux form
+         ELSE                                           ! Flux form
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
 
-                  zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))
-                  zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))
-
-                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   & 
-                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   & 
+                  IF(ln_wd) THEN
+                    zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1)
+                    zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1)
+                  ELSE
+                    zhura = hu_0(ji,jj) + zsshu_a(ji,jj)
+                    zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)
+                  END IF
+
+                  zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1))
+                  zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1))
+
+                  ua_e(ji,jj) = (                hu_e(ji,jj)  *  zun_e(ji,jj)   &
+                            &     + rdtbt * ( zhust_e(ji,jj)  *    zwx(ji,jj)   &
                             &               + zhup2_e(ji,jj)  * zu_trd(ji,jj)   &
                             &               +      hu(ji,jj)  * zu_frc(ji,jj) ) &
@@ -771 +902 @@
          !
          IF( lk_vvl ) THEN                             !* Update ocean depth (variable volume case only)
-            !                                          !  ----------------------------------------------        
-            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
-            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
+            !                                          !  ----------------------------------------------
+            IF(ln_wd) THEN
+              hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1)
+              hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1)
+            ELSE
+              hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:)
+              hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:)
+            END IF
+
             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )
             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )
@@ -781 +918 @@
          !                                                 !  -----------------------
          !
-         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries 
+         CALL lbc_lnk( ua_e  , 'U', -1._wp )               ! local domain boundaries
          CALL lbc_lnk( va_e  , 'V', -1._wp )
 
-#if defined key_bdy  
+#if defined key_bdy
                                                            ! open boundaries
          IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
 #endif
-#if defined key_agrif                                                           
+#if defined key_agrif
          IF( .NOT.Agrif_Root() )  CALL agrif_dyn_ts( jn )  ! Agrif
 #endif
@@ -807 +944 @@
          !                                             !* Sum over whole bt loop
          !                                             !  ----------------------
-         za1 = wgtbtp1(jn)                                    
+         za1 = wgtbtp1(jn)
          IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN    ! Sum velocities
-            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) 
-            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) 
-         ELSE                                ! Sum transports
+            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)
+            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)
+         ELSE                                              ! Sum transports
             ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:)
             va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:)
          ENDIF
-         !                                   ! Sum sea level
+         !                                                 ! Sum sea level
          ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
          !                                                 ! ==================== !
@@ -841 +978 @@
       !
       ! Set advection velocity correction:
-      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN     
+      IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN
          un_adv(:,:) = zu_sum(:,:)*hur(:,:)
          vn_adv(:,:) = zv_sum(:,:)*hvr(:,:)
@@ -894 +1031 @@
       !
       !
-#endif      
+#endif
       !
       !                                   !* write time-spliting arrays in the restart
@@ -905 +1042 @@
       CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a                                   )
       CALL wrk_dealloc( jpi, jpj, zhf )
+      IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy )
       !
       IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts')
@@ -918 +1056 @@
       LOGICAL, INTENT(in) ::   ll_av      ! temporal averaging=.true.
       LOGICAL, INTENT(in) ::   ll_fw      ! forward time splitting =.true.
-      INTEGER, INTENT(inout) :: jpit      ! cycle length    
+      INTEGER, INTENT(inout) :: jpit      ! cycle length
       REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) ::   zwgt1, & ! Primary weights
                                                          zwgt2    ! Secondary weights
-      
+
       INTEGER ::  jic, jn, ji                      ! temporary integers
       REAL(wp) :: za1, za2
@@ -930 +1068 @@
 
       ! Set time index when averaged value is requested
-      IF (ll_fw) THEN 
+      IF (ll_fw) THEN
          jic = nn_baro
       ELSE
@@ -938 +1076 @@
       ! Set primary weights:
       IF (ll_av) THEN
-           ! Define simple boxcar window for primary weights 
-           ! (width = nn_baro, centered around jic)     
+           ! Define simple boxcar window for primary weights
+           ! (width = nn_baro, centered around jic)
          SELECT CASE ( nn_bt_flt )
               CASE( 0 )  ! No averaging
@@ -947 +1085 @@
               CASE( 1 )  ! Boxcar, width = nn_baro
                  DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    za1 = ABS(float(jn-jic))/float(nn_baro)
                     IF (za1 < 0.5_wp) THEN
                       zwgt1(jn) = 1._wp
@@ -956 +1094 @@
               CASE( 2 )  ! Boxcar, width = 2 * nn_baro
                  DO jn = 1, 3*nn_baro
-                    za1 = ABS(float(jn-jic))/float(nn_baro) 
+                    za1 = ABS(float(jn-jic))/float(nn_baro)
                     IF (za1 < 1._wp) THEN
                       zwgt1(jn) = 1._wp
@@ -969 +1107 @@
          jpit = jic
       ENDIF
-    
+
       ! Set secondary weights
       DO jn = 1, jpit
@@ -999 +1137 @@
       !
       IF( TRIM(cdrw) == 'READ' ) THEN
-         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )   
-         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) ) 
+         CALL iom_get( numror, jpdom_autoglo, 'ub2_b'  , ub2_b  (:,:) )
+         CALL iom_get( numror, jpdom_autoglo, 'vb2_b'  , vb2_b  (:,:) )
          IF( .NOT.ln_bt_av ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )   
-            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'sshbb_e'  , sshbb_e(:,:) )
+            CALL iom_get( numror, jpdom_autoglo, 'ubb_e'    ,   ubb_e(:,:) )
             CALL iom_get( numror, jpdom_autoglo, 'vbb_e'    ,   vbb_e(:,:) )
-            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) ) 
-            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'sshb_e'   ,  sshb_e(:,:) )
+            CALL iom_get( numror, jpdom_autoglo, 'ub_e'     ,    ub_e(:,:) )
             CALL iom_get( numror, jpdom_autoglo, 'vb_e'     ,    vb_e(:,:) )
          ENDIF
@@ -1012 +1150 @@
          ! Read time integrated fluxes
          IF ( .NOT.Agrif_Root() ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )   
+            CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b'  , ub2_i_b(:,:) )
             CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b'  , vb2_i_b(:,:) )
          ENDIF
@@ -1022 +1160 @@
          !
          IF (.NOT.ln_bt_av) THEN
-            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) ) 
+            CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e'  , sshbb_e(:,:) )
             CALL iom_rstput( kt, nitrst, numrow, 'ubb_e'    ,   ubb_e(:,:) )
             CALL iom_rstput( kt, nitrst, numrow, 'vbb_e'    ,   vbb_e(:,:) )
@@ -1070 +1208 @@
       CALL wrk_alloc( jpi, jpj, zcu )
       !
-      IF (lk_vvl) THEN 
+      IF (lk_vvl) THEN
          DO jj = 1, jpj
             DO ji =1, jpi
@@ -1093 +1231 @@
       ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
       IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)
-      
+
       rdtbt = rdt / FLOAT(nn_baro)
       zcmax = zcmax * rdtbt
-                     ! Print results
+                                                        ! Print results
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface'
@@ -1129 +1267 @@
            CASE( 0 )  ;   IF(lwp) WRITE(numout,*) '           Dirac'
            CASE( 1 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_baro'
-           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro' 
+           CASE( 2 )  ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_baro'
            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' )
       END SELECT
@@ -1142 +1280 @@
       ENDIF
       IF ( zcmax>0.9_wp ) THEN
-         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )          
+         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' )
       ENDIF
       !
@@ -1165 +1303 @@
       CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
       WRITE(*,*) 'ts_rst    : You should not have seen this print! error?', kt, cdrw
-   END SUBROUTINE ts_rst  
+   END SUBROUTINE ts_rst
    SUBROUTINE dyn_spg_ts_init( kt )       ! Empty routine
       INTEGER         , INTENT(in) ::   kt         ! ocean time-step
@@ -1171 +1309 @@
    END SUBROUTINE dyn_spg_ts_init
 #endif
-   
+
    !!======================================================================
 END MODULE dynspg_ts

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

@@ -9 +9 @@
    !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
    !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
+   !!            3.?  !  2014-09  (H. Liu) add wetting and drying  
    !!----------------------------------------------------------------------
 
@@ -42 +43 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE wadlmt          ! Wetting/Drying flux limting
 
    IMPLICIT NONE
@@ -94 +96 @@
       ENDIF
       !
-      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
-      !
       z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
       IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
+      !
+      z1_rau0 = 0.5_wp * r1_rau0
+      !
+      IF(ln_wd) CALL wad_lmt(sshb, z1_rau0 * (emp_b(:,:) + emp(:,:)), z2dt)
+
+      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
+      !
 
       !                                           !------------------------------!
@@ -110 +117 @@
       ! 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)
 

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -30 +30 @@
    !!            3.4  ! 2011-11  (C. Harris) decomposition changes for running with CICE
    !!                 ! 2012-05  (C. Calone, J. Simeon, G. Madec, C. Ethe) Add grid coarsening 
+   !!            3.6.?! 2014-09  (H. Liu) Add Wetting and Drying 
    !!----------------------------------------------------------------------
 
@@ -233 +234 @@
       NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, &
          &             jpizoom, jpjzoom, jperio
+      NAMELIST/namwad/ ln_wd, rn_wdmin1, rn_wdmin2, rn_wdld, nn_wdit
       !!----------------------------------------------------------------------
       !
@@ -257 +259 @@
       READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 )
 904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )   
+
+      REWIND( numnam_ref )              ! Namelist namwad in reference namelist : Parameters for Wetting/Drying
+      READ  ( numnam_ref, namwad, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namwad in configuration namelist : Parameters for Wetting/Drying
+      READ  ( numnam_cfg, namwad, IOSTAT = ios, ERR = 906)
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namwad in configuration namelist', lwp )
+
 
 ! Force values for AGRIF zoom (cf. agrif_user.F90)

branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

@@ -116 +116 @@
 #endif
 #if defined key_top
-   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
+   USE trcstp           ! passive tracer time-stepping     (trc_stp routine)
 #endif
+   USE wadlmt           ! wetting/drying limiter           (wad_lmt routine)
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)