fix calculation of the distance between model grid and inventory sites

main/FatesInventoryInitMod.F90

@@ -77,6 +77,7 @@ module FatesInventoryInitMod
    use FatesConstantsMod,   only : fates_unset_int
    use EDCanopyStructureMod, only : canopy_summarization, canopy_structure
    use FatesRadiationMemMod, only : num_swb
+   use FatesUtilsMod,       only : GreatCircleDist
    implicit none
    private
 
@@ -105,6 +106,9 @@ module FatesInventoryInitMod
                                                            ! defined in model memory and a physical
                                                            ! site listed in the file
 
+   real(r8), parameter :: max_site_adjacency_m  = 5500._r8 ! 0.05 deg roughly equals 5.5k meters
+                                                           ! at the two tropic lines (111 km/deg)
+
    logical, parameter :: do_inventory_out = .false.
 
 
@@ -148,6 +152,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in)
       real(r8)                                     :: age_init             ! dummy value for creating a patch
       real(r8)                                     :: area_init            ! dummy value for creating a patch
       integer                                      :: s                    ! site index
+      integer                                      :: i                    ! inventory site index
       integer                                      :: ipa                  ! patch index
       integer                                      :: iv, ft, ic
       integer                                      :: total_cohorts        ! cohort counter for error checking
@@ -157,6 +162,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in)
 
       real(r8),                        allocatable :: inv_lat_list(:)      ! list of lat coords
       real(r8),                        allocatable :: inv_lon_list(:)      ! list of lon coords
+      real(r8),                        allocatable :: delta_site_list(:)   ! list of differences between model site and inv site (m)
       integer                                      :: invsite              ! index of inventory site
                                                                            ! closest to actual site
       integer                                      :: el                   ! loop counter for number of elements
@@ -210,7 +216,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in)
       allocate(inv_css_list(nfilesites))
       allocate(inv_lat_list(nfilesites))
       allocate(inv_lon_list(nfilesites))
-
+      allocate(delta_site_list(nfilesites))
 
       ! Check through the sites that are listed and do some sanity checks
       ! ------------------------------------------------------------------------------------------
@@ -231,17 +237,22 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in)
       ! For each site, identify the most proximal PSS/CSS couplet, read-in the data
       ! allocate linked lists and assign to memory
       do s = 1, nsites
-         invsite = &
-               minloc( (sites(s)%lat-inv_lat_list(:))**2.0_r8 + &
-               (sites(s)%lon-inv_lon_list(:))**2.0_r8 , dim=1)
+
+         do i = 1,nfilesites
+            ! Great circle calculates the distance in meters between two points
+            ! on the earth and also factors in the earth's curvature
+            delta_site_list(i) = &
+                 GreatCircleDist(sites(s)%lon,inv_lon_list(i),sites(s)%lat,inv_lat_list(i))
+         end do
+
+         invsite = minloc(delta_site_list(:), dim=1)
 
          ! Do a sanity check on the distance separation between physical site and model site
-         if ( sqrt( (sites(s)%lat-inv_lat_list(invsite))**2.0_r8 + &
-               (sites(s)%lon-inv_lon_list(invsite))**2.0_r8 ) > max_site_adjacency_deg ) then
+         if ( delta_site_list(invsite) > max_site_adjacency_m ) then
             write(fates_log(), *) 'Model site at lat:',sites(s)%lat,' lon:',sites(s)%lon
             write(fates_log(), *) 'has no reasonably proximal site in the inventory site list.'
             write(fates_log(), *) 'Closest is at lat:',inv_lat_list(invsite),' lon:',inv_lon_list(invsite)
-            write(fates_log(), *) 'Separation must be less than ',max_site_adjacency_deg,' degrees'
+            write(fates_log(), *) 'Separation must be less than ',max_site_adjacency_m,' meters'
             write(fates_log(), *) 'Exiting'
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if
@@ -481,7 +492,7 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in)
 
       end do
 
-      deallocate(inv_format_list, inv_pss_list, inv_css_list, inv_lat_list, inv_lon_list)
+      deallocate(inv_format_list, inv_pss_list, inv_css_list, inv_lat_list, inv_lon_list,delta_site_list)
 
       return
    end subroutine initialize_sites_by_inventory

main/FatesUtilsMod.F90

@@ -7,7 +7,6 @@ module FatesUtilsMod
   use FatesGlobals, only      : fates_log
   use FatesConstantsMod, only : nearzero
   use FatesGlobals,      only : endrun => fates_endrun
-
   use shr_log_mod , only      : errMsg => shr_log_errMsg
 
   implicit none
@@ -21,6 +20,7 @@ module FatesUtilsMod
   public :: QuadraticRootsNSWC
   public :: QuadraticRootsSridharachary
   public :: ArrayNint
+  public :: GreatCircleDist
 
   character(len=*), parameter, private :: sourcefile = &
        __FILE__
@@ -130,10 +130,15 @@ real(r8) function GreatCircleDist(slons,slonf,slats,slatf)
      real(r8)     :: x
      real(r8)     :: y
      !---------------------------------------------------------------------------------------!
-
+
+     ! ---- Make sure that longitudes are using the same convention (-180,180)
+
+     lons = modulo(slons + 180.0_r8,360.0_r8)-180.0_r8
+     lonf = modulo(slonf + 180.0_r8,360.0_r8)-180.0_r8
+
      !----- Convert the co-ordinates to double precision and to radians. --------------------!
-     lons = slons * rad_per_deg
-     lonf = slonf * rad_per_deg
+     lons = lons * rad_per_deg
+     lonf = lonf * rad_per_deg
      lats = slats * rad_per_deg
      latf = slatf * rad_per_deg
      dlon = lonf - lons

testing/great_circle/TestGreatCircle.F90

@@ -0,0 +1,34 @@
+program TestGreatCircle
+
+  use FatesConstantsMod, only : r8 => fates_r8
+  use FatesUtilsMod,     only : GreatCircleDist
+  implicit none
+
+  ! Variable declarations
+  real(r8) :: inv_lat_list(5)      ! list of lat coords
+  real(r8) :: inv_lon_list(5)      ! list of lon coords
+  real(r8) :: site_lat_list(5)      ! list of lat coords
+  real(r8) :: site_lon_list(5)      ! list of lon coords
+  real(r8) :: delta_site_list(5)
+
+  integer :: i,s
+  integer :: invsite
+
+  inv_lat_list = (/-89._r8, -20._r8, 0._r8, 60._r8, 90._r8/)
+  inv_lon_list = (/-170._r8, -20._r8, 120.5_r8, 210._r8, 90._r8/)
+
+  site_lat_list = (/-19._r8, -89._r8, 1._r8, 63._r8, 88._r8/)
+  site_lon_list = (/-21._r8, -171._r8, 118.5_r8, 214._r8, 78._r8/)
+
+  do s=1,5
+     do i =1,5
+        delta_site_list(i) = &
+             GreatCircleDist(site_lon_list(s),inv_lon_list(i),site_lat_list(s),inv_lat_list(i))
+     end do
+     invsite = minloc(delta_site_list(:), dim=1)
+     write(*,'(A,2(F6.1),A,2(F6.1))') "closest to ", site_lat_list(s),site_lon_list(s), &
+          " is: ",inv_lat_list(invsite),inv_lon_list(invsite)
+     write(*,'(A,F6.1,A)') "   with distance of:",delta_site_list(invsite)/1000._r8," km"
+  end do
+
+end program TestGreatCircle

testing/great_circle/WrapGreatCircleMod.F90

@@ -0,0 +1,45 @@
+module shr_log_mod
+  use iso_c_binding, only : c_char
+  use iso_c_binding, only : c_int
+
+  public :: shr_log_errMsg
+
+contains
+  function shr_log_errMsg(source, line) result(ans)
+    character(kind=c_char,len=*), intent(in) :: source
+    integer(c_int), intent(in) :: line
+    character(kind=c_char,len=4) :: cline  ! character version of int
+    character(kind=c_char,len=128) :: ans
+
+    write(cline,'(I4)') line
+    ans = "source: " // trim(source) // " line: "// trim(cline)
+
+  end function shr_log_errMsg
+
+end module shr_log_mod
+
+module FatesGlobals
+
+  use iso_c_binding, only : c_char
+  use iso_c_binding, only : c_int
+  use FatesConstantsMod, only : r8 => fates_r8
+
+  integer :: stdo_unit = 6
+
+contains
+
+  integer function fates_log()
+    fates_log = 6
+  end function fates_log
+
+  subroutine fates_endrun(msg)
+
+    implicit none
+    character(len=*), intent(in) :: msg    ! string to be printed
+
+    write(stdo_unit,*) msg
+
+    stop
+
+  end subroutine fates_endrun
+end module FatesGlobals

testing/great_circle/bld/README

@@ -0,0 +1 @@
+folder holder

testing/great_circle/build_gc.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+
+# Path to FATES src
+
+FC='gfortran'
+
+#F_OPTS="-fPIC -O3 -llapack"
+F_OPTS="-g -fPIC"
+F_OBJ_OPTS="-shared"
+
+FATES_PATH='../../'
+
+#F_OPTS="-fPIC -O0 -g -ffpe-trap=zero,overflow,underflow -fbacktrace -fbounds-check -Wall"
+
+MOD_FLAG="-J"
+
+rm -f bld/*.o
+rm -f bld/*.so
+rm -f bld/*.mod
+rm -f bld/*.a
+
+# Build dgesv from lapack
+${FC} ${F_OPTS} -c -I bld/ -J./bld/ -o bld/libFatesConstantsMod.so ${FATES_PATH}/main/FatesConstantsMod.F90
+${FC} ${F_OPTS} -c -I bld/ -J./bld/ -o bld/libWrapGreatCircleMod.so WrapGreatCircleMod.F90
+${FC} ${F_OPTS} -c -I bld/ -J./bld/ -o bld/libFatesUtilsMod.so ${FATES_PATH}/main/FatesUtilsMod.F90
+${FC} ${F_OPTS} -I bld/ -J./bld/ -L./bld/ -lFatesConstantsMod -lFatesUtilsMod -lWrapGreatCircleMod -o test_gc TestGreatCircle.F90
+