support C-grid connectivity in deseas

[aekiss]

fixes #37

[aekiss]

@ezhilsabareesh8 could you try this out with your C-grid workflow please?

You'll need to add --grid_type C to any invocations of deseas, fill_fraction, check_nonadvective and fix_nonadvective.

[ezhilsabareesh8]

@aekiss I tired deseas with --grid_type C, here is the output after deseas for the same input file. Marginal seas are removed in both workflows.

I am trying to generate the topography from scratch using both B grid and C grid type, will update here shortly.

[aekiss]

OK thanks @ezhilsabareesh8 - referring to this list the only differences I can see are the Malacca Strait is now open and there are a couple more wet points at the mouth of the Baltic.

Is your script on github somewhere I can take a look?

[ezhilsabareesh8]

referring to this list the only differences I can see are the Malacca Strait is now open and there are a couple more wet points at the mouth of the Baltic.

Thanks @aekiss, compared to OM2, I think the Mediterranean, red and black sea are removed in the current workflow, even when the --grid_type is C.

OM2 topo:

My script is here /g/data/tm70/ek4684/domain-tools/gen_topo.sh

[aekiss]

referring to this list the only differences I can see are the Malacca Strait is now open and there are a couple more wet points at the mouth of the Baltic.

Thanks @aekiss, compared to OM2, I think the Mediterranean, red and black sea are removed in the current workflow, even when the --grid_type is C.

I meant the only noticeable improvements in using grid_type C are Malacca Strait and the mouth of the Baltic. This suggest we'll need hand-edits prior to deseas.

[aekiss]

deseas is also generating a lot of spurious land points in the middle of the ocean, regardless of grid_type. Was this happening before my code changes?

[ezhilsabareesh8]

deseas is also generating a lot of spurious land points in the middle of the ocean, regardless of grid_type

I also noticed this. I think it is happening after code changes, it was better before

New (after code changes):

[aekiss]

do the spurious land points also appear using --grid_type B?

[ezhilsabareesh8]

do the spurious land points also appear using --grid_type B?

Yes but it is comparatively lesser than --grid_type C. Also I am getting this warning WARNING: could not number all the seas. Algorithm reached maximum number of iterations. when I run with --grid_type C

[ezhilsabareesh8]

I am trying to generate the topography from scratch using both B grid and C grid type, will update here shortly

When I generate everything from scratch it seems okay, there are no spurious land points. I did the following

# Interpolate topography on horizontal grid:
./topogtools gen_topo -i GEBCO_2024.nc -o topog_new.nc --hgrid ocean_hgrid.nc --tripolar --longitude-offset -100 

# Fill cells that have a sea area fraction smaller than 0.5:
./topogtools fill_fraction -i topog_new.nc -o topog_new_fillfraction.nc  --fraction 0.5 --grid_type B

# Remove seas:
./topogtools deseas -i topog_new_fillfraction.nc -o topog_new_fillfraction_0.5_deseas_grid_B.nc --grid_type B

--grid_type B

--grid_type c

[aekiss]

Great! With --grid_type C it's also retained the Med, Red & Baltic seas, but we still need to fix up the Black Sea, I guess by opening up the Bosphorus and Dardanelles.

Are Malacca and Sunda straits open? It's hard to tell.

[micaeljtoliveira]

Sorry for not looking at this earlier. Happy to see that this seems to work as expected. I'll do a more thorough review tomorrow.

[ezhilsabareesh8]

Malacca

Malacca and Sunda straits are open when --grid_type C

--grid_type B

--grid_type C

[aekiss]

Dardanelles is open with --grid_type C, but Bosphorus needs to be manually opened to connect the Black Sea

[aekiss]

I'm still getting this with --grid_type C (not B). It's apparently unable to remove the last 185 seas. Some sort of bug - I'll try to fix it.

Numbering seas
       # Iter   # Changes      # Seas
           1     1168182        2054
           2       39772         780
           3       14051         455
           4        7208         324
           5        6029         271
           6        4903         244
           7        4692         231
           8        4220         220
           9        3954         213
          10        3255         207
          11        3213         202
          12        2552         199
          13         661         197
          14         353         195
          15         283         193
          16         250         191
          17         243         190
          18         227         189
          19         168         188
          20         145         187
          21         128         186
          22         114         185
          23          54         185
          24          39         185
          25          39         185
...
         149          39         185
         150          39         185
WARNING: could not number all the seas. Algorithm reached maximum number of iterations.

[micaeljtoliveira]

It's apparently unable to remove the last 185 seas. Some sort of bug - I'll try to fix it.

@aekiss I'm afraid it's not a bug, but rather a "feature". The maximum number of diffusion iterations is hard-coded, most likely to avoid an infinite loop in case the code has some issue. I would suggest to increase that number to a much larger value.

EDIT: better yet, add this value to the CLI, so that the user can choose this value (and keep the default reasonable).

[micaeljtoliveira]

Oh, scrap that. I now see that the algorithm seems stuck. Yes, definitely a bug.

[aekiss]

I've pushed a commit to this PR that fixes the issue

[aekiss]

@micaeljtoliveira sorry for the churn, I think it's ready for review now.

I've streamlined the code, only requiring grid_type for deseas. The other methods obtain grid_type from the input .nc file (or default to B-grid if the file has no valid grid_type attribute, for backwards compatibility).

[micaeljtoliveira]

@aekiss Overall looks good.

That was quite a nasty little bug...

[micaeljtoliveira]

Wouldn't it be better to check if the value given in the command line is consistent with the value on the file? Otherwise, if I understand correctly, it's possible to change the grid type on file when running deseas, which is a bit weird.

[aekiss]

This is an intentional feature, otherwise there's no way to alter the grid type on file (which otherwise defaults to B)

[aekiss]

I'll add some words to that effect in the docs.

[aekiss]

The main purpose of grid_type in the file is to indicate to fill_fraction, fix_nonadvective and check_nonadvective what connectivity rules to use when checking whether the number of seas has changed, ie to indicate the rules used previously by deseas when lakes_removed was set to 'yes'.

[micaeljtoliveira]

Ah, I see. Yes, this definitely needs to be documented, both in the docs and in the code.

[micaeljtoliveira]

For cases like this one, I strongly prefer to use the select case construct instead of an if statement, as it's easier to extend and to catch possible mistakes (e.g., right now, if one would have this%grid_type == 'Z', the code would run without any error).

[aekiss]

thanks, good idea