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

support C-grid connectivity in deseas #39

Open
wants to merge 8 commits into
base: master
Choose a base branch
from
Open

support C-grid connectivity in deseas #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
1fb330c
add test input file
aekiss Oct 21, 2024
a2dc32b
optionally support c-grid connectivity in deseas - https://github.com…
aekiss Oct 21, 2024
f2dfe61
update README
aekiss Oct 21, 2024
43b590d
fix bug in topography_number_seas - se https://github.com/COSIMA/doma…
aekiss Oct 23, 2024
690da6b
really fix https://github.com/COSIMA/domain-tools/pull/39#issuecommen…
aekiss Oct 23, 2024
806f4df
improve test_topo.png
aekiss Oct 23, 2024
c176f26
use grid_type from file; remove --grid_type arg from fill_fraction, f…
aekiss Oct 24, 2024
fc7a8a1
Update src/topography.f90
aekiss Oct 25, 2024
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
30 changes: 21 additions & 9 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -45,11 +45,14 @@ Options
### deseas

```
usage: topogtools deseas --input <input_file> --output <output_file>
usage: topogtools deseas --input <input_file> --output <output_file> [--grid_type <type>]
```

Remove enclosed seas from <input_file> and writes the result to <output_file>.

Options
* `--grid_type <type>` Arakawa type of horizontal grid ('B' or 'C'; default is 'B')

### min_max_depth

```
Expand Down Expand Up @@ -80,13 +83,13 @@ to zero. Can produce non-advective cells and/or new seas.

```
usage: topogtools check_nonadvective --input <input_file>
[--vgrid <vgrid> --vgrid_type <type>
--potholes --coastal-cells]
[--vgrid <vgrid> --vgrid_type <type>
--potholes --coastal-cells]
```

Check for non-advective cells. There are two types of checks available: potholes
and non-advective coastal cells. Checking for non-advective coastal cells should
only be needed when using a B-grid.
Check topography for non-advective cells. There are two types of checks
available: potholes and non-advective coastal cells. B-grid connectivity rules
are assumed. Aborts if input_file is not on a B-grid.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
Expand All @@ -102,9 +105,9 @@ usage: topogtools fix_nonadvective --input <input_file> --output <output_file>
--potholes --coastal-cells]
```

Fix non-advective cells. There are two types of fixes available: potholes and
non-advective coastal cells. Fixes to non-advective coastal cells should only be
needed when using a B-grid.
Fix non-advective cells. There are two types of checks available: potholes and
non-advective coastal cells. B-grid connectivity rules are assumed. Aborts if
input_file is not on a B-grid.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
Expand Down Expand Up @@ -138,6 +141,15 @@ Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')


## test/png2nc.py

```
usage: png2nc.py
```

Converts `test_topo.png` to `test_topo.nc` for use as a test input file for `topogtools deseas`.


# Building and Installation

## General Instructions
Expand Down
99 changes: 70 additions & 29 deletions src/topography.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,8 @@ module topography
integer(int32) :: nyt = 0
! Depth variable and attributes
real(real32), allocatable :: depth(:,:)
character(len=3) :: lakes_removed = "no "
character(len=1) :: grid_type = 'B'
character(len=3) :: lakes_removed = 'no '
real(real32) :: min_depth = -1.0
integer :: min_level = 0
real(real32) :: max_depth = -1.0
Expand Down Expand Up @@ -44,9 +45,11 @@ module topography
contains

!-------------------------------------------------------------------------
type(topography_t) function topography_constructor(filename) result(topog)
type(topography_t) function topography_constructor(filename, grid_type) result(topog)
character(len=*), intent(in) :: filename
character(len=1), intent(in), optional :: grid_type

character(len=1):: file_grid_type
integer(int32) :: ncid, depth_id, frac_id, geolon_id, geolat_id, dids(2), history_len ! NetCDF ids

write(output_unit,'(3a)') "Reading topography from file '", trim(filename), "'"
Expand All @@ -72,6 +75,19 @@ type(topography_t) function topography_constructor(filename) result(topog)
call handle_error(nf90_get_att(ncid, depth_id, 'maximum_depth', topog%max_depth), isfatal=.false.)
call handle_error(nf90_get_att(ncid, depth_id, 'nonadvective_cells_removed', topog%nonadvective_cells_removed), isfatal=.false.)

if (present(grid_type)) then
topog%grid_type = grid_type ! grid_type arg overrides value in file
Copy link
Contributor

Choose a reason for hiding this comment

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

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

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'.

Copy link
Contributor

Choose a reason for hiding this comment

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

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

else
call handle_error(nf90_get_att(ncid, depth_id, 'grid_type', file_grid_type), isfatal=.false.)
if ( file_grid_type == 'B' .or. file_grid_type == 'C' ) then
topog%grid_type = file_grid_type
end if
end if
write(output_unit,*) " grid_type = ", topog%grid_type
if (all(topog%grid_type /= ['B', 'C'])) then
call handle_error(nf90_einval, .true., "grid_type must be B or C")
end if

! Get sea area fraction
call handle_error(nf90_inq_varid(ncid, 'sea_area_fraction', frac_id))
allocate(topog%frac(topog%nxt, topog%nyt))
Expand Down Expand Up @@ -112,6 +128,7 @@ subroutine topography_copy(topog_out, topog_in)
topog_out%min_level = topog_in%min_level
topog_out%max_depth = topog_in%max_depth
topog_out%nonadvective_cells_removed = topog_in%nonadvective_cells_removed
topog_out%grid_type = topog_in%grid_type

! Sea area fraction
allocate(topog_out%frac, source=topog_in%frac)
Expand Down Expand Up @@ -150,6 +167,7 @@ subroutine topography_write(this, filename)
call handle_error(nf90_def_var_fill(ncid, depth_id, 0, MISSING_VALUE))
call handle_error(nf90_put_att(ncid, depth_id, 'long_name', 'depth'))
call handle_error(nf90_put_att(ncid, depth_id, 'units', 'm'))
call handle_error(nf90_put_att(ncid, depth_id, 'grid_type', this%grid_type))
call handle_error(nf90_put_att(ncid, depth_id, 'lakes_removed', this%lakes_removed))
if (this%min_depth > 0.0) then
call handle_error(nf90_put_att(ncid, depth_id, 'minimum_depth', this%min_depth))
Expand Down Expand Up @@ -249,7 +267,6 @@ subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
silent_ = .false.
end if

! Do
land = this%nxt + this%nyt + 1
sea = land
do j = 1, this%nyt
Expand Down Expand Up @@ -313,20 +330,27 @@ subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
im = this%nxt
ip = 2
if (sea(i, j) < land .and. sea(i, j) > 0) then
sea(i,j) = min(sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
counter = counter + 1
new_sea = min(sea(im, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
end if
end if
do i = 2, this%nxt - 1
im = i - 1
ip = i + 1
if (sea(i, j) < land .and. sea(i, j) > 0) then
!get chokes
choke_east = .not. (any(sea(i:ip, jp) == land) .and. any(sea(i:ip, jm) == land))
choke_west = .not. (any(sea(im:i, jp) == land) .and. any(sea(im:i, jm) == land))
choke_south = .not. (any(sea(im, jm:j) == land) .and. any(sea(ip, jm:j) == land))
choke_north = .not. (any(sea(im, j:jp) == land) .and. any(sea(ip, j:jp) == land))
new_sea = min(minval([sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp)], &
mask=[choke_west, choke_east, choke_south, choke_north]), land)
if ( this%grid_type == 'C' ) then
Copy link
Contributor

Choose a reason for hiding this comment

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

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).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

thanks, good idea

new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
else
! get chokes, assuming B-grid connectivity rules
choke_east = .not. (any(sea(i:ip, jp) == land) .and. any(sea(i:ip, jm) == land))
choke_west = .not. (any(sea(im:i, jp) == land) .and. any(sea(im:i, jm) == land))
choke_south = .not. (any(sea(im, jm:j) == land) .and. any(sea(ip, jm:j) == land))
choke_north = .not. (any(sea(im, j:jp) == land) .and. any(sea(ip, j:jp) == land))
new_sea = min(sea(i, j), minval([sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp)], &
mask=[choke_west, choke_east, choke_south, choke_north]))
end if
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
Expand All @@ -337,8 +361,11 @@ subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
ip = 1
im = i - 1
if (sea(i, j) < land .and. sea(i, j) > 0) then
sea(i,j)=min(sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
counter = counter + 1
new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
end if
end if
end do

Expand All @@ -350,20 +377,27 @@ subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
im = this%nxt
ip = 2
if (sea(i, j) < land .and. sea(i, j) > 0) then
sea(i,j) = min(sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
counter = counter + 1
new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
end if
end if
do i = this%nxt - 1, 2, -1
im = i - 1
ip = i + 1
if (sea(i, j) < land .and. sea(i, j) > 0) then
!get chokes
choke_east = .not. (any(sea(i:ip, jp) == land) .and. any(sea(i:ip, jm) == land))
choke_west = .not. (any(sea(im:i, jp) == land) .and. any(sea(im:i, jm) == land))
choke_south = .not. (any(sea(im, jm:j) == land) .and. any(sea(ip, jm:j) == land))
choke_north = .not. (any(sea(im, j:jp) == land) .and. any(sea(ip, j:jp) == land))
new_sea = min(minval([sea(im, j), sea(ip, j), sea(i, jm), sea(i,jp)], &
mask=[choke_west, choke_east, choke_south, choke_north]), land)
if ( this%grid_type == 'C' ) then
new_sea = min(sea(i, j), sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
else
! get chokes, assuming B-grid connectivity rules
choke_east = .not. (any(sea(i:ip, jp) == land) .and. any(sea(i:ip, jm) == land))
choke_west = .not. (any(sea(im:i, jp) == land) .and. any(sea(im:i, jm) == land))
choke_south = .not. (any(sea(im, jm:j) == land) .and. any(sea(ip, jm:j) == land))
choke_north = .not. (any(sea(im, j:jp) == land) .and. any(sea(ip, j:jp) == land))
new_sea = min(sea(i, j), minval([sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp)], &
mask=[choke_west, choke_east, choke_south, choke_north]))
end if
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
Expand All @@ -374,8 +408,11 @@ subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
ip = 1
im = i - 1
if (sea(i, j) < land .and. sea(i, j) > 0) then
sea(i,j) = min(sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
counter = counter + 1
new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
end if
end if
end do

Expand Down Expand Up @@ -427,7 +464,7 @@ subroutine topography_deseas(this)
this%depth = MISSING_VALUE
this%frac = MISSING_VALUE
end where
this%lakes_removed = "yes"
this%lakes_removed = 'yes'

deallocate(sea)

Expand Down Expand Up @@ -484,7 +521,7 @@ subroutine topography_fill_fraction(this, sea_area_fraction)
call this%number_seas(number_of_seas = nseas, silent=.true.)
if (nseas > 1) then
write(output_unit,'(a)') "WARNING: new seas have been created. To fix, rerun deseas again."
this%lakes_removed = 'no'
this%lakes_removed = 'no '
end if
end if
end if
Expand All @@ -508,6 +545,10 @@ subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coast
integer(int32) :: im, ip, jm, jp
integer(int32) :: nseas

if ( this%grid_type /= 'B' ) then
call handle_error(nf90_einval, .true., "nonadvective: grid_type must be B")
end if

vgrid = vgrid_t(vgrid_file, vgrid_type)
write(output_unit,*) 'Zeta dimensions', 2*vgrid%nlevels + 1, vgrid%nlevels
allocate(zw(0:vgrid%nlevels))
Expand Down Expand Up @@ -618,11 +659,11 @@ subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coast
if (fix .and. (coastal .or. potholes)) then
this%nonadvective_cells_removed = 'yes'
if (changes_made .and. this%lakes_removed == 'yes') then
! Check if new lakes were created new lakes
! Check if new lakes were created
call this%number_seas(number_of_seas = nseas, silent=.true.)
if (nseas > 1) then
write(output_unit,'(a)') "WARNING: new seas have been created. To fix, rerun deseas again."
this%lakes_removed = 'no'
this%lakes_removed = 'no '
end if
end if
end if
Expand Down
26 changes: 15 additions & 11 deletions src/topogtools.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ program topogtools
'usage: topogtools [--help] <command> [<args>] ', &
' ', &
'Collection of tools to edit and manipulate ocean model topographies. ', &
'See ''topogtools --help <command>'' to read about a specific subcommand. ', &
'See ''topogtools --help <command>'' to read about a specific subcommand. ', &
' ', &
'Available commands: ', &
' gen_topo - Generate a new topography file from a bathymetry ', &
Expand All @@ -46,12 +46,16 @@ program topogtools
' tripolar grid ', &
' --longitude-offset <value> offset (in degrees) between the central ', &
' longitude of the ocean horizontal grid and of ', &
' the bathymetry grid (default ''0.0'') ', &
' the bathymetry grid (default ''0.0'') ', &
'']
help_deseas = [character(len=80) :: &
'usage: topogtools deseas --input <input_file> --output <output_file> ', &
' [--grid_type <type>] ', &
' ', &
'Remove enclosed seas from <input_file> and writes the result to <output_file>. ', &
' ', &
'Options ', &
' --grid_type <type> Arakawa type of horizontal grid (''B'' or ''C''; default ''B'')', &
'']
help_min_max_depth = [character(len=80) :: &
'usage: topogtools min_max_depth --input <input_file> --output <output_file> ', &
Expand All @@ -63,8 +67,8 @@ program topogtools
'Can produce non-advective cells. ', &
' ', &
'Options ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
'']
help_fill_fraction = [character(len=80) :: &
'usage: topogtools fill_fraction --input <input_file> --output <output_file> ', &
Expand All @@ -83,8 +87,8 @@ program topogtools
'needed when using a B-grid. ', &
' ', &
'Options ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
' --potholes fix potholes ', &
' --coastal-cells fix non-advective coastal cells ', &
'']
Expand All @@ -98,10 +102,10 @@ program topogtools
'only be needed when using a B-grid. ', &
' ', &
'Options ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --vgrid_type <type> can be ''mom5'' or ''mom6'' (default ''mom5'') ', &
' --potholes check for potholes ', &
' --coastal-cells check for non-advective coastal cells ', &
' --coastal-cells check for non-advective coastal cells (for B grid) ', &
'']
help_mask = [character(len=80) :: &
'usage: topogtools mask --input <input_file> --output <output_file> ', &
Expand All @@ -116,7 +120,7 @@ program topogtools
call set_args('--input:i "unset" --output:o "unset" --hgrid "ocean_hgrid.nc" --tripolar F --longitude-offset 0.0', &
help_gen_topo, version_text)
case ('deseas')
call set_args('--input:i "unset" --output:o "unset"', help_deseas, version_text)
call set_args('--input:i "unset" --output:o "unset" --grid_type "B"', help_deseas, version_text)
case ('min_max_depth')
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --level 0', &
help_min_max_depth, version_text)
Expand Down Expand Up @@ -174,7 +178,7 @@ program topogtools
call gen_topo(file_in, file_out, hgrid, lget('tripolar'), rget('longitude-offset'))

case ('deseas')
topog = topography_t(file_in)
topog = topography_t(file_in, grid_type=sget('grid_type'))
call topog%deseas()
call topog%update_history(get_mycommand())
call topog%write(file_out)
Expand Down
2 changes: 1 addition & 1 deletion src/utils.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ subroutine handle_error(error_flag, isfatal, err_string)
if (present(isfatal)) fatal = isfatal
if (error_flag /= nf90_noerr) then
if (fatal) then
write(error_unit,'(2a)') 'FATAL ERROR:', nf90_strerror(error_flag)
write(error_unit,'(2a)') 'FATAL ERROR: ', nf90_strerror(error_flag)
if (present(err_string)) write(error_unit,'(a)') trim(err_string)
error stop
end if
Expand Down
21 changes: 21 additions & 0 deletions test/png2nc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
#!/usr/bin/env python3
"""
Convert test_topo.png to test_topo.nc to use for testing deseas.
"""

from PIL import Image
import numpy as np
import xarray as xr

image = Image.open ("test_topo.png")
data = np.flipud(np.array(image)).astype("float")

coords = {"ny": range(0, data.shape[0]),
"nx": range(0, data.shape[1])}
da = xr.DataArray(data, dims=[ k for k in coords ], coords=coords, name="depth")
ds = xr.Dataset(data_vars={"depth": da,
"sea_area_fraction": da, # dummy data
"geolon_t": da, # dummy data
"geolat_t": da # dummy data
})
ds.to_netcdf("test_topo.nc")
Binary file added test/test_topo.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading