Issue with Bougeault and Andre length scale #56

Hello,
I am running into issues with the TKE closure with the Bougeault and Andre length scale (I believe this is the default scheme that is also used in NEMO).

For the wind driven entrainement test case with

turbulence:
  turb_method: first_order
  tke_method: 2
  len_scale_method: Bougeault_Andre

I have infinite values of dissipation epsilon. I looked at the routine potentialml and it seems that indeed we have ld(1) = 0 (at the lowest level the downward length scale is zero). This causes the length scale L(1) to be also zero and dissipation to be infinite. I suspect this is due to the goto 610 which does not allow to increment ld as it should (a couple of lines below, we have ld(i)=ld(i)+h(j) which is never realized at the bottom).

a quick fix is to add a minimum value for ld at the end of this loop: something like this fixes the problem:

610   continue
!     if (ld(i).gt.db(i)) ld(i)=db(i) !ld limited by distance to bottom
      if (ld(i).le.0.0) ld(i) = db(i)

This may not be the correct fix. I am happy to implement another suggestion if you tell me a better way.

Thank you