add area to Geometry

[travissluka]

BUMP is no longer working for us after a recent update to SABER. It appears that it now expects the area field to be provided (actually it seems that it was expecting area before, but was somehow working even though we weren't providing it.)

To confirm this, I tested by adding the following quick hack to the geometry constructor

atlas::Field area = atlasFunctionSpace_->createField<double>(
                   atlas::option::levels(1) |
                   atlas::option::name("area"));
make_view<double, 2>(area).assign(110000.0 * 110000.0);
atlasFieldSet_->add(area);

And the staticbinit.x is able to run again, though with slightly different answers.

We obviously will want to fill in the area with the actual gridcell area. I'm assuming that ATLAS should have a convenient way of calculating and returning the grid cell area, but I wasn't able to find such a method at first glance, perhaps you'll have better luck @loganchen39 ?

This will have to be fixed before the other pending PR

[loganchen39]

@travissluka Understandable, I'll work on this.

[loganchen39]

@travissluka The Atlas did not provide such method. Below is the email from Willem.

Hi Ligang, Atlas did not provide a gridcell area because everybody has his own opinion on this, and several answers are possible:

• Are the area's computed from an imaginary control volume around each grid point, or do we form elements where grid points are vertices?
• Are area's computed in radians^2 or metres^2? So for the moment it is best to create a simple method on your end that computes the areas. Best, Willem

[loganchen39]

@travissluka In order to compute the grid cell area for each grid point (with lat: -89.5, -88.5, ..., 89.5; and lon: -179.5, -178.5, ..., 179.5), what is the exact definition of the grid cell here? My initial thought is for each grid point, e.g. (85.5, -177.5), the grid cell is a 1x1 degree "square" centered at the grid point, e.g. the "square" with bottom-left (85, -176) and top-right (86, -175), and we need to compute the area of this grid cell for each grid point. Is my understanding right?

[travissluka]

I agree with you that our idea of a cell is the 1x1 grid with the given lat/lon point being the center.

I talked to the JEDI team this morning about this though, and they said that, yes, area is now required, but it does not have to be an exact area, just an approximation. I don't know how much work it would be to get the dx and dy from atlas to properly calculate the grid cell area. It might be more work than it's worth at this point. Assuming we're going to be using regular lat/lon grid for the foreseeable future, adding this chunk of code to the Geometry constructor seems to work.

    double dx = 2. * M_PI * atlas::util::DatumIFS::radius()
        / atlasFunctionSpace()->grid().nxmax();
    auto lonlat_data = make_view<double, 2>(atlasFunctionSpace_->lonlat());
    atlas::Field area = atlasFunctionSpace_->createField<double>(
        atlas::option::levels(1) | atlas::option::name("area"));
    auto area_data = make_view<double, 2>(area);
    for (int i=0; i < atlasFunctionSpace_->size(); i++) {
      area_data(i, 0) = dx*dx*cos(lonlat_data(i, 1)*M_PI/180.);
    }
    atlasFieldSet_->add(area);

This only works however for a global regular latlon grid. We might want to use other grid types, including regional grids, in the future, but I think this would work for now and we can fix it properly if we need to.

[loganchen39]

@travissluka added your code and updated the branch, please merge for now.

[travissluka]

Closed by #51