I've been working again on grid areas for PIOMAS for the last few days. Your calculation (if you remember) was
(HTN+HUS)/2 * (HTE+HUW)/2
Which gives the closest to the main PIOMAS series you can get with the data in grid.dat.pop.
I've not been able to improve on that so have contacted Dr Zhang. It turns out that HUS and HUW are vector quantities, so instead of using those figure from grid.dat.pop we should have calculated by obtaining the other edges (HUS and HUW) from adjacent grid cells. That's the method Dr Zhang uses.
So with an array Grid(1 to 120, 1 to 360,XXX), where XXX is for example HTN or HTE the grid cell area is calculated as:
(Grid(n,m,HTN)+Grid(n-1,m,HTN))/2 * (Grid(n,m,HTE)+Grid(n,m-1,HTE))/2
Using this the error drops to about 0.01%, including the landmask in calculations drops the error far lower (can't remember the figure but it's orders of magnitude - Excel tied up right now checking results).
I'll post on my blog and include a link to the revised grid cell area mask in the next couple of days. But I'm going to try to remove some more of the error first.