Regridding in GEOS-Chem
On this page, we shall describe how GEOS-Chem regrids data from one resolution to another.
The MAP_A2A algorithm
GEOS-Chem uses the regridding software MAP_A2A from S-J Lin. MAP_A2A is an area-preserving mapping scheme. For example, if you have a quantity such as kg/m2/s or W/m2, MAP_A2A will multiply by the area on the input grid, then regrid, and divide by the area on the output grid, such that the total quantity is preserved.
MAP_A2A was originally used to regrid GMAO met fields
Bob Yantosca obtained the MAP_A2A regridding package from S-J Lin. We have since used MAP_A2A to regrid all of the GMAO met data products from GEOS-3 to present.
Claire Carouge tested the MAP_A2A algorithm to ensure that it would work properly when we regrid the "raw" GMAO met field data. Her comments follow:
Claire Carouge wrote:
There are two types of physical quantities:
- Intensive
- A quantity whose value doesn't change with the grid cell size.
- Extensive
- A quantity whose value changes with the grid cell size.
For example, if you have a uniformly distributed atmosphere and a grid and you increase the resolution by 2 (grid cells are 4 times smaller), then temperature, velocities, concentrations, pressure won't change at each grid cell but the mass will be smaller for each grid cell. So mass is an extensive variable and the others are intensives.
The MAP_A2A algorithm is set up to regrid extensive quantities. So in order for us to use it to regrid winds, we must first multiply the winds by the pressure on the input grid in order to create an extensive uantity (e.g. a "mass flux"). Then we must divide by the pressure on the output grid to convert back to a wind.
In other regridding routines, (e.g. NCREGRID), you can specify if the quantity you are regridding is intensive or extensive, and it will do the regridding accordingly. I've been looking at the MAP_A2A regridding algorithm to answer two different questions:
Why does it look so complicated? Apparently we are using a more elaborate regridding algorithm than NCREGRID. It probably has some nice qualities that NCREGRID hasn't.
Is the treatment of the poles coherent with the TPCORE advection scheme? In the regridding algorithm, there is some special treatment on the poles. This is due to the fact, the algorithm is based on the calculation of the slopes between neighbor grid cells. So we need different values from neighbor grid cells. The problem with the poles is to define the grid cells j-1 (resp. j+1) at the South Pole (resp. North Pole). The method used in the regridding algorithm is the same than the one used in TPCORE: to access the grid cell j-1 at South Pole you need to go southward. If you start from the cell, i=1 then you arrive in the cell i=im/2+1 (sphere). So in the code the values are "crossed" at the poles. Also, they don't take the j=1 grid cells but the j=2 because the poles are supposed to be only one circular grid cell, so the points (i=1, j=1) and (i=im/2+1, j=1) are supposed to be the same grid cell. So the cell j-1 is the cell (i=im/2+1, j=2) ...
In addition, the wind values for the pole grid cells are averaged at the end. This is also coherent with tpcore as we consider the poles are one circular grid cell.So I think everything is done correctly in the regridding. The differences we see (between NCREGRID and MAP_A2A) are explained by the differences in algorithms used but we can't say one is better than the other.
--Bob Y. 16:01, 28 August 2012 (EDT)
Modification of MAP_A2A for use within GEOS-Chem
Matt Cooper has replaced the previous GEOS-Chem regridding routines with an algorithm adapted from the MAP_A2A regridding package (developed by S.-J. Lin and refined by Bob Yantosca). The previous GEOS-Chem regridding routines sometimes involved two separate regridding processes, passing through the GEOS 1° x 1° grid, which could lead to loss of information. The MAP_A2A algorithm regrids emissions from any arbitrary horizontal grid to the current model resolution.
MAP_A2A was introduced into GEOS-Chem v9-01-03 (and tested with 1-month benchmark v9-01-03k, which was approved on 27 Apr 2012.
Prasad Kasibhatla wrote:
I discovered that the MAP_A2A gridding algorithm was smearing emissions incorrectly (though not in a major way) - seems to arise from mapping based on some sort of piecewise polynomial interpolation between grids, rather than simply using area-overlap fractions to apportion grid-average fluxes from old to new grid.
I put in a pretty simple fix. I basically modified the XMAP and YMAP routines and got rid of all the other subroutine calls - the calling sequence is the same so I suspect it should be pretty easy for you to modify MAP_A2A.
A couple of words of caution:
My fix only works for regular lat-lon grids because it uses Δ( longitude ) and Δ( SIN(latitude) ) to do the regridding
My fix is exact for cases where you are regridding grid-average quantities (eg gridded emission fields), but I would do it differenty if I were regridding quantities that implicitly had sub-grid variation (eg if one assumes that a variable represents the center of a grid box and that there is a linear variation from one grid box center to another.
You might also be interested in how I set it up so it handles nested grids correctly. I created new functions called GET_IIIPAR and GET_JJJPAR (in GeosUtil/global_grid_mod.F90 that I use in regridding on a global grid, and then use GET_XOFFSET and GET_YOFFSET to cut out the nested grid portion.
--Bob Y. 15:50, 28 August 2012 (EDT)
Previous issues that have are now resolved
Bug fixes for the MAP_A2A regridding algorithm
This update was tested in the 1-month benchmark simulation v9-01-03n and approved on 08 Jun 2012.
Lee Murray found a typo in the latlon_geos1x1.txt file used for the MAP_A2A regridding algorithm.
Lee Murray wrote:
- There is definitely an error in the erroneous missing of 0.5°W in the latlon_geos1x1.txt file. This led to inherent inconsistencies in how latlon_geos1x1.txt and latlon_generic.txt were indexed by GEOS-Chem, since latlon_geos1x1 made up for the missing longitude edge by doubling up both 180.5°W/179.5°E and 179.5°W/180.5°E. latlon_generic.txt only repeats 180°W/180°E.
Several bugs were also fixed in the MAP_A2A regridding code to ensure consistency in the use of IM and JM, the number of longitude and latitude centers on the input grid.
Prasad Kasibhatla wrote:
- I think we are all on the same page, but just to make sure here is what I suggest we agree on:
- (1) That the calling sequence for DO_REGRID_A2A be:
CALL DO_REGRID_A2A(FILENAME,IM,JM,INGRID,OUTGRID,PERAREA)
- where
- IM = number of grid boxes in E-W direction (not number of lon edges)
- JM = number of grid boxes in N-S direction (not number of lat edges)
- so that the dimension of INGRID is (IM,JM) and IM+1 and JM+1 are number of lon edges and number of lat edges, respectively.
- (2) That the calling sequence for MAP_A2A be
CALL MAP_A2A( IM, JM, INLON, INSIN, INGRID, &
IIPAR, JJPAR, LON2, SIN2, OUTGRID, 0, 0)
- where INLON, INSIN, and INGRID are dimensioned IM+1, JM+1, and (IM,JM), respectively, and LON2, SIN2, and OUTGRID are dimensioned IIPAR+1, JJPAR+1, and (IIPAR, JJPAR), respectively.
- So, for example, when the input grid is a generic 1x1 grid, and the output grid is the geos 1x1 grid:
IM=360
JM=180
INLON(1)=-180
INLON(361)=180
INSIN(1)=-1
INSIN(181)=1
LON2(1)=-180.5
LON2(361)=179.5
SIN2(1)=-1
SIN2(181)=1
One last suggestion - if possible, change variable names LON2 and SIN2 to OUTLON and OUTSIN.
--Melissa Payer 11:44, 8 June 2012 (EDT)
Bug fixes for regional emission masks
This update was tested in the 1-month benchmark simulation v9-01-03q.
Lee Murray wrote:
- I believe that there is a bug in the masks for the regional overwrites. I added the following code to the end of READ_EUROPE_MASK in emep_mod.F:
! Plot map
DO J=JJPAR,1,-1
WRITE(6,'(10F5.1)') EUROPE_MASK(30:39,J)
ENDDO
- to output a map of EUROPE_MASK values, after it was re-gridded from 1x1 to the GEOS5 4x5 grid.
- READ_EUROPE_MASK: Reading /Data/input/geos/GEOS_1x1/EMEP_200911/EMEP_mask.geos.1x1
... etc ...
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 7.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 12.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 10.5 19.5 20.0 20.0 20.0 20.0 20.0 20.0 20.0
0.0 0.0 2.0 3.5 9.5 15.5 20.0 20.0 20.0 20.0
0.0 0.0 0.0 0.0 0.0 0.0 6.5 17.5 20.0 20.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 9.5 12.5
... etc ...
- You can see that many of the values exceed 1.0, so in many grid boxes, the emission rates for Europe would be multiplied by a factor of 20 when the mask is applied. I have also tested and get similar results for the Canada and US masks. I imagine the same holds for China and Mexico.
- We can get masks in the correct range of 0 to 1 if we use
CALL DO_REGRID_A2A( LLFILENAME, I1x1, J1x1,
& INGRID, OUTGRID, PERAREA=0 )
EUROPE_MASK = OUTGRID
- instead of
CALL DO_REGRID_A2A( LLFILENAME, I1x1, J1x1,
& INGRID, OUTGRID, PERAREA=1 )
EUROPE_MASK = OUTGRID
- I think that this is because although the masks themselves are unitless, they are actually m2 region / m2 box, they are already in a "per unit area" and don't need to be converted in DO_REGRID_A2A.
- Also, I suspect that once they are re-gridded, all values greater than 0 need to be set to 1 to be consistent with how the masks are applied to the regridded emissions. But someone who is more familiar with the emissions should verify that.
--Bob Y. 15:14, 20 December 2012 (EST)
Updates to improve performance
This update was tested in the 1-month benchmark simulation v9-01-03r and approved on 28 Aug 2012.
In GEOS-Chem v9-01-03, we have added the following updates to the MAP_A2A regridding scheme:
The driver routine DO_REGRID_A2A now reads the grid input information (i.e. longitude edges and sines of latitude edges) from netCDF files. This replaces the original implementation, which read from ASCII files.
Important DO loops in routines DO_REGRID_A2A, XMAP and YMAP (all in GeosUtil/regrid_a2a_mod.F90) have now been parallelized with OpenMP directives.
The QSUM and SUM variables in routines XMAP and YMAP have been changed from REAL*4 to REAL*8. This eliminates some minor numerical noise past the 6th decimal place.
--Bob Y. 15:33, 28 August 2012 (EDT)
Fix for regridding error in offline CO2 simulation
We have corrected this issue as a post-release patch on 22 Sep 2012. This issue does not affect the full-chemistry simulations, only the CO2 simulation.
Ray Nassar wrote:
- There is a bug in v9-01-03 for the CO2 simulation that causes a crash. It relates to the MAP_A2A regridding and the fact that the net terrestrial exchange file is named with geos.1x1 (360x181) although it is actually generic 1x1 (360x180) and it is trying to be read as geos.1x1.
Please see this wiki post on our CO2 simulation wiki page for a full description of the issue.
--Bob Y. 15:21, 20 December 2012 (EST)
Outstanding issues that have not been resolved
Bug in regridding of anthropogenic emissions
NOTE: This issue was resolved in GEOS-Chem v9-02e, which was approved on 09 Jan 2013.
It appears that there is a bug in the regridding of the NEI2005 anthropogenic emissions and the David Streets Anthropogenic Emissions. In some instances, the regridding is assuming that the emissions data have units of per unit area (i.e. per m2 or cm2) when they do not (or vice versa). We will rectify this situation shortly.
Christoph Keller wrote:
- Jenny Fisher asked me why the PERAREA flag is switched on in all the READ_STREETS calls in GeosCore/streets_anthro_mod.F except for NOx and CO? Those flags determine whether the fields are normalized by their area during the regridding or not, and in my opinion the flag should be switched on (I.e. PERAREA=1) for all compounds since all emissions are in kg/yr.
- I'm concerned that we have the wrong regridding for STREETS NOx and CO emissions. E.g. On line 686 we have:
CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
& TAU, TEMP, PERAREA=0 )
- However, these are emissions in kg/yr, so I think we should have:
CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
& TAU, TEMP, PERAREA=1 )
- This only affects CO and NOx, the other species use the PERAREA=1 flag.
- I have made a test run where I wrote out the total Streets emissions before and after regridding, this is what I get:
Original version (PERAREA = 0):
- EMISS_STREETS_ANTHRO: Reading /as/data/geos/GEOS_1x1/Streets_200812/Streets_NOx_ind_2006.generic.1x1
STREETS sum before regridding (tracer 1): 8483151381.67771
STREETS sum after regridding (tracer 1): 423913738.256104
With PERAREA = 1:
- EMISS_STREETS_ANTHRO: Reading /as/data/geos/GEOS_1x1/Streets_200812/Streets_NOx_ind_2006.generic.1x1
STREETS sum before regridding (tracer 1): 8483151381.67771
STREETS sum after regridding (tracer 1): 8483151381.67771
- I would say the latter makes much more sense! I think this is a bug which has been introduced with the new regridding package in v9-01-03, so this is only very recent.
- The PERAREA flag can be a bit confusing, though, as PERAREA can be interpreted as "is the field I want to regrid in units of per area", whereas in fact it means "do you want to convert to per area before regridding".
Melissa Payer wrote:
- After looking carefully through GeosCore/streets_anthro_mod.F, we found there was a second error in the code for regridding of the monthly variability factors which essentially negates the bug in NOx emissions.
- For regridding of NOx emissions, we call routine READ_STREETS with:
! File name
IF ( IS_2006 ) THEN
FILENAME = TRIM( STREETS_DIR ) // 'Streets_NOx_'//
& SOURCES( NSOURCE ) // '_2006.generic.1x1'
TAU = TAU2006
ELSE
FILENAME = TRIM( STREETS_DIR ) //
& 'Streets_NOx_FF_2000.generic.1x1'
TAU = TAU2000
ENDIF
! Read data
CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
& TAU, TEMP, PERAREA=0 )
- As Christoph pointed out, this should be:
! Read data
CALL READ_STREETS( FILENAME, 'ANTHSRCE', 1,
& TAU, TEMP, PERAREA=1 )
- We then apply seasonal variation to both NOx and CO:
! Seasonal Variation for NOx
!--------------------------
! Monthly Variability for any year. Variability has
! been obtained from 2004 1x1 data, for two cases:
! FF seasonality for 2000, and TOTAL (BF+FF) for 2006
! inventories (phs, 12/2/08)
IF ( IS_2006 ) THEN
FILENAME = TRIM( STREETS_DIR ) //
& 'Streets_2004_NOx_MonthFctr_total.generic.1x1'
ELSE
! redefine the entire path here
FILENAME = TRIM( DATA_DIR_1x1 ) // 'Streets_200812/'
& // 'Streets_2004_NOx_MonthFctr_FF.generic.1x1'
ENDIF
CALL READ_STREETS( FILENAME, 'RATIO-2D', 71,
$ TAUMONTH_2004, TEMP, PERAREA=1 )
NOX = NOX * TEMP
- However, since the monthly scale factors are unitless, the regridding routine should actually be called with
CALL READ_STREETS( FILENAME, 'RATIO-2D', 71,
$ TAUMONTH_2004, TEMP, PERAREA=0 )
- When we fix both of these bugs, they essentially cancel one another in the NOx totals. This is the reason that we did not see the bug in the 1-month benchmark output. The attached plots show that with both bug fixes applied surface NOx changes by a maximum of +/- 0.0522 ppb over Asia.
- If we only fix the bug in regridding of NOx emissions, we get an increase of up to 300 ppb over China, which is clearly unrealistic. It is therefore obvious that we need to fix both the regridding of NOx emissions and of monthly scale factors.
- These bugs go as far back as GEOS-Chem v9-01-03k, when the MAP_A2A regridding algorithm was implemented. However, we did not pick up on it because there was very little change in NOx emissions due to the cancelling effect of the two bugs. The source of the error likely stems from the confusion of the argument PERAREA.
- PERAREA = 1 is used for units of mass (kg/yr, etc.). The array will be divided by unit area before regridding.
- PERAREA = 0 is used for units of concentration (unitless, v/v, molec/cm2/s, etc).
- We will fix the bugs in GeosCore/streets_anthro_mod.F in the next 1-month benchmark (GEOS-Chem v9-02e). I will also take another thorough look through all of the regridding calls in the code to make sure that we consistently and properly call the regridding routines. This will likely involve renaming the argument PERAREA throughout the code to avoid confusion as to what it represents.
--Bob Y. 16:11, 19 December 2012 (EST)
Additional bug fixes for MAP_A2A regridding algorithm
NOTE: This issue was resolved in GEOS-Chem v9-02e, which was approved on 09 Jan 2013.
The source code was carefully combed to ensure that PERAREA is assigned the proper value when calling the MAP_A2A regridding routines. The following bugs were found and corrected:
- In subroutine SEASCL_EDGAR_SHIP_SO2 (edgar_mod.F), PERAREA should be 0 because scale factors are unitless.
- In subroutine GET_AEF_05x0666 (megan_mod.F), PERAREA should be 0 for all calls to DO_REGRID_A2A because AEFs are in μg C/m2/hr and do not need to be converted to units per area. This fix ensures that the nested grid code is consistent with the code in GET_AEF for global grids.
- In subroutine GET_NEI99_WKSCALE_05x0666 (nei2005_anthro_mod.F), PERAREA should be 0 because scale factors are unitless.
- In subroutine READ_BIOFUEL_SO2 (sulfate_mod.F), PERAREA should be 1 for historical emissions because data are in kg/yr.
To avoid confusion in assigning a value to PERAREA, we have renamed the parameter to IS_MASS throughout the standard source code, where:
- IS_MASS = 1 is used for units of mass (kg/yr, etc.). The array will be divided by unit area before regridding.
- IS_MASS = 0 is used for units of concentration (unitless, v/v, molec/cm2/s, etc).
--Melissa Payer 14:44, 7 January 2013 (EST)
Bug in grid_mod.F90
Christoph Keller wrote:
- I found a small bug in grid_mod.F: The sine of the last latitude edge (i.e. the north pole) is not calculated, which can cause problems for the regridding.
- I just added the following lines to the (last) latitude edge calculation (COMPUTE_GRID in grid_mod.F, L235ff):
!%%%%%% LATITUDE EDGES (the last edge) %%%%%%%%%%%%%%%%%%%%%%%%
! Test for North Pole
IF ( J2 == JNP ) THEN
! North pole case (global grids only)
DO I = I1, I2
YEDGE (I,JNP+1,L) = +90d0
YEDGE_R(I,JNP+1,L) = YEDGE(I,JNP+1,L) * PI_180
! Also compute sine of last latitude edge! (ckeller, 02/13/12)
YSIN(I,J2+1,L) = SIN ( YEDGE_R(I,J2+1,L) )
ENDDO
ELSE
! No north pole (nested grids only)
DO I = I1, I2
YEDGE (I,J2+1,L) = YEDGE(I,J2,L ) + DLAT(I,J2,L)
YEDGE_R(I,J2+1,L) = YEDGE(I,J2+1,L) * PI_180
! Also compute sine of last latitude edge! (ckeller, 02/13/12)
YSIN(I,J2+1,L) = SIN ( YEDGE_R(I,J2+1,L) )
ENDDO
ENDIF
- Looks like this fixes the problem.
--Melissa Payer 16:41, 13 February 2013 (EST)