#139closed - Defect (None)
Bad do-loop in richardson.f90
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...