From 1cda195579a9637050c82b78cb853acd67230cc1 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Mon, 27 Oct 2025 22:02:32 +1100 Subject: [PATCH 1/4] Fix NN search in cable_parameters.F90 --- src/offline/cable_parameters.F90 | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 44ebc6b9d..e7d1c66e2 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -950,8 +950,8 @@ 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. @@ -963,8 +963,11 @@ SUBROUTINE get_land_index(nlon, nlat) 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 @@ -1012,7 +1015,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 @@ -1027,12 +1030,15 @@ SUBROUTINE countPatch(nlon, nlat, npatch) landpt(:)%ilat = -999 ncount = 0 DO kk = 1, mland - distance = 300.0 ! initialise, units are degrees + distance = 300.0 ** 2 ! initialise, units are degrees 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 From b386bb1f6b4c25c0df803fbcd5a17af9e5216e10 Mon Sep 17 00:00:00 2001 From: Whyborn Date: Tue, 28 Oct 2025 08:53:13 +1100 Subject: [PATCH 2/4] Remove superfluous ncount --- src/offline/cable_parameters.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index e7d1c66e2..7bd3a096e 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -957,7 +957,6 @@ SUBROUTINE get_land_index(nlon, nlat) ! and longitude(:) has already been converted to -180 to 180 for CABLE. landpt(:)%ilon = -999 landpt(:)%ilat = -999 - ncount = 0 DO kk = 1, mland distance = 5300.0 ! initialise, units are degrees DO jj = 1, nlat From cd52d20cf28d7f70cd270c98828eb4cfad77e15d Mon Sep 17 00:00:00 2001 From: Whyborn Date: Mon, 17 Nov 2025 10:59:08 +1100 Subject: [PATCH 3/4] Set max search distance to 12 degrees --- src/offline/cable_parameters.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 7bd3a096e..404b2b427 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -958,7 +958,7 @@ SUBROUTINE get_land_index(nlon, nlat) landpt(:)%ilon = -999 landpt(:)%ilat = -999 DO kk = 1, mland - distance = 5300.0 ! initialise, units are degrees + distance = 144.0 ! initialise, units are degrees^2 DO jj = 1, nlat DO ii = 1, nlon IF (inVeg(ii,jj, 1) > 0) THEN @@ -1029,7 +1029,7 @@ SUBROUTINE countPatch(nlon, nlat, npatch) landpt(:)%ilat = -999 ncount = 0 DO kk = 1, mland - distance = 300.0 ** 2 ! initialise, units are degrees + distance = 144.0 ! initialise, units are degrees^2 DO jj = 1, nlat DO ii = 1, nlon IF (inVeg(ii,jj, 1) > 0) THEN From e8f58ba016194020a164ba8ff3bea7b8c1e6b7bd Mon Sep 17 00:00:00 2001 From: Whyborn Date: Mon, 17 Nov 2025 13:36:22 +1100 Subject: [PATCH 4/4] Add comments on why its not a true NN search and why 144 was used as initial distance^2 --- src/offline/cable_parameters.F90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/offline/cable_parameters.F90 b/src/offline/cable_parameters.F90 index 404b2b427..f3db9f2df 100644 --- a/src/offline/cable_parameters.F90 +++ b/src/offline/cable_parameters.F90 @@ -957,7 +957,13 @@ SUBROUTINE get_land_index(nlon, nlat) ! and longitude(:) has already been converted to -180 to 180 for CABLE. landpt(:)%ilon = -999 landpt(:)%ilat = -999 + + ! 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 + ! 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 @@ -1028,7 +1034,13 @@ 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 + ! 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