Changeset 15422 for NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90

Timestamp:

2021-10-21T11:19:25+02:00 (3 years ago)

Author:

jcastill

Message:

Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway

File:

• NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90 (modified) (7 diffs)

NEMO/branches/UKMO/r14075_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
@@ -58 +59 @@
    PUBLIC    dyn_nxt   ! routine called by step.F90
 
+   !! Substitution 
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -182 +185 @@
             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  
@@ -214 +221 @@
                END DO
             END DO
+            ! limit velocities  
+            IF (ln_ulimit) THEN  
+               call dyn_limit_velocity (kt)  
+            ENDIF
             !                             ! ================!
          ELSE                             ! Variable volume !
@@ -263 +274 @@
                   END DO
                END DO
+               ! limit velocities  
+               IF (ln_ulimit) THEN  
+                  call dyn_limit_velocity (kt)  
+               ENDIF
                !
             ELSE                          ! Asselin filter applied on thickness weighted velocity
@@ -290 +305 @@
                   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)
@@ -401 +420 @@
    END SUBROUTINE dyn_nxt
 
+   SUBROUTINE dyn_limit_velocity (kt)  
+      ! limits maximum values 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 = 2, jpjm1  
+            DO ji = fs_2, fs_jpim1   ! vect. opt. 
+               ! 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)  
+               ! limit currents  
+               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
+            ENDDO  
+         ENDDO  
+      ENDDO  
+       CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U',  -1., vn(:,:,:), 'V',  -1. ) 
+   
+    END SUBROUTINE dyn_limit_velocity
+
    !!=========================================================================
 END MODULE dynnxt