Changeset 11293

Timestamp:

2019-07-18T15:21:40+02:00 (6 years ago)

Author:

jchanut

Message:

#2292, Fixes courant number computation in adaptive vertical advection

File:

• NEMO/trunk/src/OCE/DYN/sshwzv.F90 (modified) (5 diffs)

NEMO/trunk/src/OCE/DYN/sshwzv.F90

@@ -284 +284 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
+      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars
       REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
       REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
@@ -297 +297 @@
          IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
          IF(lwp) WRITE(numout,*) '~~~~~ '
+         wi(:,:,:) = 0._wp
       ENDIF
       !
@@ -303 +304 @@
          DO jj = 2, jpjm1
             DO ji = 2, fs_jpim1   ! vector opt.
-               z1_e3w = 1._wp / e3w_n(ji,jj,jk)
-               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &  ! 2*rdt and not r2dt (for restartability)
-                  &                             + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
-                  &                                 MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                  &
-                  &                             + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
-                  &                                 MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
-                  &                               * r1_e1e2t(ji,jj)                                                  &
-                  &                             ) * z1_e3w
+               z1_e3t = 1._wp / e3t_n(ji,jj,jk)
+               Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )   &  ! 2*rdt and not r2dt (for restartability)
+                  &                             + ( MAX( e2u(ji  ,jj)*e3u_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
+                  &                                 MIN( e2u(ji-1,jj)*e3u_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                 &
+                  &                             + ( MAX( e1v(ji,jj  )*e3v_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
+                  &                                 MIN( e1v(ji,jj-1)*e3v_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
+                  &                               * r1_e1e2t(ji,jj)                                                 &
+                  &                             ) * z1_e3t
             END DO
          END DO
@@ -320 +321 @@
       !
       IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
-         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
+         DO jk = jpkm1, 2, -1                           ! or scan Courant criterion and partition
             DO jj = 1, jpj                              ! w where necessary
                DO ji = 1, jpi
                   !
-                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
+                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
+! alt:
+!                  IF ( wn(ji,jj,jk) > 0._wp ) THEN 
+!                     zCu =  Cu_adv(ji,jj,jk) 
+!                  ELSE
+!                     zCu =  Cu_adv(ji,jj,jk-1)
+!                  ENDIF 
                   !
                   IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
@@ -343 +350 @@
             END DO
          END DO
+         Cu_adv(:,:,1) = 0._wp 
       ELSE
          ! Fully explicit everywhere