From 8c7fd1d1b6580e47d5911f10bb3e3734fc5a2de1 Mon Sep 17 00:00:00 2001 From: Micael Oliveira Date: Wed, 18 Sep 2024 09:13:01 +1000 Subject: [PATCH] Allow for MOM6 vertical grids through a new command line option. These grids have only n+1 cell boundaries instead of the 2n+1 interleaved cell boundaries and centres of the MOM5 grids. --- README.md | 31 ++++++++++++------- src/float_vgrid.f90 | 4 +-- src/topography.f90 | 12 ++++---- src/topogtools.f90 | 31 +++++++++++-------- src/vgrid.f90 | 73 ++++++++++++++++++++++++++++++++++----------- 5 files changed, 103 insertions(+), 48 deletions(-) diff --git a/README.md b/README.md index 97bd7ab..16eb61c 100644 --- a/README.md +++ b/README.md @@ -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 `` is assumed to be a MOM5 vertical grid file (with $2n+1$ values for an $n$-level model). Using a MOM6 `` 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) @@ -51,7 +54,8 @@ Remove enclosed seas from and writes the result to . ``` usage: topogtools min_max_depth --input --output - --level [--vgrid ] + --level + [--vgrid --vgrid_type ] ``` Set minimum depth to the depth at a specified level and set maximum depth to @@ -59,7 +63,8 @@ deepest in ``. `` is the minimum number of depth levels (e.g. 4). Can produce non-advective cells. Options - * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') + * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') + * `--vgrid_type ` can be mom5 or mom6 (default 'mom5') ### fill_fraction @@ -75,7 +80,8 @@ to zero. Can produce non-advective cells and/or new seas. ``` usage: topogtools check_nonadvective --input - [--vgrid --potholes --coastal-cells] + [--vgrid --vgrid_type + --potholes --coastal-cells] ``` Check for non-advective cells. There are two types of checks available: potholes @@ -83,15 +89,17 @@ and non-advective coastal cells. Checking for non-advective coastal cells should only be needed when using a B-grid. Options - * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') - * `--potholes` check for potholes - * `--coastal-cells` check for non-advective coastal cells + * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') + * `--vgrid_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 --output - [--vgrid --potholes --coastal-cells] + [--vgrid --vgrid_type + --potholes --coastal-cells] ``` Fix non-advective cells. There are two types of fixes available: potholes and @@ -99,9 +107,10 @@ non-advective coastal cells. Fixes to non-advective coastal cells should only be needed when using a B-grid. Options - * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') - * `--potholes` fix potholes - * `--coastal-cells` fix non-advective coastal cells + * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') + * `--vgrid_type ` can be mom5 or mom6 (default 'mom5') + * `--potholes` fix potholes + * `--coastal-cells` fix non-advective coastal cells ### mask diff --git a/src/float_vgrid.f90 b/src/float_vgrid.f90 index f4dd3b9..793d151 100644 --- a/src/float_vgrid.f90 +++ b/src/float_vgrid.f90 @@ -18,12 +18,12 @@ program float_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) - vgrid = vgrid_t(file) + vgrid = vgrid_t(file, sget('vgrid_type')) call vgrid%float() call vgrid%update_history(get_mycommand()) call vgrid%write(file) diff --git a/src/topography.f90 b/src/topography.f90 index 4545677..c562b57 100644 --- a/src/topography.f90 +++ b/src/topography.f90 @@ -434,9 +434,9 @@ subroutine topography_deseas(this) end subroutine topography_deseas !------------------------------------------------------------------------- - subroutine topography_min_max_depth(this, vgrid_file, level) + subroutine topography_min_max_depth(this, vgrid_file, vgrid_type, level) class(topography_t), intent(inout) :: this - character(len=*), intent(in) :: vgrid_file + character(len=*), intent(in) :: vgrid_file, vgrid_type integer, intent(in) :: level integer(int32) :: i,j @@ -450,7 +450,7 @@ subroutine topography_min_max_depth(this, vgrid_file, level) this%min_level = level - vgrid = vgrid_t(vgrid_file) + vgrid = vgrid_t(vgrid_file, vgrid_type) this%min_depth = vgrid%zeta(this%min_level) this%max_depth = vgrid%zeta(vgrid%nlevels) @@ -498,9 +498,9 @@ subroutine topography_fill_fraction(this, sea_area_fraction) end subroutine topography_fill_fraction !------------------------------------------------------------------------- - subroutine topography_nonadvective(this, vgrid_file, 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_file + character(len=*), intent(in) :: vgrid_file, vgrid_type logical, intent(in) :: potholes, coastal, fix real(real32), allocatable :: depth_halo(:,:) @@ -516,7 +516,7 @@ subroutine topography_nonadvective(this, vgrid_file, potholes, coastal, fix) integer(int32) :: im, ip, jm, jp integer(int32) :: nseas - vgrid = vgrid_t(vgrid_file) + 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) diff --git a/src/topogtools.f90 b/src/topogtools.f90 index de8c075..dbc4313 100644 --- a/src/topogtools.f90 +++ b/src/topogtools.f90 @@ -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 @@ -55,14 +55,16 @@ program topogtools ''] help_min_max_depth = [character(len=80) :: & 'usage: topogtools min_max_depth --input --output ', & - ' --level [--vgrid ] ', & + ' --level ', & + ' [--vgrid --vgrid_type ] ', & ' ', & 'Set minimum depth to the depth at a specified level and set maximum depth to ', & 'deepest in . is the minimum number of depth levels (e.g. 4). ', & 'Can produce non-advective cells. ', & ' ', & 'Options ', & - ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & ''] help_fill_fraction = [character(len=80) :: & 'usage: topogtools fill_fraction --input --output ', & @@ -73,7 +75,8 @@ program topogtools ''] help_fix_nonadvective = [character(len=80) :: & 'usage: topogtools fix_nonadvective --input --output ', & - ' [--vgrid --potholes --coastal_cells] ', & + ' [--vgrid --vgrid_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', & @@ -81,12 +84,14 @@ program topogtools ' ', & 'Options ', & ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_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 ', & - ' [--vgrid --potholes --coastal_cells] ', & + ' [--vgrid --vgrid_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', & @@ -94,6 +99,7 @@ program topogtools ' ', & 'Options ', & ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & + ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & ' --potholes check for potholes ', & ' --coastal-cells check for non-advective coastal cells ', & ''] @@ -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) @@ -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) @@ -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) diff --git a/src/vgrid.f90 b/src/vgrid.f90 index a33d75a..8627706 100644 --- a/src/vgrid.f90 +++ b/src/vgrid.f90 @@ -9,6 +9,7 @@ module vgrid type vgrid_t ! Vertical grid variable + character(len=:), allocatable :: type integer :: nlevels = 0 real(real64), allocatable :: zeta(:) real(real64), allocatable :: zeta_super(:) @@ -31,14 +32,22 @@ module vgrid contains !------------------------------------------------------------------------- - type(vgrid_t) function vgrid_constructor(filename) result(vgrid) - character(len=*), intent(in) :: filename + type(vgrid_t) function vgrid_constructor(filename, type) result(vgrid) + character(len=*), intent(in) :: filename, type integer(int32) :: ncid, zeta_id, did(1), zeta_len, author_len, history_len ! NetCDF ids and dims + real(real64), allocatable :: zeta_var(:) write(output_unit,'(3a)') "Reading vgrid from file '", trim(filename), "'" vgrid%original_file = filename + vgrid%type = type + select case (type) + case ("mom5", "mom6") + case default + write(error_unit,'(3a)') "ERROR: '", type, "' is not a valid vertical grid type." + error stop + end select ! Open file call handle_error(nf90_open(trim(filename), nf90_nowrite, ncid)) @@ -46,14 +55,11 @@ type(vgrid_t) function vgrid_constructor(filename) result(vgrid) ! Get dimension call handle_error(nf90_inq_dimid(ncid, 'nzv', did(1))) call handle_error(nf90_inquire_dimension(ncid, did(1), len=zeta_len)) - vgrid%nlevels = zeta_len/2 ! Get zeta - allocate(vgrid%zeta_super(zeta_len)) - allocate(vgrid%zeta(0:vgrid%nlevels)) + allocate(zeta_var(zeta_len)) call handle_error(nf90_inq_varid(ncid, 'zeta', zeta_id)) - call handle_error(nf90_get_var(ncid, zeta_id, vgrid%zeta_super)) - vgrid%zeta = vgrid%zeta_super(1:zeta_len:2) + call handle_error(nf90_get_var(ncid, zeta_id, zeta_var)) if (nf90_inquire_attribute(ncid, zeta_id, 'author', len=author_len) == nf90_noerr) then allocate(character(len=author_len) :: vgrid%author) call handle_error(nf90_get_att(ncid, zeta_id, 'author', vgrid%author)) @@ -65,6 +71,21 @@ type(vgrid_t) function vgrid_constructor(filename) result(vgrid) call handle_error(nf90_get_att(ncid, nf90_global, 'history', vgrid%history)) end if + ! Handle the different types of grids + select case (type) + case ("mom5") + vgrid%nlevels = zeta_len/2 + allocate(vgrid%zeta_super(zeta_len)) + allocate(vgrid%zeta(0:vgrid%nlevels)) + vgrid%zeta_super = zeta_var + vgrid%zeta = zeta_var(1:zeta_len:2) + + case ("mom6") + vgrid%nlevels = zeta_len - 1 + allocate(vgrid%zeta(0:vgrid%nlevels)) + vgrid%zeta = zeta_var + end select + ! Close file call handle_error(nf90_close(ncid)) @@ -75,12 +96,16 @@ subroutine vgrid_copy(vgrid_out, vgrid_in) class(vgrid_t), intent(out) :: vgrid_out class(vgrid_t), intent(in) :: vgrid_in + vgrid_out%type = vgrid_in%type + ! Dimension vgrid_out%nlevels = vgrid_in%nlevels ! Zeta variable and attributes allocate(vgrid_out%zeta, source=vgrid_in%zeta) - allocate(vgrid_out%zeta_super, source=vgrid_in%zeta_super) + if (allocated(vgrid_in%zeta_super)) then + allocate(vgrid_out%zeta_super, source=vgrid_in%zeta_super) + end if if (allocated(vgrid_in%author)) then vgrid_out%author = vgrid_in%author end if @@ -95,18 +120,30 @@ end subroutine vgrid_copy !------------------------------------------------------------------------- subroutine vgrid_write(this, filename) - class(vgrid_t), intent(in) :: this + class(vgrid_t), target, intent(in) :: this character(len=*), intent(in) :: filename - integer(int32) :: ncid, zeta_id, did(1) ! NetCDF ids + integer(int32) :: ncid, zeta_id, zeta_len, did(1) ! NetCDF ids + real(real64), pointer :: zeta_var(:) write(output_unit,'(3a)') "Writing vgrid to file '", trim(filename), "'" + ! Handle the different types of grids + select case (this%type) + case ("mom5") + zeta_len = 2*this%nlevels + 1 + zeta_var(1:zeta_len) => this%zeta_super(1:zeta_len) + + case ("mom6") + zeta_len = this%nlevels + 1 + zeta_var(1:zeta_len) => this%zeta(0:zeta_len - 1) + end select + ! Open file call handle_error(nf90_create(trim(filename), ior(nf90_netcdf4, nf90_clobber), ncid)) ! Write dimension - call handle_error(nf90_def_dim(ncid, 'nzv', 2*this%nlevels+1, did(1))) + call handle_error(nf90_def_dim(ncid, 'nzv', zeta_len, did(1))) ! Write zeta call handle_error(nf90_def_var(ncid, 'zeta', nf90_double, did, zeta_id)) @@ -114,7 +151,7 @@ subroutine vgrid_write(this, filename) call handle_error(nf90_put_att(ncid, zeta_id, 'standard_name', 'vertical_grid_vertex')) call handle_error(nf90_put_att(ncid, zeta_id, 'long_name', 'vgrid')) call handle_error(nf90_put_att(ncid, zeta_id, 'author', trim(this%author))) - call handle_error(nf90_put_var(ncid, zeta_id, this%zeta_super)) + call handle_error(nf90_put_var(ncid, zeta_id, zeta_var)) ! Write global attributes call handle_error(nf90_put_att(ncid, nf90_global, 'original_file', trim(this%original_file))) @@ -151,13 +188,15 @@ subroutine vgrid_float(this) real(real32), allocatable :: zeta_float(:) - allocate(zeta_float(2*this%nlevels+1)) - zeta_float = this%zeta_super - this%zeta_super = zeta_float - deallocate(zeta_float) + if (allocated(this%zeta_super)) then + allocate(zeta_float(2*this%nlevels+1)) + zeta_float = real(this%zeta_super, real32) + this%zeta_super = zeta_float + deallocate(zeta_float) + end if allocate(zeta_float(0:this%nlevels)) - zeta_float = this%zeta + zeta_float = real(this%zeta, real32) this%zeta = zeta_float deallocate(zeta_float)