Changeset 4375 for branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

Timestamp:

2014-01-28T14:55:35+01:00 (10 years ago)

Author:

hliu

Message:

updated gravity filters

File:

• branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (4 diffs)

branches/2013/dev_r4050_NOC_WaD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

@@ -76 +76 @@
       !!             - Save the trend (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt                               ! ocean time-step index
+      INTEGER             ::   ji, jj, jk                       ! dummy loop indices
+      REAL                ::   zcpx, zcpy                       ! gravity filter for W/D
+      REAL                ::   zij, zim1j, zijm1, zim1jm1   ! local variables
+      REAL                ::   dij, dim1j, dijm1, dim1jm1   ! local variables
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
       !!----------------------------------------------------------------------
@@ -82 +86 @@
       IF( nn_timing == 1 )  CALL timing_start('dyn_hpg')
       !
-      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
+
+      IF( l_trddyn .OR. ln_wad ) THEN               ! Temporary saving of ua and va trends (l_trddyn)
          CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
+
       !
       SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
@@ -95 +101 @@
       CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
       END SELECT
+
+      !! Gravity filter for W/D
+      IF(ln_wad) THEN
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zij     =  sshn(ji   , jj   )
+                  zim1j   =  sshn(ji- 1, jj   )
+                  zijm1   =  sshn(ji   , jj- 1)
+                  zim1jm1 =  sshn(ji-1 , jj- 1)
+                  dij     = bathy(ji   , jj   )
+                  dim1j   = bathy(ji- 1, jj   )
+                  dijm1   = bathy(ji   , jj- 1)
+                  dim1jm1 = bathy(ji-1 , jj- 1)
+
+                  zcpx = 0.5_wp + sign(0.5_wp, min(zij,zim1j) - max(-dij,-dim1j))
+                  zcpy = 0.5_wp + sign(0.5_wp, min(zij,zijm1) - max(-dij,-dijm1))
+                  
+                  ua(ji,jj,jk) = ztrdu(ji,jj,jk) + zcpx * ( ua(ji,jj,jk) - ztrdu(ji,jj,jk))
+                  va(ji,jj,jk) = ztrdv(ji,jj,jk) + zcpy * ( va(ji,jj,jk) - ztrdv(ji,jj,jk))
+               END DO
+            END DO
+         END DO
+      END IF
       !
       IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
@@ -100 +130 @@
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
          CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
+      ENDIF
+
+      IF( l_trddyn .OR. ln_wad ) THEN 
          CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
       ENDIF
+
       !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &