Bad do-loop in richardson.f90

[adingwell]

When determining the levels to use for the Richardson number, the following loop is executed:

  ! Integrate z up to one level above zt
  !*************************************

  do k=2,nuvz
    ...
    if (ri.gt.ric .and. thetaold.lt.theta) goto 20
    ...
  end do
  
  20   continue

The final k-value is then used for the main calculation in a second loop:

  do i=1,20
    zl=zold+real(i)/20.*(z-zold)
    ul=ulev(k-1)+real(i)/20.*(ulev(k)-ulev(k-1))
    vl=vlev(k-1)+real(i)/20.*(vlev(k)-vlev(k-1))
    thetal=thetaold+real(i)/20.*(theta-thetaold)
    rhl=rhold+real(i)/20.*(rh-rhold)
    ril=ga/thetaref*(thetal-thetaref)*(zl-zref)/ &
      max(((ul-ulev(2))**2+(vl-vlev(2))**2+b*ust**2),0.1)
    zl2=zl
    theta2=thetal
    if (ril.gt.ric) goto 25
    zl1=zl
    theta1=thetal
  end do

When a level is caught by the if-statement in the first loop, everything works as expected. However, if the expression in that if-statement is never satisfied, the loop exits by reaching the upper limit (k=nuvz). In standard Fortran this means that the final value will be k=nuvz+1, which means that ulev and vlev will exceed their bounds in the second loop!

This problem might go undetected without strict compiler flags, but it's probably causing some strange behaviour even when the error is not raised.

A quick work-around to this problem is to add "k=k-1" between "end do" and "20 continue". I'm not sure if this is a good solution from a physical point of view...

[espen]

Fixed for version 10 (git branch "dev")

eso@…

[pesei]

Thanks, should also be fixed in Fp9.2

[pesei]

  ! fix bug ticket:139 as follows:
  print*,'Warning - in richardson.f90, no Ri_c found'
  k = nuvz

added for Flexpart_9.0.3