error when computing h_rnf in sbcrnf

There is a bug in sbcrnf.F90 when computing h_rnf and phdivn in non-linear free surface:

            DO_2D( 1, 2, 1, 2 )         ! update the depth over which runoffs are distributed
               h_rnf(ji,jj) = 0._wp
               DO jk = 1, nk_rnf(ji,jj)
                  h_rnf(ji,jj)     = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm) ! recalculates h_rnf to be the depth in metres to the bottom of the relevant grid box
#if defined key_RK3
                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj)                    ! RK3: rnf forcing at n+1/2
#else
                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)   ! MLF: rnf forcing at Kmm (n)
#endif
               END DO
            END_2D

This loop was originally (in 4.2) in 2 parts:

            DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )         ! update the depth over which runoffs are distributed
               h_rnf(ji,jj) = 0._wp
               DO jk = 1, nk_rnf(ji,jj)                             ! recalculates h_rnf to be the depth in metres
                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box
               END DO
            END_2D
            DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )         ! apply the runoff input flow
               DO jk = 1, nk_rnf(ji,jj)
                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj)
               END DO
            END_2D

It must be kept in 2 parts whereas h_rnf and phdivn are wrong!