Opened 5 years ago

Last modified 18 months ago

#48 assigned Defect

Crash with grib1 GFS/FNL

Reported by: bzhang Owned by: ignacio
Priority: major Milestone: FLEXPART 9.2
Component: FP input data Version: flexpart 8.2-2
Keywords: Cc:


Hi all,

I am trying to create a transport climatology for an observatory station. I collected some GFS wind fields of 2005 and 2006 which are in old grib1 format. Then I got error message generated by "readoutgrid" like this:

Mother domain:

Longitude range: 0.00 360.00 Grid distance: 1.00
Latitude range: -90.00 90.00 Grid distance: 1.00


Then I printed some related variable:
write(*,*) outlon0,outlat0,eps,dx,dy,nxmin1,nymin1
write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout

-179.0000 0.000000 9.9999997E-05 1.000000 1.000000 360 180

360.0000 90.00000 0.000000 -90.00000 181.0000 90.00000 1.000000 1.000000

In OUTGRID, I normally set the longitude domain from -179 to 181, which runs flawlessly with grib2 GFS data after 2007. Same settings cannot get through with grib1 GFS data. I personally remember fields are set in a different sequence in grib1, which could be the potential reason for this. I am current looking for the files set the domain. However, I have never digged too deep into FLEXPART codes. I hope someone can point me to the right files or provide solutions/suggestions. I will also post any progress I make in the future. Thanks.

Change History (8)

comment:1 Changed 5 years ago by kji


Your problem may be due to the xaux1, xaux2 checks in "gridcheck_gfs", ln 237 for the .f90 version (if you're using the older 82-2 version, it may be a few lines below at around 242).

The if statement goes:

if((abs(xaux1).lt.eps).and.( then

I would try printing out the xaux1 and xaux2 values of your GRIB1 files to make sure this IF statement is satisfied. I've seen GRIB1 files where the xaux2 would be "-1" (as in, 360-1, I guess) instead of "359", in which case the statement won't satisfy.

I would note that this hint was given to me a few years back by John Miller at NOAA.

Hope this helps!

comment:2 Changed 5 years ago by bzhang

Thanks a lot!

The problem is exactly the same as you said.

I added a line 'if (xaux2.eq.-1) xaux2 = 359' before the IF statement, which I guess would survive for both grib1 and grib2 cases.

I will post how it goes when I see the results on map.


comment:3 Changed 5 years ago by pesei

Hello Zhang, is the problem solved?

comment:4 Changed 5 years ago by bzhang

The problem is solved. Thanks a lot!

comment:5 Changed 5 years ago by pesei

  • Milestone set to FLEXPART 9.2

I am adding this ticket to the FLEXPART 9.2 milestone. Please, can NILU check whether the fix proposed here needs to be integrated into the next release.

comment:6 Changed 5 years ago by ignacio

  • Owner changed from somebody to ignacio
  • Status changed from new to assigned

The fix will be available in version 9.2

comment:7 Changed 4 years ago by pesei

Fix is not implemented yet (9.2beta or even changeset:31)

comment:8 Changed 18 months ago by donaldjmorton

I am here three years later, working on FPv10.1 and, yet again, this problem has appeared. I will try to get NILU to put this change in their version of the code as they and the community consider adopting some of our other CTBTO-inspired changes in the near future.

Note that this fix was applied to FPv9.3 about 14 months ago.

For now, in FPv10.1, I've done roughly what bzhang did four years ago. In gridcheck_gfs.f90 I inserted the following code just after xaux1in, etc. are read from the GRIB file:

  call grib_get_real8(igrib,'latitudeOfFirstGridPointInDegrees', &
  call grib_check(iret,gribFunction,gribErrorMsg)

!!!!!!!!  DJM ARTIFICIAL CHANGE FOR NCEP GRIB1 - change value from -1 to 359
!!!!!!!!  See CTBTO project Ticket #112
!!!!!!!!  Also see mainstream Ticket #48
  if (xaux2in .lt. 0) xaux2in = 359.0

Note: See TracTickets for help on using tickets.
hosted by ZAMG