Opened 4 years ago

Last modified 16 months ago

#139 accepted Defect

Bad do-loop in richardson.f90

Reported by: adingwell Owned by: pesei
Priority: minor Milestone: FLEXPART 9.2
Component: FP coding/compilation Version: FLEXPART 9.0.2
Keywords: Cc:

Description

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...

Change History (3)

comment:1 Changed 4 years ago by espen

Fixed for version 10 (git branch "dev")

eso@…

comment:2 Changed 3 years ago by pesei

  • Milestone set to FLEXPART 9.2
  • Owner set to pesei
  • Status changed from new to accepted

Thanks, should also be fixed in Fp9.2

comment:3 Changed 16 months ago by pesei

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

fixed in Flexpart_9.0.3

Last edited 16 months ago by pesei (previous) (diff)
Note: See TracTickets for help on using tickets.
hosted by ZAMG