PARDSP causing NaN particle trajectories in PARDUMPs
Posted: 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
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