The estimation of the domain size in the interpolation has a small bug

[jpolcher]

The routine which determines the region over which we have to interpolate (aggregate_vec and aggregate_2d), has a small bug in estimating the maximal size of the domain. It is less critical in the second case.

The problem is in the calculation of domain_maxlon and domain_minlon where we need to add or subtract half of the grid-size (resolution). The coslat used in this computation should be done with the largest latitude. So the code should look something like this :

ff = MAXLOC(ABS(lalo(:,1)))
ll = MINLOC(lalo(:,2))
coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
domain_minlon = lalo(ll(1),2) - resolution(ll(1),1)/(2.0*coslat)
ll = MAXLOC(lalo(:,2))
coslat = MAX(COS(lalo(ff(1),1) * pi/180. ), mincos )*pi/180. * R_Earth
domain_maxlon = lalo(ll(1),2) + resolution(ll(1),1)/(2.0*coslat)

This only affects people working on certain domains. While we are at it we could also fix it for aggregate_2d.