@@ -44,23 +44,31 @@ subroutine haversine_distmat(xa, xb, dist)
4444 real (kind= 8 ), parameter :: erad = 6371315.0
4545 real (kind= 8 ), dimension (:,:), intent (inout ) :: dist
4646 external :: haversine
47-
47+ ! write (*,*) 'aaa'
4848 d2r = atan2 (0 .,- 1 .)/ 180 . ! atan2(0.,-1.) == pi
4949 m = size (xa, 1 )
5050 n = size (xb, 1 )
51-
51+ allocate (xa1(m))
52+ allocate (xa2(m))
53+ allocate (xb1(n))
54+ allocate (xb2(n))
5255 xa1 = xa(:,1 )
5356 xa2 = xa(:,2 )
5457 xb1 = xb(:,1 )
5558 xb2 = xb(:,2 )
59+ ! write (*,*) 'bbb', xa1
5660
5761! 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
6366 enddo
67+ enddo
68+ deallocate (xa1)
69+ deallocate (xa2)
70+ deallocate (xb1)
71+ deallocate (xb2)
6472 dist = dist * erad
6573 end subroutine haversine_distmat
6674
@@ -81,10 +89,10 @@ subroutine haversine_distvec(lon1, lat1, lon2, lat2, dist)
8189 m = size (lon1)
8290
8391! 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
8896 dist = dist * erad
8997
9098 end subroutine haversine_distvec
0 commit comments