diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 6673f4b819..32f5c6714b 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/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 diff --git a/main/FatesUtilsMod.F90 b/main/FatesUtilsMod.F90 index 03537bd226..1bd644bf4c 100644 --- a/main/FatesUtilsMod.F90 +++ b/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 diff --git a/testing/great_circle/TestGreatCircle.F90 b/testing/great_circle/TestGreatCircle.F90 new file mode 100644 index 0000000000..87316d84fa --- /dev/null +++ b/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 diff --git a/testing/great_circle/WrapGreatCircleMod.F90 b/testing/great_circle/WrapGreatCircleMod.F90 new file mode 100644 index 0000000000..2728111c34 --- /dev/null +++ b/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 diff --git a/testing/great_circle/bld/README b/testing/great_circle/bld/README new file mode 100644 index 0000000000..5954904f79 --- /dev/null +++ b/testing/great_circle/bld/README @@ -0,0 +1 @@ +folder holder diff --git a/testing/great_circle/build_gc.sh b/testing/great_circle/build_gc.sh new file mode 100755 index 0000000000..c3c065c7b5 --- /dev/null +++ b/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 +