Changeset 12453 for NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90

Timestamp:

2020-02-25T16:29:34+01:00 (4 years ago)

Author:

jcastill

Message:

First implementation of the branch - compiling after merge

File:

• NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90 (modified) (6 diffs)

NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90

@@ -33 +33 @@
    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
@@ -178 +179 @@
             vn(:,:,jk) = va(:,:,jk)
          END DO
+         ! limit velocities 
+         IF (ln_ulimit) THEN 
+            call dyn_limit_velocity (kt) 
+         ENDIF
          IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
 !!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a  
@@ -210 +215 @@
                END DO
             END DO
+            ! limit velocities 
+            IF (ln_ulimit) THEN 
+               call dyn_limit_velocity (kt) 
+            ENDIF
             !                             ! ================!
          ELSE                             ! Variable volume !
@@ -275 +284 @@
                   END DO
                END DO
+               ! limit velocities 
+               IF (ln_ulimit) THEN 
+                  call dyn_limit_velocity (kt) 
+               ENDIF
                !
             ELSE                          ! Asselin filter applied on thickness weighted velocity
@@ -302 +315 @@
                   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)
@@ -377 +394 @@
    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