Changeset 6603 for branches/2016/dev_r6393_NOC_WAD

Timestamp:

2016-05-23T15:54:28+02:00 (8 years ago)

Author:

hliu

Message:

wet/dry for dynspg_ts.F90

File:

• branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (8 diffs)

branches/2016/dev_r6393_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -157 +157 @@
       REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy                 ! Wetting/Dying gravity filter coef.
       REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1           ! Wetting/Dying velocity filter coef.
+      REAL(wp), POINTER, DIMENSION(:,:) :: sshai                      ! Wetting/Dying velocity filter coef.
       !!----------------------------------------------------------------------
       !
@@ -168 +169 @@
       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, wduflt1, wdvflt1 )
+      IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1,sshai )
       !
       zmdi=1.e+20                               !  missing data indicator for masking
@@ -596 +597 @@
       !
       ! Initialize sums:
+      IF(ln_wd) sshai(:,:) = ssha(:,:)
       ua_b  (:,:) = 0._wp       ! After barotropic velocities (or transport if flux form)          
       va_b  (:,:) = 0._wp
@@ -687 +689 @@
          ENDIF
 #endif
-         IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
+         !IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)
          !
          ! Sum over sub-time-steps to compute advective velocities
@@ -701 +703 @@
             END DO
          END DO
-         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
-         IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 
+
+         IF(ln_wd) THEN
+           ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) * wdmask(:,:) 
+           &                + (1._wd - wdmask(:,:)) * (sshai(:,:) - sshn_e(:,:))) * ssmask(:,:)
+         ELSE
+           ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:)
+         END IF
+         
          CALL lbc_lnk( ssha_e, 'T',  1._wp )
 
@@ -888 +896 @@
                  zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj)
                  zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj)
-                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 
-                 zwy(ji,jj) = zv_spg * zcpy(ji,jj)
+                 zwx(ji,jj) = zu_spg * zcpx(ji,jj) * wdmask(ji,jj) * wdmask(ji+1, jj) 
+                 zwy(ji,jj) = zv_spg * zcpy(ji,jj) * wdmask(ji,jj) * wdmask(ji, jj+1)
               END DO
            END DO
@@ -1000 +1008 @@
          !                                   ! Sum sea level
          ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)
+         IF(ln_wd) ssha(:,:) = sshai(:,:)
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
@@ -1024 +1033 @@
       !
       ! Update barotropic trend:
-      IF( ln_dynadv_vec .OR. ln_linssh ) THEN
-         DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
-            va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
-         END DO
+      IF(ln_wd) THEN
+        IF( ln_dynadv_vec .OR. ln_linssh ) THEN
+           DO jk=1,jpkm1
+              ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b * wdmask(:,:)
+              va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b * wdmask(:,:)
+           END DO
+        ELSE
+           ! At this stage, ssha has been corrected: compute new depths at velocity points
+           DO jj = 1, jpjm1
+              DO ji = 1, jpim1      ! NO Vector Opt.
+                 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
+                    &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
+                    &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
+                 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
+                    &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
+                    &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
+              END DO
+           END DO
+           CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
+           !
+           DO jk=1,jpkm1
+              ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b * wdmask(:,:)
+              va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b * wdmask(:,:)
+           END DO
+           ! Save barotropic velocities not transport:
+           ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
+           va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
+        ENDIF
       ELSE
-         ! At this stage, ssha has been corrected: compute new depths at velocity points
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
-                  &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
-                  &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
-               zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
-                  &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
-                  &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
-         CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
-         !
-         DO jk=1,jpkm1
-            ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
-            va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
-         END DO
-         ! Save barotropic velocities not transport:
-         ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
-         va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
-      ENDIF
+        IF( ln_dynadv_vec .OR. ln_linssh ) THEN
+           DO jk=1,jpkm1
+              ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b
+              va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b
+           END DO
+        ELSE
+           ! At this stage, ssha has been corrected: compute new depths at velocity points
+           DO jj = 1, jpjm1
+              DO ji = 1, jpim1      ! NO Vector Opt.
+                 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj) &
+                    &              * ( e1e2t(ji  ,jj) * ssha(ji  ,jj)    &
+                    &              +   e1e2t(ji+1,jj) * ssha(ji+1,jj) )
+                 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj) &
+                    &              * ( e1e2t(ji,jj  ) * ssha(ji,jj  )    &
+                    &              +   e1e2t(ji,jj+1) * ssha(ji,jj+1) )
+              END DO
+           END DO
+           CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
+           !
+           DO jk=1,jpkm1
+              ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b
+              va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b
+           END DO
+           ! Save barotropic velocities not transport:
+           ua_b(:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
+           va_b(:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
+        ENDIF
+
+      END IF
       !
       DO jk = 1, jpkm1