Bug in dyn_zdf when using vector formulation

Context

• Branches impacted: main, branch_5.0

Analysis

There is a possible bug in dyn_zdf when using the vector formulation.

The surface boundary conditions for the u and v components are given as:

! u direction
DO_1Di( 0, 0 )   !* Surface boundary conditions
   zwi(ji,1) = 0._wp
   zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji  ,jj,2) )   &
   &           / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
   !                    ! w at u-point
   zWu  = zDt_4 * ( e1e2t(ji,jj)*wi(ji,jj,2) + e1e2t(ji+1,jj)*wi(ji+1,jj,2) ) * r1_e1e2u(ji,jj)
   zzWp_e3 = MAX( zWu, 0._wp ) / e3uw(ji,jj,2,Kaa)
   zws(ji,1) = zzws - zzWp_e3
   zwd(ji,1) = 1._wp - zzws + zzWp_e3
END_1D  

! v direction
DO_1Di( 0, 0 )   !* Surface boundary conditions
   zwi(ji,1) = 0._wp
   zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) )    &
   &           / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
   !                    ! w at u-point
   zWv = zDt_4 * ( e1e2t(ji,jj)*wi(ji  ,jj,2) +  e1e2t(ji,jj+1)*wi(ji,jj+1,2) )
   zzWp_e3 = MAX( zWv, 0._wp ) / e3vw(ji,jj,2,Kaa)
   zws(ji,1) = zzws - zzWn_e3
   zwd(ji,1) = 1._wp - zzws + zzWp_e3
END_1D            

The two components have inconsistent implementations. For the v component, r1_e1e2v is not used in the expression for zWv and zws uses zzWn_e3, which is not defined in the loop (it should probably be zzWp_e3):

- zWu = zDt_4 * ( e1e2t(ji,jj)*wi(ji,jj,2) + e1e2t(ji+1,jj)*wi(ji+1,jj,2) ) * r1_e1e2u(ji,jj)
+ zWv = zDt_4 * ( e1e2t(ji,jj)*wi(ji  ,jj,2) +  e1e2t(ji,jj+1)*wi(ji,jj+1,2) )
- zws(ji,1) = zzws - zzWp_e3
+ zws(ji,1) = zzws - zzWn_e3

Because zzWn_e3 is not defined in the loop for the v component, the definition from the previous loop is used instead (which handles internal depths):

DO_2Dik( 0, 0,    1, jpkm1, 1 )
   zzwi = - zDt_2 * ( avm(ji,jj+1,jk  ) + avm(ji,jj,jk  ) )    &
   &           / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk  ,Kmm) ) * wvmask(ji,jj,jk  )
   zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) )    &
   &           / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
   !                    ! w at u-point
   zWv  = zDt_4 * (  ( e1e2t(ji,jj)*wi(ji,jj,jk  ) + e1e2t(ji,jj+1)*wi(ji,jj+1,jk  ) )  &
   &            + ( e1e2t(ji,jj)*wi(ji,jj,jk+1) + e1e2t(ji,jj+1)*wi(ji,jj+1,jk+1) )  ) * r1_e1e2v(ji,jj)
   zzWn_e3 = MIN( zWv, 0._wp ) / e3vw(ji,jj,jk  ,Kaa)
   zzWp_e3 = MAX( zWv, 0._wp ) / e3vw(ji,jj,jk+1,Kaa)
   zwi(ji,jk) = zzwi + zzWn_e3
   zws(ji,jk) = zzws - zzWp_e3
   zwd(ji,jk) = 1._wp - zzwi - zzws - zzWn_e3 + zzWp_e3
END_2D

I'm not familiar with the underlying algorithms, but this doesn't seem correct to me.

Fix

Unsure- I'm not familiar with the science behind this code.