COSIMA / regional-mom6

Automatic generation of regional configurations of the Modular Ocean Model 6 (MOM6) in Python
https://regional-mom6.readthedocs.io/en/latest
MIT License
24 stars 13 forks source link

Clean up construction of vertical grid + include surface `zi = 0` #107

Closed navidcy closed 8 months ago

navidcy commented 9 months ago

I don't understand why we use self.vlayers + 1 as input in the dz_hyperbolictan here:

https://github.com/COSIMA/regional-mom6/blob/3e9844286db015e17940bbbacc2c35f7ddc82e23/regional_mom6/regional_mom6.py#L705-L713

The dz_hyperbolictan method actually spits out the layer thicknesses, thus it will spit out an array of length self.vlayers if we provide self.vlayers as its first argument. Why do we give +1?

I think the code in the above-mentioned lines should be replaced with:

        thicknesses = dz_hyperbolictan(self.vlayers, self.dz_ratio, self.depth)

        zi = np.cumsum(thicknesses)
        zi = np.insert(zi, 0, 0.0)

        vcoord = xr.Dataset(
            {
                "zi": ("zi", zi),
                "zl": ("zl", zi[0:-1] + thickness/2),
            }
        )

According to what I suggest, thicknesses is an array of size self.vlayers. Then zi are the interfaces height and is an array of size self.vlayers + 1. Last, zl, are the heights of the midpoints of each layer and is an array of size self.vlayers.

navidcy commented 8 months ago

Turns out the current implementation for constructing zi precludes the surface zi = 0. E.g.,

>>> import numpy as np
>>> thicknesses = np.array([10, 20, 30])
>>> thicknesses
array([10, 20, 30])
>>> np.cumsum(thicknesses)
array([10, 30, 60])

See also https://github.com/COSIMA/regional-mom6/pull/106#issuecomment-1956101121