diff --git a/CMakeLists.txt b/CMakeLists.txt index 9e811d7..b7dd3a8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 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 32e4bbc..793d151 100644 --- a/src/float_vgrid.f90 +++ b/src/float_vgrid.f90 @@ -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) :: & @@ -23,17 +18,14 @@ 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) - 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 diff --git a/src/topography.f90 b/src/topography.f90 index 7663bc6..c562b57 100644 --- a/src/topography.f90 +++ b/src/topography.f90 @@ -2,6 +2,7 @@ module topography use iso_fortran_env use netcdf use utils + use vgrid implicit none type topography_t @@ -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 @@ -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' @@ -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. @@ -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)) @@ -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 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 new file mode 100644 index 0000000..9ade85e --- /dev/null +++ b/src/vgrid.f90 @@ -0,0 +1,210 @@ +module vgrid + use iso_fortran_env + use netcdf + use utils + implicit none + + private + public :: vgrid_t + + type vgrid_t + ! Vertical grid variable + character(len=:), allocatable :: type + integer :: nlevels = 0 + real(real64), allocatable :: zeta(:) + real(real64), allocatable :: zeta_super(:) + character(len=:), allocatable :: author + ! Global attributes + character(len=:), allocatable :: original_file + character(len=:), allocatable :: history + contains + procedure :: copy => vgrid_copy + generic :: assignment(=) => copy + procedure :: write => vgrid_write + procedure :: update_history => vgrid_update_history + procedure :: float => vgrid_float + end type vgrid_t + + interface vgrid_t + module procedure vgrid_constructor + end interface vgrid_t + +contains + + !------------------------------------------------------------------------- + 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)) + + ! Get dimension + call handle_error(nf90_inq_dimid(ncid, 'nzv', did(1))) + call handle_error(nf90_inquire_dimension(ncid, did(1), len=zeta_len)) + + ! Get zeta + allocate(zeta_var(zeta_len)) + call handle_error(nf90_inq_varid(ncid, 'zeta', zeta_id)) + 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)) + end if + + ! History (might not be present) + if (nf90_inquire_attribute(ncid, nf90_global, 'history', len=history_len) == nf90_noerr) then + allocate(character(len=history_len) :: vgrid%history) + 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") + if (mod(zeta_len, 2) == 0) then + write(error_unit,'(a)') "ERROR: MOM5 vertical grid has an even number of points, which should never happen." + error stop + end if + 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)) + + end function vgrid_constructor + + !------------------------------------------------------------------------- + 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) + 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 + + ! Global attributes + vgrid_out%original_file = vgrid_in%original_file + if (allocated(vgrid_in%history)) then + vgrid_out%history = vgrid_in%history + end if + + end subroutine vgrid_copy + + !------------------------------------------------------------------------- + subroutine vgrid_write(this, filename) + class(vgrid_t), target, intent(in) :: this + character(len=*), intent(in) :: filename + + 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', zeta_len, did(1))) + + ! Write zeta + call handle_error(nf90_def_var(ncid, 'zeta', nf90_double, did, zeta_id)) + call handle_error(nf90_put_att(ncid, zeta_id, 'units', 'meters')) + 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, zeta_var)) + + ! Write global attributes + call handle_error(nf90_put_att(ncid, nf90_global, 'original_file', trim(this%original_file))) + call handle_error(nf90_put_att(ncid, nf90_global, 'history', trim(this%history))) + + ! Close file + call handle_error(nf90_enddef(ncid)) + call handle_error(nf90_close(ncid)) + + end subroutine vgrid_write + + !------------------------------------------------------------------------- + subroutine vgrid_update_history(this, command) + use iso_c_binding + class(vgrid_t), intent(inout) :: this + character(len=*), intent(in) :: command + + character(len=:), allocatable :: new_history, old_history + + new_history = date_time() // ": " // trim(command) + if (allocated(this%history)) then + old_history = this%history + deallocate(this%history) + this%history = trim(new_history) // C_NEW_LINE // old_history + else + this%history = trim(new_history) + end if + + end subroutine vgrid_update_history + + !------------------------------------------------------------------------- + subroutine vgrid_float(this) + class(vgrid_t), intent(inout) :: this + + real(real32), allocatable :: 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 = real(this%zeta, real32) + this%zeta = zeta_float + deallocate(zeta_float) + + end subroutine vgrid_float + + +end module vgrid