diff --git a/src/topography.f90 b/src/topography.f90 index 6f722ec..7663bc6 100644 --- a/src/topography.f90 +++ b/src/topography.f90 @@ -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 @@ -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 @@ -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 @@ -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))) @@ -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 !------------------------------------------------------------------------- @@ -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 @@ -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 !------------------------------------------------------------------------- @@ -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)) @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/topogtools.f90 b/src/topogtools.f90 index ee238d9..de8c075 100644 --- a/src/topogtools.f90 +++ b/src/topogtools.f90 @@ -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)