add_t_cell_cutoff_function

README.md

@@ -137,6 +137,13 @@ double-precision topography file.
 Options
   * `--vgrid <vgrid>`  vertical grid (default 'ocean_vgrid.nc')
 
+## cut_off_T_cells
+
+```
+usage: topogtools cut_off_T_cells --input <input_file> --output <output_file> --hgrid <grid> --cutoff <cutoff_value>       
+```
+
+Cut off T cells with size smaller than <cutoff_value>. Cut off should be in kilometers
 
 # 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 :: cut_off_T_cells => topography_cut_off_T_cells
   end type topography_t
 
   interface topography_t
@@ -464,6 +465,55 @@ subroutine topography_min_max_depth(this, vgrid_file, vgrid_type, level)
 
   end subroutine topography_min_max_depth
 
+ subroutine topography_cut_off_T_cells(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
+    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
+    real(real64), allocatable :: dy_t(:,:)       ! To store dy_t (new array)
+
+    ! Read hgrid to get dy
+    print*, 'Attempting to open:', 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=ny_len))
+    call handle_error(nf90_inquire_dimension(ncid_hgrid, dids_dy(2), len=nxp_len))
+
+    ! Allocate memory for dy based on its dimensions
+    allocate(dy(ny_len, nxp_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 dy_t based on dy
+    ! dy_t = dy[::2, 1::2] + dy[1::2, 1::2]
+    ! This means dy_t will have half the number of rows and columns as dy
+    allocate(dy_t(int(ny_len / 2), int((nxp_len - 1) / 2)))
+
+    do i = 1, int(ny_len / 2)
+      do j = 1, int((nxp_len - 1) / 2)
+          dy_t(i, j) = dy(2 * i - 1, 2 * j) + dy(2 * i, 2 * j)
+      end do
+    end do
+
+    ! Apply cutoff to depth based on the provided T-cell cutoff value in kilometers
+    do i = 1, int(ny_len / 2)
+      do j = 1, int((nxp_len - 1) / 2)
+            if (dy_t(i, j) < cutoff*1000.0) then  !Input cutoff in Kilometers covert it to 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_cut_off_T_cells
+
   !-------------------------------------------------------------------------
   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                                                          ', &
+    '  cut_off_T_cells - Cut off T cells below a certain cell size                   ', &
     '']
   help_gen_topo = [character(len=80) :: &
     'usage: topogtools gen_topo --input <input_file> --output <output_file>          ', &
@@ -109,6 +111,14 @@ program topogtools
     'Creates a land mask from a topography.                                          ', &
     '']
 
+  help_cutoff = [character(len=80) :: &
+    'usage: topogtools cut_off_T_cells --input <input_file> --output <output_file>   ', &
+    '                                  --hgrid <grid> --cutoff <cutoff_value>        ', &
+    '                                                                                ', &
+    'Cut off T cells with size smaller than <cutoff_value>. Writes the               ', &
+    'result to <output_file>.                                                        ', &
+    'Cut off should be in kilometers']
+
   ! Read command line
   name = get_subcommand()
   select case (name)
@@ -130,6 +140,8 @@ program topogtools
       help_check_nonadvective, version_text)
   case ('mask')
     call set_args('--input:i "unset" --output:o "unset"', help_mask, version_text)
+  case ('cut_off_T_cells')
+    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 +222,24 @@ program topogtools
     topog = topography_t(file_in)
     call topog%mask(file_out)
 
+  case ('cut_off_T_cells')
+    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%cut_off_T_cells(hgrid,cutoff)
+    call topog%update_history(get_mycommand())
+    call topog%write(file_out)
+
   end select
 
 end program topogtools