@@ -44,23 +44,31 @@ subroutine haversine_distmat(xa, xb, dist)
44
44
real (kind= 8 ), parameter :: erad = 6371315.0
45
45
real (kind= 8 ), dimension (:,:), intent (inout ) :: dist
46
46
external :: haversine
47
-
47
+ ! write (*,*) 'aaa'
48
48
d2r = atan2 (0 .,- 1 .)/ 180 . ! atan2(0.,-1.) == pi
49
49
m = size (xa, 1 )
50
50
n = size (xb, 1 )
51
-
51
+ allocate (xa1(m))
52
+ allocate (xa2(m))
53
+ allocate (xb1(n))
54
+ allocate (xb2(n))
52
55
xa1 = xa(:,1 )
53
56
xa2 = xa(:,2 )
54
57
xb1 = xb(:,1 )
55
58
xb2 = xb(:,2 )
59
+ ! write (*,*) 'bbb', xa1
56
60
57
61
! Loop over empty dist matrix and fill
58
- do j = 1 , m
59
- do i = 1 , n
60
- call haversine(xa1(j), xa2(j), xb1(i), xb2(i), d2r, thedist)
61
- dist(j,i) = thedist
62
- enddo
62
+ do j = 1 , m
63
+ do i = 1 , n
64
+ call haversine(xa1(j), xa2(j), xb1(i), xb2(i), d2r, thedist)
65
+ dist(j,i) = thedist
63
66
enddo
67
+ enddo
68
+ deallocate (xa1)
69
+ deallocate (xa2)
70
+ deallocate (xb1)
71
+ deallocate (xb2)
64
72
dist = dist * erad
65
73
end subroutine haversine_distmat
66
74
@@ -81,10 +89,10 @@ subroutine haversine_distvec(lon1, lat1, lon2, lat2, dist)
81
89
m = size (lon1)
82
90
83
91
! Loop over empty dist matrix and fill
84
- do i = 1 , m
85
- call haversine(lon1(i), lat1(i), lon2(i), lat2(i), d2r, thedist)
86
- dist(i) = thedist
87
- enddo
92
+ do i = 1 , m
93
+ call haversine(lon1(i), lat1(i), lon2(i), lat2(i), d2r, thedist)
94
+ dist(i) = thedist
95
+ enddo
88
96
dist = dist * erad
89
97
90
98
end subroutine haversine_distvec
0 commit comments