With key_vvl only:

• In dynspg_flt.F90, next transport is computed with now vertical scale factors so that:

spgu(ji,jj) = spgu(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk)

spgv(ji,jj) = spgv(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk)

it should rather be:

spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)

spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)

• In bdydyn.F90:

Barotropic velocity is computed (then removed) from ua/va arrays again with "now" vertical scale factors:

pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
...
pu2d(:,:) = pu2d(:,:) * phur(:,:)
pv2d(:,:) = pv2d(:,:) * phvr(:,:)

again fse3u_a/fse3v_a arrays should be used, as well as inverse depths at next time step.

in bdyvol.F90 the sum of the forcing contribution is performed using e1v * e2u whereas the surface of a box at T-point is simply e1t*e2t. Therefore the folling lines:

z_cflxemp = 0.e0 DO jj = 1, jpj

DO ji = 1, jpi

z_cflxemp = z_cflxemp + emp(ji,jj) * bdytmask(ji,jj) * tmask_i(ji,jj) / rau0 * e1v(ji,jj) * e2u(ji,jj)

END DO

END DO

should be changed into

z_cflxemp = 0.e0 DO jj = 1, jpj

DO ji = 1, jpi

z_cflxemp = z_cflxemp + emp(ji,jj) * bdytmask(ji,jj) * tmask_i(ji,jj) / rauw * e1t(ji,jj) * e2t(ji,jj)

END DO

END DO

or in a more compact form (like the one used in obcvol.F90)

z_cflxemp = SUM ( emp(:,:) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) / rauw )

by the way, introducing a scalar equal to 1/rauw will reduce the cpu time (less division):

REAL(wp)

zraur ! temporary scalar

...

zraur = 1.e0 / rauw z_cflxemp = SUM ( emp(:,:) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) * zraur )

gurvan

1. In VVL case, temporary array for vertical scale factors are swapped after the divergence computation, so the 3d computation gets wrong.
   
2. For the same raison, the 1rst level vertical scale factor used in sea-ice is also wrong. Furthermore, we should used an average scale factor on nn_sbc.
   

Fix:

1. Move swap of arrays at the beginning of ssh_wzv and call divcur just after inside ssh_wzv. An improvement could be to call cla_div inside div_cur
   
2. Define in domzgr_substitute, fse3t_m corresponding to averaged sshm
   

Thanks to Dave and Rachel for addressing this.