PARDSP causing NaN particle trajectories in PARDUMPs

Post any defects you find in the HYSPLIT software here. The HYSPLIT Developers carefully monitor this list and will work diligently to repair any reported problems. When posting a bug report, please specify both the HYSPLIT version and operating system you are using.
Post Reply
wsr
Posts: 9
Joined: August 23rd, 2017, 12:08 am
Registered HYSPLIT User: No

PARDSP causing NaN particle trajectories in PARDUMPs

Post by wsr » May 23rd, 2019, 4:18 pm

ISSUE: NaNs generated for Hysplit dispersion particle trajectories

CONTEXT:

When running the Hysplit dispersion model, some of the PARDUMP particle elevations showed up as NaNs. In tracking this down, it appears that elements of the ZPOS vector were NaN's and those NaNs originate in the pardsp.f routine. In particular, the NaNs result from the computation of WRATIO. WRATIO is passed in as an argument to the PARDSP routine and also is modified by PARDSP (its INTENT is (INOUT)).

DETAILED DESCRIPTION:

The NaNs arise from the computation of WRATIO the line 267:

267 WRATIO = RAUTO*WRATIO+SQRT(1.0-RAUTO*RAOUT)*WW/SIGW + VSCALE*(1.0-RAUTO)*SIGDZ

whenever WW and SIGW are both 0. Earlier in the routine, SIGW is computed from the following three statements in lines 240, 241 and 244

240 FACT = MAX(0.0,MIN(1.0,(ZSG(KB)-ZPOS)/(ZSG(KB)-ZSG(KT))
241 VM = (WMIX(KT)-WMIX(KB))*FACT + WMIX(KB)
244 SIGW = SQRT(VM)

Apparently it is possible for SIGW to be 0, for example if WMIX(KT)=WMIX(KB)=0, or if WMIX(KB) = 0 and FACT = 0 which is possible if ZPOS < ZSG(KB) < ZSG(KT)
or ZSG(KT) < ZSG(KB) < ZPOS.

WW is computed as SIGW times a random_number (lines 246 to 251), where the random number comes from calling GASDEV inside of the call to PARVAR or if KRAND=1 from an element of the rannumb array.

In either case, whenever SIGW is 0, WW is also 0. WW is not used in any other spot in PARDSP and since it is not in a common block nor passed as an argument, and generates no other side effects.

When WRATIO is a NaN (SIGW = 0 and WW = 0), then in line 269:

269 WPRIME = WRATIO*(SIGW+(WPRIME*DELT*SIGDZ))

WPRIME becomes a NaN regardless of the values of SIGW, DELT and SGIDZ. This in turn causes the calculation of ZPOS in the line 272,

272 ZPOS=ZPOS-(WPRIME/(ZMDL-ZSF))*DELT

to yield a NaN. This seems to be a bug, especially since WW = random_number*SIGW, so the ratio WW/SIGW in the calculation of WRATIO should always yield just the random number.

PROPOSED FIX:

1. Add declaration

REAL :: GASDEV

to the declarations section of pardsp.f

2. Replace the lines 246 to 251

IF (KRAND.EQ.1) THEN
if ((nrand+knum) .GT. maxrand)nrand = 1
WW = rannmb(nrand+knum)*SIGW
ELSE
CALL PARVAR(SIGW,WW,ISEED)
ENDIF

with

IF (KRAND.EQ.1) THEN
if ((nrand+knum) .GT. maxrand) nrand = 1
WW = rannmb(nrand+knum)
ELSE
WW=GASDEV(iseed)
ENDIF

Note that the WW=GASDEV(iseed) is taken from the single line of code in the parvar routine, avoiding the multiply by SIGW.

3. Replace the lines 267 and 268:

WRATIO=RAUTO*WRATIO+SQRT(1.0-RAUTO*RAUTO)*WW/SIGW &
+VSCALE*(1.0-RAUTO)*SIGDZ

with the lines

WRATIO=RAUTO*WRATIO+SQRT(1.0-RAUTO*RAUTO)*WW &
+VSCALE*(1.0-RAUTO)*SIGDZ

alicec
Posts: 158
Joined: February 8th, 2016, 12:56 pm
Registered HYSPLIT User: Yes

Re: PARDSP causing NaN particle trajectories in PARDUMPs

Post by alicec » May 30th, 2019, 12:58 pm

Could you share with us the setup you are using so we can reproduce this? In particular we would like know the value of KBLT in the SETUP.CFG file and the meteorological data you are using. Complete CONTROL and SETUP.CFG files would be helpful.

The SIGW in the pardsp routine is the width of the turbulent velocity distribution.
The WW is thus random number drawn from a distribution with that width.
WPRIME is then the turbulent velocity component. See the following reference for how wprime is calculated.
Chock and Winkler (1994)
J. Geophys. Res., Vol 99, No D1, Pages 1033-1041

A zero SIGW would mean no turbulence, so the correct thing to do would be to simply make WPRIME zero in that case, and the zpos would not be affected by turbulent dispersion and the computational particle would simply be advected by the mean wind.

However, for most options, we make the assumption that there is always some turbulent mixing and set a minimum value for the wmix (from which sigw is calculated) in the routines that calculate the mixing.

We will have to reproduce what you are doing and understand which options are producing the 0 sigma values so we can create an appropriate fix. It would be helpful if you could send us information on your run - including SETUP and CONTROL files - particularly which mixing options (value of KBLT) you are using in the SETUP file and what meteorological data.

Post Reply