Hello,
A geographical file (and the WRF output) includes the following attributes:
:DYN_OPT = 2 ;
:CEN_LAT = 47.58337f ;
:CEN_LON = -122.3172f ;
:TRUELAT1 = 40.f ;
:TRUELAT2 = 50.f ;
:MOAD_CEN_LAT = 45.25f ;
:STAND_LON = -120.f ;
:POLE_LAT = 90.f ;
:POLE_LON = 0.f ;
:corner_lats = 45.43848f, 49.64059f, 49.7053f, 45.49858f, 45.43816f, 49.64025f, 49.70536f, 45.49863f, 45.43172f, 49.64734f, 49.71205f, 45.4918f, 45.4314f, 49.64699f, 49.7121f, 45.49186f ;
:corner_lons = -123.7844f, -124.0875f, -120.7305f, -120.6763f, -123.7941f, -124.0979f, -120.7201f, -120.6666f, -123.784f, -124.088f, -120.7306f, -120.6762f, -123.7936f, -124.0984f, -120.7202f, -120.6665f ;
:MAP_PROJ = 1 ;
This corresponds to a nested configuration where the parent domain is centered at 45.25,-120, and the child domain at 47.58337; -122.3172. Maybe I’m interpreting this wrong, but shouldn’t MCIP locate the center of the CMAQ domain at 47.58337; -122.3172? It’s actually locating it at 47.58337; -120.0 as the former value is set by WRF_LC_REF_LAT. The latter is set at setgriddefs with the variable met_proj_lon (which is read from STAND_LON) rather than met_cen_lon (read from CEN_LON).
I was able to tweak this. However, the next problem is that MCIP set xorig & yorig to some odd values (maybe a combination between parent and chil coordinates?). Shouldn’t these values be equal to -nxdx0.5 and -nydy0.5 even for a nested domain?
Thanks.
See On The Definition of Horizontal and Vertical Grids and Coordinates for Models-3 (used by MCIP) which substantially pre-dates WRF, and which normally leads to coordinate and grid definitions from which one can use simple arithmetic to answer questions like “Is this grid a proper nest of that one?”.
The usual (and, IMNHO perverse) WRF conventions require sophisticated spherical-trigonometry calculations to make such a determination. And which make a mismash of proper mathematical/geophysical-structure analysis of the relationship between grids and coordinate systems.
Hi @a_fernandez!
Thank you for your questions.
As suggested by @cjcoats, translating the coordinate descriptions from WRF into those used by CMAQ is not straightforward. The process can be very messy because there is considerable flexibility in the way WRF domains and nests can be defined. (By contrast, the MM5 domain definitions were much more constrained, and those constraints aligned nicely with the expectations of the geographical definitions used by CMAQ.) The flexibilities in WRF domain definitions sometimes require that the translations for CMAQ break the “domain family” paradigm that the GRIDDESC file is designed to use.
For CMAQ, the domain locations within a given projection are defined using a reference point (CEN_LON and CEN_LAT) and an easting and northing (XORIG and YORIG) from that reference point to get the coordinates of (1,1) in a 1-based count of the columns and rows in the domain. Within the “domain family” system, the standard longitude in a projected domain was required to bisect the “mother of all domains” (MOAD). In WRF, the MOAD can be shifted so that the standard longitude does not bisect that domain. However, in Lambert conformal projection, XORIG and YORIG must be with respect to a point along the standard longitude. (In practice, you want this point to fall somewhere within your domain.) This is why CEN_LON is set to -120.0 (your standard longitude) in MCIP.
The setting for CEN_LAT is somewhat arbitrary, though (as I mentioned) you want it to fall somewhere in your domain. In a highly constrained case where the standard longitude bisects the domain, you could select the CEN_LAT to be the true center of the domain. In that case, the XORIG and YORIG for (1,1) would have to be factors of 0.5 delta-x (depending on whether there are odd or even number of cells in each dimension). Again, the flexibility in defining the WRF domains is such that the standard longitude may not necessarily (and often does not) fall along a column that is a multiple of 0.5 delta-x. Therefore, the values of XORIG and YORIG are less likely to be “neat” or “round” for WRF domains.
Hope this helps.
--Tanya
Hi Tanya and @cjcoats.
Thank you for the feedback. I conceptually understand what you’re saying (would need to go over the math for full comprehension). However, I need to find a practical solution and have a follow-up question. The cell centers of WRF and CMAQ (for a general domain) share latitudes and longitudes; from this point, MCIP performs whatever calculation it needs to do based on the staggered WRF grid. What it’s bugging me (based on your input) is whether MCIP will keep this condition for the nested domain. If that condition still holds true, why not treat the nested domain like just a regular domain (the CMAQ simulations are completely separate items, and BCON & ICON will take care of computing the initial/boundary conditions)? If the (WRF & CMAQ) cell centers for nested domains do not share lat & lon, I guess that MCIP will account for this.
Sorry for the rambling but I’d appreciate if you could clarify the relationship of the cell centers (WRF and CMAQ) for a nested domain.
Thank you,
Arturo
OK. Let’s talk about the most common situation – grids defined on Lambert conformal conic map projections.
What does it mean to define such a map projection? – first, select a cone with central axis the same as the N-S Earth polar axis, that intersects the Earth – at two latitudes ALPHA and BETA. On this cone, select a “true north line” at longitude GAMMA, slice the cone at diametrically-opposite longitude -GAMMA, and unroll it. This gives what we mathematicians call an “affine plane” – it doesn’t have coordinates yet. To get coordinates, you have to pick a coordinate-origin (where XCENT=YCENT=0).
Then what does it mean to have a (regular) grid in this coordinate system ? n Select a starting corner (XORIG,YORIG), a cell-size (DX,DY), and dimensions NCOLS,NROWS, and count off the cells… That’s what the mathematician-definition (and the CMAQ/Models-3 system) does.
Suppose you have two grids with different coordinate systems (but using the same cone). The only way you can tell whether one nests into the other (“everything lines up, and the cells of one are multiples of the cells of the other”) is to re-describe these grids in the same coordinate system. And someone needs to do this re-description: be glad Tanya has done this for you. Without that common description, there’s no way you can even think about boundary conditions, for example.
And that’s a messy piece of spherical trigonometry (or worse, a messy elliptic-function calculation). And it is an important question – there is a reason CMAQ-related codes check grid consistency on startup: people do make mistakes about this and the code needs to catch these mistakes.
In a common coordinate system, the nest-check is easy, and is purely normal arithmetic:
- Is the ratio of the cell-sizes a (small) integer N
- Is the (vector) distance from the nest starting-corner to the coarse starting-corner an integer multiple of the nest cell-size?
- Are the nest-grid starting-corner and ending-corner (XORIG+NCOLS*DX,YORIG+NROWS*DY) inside the coverage of the coarse grid?
Hello,
I probably did something wrong on Thursday as I tried again yesterday and the cell centers for the nested WRF domain and for the mesh generated by MCIP now overlap. Just for clarification, I can confirm that it’s a Lambert projection (“+proj=lcc +lat_1=40.0 +lat_2=50.0 +lat_0=45.25 +lon0=-120.0 +x_0=0.0 +y_0=0.0 +datum=WGS84”) and that BTRIM is set to 5. I understand that MCIP is double-checking the fitness of the grid, but I also need the parameters that MCIP is using for other tasks. The problem is twofold:
1 - How it calculates the latitudes and longitudes from a nested domain, which is now clear (the very different behaviour of latitudes and longitudes tripped me).
2 - The calculation of XORIG and YORIG, which has brought a different question. According to the Models-3 document, X_ORIG is the X coordinate of the grid origin (lower left corner of the cell at column=row=1), given in map projection units (meters, except in Lat-Lon coordinate systems) and this should be true for the parent domain. However, and for the nested domain, the south-west corner of the 1st-row & 1st column cell is located at 123.6832W-45.51632N, whose coordinates in the XY plane are -286605.78473984485, 36010.18388204688. If one computes the distance between this corner and the position given by the longitude and latitude written at GRIDDESC (120.0W,47.58337N), the distance or relative position is equal to -286605.78473984485 -222470.97698825411. This agrees (pretty closely) with the values generated by MCIP, which are -285750, -222520. As far the rounding difference, my numbers come form OSGeo (or alternatively pyproj) so they should be pretty accurate.
Thanks.
Hmmm…
One issue which you almost-certainly did not consider: WRF “believes in” a spherical Earth – and all of its calculations–dynamics, microphysics, surface exchanges, land-cover attributes, etc. – are done in terms of that, whereas you most probably did your calculations in terms of the WGS84 spheroid you mentioned. That will cause placement-differences on the order of 20 KM at continental scales, and well may be the cause of the discrepancy you see. In NCAR’s defense, doing dynamics on an ellipsoid is much more difficult than dynamics on a sphere (the whole concepts of “map scale factors” and conformal map projections break down, and …)
Many years ago, a number of us in the EPA air quality community tried to get NCAR to deal with things differently, in at least a couple of ways:
- do the spherical-approximation dynamics (as an explicit approximation), but use a proper geodetic spheroid like WGS84 for everything else (land cover, surface exchanges, etc.);
- do geographic calculations in double precision, since single precision is insufficiently accurate for grid related calculations at urban scales.
As you can see, they didn’t accept our proposal ;-(
Hi @cjcoats,
Firstly my apologies for the delay as I’ve been out of the office all week long due to some situation. Yes,I had forgotten that WRF considers a perfect sphere (and should have remembered it because I’ve been tinkering with MPAS and it’s the same situation). Anyway, changing the Lambert projection to “+proj=lcc +lat_1=40.0 +lat_2=50.0 +lat_0=45.25 +lon0=-120.0 +R=6370000 +units=m +no_defs” returns -286605.8 -222509.9, were the X-distance has not changed (makes sense) and is closer to the MCIP calculation.
Thanks.