add_t_cell_cutoff_function

README.md

@@ -116,15 +116,10 @@ Options
 
 ```
 usage: topogtools mask  --input <input_file> --output <output_file>
-                        [--fraction <frac>]
 ```
 
 Creates a land mask from a topography.
 
-Options
-  * `--fraction <frac>`  cells with a fraction of sea area smaller than `<frac>` will be set as land (default '0.0')
-
-
 ## float_vgrid
 
 ```
@@ -137,6 +132,16 @@ double-precision topography file.
 Options
   * `--vgrid <vgrid>`  vertical grid (default 'ocean_vgrid.nc')
 
+## min_dy
+
+```
+usage: topogtools min_dy --input <input_file> --output <output_file> --hgrid <grid> --cutoff <cutoff_value>       
+```
+
+Convert ocean cells into land if their y size is smaller than <cutoff_value> metres.
+
+Options
+  * `--hgrid <hgrid>`  vertical grid (default 'ocean_hgrid.nc')
 
 # Building and Installation
 

src/topography.f90

@@ -35,6 +35,7 @@ module topography
     procedure :: nonadvective => topography_nonadvective
     procedure :: min_max_depth => topography_min_max_depth
     procedure :: mask => topography_mask
+    procedure :: min_dy => topography_min_dy
   end type topography_t
 
   interface topography_t
@@ -464,6 +465,45 @@ subroutine topography_min_max_depth(this, vgrid_file, vgrid_type, level)
 
   end subroutine topography_min_max_depth
 
+ subroutine topography_min_dy(this, hgrid, cutoff)
+    class(topography_t), intent(inout) :: this
+    character(len=*), intent(in) :: hgrid
+    real(real64), intent(in) :: cutoff
+
+    integer(int32) :: i, j
+    integer(int32) :: ncid_hgrid, dy_id          ! NetCDF ids for hgrid and dy
+    integer(int32) :: dids_dy(2)                 ! NetCDF ids for dimensions
+    integer(int32) :: ny_len, nxp_len, nx_len    ! dimensions for hgrid
+    real(real64), allocatable :: dy(:,:)         ! To store dy variable from hgrid
+
+    ! Read hgrid to get dy
+    write(output_unit,'(3a)') "Attempting to open file '", trim(hgrid), "'"
+    call handle_error(nf90_open(trim(hgrid), nf90_nowrite, ncid_hgrid))
+    call handle_error(nf90_inq_varid(ncid_hgrid, 'dy', dy_id))
+    call handle_error(nf90_inquire_variable(ncid_hgrid, dy_id, dimids=dids_dy))
+    call handle_error(nf90_inquire_dimension(ncid_hgrid, dids_dy(1), len=nxp_len))
+    call handle_error(nf90_inquire_dimension(ncid_hgrid, dids_dy(2), len=ny_len))
+
+    ! Allocate memory for dy based on its dimensions
+    allocate(dy(nxp_len, ny_len))
+
+    ! Read the dy variable from hgrid
+    call handle_error(nf90_get_var(ncid_hgrid, dy_id, dy))
+    call handle_error(nf90_close(ncid_hgrid))
+
+    ! Calculate T cell size based on dy
+    ! For each point, the T cell size is a sum of dy(2*i-1, 2*j) and dy(2*i, 2*j)    
+    ! Apply cutoff to depth based on the provided T-cell cutoff value in meters
+    do j = 1, ny_len / 2
+      do i = 1, (nxp_len - 1) / 2
+          if (dy(2 * i - 1, 2 * j) + dy(2 * i, 2 * j) < cutoff) then  !Input cutoff in meters
+            this%depth(i, j) = MISSING_VALUE  ! Set values below cutoff to zero or another value as needed
+          end if
+        end do
+    end do   
+
+end subroutine topography_min_dy
+
   !-------------------------------------------------------------------------
   subroutine topography_fill_fraction(this, sea_area_fraction)
     class(topography_t), intent(inout) :: this

src/topogtools.f90

@@ -8,12 +8,13 @@ program topogtools
   character(len=5), PARAMETER :: VERSION = "1.0.0"
 
   character(len=:), allocatable :: name
-  character(len=:), allocatable :: help_general(:), help_gen_topo(:), help_deseas(:), help_min_max_depth(:)
+  character(len=:), allocatable :: help_general(:), help_gen_topo(:), help_deseas(:), help_min_max_depth(:), help_cutoff(:)
   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
   type(topography_t) :: topog
   real(real32) :: sea_area_fraction
+  real(real64) :: cutoff
   integer :: ii
 
   version_text(1) = 'topogtools version '//VERSION
@@ -33,6 +34,7 @@ program topogtools
     '  check_nonadvective - Check for non-advective cells                            ', &
     '  fix_nonadvective - Fix non-advective cells                                    ', &
     '  mask - Generate mask                                                          ', &
+    '  min_dy - Set small ocean cells to land                                        ', &
     '']
   help_gen_topo = [character(len=80) :: &
     'usage: topogtools gen_topo --input <input_file> --output <output_file>          ', &
@@ -109,6 +111,18 @@ program topogtools
     'Creates a land mask from a topography.                                          ', &
     '']
 
+  help_cutoff = [character(len=80) :: &
+    'usage: topogtools min_dy --input <input_file> --output <output_file>   ', &
+    '                                  --cutoff <cutoff_value>        ', &
+    '                                  [--hgrid <hgrid_file>]        ', &
+    '                                                                                ', &
+    'Convert ocean cells into land if their y size is less than <cutoff_value> in    ', &
+    'kilometres. Writes the result to <output_file>.                                 ', &
+    '                                                                                ', &
+    'Options                                                                         ', &
+    '   --hgrid <hgrid_file>     horizontal supergrid (default ''ocean_hgrid.nc'')   ', &
+    '']
+
   ! Read command line
   name = get_subcommand()
   select case (name)
@@ -130,6 +144,8 @@ program topogtools
       help_check_nonadvective, version_text)
   case ('mask')
     call set_args('--input:i "unset" --output:o "unset"', help_mask, version_text)
+  case ('min_dy')
+    call set_args('--input:i "unset" --output:o "unset" --hgrid "ocean_hgrid.nc" --cutoff 0.0', help_cutoff, version_text)
   case ('')
     ! Print help in case the user specified the --help flag
     call set_args(' ', help_general, version_text)
@@ -210,6 +226,24 @@ program topogtools
     topog = topography_t(file_in)
     call topog%mask(file_out)
 
+  case ('min_dy')
+    hgrid = sget('hgrid')
+    call check_file_exist(hgrid)
+    cutoff = rget('cutoff')
+    if (cutoff <= 0.0) then
+      write(error_unit,'(a)') "ERROR: cutoff value must be larger than 0 "
+      error stop
+    end if
+    file_out = sget('output')
+    if (file_out == 'unset') then
+      write(error_unit,'(a)') 'ERROR: no output file specified'
+      error stop
+    end if  
+    topog = topography_t(file_in)
+    call topog%min_dy(hgrid, cutoff)
+    call topog%update_history(get_mycommand())
+    call topog%write(file_out)
+
   end select
 
 end program topogtools