From cc75731bbb2450609266d3a4530cd3c02a78f805 Mon Sep 17 00:00:00 2001 From: ezhilsabareesh8 Date: Thu, 24 Oct 2024 12:40:11 +1100 Subject: [PATCH] add_t_cell_cutoff_function --- README.md | 7 +++++++ src/topography.f90 | 50 ++++++++++++++++++++++++++++++++++++++++++++++ src/topogtools.f90 | 32 ++++++++++++++++++++++++++++- 3 files changed, 88 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 16eb61c..71da438 100644 --- a/README.md +++ b/README.md @@ -137,6 +137,13 @@ double-precision topography file. Options * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') +## cut_off_T_cells + +``` +usage: topogtools cut_off_T_cells --input --output --hgrid --cutoff +``` + +Cut off T cells with size smaller than . Cut off should be in kilometers # Building and Installation diff --git a/src/topography.f90 b/src/topography.f90 index 8651e21..b411a1c 100644 --- a/src/topography.f90 +++ b/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 diff --git a/src/topogtools.f90 b/src/topogtools.f90 index dbc4313..c5ecf7c 100644 --- a/src/topogtools.f90 +++ b/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 --output ', & @@ -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 --output ', & + ' --hgrid --cutoff ', & + ' ', & + 'Cut off T cells with size smaller than . Writes the ', & + 'result to . ', & + '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