Skip to content

Commit

Permalink
Merge pull request #30 from COSIMA/improve_workflow_and_output
Browse files Browse the repository at this point in the history 
Improve workflow and output
  • Loading branch information
@micaeljtoliveira
micaeljtoliveira authored Feb 2, 2023
2 parents c6536ae + e51682d commit d4dcf41
Show file tree
Hide file tree
Showing 2 changed files with 102 additions and 41 deletions.
139 changes: 100 additions & 39 deletions src/topography.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ module topography
procedure :: copy => topography_copy
generic :: assignment(=) => copy
procedure :: update_history => topography_update_history
procedure :: number_seas => topography_number_seas
procedure :: deseas => topography_deseas
procedure :: fill_fraction => topography_fill_fraction
procedure :: nonadvective => topography_nonadvective
Expand Down Expand Up @@ -220,23 +221,34 @@ subroutine topography_update_history(this, command)
end subroutine topography_update_history

!-------------------------------------------------------------------------
subroutine topography_deseas(this)
class(topography_t), intent(inout) :: this
subroutine topography_number_seas(this, sea_number, number_of_seas, silent)
class(topography_t), intent(in) :: this
integer(int16), intent(out), target, optional :: sea_number(:,:)
integer(int32), intent(out), optional :: number_of_seas
logical, intent(in), optional :: silent

integer(int32) :: i, j, counter, its, its1, its2, sea_num, iblock, jblock
integer(int32) :: i, j, counter, its, its1, its2, sea_num
integer(int32) :: im, ip, jm, jp, land

integer(int32) :: ncid, sea_id, dids(2) ! NetCDF ids

integer(int16) :: new_sea
integer(int16), allocatable :: sea(:,:)
integer(int16), pointer :: sea(:,:)

logical :: choke_west, choke_east, choke_north, choke_south
logical :: silent_, choke_west, choke_east, choke_north, choke_south

allocate(sea(this%nxt, this%nyt))
integer(int16), parameter :: MAX_ITER = 150

if (present(sea_number)) then
sea => sea_number
else
allocate(sea(this%nxt, this%nyt))
end if
if (present(silent)) then
silent_ = silent
else
silent_ = .false.
end if

! Do
write(output_unit,'(a)') "Removing seas"
land = this%nxt + this%nyt + 1
sea = land
do j = 1, this%nyt
Expand All @@ -246,7 +258,9 @@ subroutine topography_deseas(this)
if (all(this%depth(:, j) > 0.0)) sea(:, j) = 0 ! Southern Ocean all water across
end do

do its = 1, 150 ! Only need high number after massive editing session with fjords. Normally 10 or so sweeps works.
if (.not. silent_) write(output_unit,'(a)') " # Iter # Changes # Seas"

do its = 1, MAX_ITER ! Only need high number after massive editing session with fjords. Normally 10 or so sweeps works.
counter = 0
sea_num = 0

Expand Down Expand Up @@ -363,24 +377,40 @@ subroutine topography_deseas(this)
counter = counter + 1
end if
end do
write(output_unit,*) counter, sea_num

if (.not. silent_) write(output_unit,*) its, counter, sea_num + 1

! If we only have one sea or no changes are made we are finished.
if (counter == 0 .or. sea_num == 1) exit
if (counter == 0 .or. sea_num + 1 == 1) exit
if (its == MAX_ITER) then
write(output_unit, '(a)') "WARNING: could not number all the seas. Algorithm reached maximum number of iterations."
end if
end do

write(output_unit,'(a)') "Done"
if (present(sea_number)) then
nullify(sea)
else
deallocate(sea)
end if

! Write out new topography
do j = 1, this%nyt
do i = 1, this%nxt
if (sea(i, j) > 0) then
this%depth(i, j) = MISSING_VALUE
this%frac(i, j) = MISSING_VALUE
end if
end do
end do
this%lakes_removed = "yes"
if (present(number_of_seas)) then
number_of_seas = sea_num + 1
end if

end subroutine topography_number_seas

!-------------------------------------------------------------------------
subroutine topography_deseas(this)
class(topography_t), intent(inout) :: this

integer(int32) :: ncid, sea_id, dids(2) ! NetCDF ids
integer(int16), allocatable :: sea(:,:)

allocate(sea(this%nxt, this%nyt))

write(output_unit,'(a)') "Numbering seas"

call this%number_seas(sea_number=sea)

call handle_error(nf90_create(trim('sea_num.nc'), ior(nf90_netcdf4, nf90_clobber), ncid))
call handle_error(nf90_def_dim(ncid, 'nx', this%nxt, dids(1)))
Expand All @@ -391,6 +421,15 @@ subroutine topography_deseas(this)
call handle_error(nf90_put_var(ncid, sea_id, sea))
call handle_error(nf90_close(ncid))

write(output_unit,'(a)') "Removing seas"
where(sea > 0)
this%depth = MISSING_VALUE
this%frac = MISSING_VALUE
end where
this%lakes_removed = "yes"

deallocate(sea)

end subroutine topography_deseas

!-------------------------------------------------------------------------
Expand Down Expand Up @@ -444,6 +483,8 @@ subroutine topography_fill_fraction(this, sea_area_fraction)
class(topography_t), intent(inout) :: this
real(real32), intent(in) :: sea_area_fraction

integer(int32) :: nseas

write(output_unit,'(a,f7.2)') "Filling cells that have a sea area fraction smaller than ", sea_area_fraction

if (any(this%frac < sea_area_fraction)) then
Expand All @@ -453,13 +494,15 @@ subroutine topography_fill_fraction(this, sea_area_fraction)
end where

if (this%lakes_removed == 'yes') then
! We might have created new lakes, so update the corresponding attribute
this%lakes_removed = 'no'
! Check if new seas have been 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."
this%lakes_removed = 'no'
end if
end if
end if

write(output_unit,'(a)') "Done"

end subroutine topography_fill_fraction

!-------------------------------------------------------------------------
Expand All @@ -471,13 +514,14 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
real(real32), allocatable :: depth_halo(:,:)
integer(int32), allocatable :: num_levels(:,:)
real(real32), allocatable :: zw(:), zeta(:)
integer(int32) :: passes, i, j, k, ni, nj, nzeta, nz, its, counter
integer(int32) :: passes, i, j, k, ni, nj, nzeta, nz, its, coastal_counter, potholes_counter
integer(int32) :: ncid, vid
integer(int32) :: dids(2)
logical :: se, sw, ne, nw ! .TRUE. if C-cell centre is shallower than T cell centre.
logical :: changes_made = .false.
integer(int32) :: kse, ksw, kne, knw, kmu_max
integer(int32) :: im, ip, jm, jp
integer(int32) :: nseas

call handle_error(nf90_open(trim(vgrid), nf90_nowrite, ncid))
call handle_error(nf90_inq_varid(ncid, 'zeta', vid))
Expand All @@ -496,11 +540,15 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)

if (fix) then
passes = 20
write(output_unit,'(a)') "Fixing non-advective cells"
else
write(output_unit,'(a)') "Checking for non-advective cells"
passes = 1
end if
do its = 1, passes
counter = 0
if (fix) write(output_unit,'(a,i0)') " Pass # ", its
coastal_counter = 0
potholes_counter = 0
num_levels = 0

depth_halo = 0
Expand All @@ -522,6 +570,10 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
end do

if (coastal) then
if (.not. fix) then
write(output_unit,'(a)') " Coastal cells"
write(output_unit,'(a)') " i j Depth"
end if
do j = 2, this%nyt - 1
do i = 1, this%nxt
if (depth_halo(i, j) > 0.5) then
Expand All @@ -534,18 +586,23 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
depth_halo(i, j) = 0.0
this%frac(i, j) = 0.0
num_levels(i, j) = 0
else
write(output_unit,*) i, j, 0.0 ,' ! nonadvective'
end if
counter = counter + 1
write(output_unit,*) i, j, 0.0 ,' ! nonadvective'
coastal_counter = coastal_counter + 1
end if
end if
end do
end do

write(output_unit,*) '1', counter
write(output_unit,'(a,i0,a)') ' Found ', coastal_counter, ' non-advective coastal cells'
end if

if (potholes) then
if (.not. fix) then
write(output_unit,'(a)') " Potholes"
write(output_unit,'(a)') " i j Level Max. halo level"
end if
do j = 2, this%nyt
jm = j - 1
jp = j + 1
Expand All @@ -565,25 +622,29 @@ subroutine topography_nonadvective(this, vgrid, potholes, coastal, fix)
num_levels(i, j) = kmu_max
depth_halo(i, j) = zw(kmu_max)
else
write(output_unit,*) i, j, ' nonadvective, Deep', num_levels(i, j), kmu_max
write(output_unit,*) i, j, num_levels(i, j), kmu_max
end if
counter = counter + 1
potholes_counter = potholes_counter + 1
end if
end if
end do
end do
write(output_unit,*) counter
write(output_unit,'(a,i0,a)') ' Found ', potholes_counter, ' non-advective potholes'
end if
if (counter > 0 .and. fix) changes_made = .true.
if ((coastal_counter > 0 .or. potholes_counter > 0) .and. fix) changes_made = .true.
this%depth = depth_halo(1:this%nxt, 1:this%nyt)
if (counter == 0 .and. fix) exit
if (coastal_counter == 0 .and. potholes_counter == 0 .and. fix) exit
end do

if (fix .and. (coastal .or. potholes)) then
this%nonadvective_cells_removed = 'yes'
if (changes_made .and. this%lakes_removed == 'yes') then
! We might have created new lakes, so rerun deseas
call this%deseas()
! Check if new lakes were created new lakes
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."
this%lakes_removed = 'no'
end if
end if
end if

Expand Down
4 changes: 2 additions & 2 deletions src/topogtools.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -116,10 +116,10 @@ program topogtools
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 T --coastal-cells F', &
call set_args('--input:i "unset" --output:o "unset" --vgrid "ocean_vgrid.nc" --potholes F --coastal-cells F', &
help_fix_nonadvective, version_text)
case ('check_nonadvective')
call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --potholes T --coastal-cells F', &
call set_args('--input:i "unset" --vgrid "ocean_vgrid.nc" --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)
Expand Down

0 comments on commit d4dcf41

Please sign in to comment.