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

Improve handling of vertical grid, adding support for MOM6 grids #35

Merged
merged 5 commits into from
Sep 18, 2024
Merged

Improve handling of vertical grid, adding support for MOM6 grids #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1ef3716
Add a class to handle vertical grids.
micaeljtoliveira Sep 17, 2024
3da1007
Store super and normal vertical grid in vgrid type, instead of just s…
micaeljtoliveira Sep 17, 2024
c5919ab
Use new vgrid class in topography tools.
micaeljtoliveira Sep 17, 2024
66db2c7
Allow for MOM6 vertical grids through a new command line option. Thes…
micaeljtoliveira Sep 17, 2024
f651e5a
Check if number of points in MOM5 vertical grid is odd, stop otherwise.
micaeljtoliveira Sep 17, 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
1 change: 1 addition & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ add_library(M_CLI2

add_library(topography
src/utils.f90
src/vgrid.f90
src/topography.f90
src/gen_topo.f90
src/kdtree2.f90
Expand Down
31 changes: 20 additions & 11 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,10 @@ Code and tools to edit and manipulate ocean model grids and topographies.

Below is a list of included tools and short documentation for each.

**Note:** in all cases `<vgrid>` is assumed to be a MOM5 vertical grid file (with $2n+1$ values for an $n$-level model). Using a MOM6 `<vgrid>` file (with $n+1$ values) will produce incorrect results.
**Note:** these tools support two types of vertical grids: MOM5 grids (with
$2n+1$ values for an $n$-level model) and MOM6 grids (with $n+1$ values). It is
important to select the correct type or the tools will produce incorrect
results.

## topogtools (Russ' Fortran tools)

Expand Down Expand Up @@ -51,15 +54,17 @@ Remove enclosed seas from <input_file> and writes the result to <output_file>.

```
usage: topogtools min_max_depth --input <input_file> --output <output_file>
--level <level> [--vgrid <vgrid>]
--level <level>
[--vgrid <vgrid> --vgrid_type <type>]
```

Set minimum depth to the depth at a specified level and set maximum depth to
deepest in `<vgrid>`. `<level>` is the minimum number of depth levels (e.g. 4).
Can produce non-advective cells.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
* `--vgrid_type <type>` can be mom5 or mom6 (default 'mom5')

### fill_fraction

Expand All @@ -75,33 +80,37 @@ to zero. Can produce non-advective cells and/or new seas.

```
usage: topogtools check_nonadvective --input <input_file>
[--vgrid <vgrid> --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.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
* `--potholes` check for potholes
* `--coastal-cells` check for non-advective coastal cells
* `--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

### fix_nonadvective

```
usage: topogtools fix_nonadvective --input <input_file> --output <output_file>
[--vgrid <vgrid> --potholes --coastal-cells]
[--vgrid <vgrid> --vgrid_type <type>
--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.

Options
* `--vgrid <vgrid>` vertical grid (default 'ocean_vgrid.nc')
* `--potholes` fix potholes
* `--coastal-cells` fix non-advective coastal cells
* `--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

### mask

Expand Down
30 changes: 11 additions & 19 deletions src/float_vgrid.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,16 +1,11 @@
program float_vgrid
! Zeta is double precision. Convert to single and rewrite. This stops small floating point errors.
use netcdf
use, intrinsic :: iso_fortran_env
use utils
use M_CLI2
use vgrid
use utils
implicit none

real(real32), allocatable :: zeta_float(:)
real(real64), allocatable :: zeta_dp(:)
integer :: ierr, nzeta
integer :: ncid, vid
integer :: dids(1)
type(vgrid_t) :: vgrid
character(len=:), allocatable :: help_text(:), file

help_text = [character(len=80) :: &
Expand All @@ -23,17 +18,14 @@ program float_vgrid
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
'']

call set_args('--vgrid "ocean_vgrid.nc"', help_text)
call set_args('--vgrid "ocean_vgrid.nc" --vgrid_type "mom5"', help_text)

file = sget('vgrid')
call check_file_exist(file)

call handle_error(nf90_open(sget('vgrid'), nf90_write,ncid))
call handle_error(nf90_inq_varid(ncid,'zeta', vid))
call handle_error(nf90_inquire_variable(ncid, vid, dimids=dids))
call handle_error(nf90_inquire_dimension(ncid, dids(1), len=nzeta))
allocate(zeta_dp(nzeta), zeta_float(nzeta))
call handle_error(nf90_get_var(ncid, vid, zeta_dp))
zeta_float = zeta_dp
zeta_dp = zeta_float
call handle_error(nf90_put_var(ncid, vid, zeta_dp))
call handle_error(nf90_close(ncid))
vgrid = vgrid_t(file, sget('vgrid_type'))
call vgrid%float()
call vgrid%update_history(get_mycommand())
call vgrid%write(file)

end program float_vgrid
47 changes: 17 additions & 30 deletions src/topography.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ module topography
use iso_fortran_env
use netcdf
use utils
use vgrid
implicit none

type topography_t
Expand Down Expand Up @@ -433,9 +434,9 @@ subroutine topography_deseas(this)
end subroutine topography_deseas

!-------------------------------------------------------------------------
subroutine topography_min_max_depth(this, vgrid, level)
subroutine topography_min_max_depth(this, vgrid_file, vgrid_type, level)
class(topography_t), intent(inout) :: this
character(len=*), intent(in) :: vgrid
character(len=*), intent(in) :: vgrid_file, vgrid_type
integer, intent(in) :: level

integer(int32) :: i,j
Expand All @@ -445,22 +446,13 @@ subroutine topography_min_max_depth(this, vgrid, level)
integer(int32) :: dids_lev(1) ! NetCDF ids
integer(int32) :: zlen ! length of zeta array

real(real64) :: zeta
real(real64), allocatable :: zeta_arr(:)
type(vgrid_t) :: vgrid

this%min_level = level

call handle_error(nf90_open(trim(vgrid), nf90_nowrite, ncid_lev))
call handle_error(nf90_inq_varid(ncid_lev, 'zeta', lev_id))
call handle_error(nf90_get_var(ncid_lev, lev_id, zeta, start=[2*this%min_level+1]))
this%min_depth = zeta

call handle_error(nf90_inquire_variable(ncid_lev, lev_id, dimids=dids_lev))
call handle_error(nf90_inquire_dimension(ncid_lev, dids_lev(1), len=zlen))
call handle_error(nf90_get_var(ncid_lev, lev_id, zeta, start=[zlen]))
this%max_depth = zeta

call handle_error(nf90_close(ncid_lev))
vgrid = vgrid_t(vgrid_file, vgrid_type)
this%min_depth = vgrid%zeta(this%min_level)
this%max_depth = vgrid%zeta(vgrid%nlevels)

write(output_unit,'(a,f7.2,a)') 'Setting minimum depth to ', this%min_depth, ' m'
write(output_unit,'(a,f7.2,a)') 'Setting maximum depth to ', this%max_depth, ' m'
Expand Down Expand Up @@ -506,15 +498,16 @@ subroutine topography_fill_fraction(this, sea_area_fraction)
end subroutine topography_fill_fraction

!-------------------------------------------------------------------------
subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coastal, fix)
class(topography_t), intent(inout) :: this
character(len=*), intent(in) :: vgrid
character(len=*), intent(in) :: vgrid_file, vgrid_type
logical, intent(in) :: potholes, coastal, fix

real(real32), allocatable :: depth_halo(:,:)
type(vgrid_t) :: vgrid
integer(int32), allocatable :: num_levels(:,:)
real(real32), allocatable :: zw(:), zeta(:)
integer(int32) :: passes, i, j, k, ni, nj, nzeta, nz, its, coastal_counter, potholes_counter
real(real32), allocatable :: zw(:)
integer(int32) :: passes, i, j, k, ni, nj, its, coastal_counter, potholes_counter
integer(int32) :: ncid, vid
integer(int32) :: dids(2)
logical :: se, sw, ne, nw ! .TRUE. if C-cell centre is shallower than T cell centre.
Expand All @@ -523,16 +516,10 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
integer(int32) :: im, ip, jm, jp
integer(int32) :: nseas

call handle_error(nf90_open(trim(vgrid), nf90_nowrite, ncid))
call handle_error(nf90_inq_varid(ncid, 'zeta', vid))
call handle_error(nf90_inquire_variable(ncid, vid, dimids=dids))
call handle_error(nf90_inquire_dimension(ncid, dids(1), len=nzeta))
nz = nzeta/2
write(output_unit,*) 'Zeta dimensions', nzeta, nz
allocate(zeta(nzeta), zw(0:nz))
call handle_error(nf90_get_var(ncid, vid, zeta))
call handle_error(nf90_close(ncid))
zw(:) = zeta(1:nzeta:2)
vgrid = vgrid_t(vgrid_file, vgrid_type)
write(output_unit,*) 'Zeta dimensions', 2*vgrid%nlevels + 1, vgrid%nlevels
allocate(zw(0:vgrid%nlevels))
zw = real(vgrid%zeta)

write(output_unit,*) 'depth dimensions', this%nxt, this%nyt
allocate(depth_halo(0:this%nxt+1, this%nyt+1))
Expand All @@ -559,7 +546,7 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
do j = 1, this%nyt + 1
do i = 0, this%nxt + 1
if (depth_halo(i, j) > 0.0) then
kloop: do k = 2, nz
kloop: do k = 2, vgrid%nlevels
if (zw(k) >= depth_halo(i, j)) then
num_levels(i, j) = k
exit kloop
Expand Down
31 changes: 19 additions & 12 deletions src/topogtools.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ program topogtools
character(len=:), allocatable :: help_general(:), help_gen_topo(:), help_deseas(:), help_min_max_depth(:)
character(len=:), allocatable :: help_fill_fraction(:), help_fix_nonadvective(:), help_check_nonadvective(:), help_mask(:)
character(len=80) :: version_text(1)
character(len=:), allocatable :: file_in, file_out, hgrid, vgrid, grid_type
character(len=:), allocatable :: file_in, file_out, hgrid, vgrid
type(topography_t) :: topog
real(real32) :: sea_area_fraction
integer :: ii
Expand Down Expand Up @@ -55,14 +55,16 @@ program topogtools
'']
help_min_max_depth = [character(len=80) :: &
'usage: topogtools min_max_depth --input <input_file> --output <output_file> ', &
' --level <level> [--vgrid <vgrid>] ', &
' --level <level> ', &
' [--vgrid <vgrid> --vgrid_type <type>] ', &
' ', &
'Set minimum depth to the depth at a specified level and set maximum depth to ', &
'deepest in <vgrid>. <level> is the minimum number of depth levels (e.g. 4). ', &
'Can produce non-advective cells. ', &
' ', &
'Options ', &
' --vgrid <vgrid> vertical grid (default ''ocean_vgrid.nc'') ', &
' --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 @@ -73,27 +75,31 @@ program topogtools
'']
help_fix_nonadvective = [character(len=80) :: &
'usage: topogtools fix_nonadvective --input <input_file> --output <output_file> ', &
' [--vgrid <vgrid> --potholes --coastal_cells] ', &
' [--vgrid <vgrid> --vgrid_type <type> ', &
' --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. ', &
' ', &
'Options ', &
' --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 ', &
'']
help_check_nonadvective = [character(len=80) :: &
'usage: topogtools check_nonadvective --input <input_file> ', &
' [--vgrid <vgrid> --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. ', &
' ', &
'Options ', &
' --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 ', &
'']
Expand All @@ -112,14 +118,15 @@ program topogtools
case ('deseas')
call set_args('--input:i "unset" --output:o "unset"', help_deseas, version_text)
case ('min_max_depth')
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --level 0', help_min_max_depth, version_text)
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --level 0', &
help_min_max_depth, version_text)
case('fill_fraction')
call set_args('--input:i "unset" --output:o "unset" --fraction 0.0', help_fill_fraction, version_text)
case ('fix_nonadvective')
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --potholes F --coastal-cells F', &
help_fix_nonadvective, version_text)
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --potholes F &
&--coastal-cells F', help_fix_nonadvective, version_text)
case ('check_nonadvective')
call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --potholes F --coastal-cells F', &
call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --potholes F --coastal-cells F', &
help_check_nonadvective, version_text)
case ('mask')
call set_args('--input:i "unset" --output:o "unset"', help_mask, version_text)
Expand Down Expand Up @@ -174,7 +181,7 @@ program topogtools

case ('min_max_depth')
topog = topography_t(file_in)
call topog%min_max_depth(vgrid, iget('level'))
call topog%min_max_depth(vgrid, sget('vgrid_type'), iget('level'))
call topog%update_history(get_mycommand())
call topog%write(file_out)

Expand All @@ -191,13 +198,13 @@ program topogtools

case ('fix_nonadvective')
topog = topography_t(file_in)
call topog%nonadvective(vgrid, lget('potholes'), lget('coastal-cells'), fix=.true.)
call topog%nonadvective(vgrid, sget('vgrid_type'), lget('potholes'), lget('coastal-cells'), fix=.true.)
call topog%update_history(get_mycommand())
call topog%write(file_out)

case ('check_nonadvective')
topog = topography_t(file_in)
call topog%nonadvective(vgrid, lget('potholes'), lget('coastal-cells'), fix=.false.)
call topog%nonadvective(vgrid, sget('vgrid_type'), lget('potholes'), lget('coastal-cells'), fix=.false.)

case ('mask')
topog = topography_t(file_in)
Expand Down
Loading
Loading