The ice sheet flow equation is not properly discretized to conserved mass. Admittedly, this is a small error compared to other aspects of the EBM, but could still be improved.
-
First, the array daYBoundary should probably be renamed for clarity. Something like daDeltaYBoundary. We'll also need to introduce a second array, daDeltaYCenter. daDeltaYBoundary will be Delta y, centered on the boundaries, i.e., the difference in y coordinate between cell centers. daDeltaYCenter will then be centered on the cell centers, i.e, the difference in y coordinate between cell boundaries.
-
Second, the mass conservation part. Everywhere that terms involving Delta y^2 appear, these need to use a boundary value and a centered value. Terms that look like:
(body[iBody].daIceSheetDiff[iLat] / (body[iBody].daYBoundary[iLat] * body[iBody].daYBoundary[iLat])
should instead look like:
(body[iBody].daIceSheetDiff[iLat] / (body[iBody].daDeltaYBoundary[iLat] * body[iBody].daDeltaYCenter[iLat])
etc. I think it will be relatively straightforward to replace the second reference to the boundary value with the center value in each term. This goes for every value of daIceSheetMat and daIcePropsTmp.
The ice sheet flow equation is not properly discretized to conserved mass. Admittedly, this is a small error compared to other aspects of the EBM, but could still be improved.
First, the array daYBoundary should probably be renamed for clarity. Something like daDeltaYBoundary. We'll also need to introduce a second array, daDeltaYCenter. daDeltaYBoundary will be Delta y, centered on the boundaries, i.e., the difference in y coordinate between cell centers. daDeltaYCenter will then be centered on the cell centers, i.e, the difference in y coordinate between cell boundaries.
Second, the mass conservation part. Everywhere that terms involving Delta y^2 appear, these need to use a boundary value and a centered value. Terms that look like:
(body[iBody].daIceSheetDiff[iLat] / (body[iBody].daYBoundary[iLat] * body[iBody].daYBoundary[iLat])should instead look like:
(body[iBody].daIceSheetDiff[iLat] / (body[iBody].daDeltaYBoundary[iLat] * body[iBody].daDeltaYCenter[iLat])etc. I think it will be relatively straightforward to replace the second reference to the boundary value with the center value in each term. This goes for every value of daIceSheetMat and daIcePropsTmp.