Allow for MOM6 vertical grids through a new command line option. Thes…

…e grids have only n+1 cell boundaries instead of the 2n+1 interleaved cell boundaries and centres of the MOM5 grids.

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 `<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)
 
@@ -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
 
@@ -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
 

src/float_vgrid.f90

@@ -18,12 +18,12 @@ 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)
 
-  vgrid = vgrid_t(file)
+  vgrid = vgrid_t(file, sget('vgrid_type'))
   call vgrid%float()
   call vgrid%update_history(get_mycommand())
   call vgrid%write(file)

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)

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 <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>     ', &
@@ -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                      ', &
     '']
@@ -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)

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,29 +32,34 @@ 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))
 
     ! 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,26 +120,38 @@ 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))
     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, 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)