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

support C-grid connectivity in deseas #39

Merged
merged 12 commits into from
Nov 4, 2024
Merged

support C-grid connectivity in deseas #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 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
afe86c7
improve grid_type and other docs in README
aekiss Oct 30, 2024
72ab723
improve grid_type docs in code and README - see https://github.com/CO…
aekiss Nov 1, 2024
29cb4f5
review response - use select case (this%grid_type)
aekiss Nov 1, 2024
6c49dd9
stop if invalid grid_type in topography_number_seas
aekiss Nov 3, 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
71 changes: 58 additions & 13 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,22 @@ Available commands:
* `gen_topo` - Generate a new topography file from a bathymetry
* `deseas` - Remove enclosed seas
* `min_max_depth` - Set minimum and maximum depth
* `fill_fraction` - Set cells with unsufficient ocean fraction to land
* `check_nonadvective` - Check for non-advective cells
* `fix_nonadvective` - Fix non-advective cells
* `mask` - Generate mask

All commands other than `gen_topo` will add the following attributes to the
`depth` variable in `<output_file>`. Except where explained below, these
attributes will be copied from `<input_file>` if present, or otherwise set to
the default values in this table.

| Attribute | Default | Description |
|---|:-:|---|
| `grid_type` | 'B' | Arakawa grid type (B or C); determines advective connectivity between cells when counting seas in `deseas`, `fill_fraction`, `check_nonadvective` and `fix_nonadvective` |
| `lakes_removed` | 'no ' | Whether all isolated water bodies have been replaced by land |
| `nonadvective_cells_removed` | 'yes' | Whether `fix_nonadvective` has been applied (BUG: always 'yes') |

### gen_topo

```
Expand All @@ -45,10 +57,21 @@ 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>.
Remove enclosed seas from `<input_file>` and write the result to
`<output_file>`, using advective connectivity rules set by `--grid_type` (B or
C); if `--grid_type` is not specified, the `grid_type` attribute of `depth` in
`<input_file>` is used, defaulting to B if that attribute is absent.

Sets the `lakes_removed` attribute in `<output_file>` to 'yes'. If `--grid_type`
is specified, this sets the `grid_type` attribute in `<output_file>`.

Also creates a `sea_num.nc` file showing how the seas are numbered.

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

### min_max_depth

Expand All @@ -58,8 +81,11 @@ usage: topogtools min_max_depth --input <input_file> --output <output_file>
[--vgrid <vgrid> --vgrid_type <type>]
```

Set minimum depth to the depth at a specified level and set maximum depth to
Set minimum depth to the depth at level `<level>` and set maximum depth to
deepest in `<vgrid>`. `<level>` is the minimum number of depth levels (e.g. 4).
These values are recorded in the `minimum_levels`, `minimum_depth` and
`maximum_depth` attributes of `depth` in `<output_file>`.

Can produce non-advective cells.

Options
Expand All @@ -73,20 +99,26 @@ usage: topogtools fill_fraction --input <input_file> --output <output_file>
--fraction <frac>
```

Cells with a fraction of sea area smaller than <frac> will have their depth set
to zero. Can produce non-advective cells and/or new seas.
Cells with a fraction of sea area smaller than `<frac>` will have their depth
set to zero.

Can produce non-advective cells.

Can also produce new isolated seas - if this is the case, a warning is given and
the `lakes_removed` attribute of `depth` is set to 'no '.

### check_nonadvective

```
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 `grid_type` attribute of `depth` in `<input_file>` is
present and not 'B'.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
Expand All @@ -102,9 +134,13 @@ 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. If either is used, the `nonadvective_cells_removed`
attribute of `depth` is set to 'yes'. B-grid connectivity rules are assumed.
Aborts if `grid_type` attribute of `depth` in `<input_file>` is present and not 'B'.

Can produce new isolated seas. If this is the case, a warning is given and the
`lakes_removed` attribute of `depth` is set to 'no '.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
Expand Down Expand Up @@ -138,6 +174,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
105 changes: 76 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
micaeljtoliveira marked this conversation as resolved.
Show resolved Hide resolved
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,30 @@ 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)
select case (this%grid_type)
case ('B')
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]))
case ('C')
new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
case default
write(error_unit, '(a)') "topogtools: grid_type must be B or C"
error stop
end select
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
Expand All @@ -337,8 +364,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 +380,30 @@ 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)
select case (this%grid_type)
case ('B')
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]))
case ('C')
new_sea = min(sea(i, j), sea(im, j), sea(ip, j), sea(i, jm), sea(i, jp))
case default
write(error_unit, '(a)') "topogtools: grid_type must be B or C"
error stop
end select
if (sea(i, j) /= new_sea) then
sea(i, j) = new_sea
counter = counter + 1
Expand All @@ -374,8 +414,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 +470,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 +527,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 +551,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 +665,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
Loading
Loading