Changeset 11080 for branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

Timestamp:

2019-06-06T15:32:24+02:00 (5 years ago)

Author:

jcastill

Message:

Latest changes from Jason

File:

• branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90 (modified) (6 diffs)

branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -31 +31 @@
    USE dynadv         ! dynamics: vector invariant versus flux form
    USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
+   USE dynspg
    USE domvvl         ! variable volume
    USE bdy_oce   , ONLY: ln_bdy
@@ -179 +180 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         ! limit velocities
+         IF (ln_ulimit) THEN
+            call dyn_limit_velocity (kt)
+         ENDIF
          IF(.NOT.ln_linssh ) THEN
             DO jk = 1, jpkm1
@@ -203 +208 @@
                END DO
             END DO
+            ! limit velocities
+            IF (ln_ulimit) THEN
+               call dyn_limit_velocity (kt)
+            ENDIF
             !                             ! ================!
          ELSE                             ! Variable volume !
@@ -250 +259 @@
                   END DO
                END DO
+               ! limit velocities
+               IF (ln_ulimit) THEN
+                  call dyn_limit_velocity (kt)
+               ENDIF
                !
             ELSE                          ! Asselin filter applied on thickness weighted velocity
@@ -277 +290 @@
                   END DO
                END DO
+               ! limit velocities
+               IF (ln_ulimit) THEN
+                  call dyn_limit_velocity (kt)
+               ENDIF
                e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
                e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
@@ -353 +370 @@
    END SUBROUTINE dyn_nxt
 
+   SUBROUTINE dyn_limit_velocity (kt)
+      ! limits maxming vlaues of un and vn by volume courant number
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) :: zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm
+
+      ! limit fluxes
+      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number
+      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! U direction
+               zzu = un(ji,jj,jk)
+               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)
+               ! ips is 1 if flow is positive othersize zero
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
+               ! idev is 1 if divergent flow otherwise zero
+               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )
+               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
+               zzcn = zzcn * zcn
+               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
+               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &
+                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)
+               ! limit currents
+               un(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
+               ! V  direction
+               zzu = vn(ji,jj,jk)
+               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)
+               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )
+               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )
+               ! idev is 1 if divergent flow otherwise zero
+               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )
+               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )
+               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp
+               zzcn = zzcn * zcn
+               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
+               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &
+                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)
+               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max (zzu,zmlim) *ism
+            ENDDO
+         ENDDO
+      ENDDO
+
+    END SUBROUTINE dyn_limit_velocity
+
    !!=========================================================================
 END MODULE dynnxt