Skip to content

Fix NN search #644

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Whyborn wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Fix NN search #644

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1cda195
Fix NN search in cable_parameters.F90
Whyborn Oct 27, 2025
b386bb1
Remove superfluous ncount
Whyborn Oct 27, 2025
cd52d20
Set max search distance to 12 degrees
Whyborn Nov 16, 2025
e8f58ba
Add comments on why its not a true NN search and why 144 was used as …
Whyborn Nov 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 27 additions & 10 deletions src/offline/cable_parameters.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -950,21 +950,29 @@ SUBROUTINE get_land_index(nlon, nlat)
INTEGER, INTENT(IN) :: nlon, nlat

! local variables
REAL :: lon2, distance, newLength
INTEGER :: ii, jj, kk, tt, ncount
REAL :: tmp_lon, distance, newLength
INTEGER :: ii, jj, kk

! range of longitudes from input file (inLon) should be -180 to 180,
! and longitude(:) has already been converted to -180 to 180 for CABLE.
landpt(:)%ilon = -999
landpt(:)%ilat = -999
ncount = 0

! This search is not a nearest neighbour in the true sense, due to the latitude dependence on
! arc length. It is only a nearest neighbour in the latitute/longitude coordinate space.
DO kk = 1, mland
distance = 5300.0 ! initialise, units are degrees
! This is the maximum distance in degrees^2 for a gridinfo point to be considered valid for
! the nearest neighbour search. If no valid neighbours are found within this 12 degree radius,
! consider the search failed and crash out.
distance = 144.0 ! initialise, units are degrees^2
Copy link
Collaborator

@har917 har917 Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure on the choice of 144. here - this need to be larger (in degrees^2) than the expected largest gridcell. 144 here corresponds to a 12 degree-squared cell (so 30x15 cells in a global simulation) - so coarse but as coarse as it could be

I would have thought that 64800 would have been the safest option for this initialisation

Copy link
Contributor Author

@Whyborn Whyborn Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yea, my only reason for choosing this value is that it's what was chosen in the CABLE-POP_TRENDY branch- effectively saying if you can't find any valid source data within a 12 degree circle, crash out. Is this something desired? Or should it be if it finds a valid neighbour anywhere on the globe, then that's ok?

Copy link
Collaborator

@har917 har917 Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

okay - that makes some sense. In effect it's trying to say if there isn't a valid lat/lon within the 12-degree circle then you need to fix your gridinfo - we could encouter this type of thing in high resolution simulations where we get some new land appearing (think Pacific/southern ocean islands).

Probably just needs a comment to explain the choice of 144

Copy link
Contributor Author

@Whyborn Whyborn Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comment has been added

Copy link
Contributor Author

@Whyborn Whyborn Nov 17, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To be clear, I don't know what the best behaviour should be here. I'll try to raise this with people using the main branch of CABLE for offline simulations and get their opinions before merging this.

DO jj = 1, nlat
DO ii = 1, nlon
IF (inVeg(ii,jj, 1) > 0) THEN
newLength = SQRT((inLon(ii) - longitude(kk))**2 &
+ (inLat(jj) - latitude(kk))**2)
tmp_lon = inLon(ii) - longitude(kk)
if (tmp_lon > 180.0) then
tmp_lon = 360 - tmp_lon
end if
newLength = (tmp_lon)**2 + (inLat(jj) - latitude(kk))**2
IF (newLength < distance) THEN
distance = newLength
landpt(kk)%ilon = ii
Expand Down Expand Up @@ -1012,7 +1020,7 @@ SUBROUTINE countPatch(nlon, nlat, npatch)
INTEGER, INTENT(IN) :: nlon, nlat, npatch

! local variables
REAL :: lon2, distance, newLength
REAL :: tmp_lon, distance, newLength
INTEGER :: ii, jj, kk, tt, ncount

! mpatch setting below introduced by rk4417 - phase2
Expand All @@ -1026,13 +1034,22 @@ SUBROUTINE countPatch(nlon, nlat, npatch)
landpt(:)%ilon = -999
landpt(:)%ilat = -999
ncount = 0

! This search is not a nearest neighbour in the true sense, due to the latitude dependence on
! arc length. It is only a nearest neighbour in the latitute/longitude coordinate space.
DO kk = 1, mland
distance = 300.0 ! initialise, units are degrees
! This is the maximum distance in degrees^2 for a gridinfo point to be considered valid for
! the nearest neighbour search. If no valid neighbours are found within this 12 degree radius,
! consider the search failed and crash out.
distance = 144.0 ! initialise, units are degrees^2
DO jj = 1, nlat
DO ii = 1, nlon
IF (inVeg(ii,jj, 1) > 0) THEN
newLength = SQRT((inLon(ii) - longitude(kk))**2 &
+ (inLat(jj) - latitude(kk))**2)
tmp_lon = inLon(ii) - longitude(kk)
if (tmp_lon > 180.0) then
tmp_lon = 360.0 - tmp_lon
end if
newLength = tmp_lon**2 + (inLat(jj) - latitude(kk))**2
IF (newLength < distance) THEN
distance = newLength
landpt(kk)%ilon = ii
Expand Down