Difference between revisions of "Regridding in GEOS-Chem"

From Geos-chem
Jump to: navigation, search
(Using MAP_A2A to regrid GMAO met fields)
 
(88 intermediate revisions by 3 users not shown)
Line 5: Line 5:
 
GEOS-Chem uses the regridding software <tt>MAP_A2A</tt> 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, <tt>MAP_A2A</tt> 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.   
 
GEOS-Chem uses the regridding software <tt>MAP_A2A</tt> 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, <tt>MAP_A2A</tt> 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.   
  
=== Original implemenation: Using MAP_A2A to regrid GMAO met fields ===
+
=== MAP_A2A was originally used to regrid GMAO met fields ===
  
Claire Carouge tested the <tt>MAP_A2A</tt> algorithm to ensure that it would work properly when we regrid the "raw" GMAO met field data into a form that GEOS-Chem can read.  Her comments follow:
+
Bob Yantosca obtained the <tt>MAP_A2A</tt> regridding package from S-J Lin.  We have since used <tt>MAP_A2A</tt> to regrid all of the GMAO met data products from GEOS-3 to present.
 +
 
 +
Claire Carouge tested the <tt>MAP_A2A</tt> algorithm to ensure that it would work properly when we regrid the "raw" GMAO met field data.  Her comments follow:
  
 
'''''Claire Carouge wrote:'''''
 
'''''Claire Carouge wrote:'''''
Line 19: Line 21:
 
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.
 
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.
+
The <tt>MAP_A2A</tt> 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:
+
In other regridding routines, (e.g. <tt>NCREGRID</tt>), 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 <tt>MAP_A2A</tt> regridding algorithm to answer two different questions:
  
#<p>'''''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.</p>
+
#<p>'''''Why does it look so complicated?'''''  Apparently we are using a more elaborate regridding algorithm than <tt>NCREGRID</tt>.  It probably has some nice qualities that <tt>NCREGRID</tt> hasn't.</p>
#<p>'''''Is the treatment of the poles coherent with TPCORE?'''''  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 <tt>j-1</tt> (resp. <tt>j+1</tt>) 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 <tt>j-1</tt> at South Pole you need to go southward. If you start from the cell, <tt>i=1</tt> then you arrive in the cell <tt>i=im/2+1</tt> (sphere). So in the code the values are "crossed" at the poles.  Also, they don't take the <tt>j=1</tt> grid cells but the <tt>j=2</tt> because the poles are supposed to be only one circular grid cell, so the points <tt>(i=1, j=1)</tt> and <tt>(i=im/2+1, j=1)</tt> are supposed to be the same grid cell.  So the cell <tt>j-1</tt> is the cell <tt>(i=im/2+1, j=2)</tt> ...<br>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.</p>
+
#<p>'''''Is the treatment of the poles coherent with the <tt>TPCORE</tt> 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 <tt>j-1</tt> (resp. <tt>j+1</tt>) at the South Pole (resp. North Pole). The method used in the regridding algorithm is the same than the one used in <tt>TPCORE</tt>: to access the grid cell <tt>j-1</tt> at South Pole you need to go southward. If you start from the cell, <tt>i=1</tt> then you arrive in the cell <tt>i=im/2+1</tt> (sphere). So in the code the values are "crossed" at the poles.  Also, they don't take the <tt>j=1</tt> grid cells but the <tt>j=2</tt> because the poles are supposed to be only one circular grid cell, so the points <tt>(i=1, j=1)</tt> and <tt>(i=im/2+1, j=1)</tt> are supposed to be the same grid cell.  So the cell <tt>j-1</tt> is the cell <tt>(i=im/2+1, j=2)</tt> ...<br>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.</p>
  
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.
+
So I think everything is done correctly in the regridding.  The differences we see (between <tt>NCREGRID</tt> and <tt>MAP_A2A</tt>) are explained by the differences in algorithms used but we can't say one is better than the other.
 
</blockquote>
 
</blockquote>
  
--[[User:Bmy|Bob Y.]] 15:32, 28 August 2012 (EDT)
+
--[[User:Bmy|Bob Y.]] 16:01, 28 August 2012 (EDT)
  
=== Using MAP_A2A within GEOS-Chem ===
+
=== Modification of MAP_A2A for use within GEOS-Chem ===
  
'''[mailto:COOPERM2@dal.ca Matt Cooper]''' has replaced the previous GEOS-Chem regridding routines with an algorithm adapted from the <tt>MAP_A2A</tt> 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 [http://acmg.seas.harvard.edu/geos/doc/man/appendix_2.html#A2.5 GEOS 1&deg; x 1&deg; grid], which could lead to loss of information.  The <tt>MAP_A2A</tt> algorithm regrids emissions from any arbitrary horizontal grid to the current model resolution.
+
'''[mailto:COOPERM2@dal.ca Matt Cooper]''' has replaced the previous GEOS-Chem regridding routines with an algorithm adapted from the [[Regridding in GEOS-Chem#The_MAP_A2A_algorithm|<tt>MAP_A2A</tt> 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-Chem_horizontal_grids#Generic_1_x_1_grid|GEOS 1&deg; x 1&deg; grid]], which could lead to loss of information.  The <tt>MAP_A2A</tt> algorithm regrids emissions from any arbitrary horizontal grid to the current model resolution.
  
 
<tt>MAP_A2A</tt> was introduced into GEOS-Chem v9-01-03 (and tested with 1-month benchmark [[GEOS-Chem v9-01-03 benchmark history#v9-01-03k|'''v9-01-03k''']], which was approved on 27 Apr 2012.
 
<tt>MAP_A2A</tt> was introduced into GEOS-Chem v9-01-03 (and tested with 1-month benchmark [[GEOS-Chem v9-01-03 benchmark history#v9-01-03k|'''v9-01-03k''']], which was approved on 27 Apr 2012.
Line 42: Line 44:
 
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 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, that you should be able to implement pretty easily. I basically modified the <tt>XMAP</tt> and <tt>YMAP</tt> 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 <tt>MAP_A2A</tt>.
+
I put in a pretty simple fix. I basically modified the <tt>XMAP</tt> and <tt>YMAP</tt> 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 <tt>MAP_A2A</tt>.
  
 
A couple of words of caution:
 
A couple of words of caution:
  
#My fix only works for regular lat-lon grids because it uses delta(longitude) and delta(sin(latitude)) to do the regridding
+
#<p>My fix only works for regular lat-lon grids because it uses &Delta;( longitude ) and &Delta;( SIN(latitude) ) to do the regridding</p>
#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 was 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
+
#<p>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.</p>
#You might also be interested in how I set it up so it handles nested grids correctly. I created new functions called <tt>GET_IIIPAR</tt> and <tt>GET_JJJPAR</tt> (in <tt>GeosUtil/global_grid_mod.F90</tt> that I use in regridding on a global grid, and then use <tt>GET_XOFFSET</tt> and <tt>GET_YOFFSET</tt> to cut out the nested grid portion.
+
#<p>You might also be interested in how I set it up so it handles nested grids correctly. I created new functions called <tt>GET_IIIPAR</tt> and <tt>GET_JJJPAR</tt> (in <tt>GeosUtil/global_grid_mod.F90</tt> that I use in regridding on a global grid, and then use <tt>GET_XOFFSET</tt> and <tt>GET_YOFFSET</tt> to cut out the nested grid portion.</p>
 
</blockquote>
 
</blockquote>
  
 
--[[User:Bmy|Bob Y.]] 15:50, 28 August 2012 (EDT)
 
--[[User:Bmy|Bob Y.]] 15:50, 28 August 2012 (EDT)
 
== Important Updates ==
 
 
In [[GEOS-Chem v9-01-03]], we have added the following updates to the MAP_A2A regridding scheme:
 
 
#<p>The driver routine <tt>DO_REGRID_A2A</tt> 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.</p>
 
#<p>Important DO loops in routines <tt>DO_REGRID_A2A</tt>, <tt>XMAP</tt> and <tt>YMAP</tt> (all in <tt>GeosUtil/regrid_a2a_mod.F90</tt>) have now been parallelized with [http://www.openmp.org OpenMP directives].</p>
 
#<p>The <tt>QSUM</tt> and <tt>SUM</tt> variables in routines <tt>XMAP</tt> and <tt>YMAP</tt> have been changed from <tt>REAL*4</tt> to <tt>REAL*8</tt>.  This eliminates some minor numerical noise past the 6th decimal place.
 
 
--[[User:Bmy|Bob Y.]] 15:33, 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 [[GEOS-Chem v9-01-03 benchmark history#v9-01-03n|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.
 
 
'''''[mailto:ltmurray@seas.harvard.edu Lee Murray] wrote:'''''
 
:There is definitely an error in the erroneous missing of 0.5&deg;W in the <tt>latlon_geos1x1.txt</tt> file. This led to inherent inconsistencies in how <tt>latlon_geos1x1.txt</tt> and <tt>latlon_generic.txt</tt> were indexed by GEOS-Chem, since <tt>latlon_geos1x1</tt> made up for the missing longitude edge by doubling up both 180.5&deg;W/179.5&deg;E and 179.5&deg;W/180.5&deg;E. <tt>latlon_generic.txt</tt> only repeats 180&deg;W/180&deg;E.
 
 
Several bugs were also fixed in the MAP_A2A regridding code to ensure consistency in the use of <tt>IM</tt> and <tt>JM</tt>, the number of longitude and latitude centers on the input grid.
 
 
'''''[mailto:psk9@duke.edu 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 <tt>DO_REGRID_A2A</tt> be:
 
      CALL DO_REGRID_A2A(FILENAME,IM,JM,INGRID,OUTGRID,PERAREA)
 
:where
 
::<tt>IM</tt> = number of grid boxes in E-W direction (not number of lon edges)
 
::<tt>JM</tt> = number of grid boxes in N-S direction (not number of lat edges)
 
:so that the dimension of <tt>INGRID</tt> is <tt>(IM,JM)</tt> and <tt>IM+1</tt> and <tt>JM+1</tt> are number of lon edges and number of lat edges, respectively.
 
 
:(2) That the calling sequence for <tt>MAP_A2A</tt> be
 
      CALL MAP_A2A( IM,    JM,    INLON, INSIN, INGRID,  &
 
                    IIPAR, JJPAR, LON2,  SIN2,  OUTGRID,  0, 0)
 
:where <tt>INLON</tt>, <tt>INSIN</tt>, and <tt>INGRID</tt> are dimensioned <tt>IM+1</tt>, <tt>JM+1</tt>, and <tt>(IM,JM)</tt>, respectively, and <tt>LON2</tt>, <tt>SIN2</tt>, and <tt>OUTGRID</tt> are dimensioned <tt>IIPAR+1</tt>, <tt>JJPAR+1</tt>, and <tt>(IIPAR, JJPAR)</tt>, 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 <tt>LON2</tt> and </tt>SIN2<?tt> to <tt>OUTLON</tt> and <tt>OUTSIN</tt>.
 
 
--[[User:Melissa Payer|Melissa Payer]] 11:44, 8 June 2012 (EDT)
 
 
'''''NOTE:''''' Before the release of v9-01-03 (maybe in v9-01-03o), we will convert the <tt>latlon_geos_1x1.txt</tt> and <tt>latlon_generic.1x1.txt</tt> to netCDF for better compatibility with the [[Grid-independent GEOS-Chem]] modifications.
 
 
--[[User:Bmy|Bob Y.]] 15:39, 28 August 2012 (EDT)
 
 
== Outstanding issues==
 
 
None reported at this time.
 
 
--[[User:Bmy|Bob Y.]] 15:33, 28 August 2012 (EDT)
 

Latest revision as of 21:27, 13 July 2023

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:

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

  2. 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:

  1. My fix only works for regular lat-lon grids because it uses Δ( longitude ) and Δ( SIN(latitude) ) to do the regridding

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

  3. 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)