Opened 11 years ago
Last modified 7 years 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: |
Description
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.
Change History (8)
comment:1 Changed 11 years ago by kji
comment:2 Changed 11 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.
Bo
comment:3 Changed 10 years ago by pesei
Hello Zhang, is the problem solved?
comment:4 Changed 10 years ago by bzhang
The problem is solved. Thanks a lot!
Bo
comment:5 Changed 10 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 10 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 9 years ago by pesei
Fix is not implemented yet (9.2beta or even changeset:31)
comment:8 Changed 7 years 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', & 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 . . .
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!