Bugs with adaptive vertical advection: 3) Courant number computation

[jchanut]

Context

Adaptive (based on courant number) vertical advection, ln_zad_Aimp=.true.

Analysis

Courant number computation in wAimp (sshwzv.F90) is made at T-points. However, vertical scale factors at w interfaces are used:

              DO jk = 1, jpkm1            ! calculate Courant numbers
	         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
	            END DO
	         END DO
	      END DO

It should be:

              DO jk = 1, jpkm1            ! calculate Courant numbers
	         DO jj = 2, jpjm1
	            DO ji = 2, fs_jpim1   ! vector opt.
	               z1_e3w = 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_e3w
	            END DO
	         END DO
	      END DO

Then, the partitioning function between implicit and explicit vertical velocities is computed at w-points, however it considers the temporary Courant number at interfaces:

           zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )

which should be:

           zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )

Note that due to the later assignment of Cu_adv array in the loop, jk vertical loop needs to be reversed.

Question: Should we rather consider using Cu_adv in the upstream direction like in Shchepetkin (2015) paper ?

...

[jchanut]

In 11293:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found

[jchanut]

In 11294:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found

[jchanut]

In 11414:

Error: Failed to load processor CommitTicketReference

No macro or processor named 'CommitTicketReference' found