Computation of BVF and dRz (in rho_eos.F and step3d_t.F) are (a little bit) wrong when using SPLIT_EOS
Computation of BVF in rho_eos.F seems to have an erroneous extra 2 factor in front of "qp2":
ifdef SPLIT_EOS
dpth=z_w(i,j,N)-0.5*(z_r(i,j,k+1)+z_r(i,j,k))
cff2=( rho1(i,j,k+1)-rho1(i,j,k) ! Elementary
& +(qp1(i,j,k+1)-qp1(i,j,k)) ! adiabatic
& *dpth*(1.-2.*qp2*dpth) ! difference
& )
bvf(i,j,k)=-cff*cff2 / (z_r(i,j,k+1)-z_r(i,j,k))
Same for dRz in step3d_t.F:
ifdef SPLIT_EOS
dpth=z_w(i,j,N)-0.5*(z_r(i,j,k+1)+z_r(i,j,k))
dRz =rho1(i,j,k+1)-rho1(i,j,k)
& +(qp1(i,j,k+1)- qp1(i,j,k))
& *dpth*(1.-2.*qp2*dpth)
to be compared with the computation in prsgrd.F:
ifdef SPLIT_EOS
dpth=z_w(i,j,N)-0.5*(z_r(i,j,k+1)+z_r(i,j,k))
dR(i,k)=rho1(i,j,k+1)-rho1(i,j,k) ! Elementary
& +(qp1(i,j,k+1)-qp1(i,j,k)) ! adiabatic
& *dpth*(1.-qp2*dpth) ! difference
#########
With SPLIT_EOS, we have : rho = rho1 + qp1 * z + (-qp1 qp2) * z^2
So the elementary adiabatic difference in the vertical should be computed as:
drho/dz(adiab) = drho1/ds + dqp1/ds * z - d(qp1*qp2)/ds * z^2 = drho1/ds + dqp1/ds * z * ( 1 - qp2 z )
as qp2 is a constant here.
See equ. (6.6) in Shchepetkin & McWilliams, 2008 (https://doi.org/10.1016/S1570-8659(08)01202-0)