Issue using *.grb and *.grib2

[tjaya]

We have started working recently in the field of air pollution modelling. FLEXPART is a software we found apt for our application. There's an enquiry that we have regarding usage of input data set, kindly help us to resolve.

We have run the model for 20100109 date using both *.grb and *.grib2 formats (bothFNL data). The following are the two changes done in the files for making it run for the *.grb files.

1. For using the *.grb format, we had commented the line no. 290 [ if (xaux.eq.0.) xaux=-179.0 ] in readwind_gfs.f90.
2. We had also changed the value of OUTLONLEFT (GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID ) in OUTGRID to 0.0 (default was -179.0).

I am hereby attaching the plots generated from pflexible for your reference. The files "AIRTRACER_fp_grib1.png" and "AIRTRACER_tc_grib1.png" have been generated using *.grb format input files and the files "AIRTRACER_fp_grib2.png" and "AIRTRACER_tc_grib2.png" have been generated using *.grib2 input files.

We were expecting similar sensitivity output plots, whereas the course of travel in the plots also looks different in addition to the max. values plotted. We would like to know if there is any other change required for making the run using *.grb format input files. Since we are just beginners in this area of study, could you please help us to understand if this is acceptable and why the difference.

[donaldjmorton]

I have been facing the same question - how to make FLEXPART run with NCEP GRIB1 met data - and today sent the following email to the fp-dev group. I will be exploring this in more detail. I suspect that the change to readwinds_gfs is only a bandaid, and that one needs to go back to gridcheck_gfs for the full picture.

Here is the message I sent to fp-dev:

We have been working on replacing all hard-coded GRIB parameter codes in readwind* and gridcheck* with variable names, that then reference a variable table (Vtable) to get the parameters - much like WRF does. Good progress has been made in the ECMWF GRIB1, GRIB2 and GRIB1/2 inputs, as well as NCEP GRIB2.

I'm now fighting a bit with trying to get FLEXPART running with NCEP GRIB1 files, and I came across Ticket #134 by tjaya, and it appears they did exactly the same thing I did to get it running -

For using the *.grb format, we had commented the line no. 290 [ if (xaux.eq.0.) xaux=-179.0 ] in readwind_gfs.f90.
We had also changed the value of OUTLONLEFT (GEOGRAFICAL LONGITUDE OF LOWER LEFT CORNER OF OUTPUT GRID ) in OUTGRID to 0.0 (default was -179.0).
Like them, I didn't feel 100% comfortable with making this change, though.

I queried with Delia, and indirectly Sabine, today, wondering why a grid that has longitudes running from 0 to 360 is shifted this way, and my understanding is that this is done to match the range of the ECMWF inputs, as well as to place EU in the center, rather than the edge. All perfectly good reasons, in my opinion, and I certainly don't want to mess with that policy.

So, the question I have is, is it standard operating procedure to make code adjustments in order to get FLEXPART working with NCEP GRIB1 data? The key variable in readwind_gfs that seems to influence this is xlon0, and it looks like that is initially set in gridcheck_gfs.F90, dependent on xaux1, which comes from "longitudeOfFirstGridPointInDegrees." I still have to dig deeper, but it appears that things behave differently in GRIB1 and GRIB2 (though in both cases that first longitude in the grib files is 0.0)

Before I dive too deep into the black hole, I'm just wondering if there's anybody that can shed some history on this. Was there a time that these routines worked, without modification, on NCEP GRIB1, or were they written primarily for NCEP GRIB2?

[donaldjmorton]

Hello, I have investigated this more thoroughly, running NCEP GRIB1 and GRIB2 next to each other with the same source code. What I have found to be differentiating is that in GRIB1, the longitudeOfLastGridPointInDegrees is -1, but in GRIB2 it is 359.

This starts to become a problem in gridcheck_gfs(), where that value is read in as xaux2in.

Because of that difference, in GRIB1, the mother domain looks like this in readwind_gfs():

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

  Releasepoints allocated:            1
  Particles allocated for this run:      2000000 , released in simulation:         4000
  Allocating fields for nested and global output (x,y):          360         180
 ---- readwind_gfs() dump ----
 nxshift, dx:            0   1.00000000    
 xaux, xauxin:    0.00000000       0.0000000000000000     
 xlon0:    0.00000000    
 -----------------------------

and it looks like this in GRIB2:

Mother domain:
  Longitude range:    -179.00     181.00   Grid distance:       1.00
  Latitude range:     -90.00      90.00   Grid distance:       1.00

  Releasepoints allocated:            1
  Particles allocated for this run:      2000000 , released in simulation:         4000
  Allocating fields for nested and global output (x,y):          360         180
 ---- readwind_gfs() dump ----
 nxshift, dx:            0   1.00000000    
 xaux, xauxin:    0.00000000       0.0000000000000000     
 xlon0:   -179.000000    
 -----------------------------

The guilty party is that value of xlon0, and this value is ultimately determined in gridcheck_gfs(), as a function of xauxin2 (in a very indirect way).

I "hacked" a solution in gridcheck_gfs() right after xaux2in is read from the GRIB file. I simply added in

if (xaux2in .lt. 0) xaux2in = 359.0

I believe this has resolved the problem right at the source. Of course, you have to make sure your OUTGRID lower left lon is -179.0 (like you do in GRIB2), but it all seems to work for me.

I'm new at this, and I realize I might be missing something of historical significance, but I have a fair amount of confidence in this fix. I will add it to our CTBTO code for further testing.

[pesei]

I think this should be classified as Defect (major) rather than Support (minor), and as nobody has yet taken the ticket, I do so.