Crash with grib1 GFS/FNL

[bzhang]

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

#### FLEXPART MODEL ERROR! PART OF OUTPUT ####
#### GRID IS OUTSIDE MODEL DOMAIN. CHANGE ####
#### FILE OUTGRID IN DIRECTORY ####

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.

[kji]

Hello,

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.(xaux2.ge.359)) 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!

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

Bo

[pesei]

Hello Zhang, is the problem solved?

[bzhang]

The problem is solved. Thanks a lot!
Bo

[pesei]

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.

[ignacio]

The fix will be available in version 9.2

[pesei]

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

[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', &
       yaux2in,iret)
  call grib_check(iret,gribFunction,gribErrorMsg)


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


  xaux1=xaux1in
  xaux2=xaux2in
  yaux1=yaux1in
  yaux2=yaux2in
.
.
.