diff --git a/README.md b/README.md index 7310fba..bd9fa4a 100644 --- a/README.md +++ b/README.md @@ -73,30 +73,25 @@ Options ``` usage: topogtools fill_fraction --input --output - --fraction [--grid_type ] + --fraction ``` Cells with a fraction of sea area smaller than will have their depth set to zero. Can produce non-advective cells and/or new seas. -Options - * `--grid_type ` Arakawa type of horizontal grid ('B' or 'C'; default is 'B') - ### check_nonadvective ``` usage: topogtools check_nonadvective --input - [--grid_type - --vgrid --vgrid_type + [--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 -only be needed when using a B-grid. +Check topography for non-advective cells. There are two types of checks +available: potholes and non-advective coastal cells. B-grid connectivity rules +are assumed. Aborts if input_file is not on a B-grid. Options - * `--grid_type ` Arakawa type of horizontal grid ('B' or 'C'; default is 'B') * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') * `--vgrid_type ` can be mom5 or mom6 (default 'mom5') * `--potholes` check for potholes @@ -106,17 +101,15 @@ Options ``` usage: topogtools fix_nonadvective --input --output - [--grid_type - --vgrid --vgrid_type + [--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 -needed when using a B-grid. +Fix non-advective cells. There are two types of checks available: potholes and +non-advective coastal cells. B-grid connectivity rules are assumed. Aborts if +input_file is not on a B-grid. Options - * `--grid_type ` Arakawa type of horizontal grid ('B' or 'C'; default is 'B') * `--vgrid ` vertical grid (default 'ocean_vgrid.nc') * `--vgrid_type ` can be mom5 or mom6 (default 'mom5') * `--potholes` fix potholes diff --git a/src/topography.f90 b/src/topography.f90 index 89ecc57..a503949 100644 --- a/src/topography.f90 +++ b/src/topography.f90 @@ -49,18 +49,9 @@ type(topography_t) function topography_constructor(filename, grid_type) result(t character(len=*), intent(in) :: filename character(len=1), intent(in), optional :: grid_type + character(len=1):: file_grid_type integer(int32) :: ncid, depth_id, frac_id, geolon_id, geolat_id, dids(2), history_len ! NetCDF ids - if (present(grid_type)) then - if ( grid_type == 'B' .or. grid_type == 'C' ) then - topog%grid_type = grid_type - else - call handle_error(nf90_einval, .true., "grid_type must be B or C") - end if - else - topog%grid_type = 'B' - end if - write(output_unit,'(3a)') "Reading topography from file '", trim(filename), "'" topog%original_file = filename @@ -84,6 +75,19 @@ type(topography_t) function topography_constructor(filename, grid_type) result(t call handle_error(nf90_get_att(ncid, depth_id, 'maximum_depth', topog%max_depth), isfatal=.false.) call handle_error(nf90_get_att(ncid, depth_id, 'nonadvective_cells_removed', topog%nonadvective_cells_removed), isfatal=.false.) + if (present(grid_type)) then + topog%grid_type = grid_type ! grid_type arg overrides value in file + else + call handle_error(nf90_get_att(ncid, depth_id, 'grid_type', file_grid_type), isfatal=.false.) + if ( file_grid_type == 'B' .or. file_grid_type == 'C' ) then + topog%grid_type = file_grid_type + end if + end if + write(output_unit,*) " grid_type = ", topog%grid_type + if (.not.( topog%grid_type == 'B' .or. topog%grid_type == 'C' )) then + call handle_error(nf90_einval, .true., "grid_type must be B or C") + end if + ! Get sea area fraction call handle_error(nf90_inq_varid(ncid, 'sea_area_fraction', frac_id)) allocate(topog%frac(topog%nxt, topog%nyt)) @@ -541,6 +545,10 @@ subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coast integer(int32) :: im, ip, jm, jp integer(int32) :: nseas + if ( this%grid_type /= 'B' ) then + call handle_error(nf90_einval, .true., "nonadvective: grid_type must be B") + end if + vgrid = vgrid_t(vgrid_file, vgrid_type) write(output_unit,*) 'Zeta dimensions', 2*vgrid%nlevels + 1, vgrid%nlevels allocate(zw(0:vgrid%nlevels)) @@ -651,7 +659,7 @@ subroutine topography_nonadvective(this, vgrid_file, vgrid_type, potholes, coast if (fix .and. (coastal .or. potholes)) then this%nonadvective_cells_removed = 'yes' if (changes_made .and. this%lakes_removed == 'yes') then - ! Check if new lakes were created new lakes + ! Check if new lakes were created call this%number_seas(number_of_seas = nseas, silent=.true.) if (nseas > 1) then write(output_unit,'(a)') "WARNING: new seas have been created. To fix, rerun deseas again." diff --git a/src/topogtools.f90 b/src/topogtools.f90 index 918cb21..4a6285f 100644 --- a/src/topogtools.f90 +++ b/src/topogtools.f90 @@ -23,7 +23,7 @@ program topogtools 'usage: topogtools [--help] [] ', & ' ', & 'Collection of tools to edit and manipulate ocean model topographies. ', & - 'See ''topogtools --help '' to read about a specific subcommand. ', & + 'See ''topogtools --help '' to read about a specific subcommand. ', & ' ', & 'Available commands: ', & ' gen_topo - Generate a new topography file from a bathymetry ', & @@ -46,7 +46,7 @@ program topogtools ' tripolar grid ', & ' --longitude-offset offset (in degrees) between the central ', & ' longitude of the ocean horizontal grid and of ', & - ' the bathymetry grid (default ''0.0'') ', & + ' the bathymetry grid (default ''0.0'') ', & ''] help_deseas = [character(len=80) :: & 'usage: topogtools deseas --input --output ', & @@ -67,23 +67,19 @@ program topogtools 'Can produce non-advective cells. ', & ' ', & 'Options ', & - ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & - ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & + ' --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 ', & - ' --fraction [--grid_type ] ', & + ' --fraction ', & ' ', & 'Cells with a fraction of sea area smaller than will have their depth set ', & 'to zero. Can produce non-advective cells and/or new seas. ', & - ' ', & - 'Options ', & - ' --grid_type Arakawa type of horizontal grid (''B'' or ''C''; default ''B'')', & ''] help_fix_nonadvective = [character(len=80) :: & 'usage: topogtools fix_nonadvective --input --output ', & - ' [--grid_type ', & - ' --vgrid --vgrid_type ', & + ' [--vgrid --vgrid_type ', & ' --potholes --coastal_cells] ', & ' ', & 'Fix non-advective cells. There are two types of fixes available: potholes and ', & @@ -91,16 +87,14 @@ program topogtools 'needed when using a B-grid. ', & ' ', & 'Options ', & - ' --grid_type Arakawa type of horizontal grid (''B'' or ''C''; default ''B'')', & - ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & - ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & + ' --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 ', & - ' [--grid_type ', & - ' --vgrid --vgrid_type ', & + ' [--vgrid --vgrid_type ', & ' --potholes --coastal_cells] ', & ' ', & 'Check for non-advective cells. There are two types of checks available: potholes', & @@ -108,11 +102,10 @@ program topogtools 'only be needed when using a B-grid. ', & ' ', & 'Options ', & - ' --grid_type Arakawa type of horizontal grid (''B'' or ''C''; default ''B'')', & - ' --vgrid vertical grid (default ''ocean_vgrid.nc'') ', & - ' --vgrid_type can be ''mom5'' or ''mom6'' (default ''mom5'') ', & + ' --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 ', & + ' --coastal-cells check for non-advective coastal cells (for B grid) ', & ''] help_mask = [character(len=80) :: & 'usage: topogtools mask --input --output ', & @@ -132,12 +125,12 @@ program topogtools 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 --grid_type "B"', help_fill_fraction, version_text) + 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" --grid_type "B" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --potholes F & + 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" --grid_type "B" --vgrid "ocean_vgrid.nc" --vgrid_type "mom5" --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) @@ -202,19 +195,19 @@ program topogtools write(error_unit,'(a)') "ERROR: sea area fraction must be larger than 0 and smaller than 1" error stop end if - topog = topography_t(file_in, grid_type=sget('grid_type')) + topog = topography_t(file_in) call topog%fill_fraction(sea_area_fraction) call topog%update_history(get_mycommand()) call topog%write(file_out) case ('fix_nonadvective') - topog = topography_t(file_in, grid_type=sget('grid_type')) + topog = topography_t(file_in) 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, grid_type=sget('grid_type')) + topog = topography_t(file_in) call topog%nonadvective(vgrid, sget('vgrid_type'), lget('potholes'), lget('coastal-cells'), fix=.false.) case ('mask')