Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Do not drop indexes when computing rmax #40

Merged
merged 8 commits into from
Jun 27, 2021
Merged

Do not drop indexes when computing rmax #40

Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
34f32e4
Do not drop indexes when computing rmax
malmans2 Jun 21, 2021
9d928d9
implement #39
malmans2 Jun 21, 2021
9a382a8
replace land with 0s
malmans2 Jun 21, 2021
56f2193
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jun 21, 2021
d93b66b
handle land properly
malmans2 Jun 21, 2021
f0ae17f
Merge branch 'pyNEMO:main' into rmax_indexes
malmans2 Jun 22, 2021
944226c
use max rather than mean and use xarray max instead of np maximum
malmans2 Jun 25, 2021
6e8ed5e
same as pyroms
malmans2 Jun 26, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 17 additions & 12 deletions pydomcfg/tests/bathymetry.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,29 +147,34 @@ def _calc_rmax(depth):
Parameters
----------
depth: float
Bottom depth (units: m).
Bottom depth (units: m).

Returns
-------
rmax: float
Slope steepness value (units: None)
Slope steepness value (units: None)
"""
depth = depth.reset_index(list(depth.dims))

# Replace land with zeros
depth = depth.where(depth > 0, 0)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we really need this here? I think when bathy or envelope are passed to this function they should have already zeros for land

Copy link
Member Author

@malmans2 malmans2 Jun 21, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happen if the user bathymetry has negative values or NaN on land?
Should we remove it and clarify in the documentation that Bathymetry must be >= 0?

It's not a bad idea to get rid of this because it is potentially a very expensive but also useless operation...

Copy link
Contributor

@oceandie oceandie Jun 21, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would say yes, since I think we should write a function that put all land values to zero when interpolating the input bathymetry onto the model grid. So the model bathymetry will be positive defined by definition, with zeros for land (also I think NEMO will fail if this conditions are not met, but not sure about this last point)

Copy link
Contributor

@jdha jdha Jun 21, 2021 •
edited by malmans2
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I've been dipping in and out of this today (so no continuity in my thought process). @oceandie re NEMO will fail, NEMO as far as I understand doesn't use bathymetry field, just scale factors and wet cells [i.e. if top_level is equal to zero it's land]. So in that sense we can define land values in bathymetry as we see fit. I'll look over the python again tomorrow (and you may have already included some of this), but the key user defined vars are rn_sbot_min (for ln_sco) and rn_hmin (for ln_zco, ln_zps; which take different meanings depending on whether they are -ve or +ve: min number of levels or min depth). So when calculating rmax this user choices will affect the outcome (although, not in the case of tests).

from the Fortran (where the bathy is updated):

IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
         ELSE                          ;   ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 )  ! from a depth
         ENDIF
         zhmin = gdepw_1d(ik+1)                                                        ! minimum depth = ik+1 w-levels 
         WHERE( bathy(:,:) <= 0._wp )   ;   bathy(:,:) = 0._wp                         ! min=0     over the lands
         ELSE WHERE ( risfdep == 0._wp );   bathy(:,:) = MAX(  zhmin , bathy(:,:)  )   ! min=zhmin over the oceans
         END WHERE
         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
ENDIF

and

bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )

         DO jj = 1, jpj
            DO ji = 1, jpi
              IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
            END DO
         END DO

note there is also a rn_sbot_max

Apologies if I'm going over old ground...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

HI @jdha, yes sorry, I was still thinking about NEMO3.6 or DOMCFG, that we are actually trying to replicate :)
Regarding all the params and bits of code you mentioned, yes I started to consider them - soon I will push where I am now in sco_dev , then any feedback is more than welcome ;)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, this is helpful! (maybe open an issue about this as this comment will get kind of lost when we merge this PR?)

This PR should be just the starting point, and _calc_rmax will become the general function to compute rmax. Once we merge this PR, I think @oceandie will merge main into #33 and will move the whole _calc_rmax function into utils. I guess vcoord-specific parameters such as rn_sbot_min will be handled internally by our domzgr classes (but all classes will make use of the same _calc_rmax function).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @malmans2 and @jdha , how can I see the original code for calc_rmax by @jdha (the one using numpy) ?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR started from the xarray version. You'd have to look in James' PR (already merged). Should be this one: #17

I'd go there and look at the appropriate commit (e.g., click on commits and select the first commit).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks!


both_rmax = []
for dim in depth.dims:

# (H[0] - H[1]) / (H[0] + H[1])
depth_diff = depth.diff(dim)
depth_rolling_sum = depth.rolling({dim: 2}).sum().dropna(dim)
rmax = depth_diff / depth_rolling_sum
rolled = depth.rolling({dim: 2}).construct("tmp_dim")

# |(H[0] - H[1])| / (H[0] + H[1])
# First value is NaN
diff = rolled.diff("tmp_dim").squeeze("tmp_dim")
rmax = np.abs(diff) / rolled.sum("tmp_dim")

# (R[0] + R[1]) / 2
rmax = rmax.rolling({dim: 2}).mean().dropna(dim)
# (rmax[0] + rmax[1]) / 2
# First two values are NaN
rmax = rmax.rolling({dim: 2}).mean(skipna=True)
malmans2 marked this conversation as resolved.
Show resolved Hide resolved

# Fill first row and column
rmax = rmax.pad({dim: (1, 1)}, constant_values=0)
# First and last values are zero
rmax = rmax.shift({dim: -1}).fillna(0)

both_rmax.append(np.abs(rmax))
both_rmax.append(rmax)

return np.maximum(*both_rmax)