1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
!-----------------------------------------------------------------------
function distance2(rlat1,rlon1,rlat2,rlon2)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
!
! SUBPROGRAM: GCDIST COMPUTE GREAT CIRCLE DISTANCE
! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM COMPUTES GREAT CIRCLE DISTANCE
! BETWEEN TWO POINTS ON THE EARTH. COORDINATES ARE GIVEN IN RADIANS!
!
! PROGRAM HISTORY LOG:
! 96-04-10 IREDELL
!
! USAGE: ...GCDIST(RLAT1,RLON1,RLAT2,RLON2)
!
! INPUT ARGUMENT LIST:
!rlat1 - REAL LATITUDE OF POINT 1 IN RADIANS
!rlon1 - REAL LONGITUDE OF POINT 1 IN RADIANS
!rlat2 - REAL LATITUDE OF POINT 2 IN RADIANS
!rlon2 - REAL LONGITUDE OF POINT 2 IN RADIANS
!
! OUTPUT ARGUMENT LIST:
!distance2 - REAL GREAT CIRCLE DISTANCE IN KM
!
! ATTRIBUTES:
! LANGUAGE: Fortran 90
!
!$$$
use par_mod, only: dp
implicit none
real :: rlat1,rlon1,rlat2,rlon2,distance2
real(kind=dp) :: clat1,clat2,slat1,slat2,cdlon,crd
real(kind=dp),parameter :: rerth=6.3712e6_dp
real(kind=dp),parameter :: pi=3.14159265358979_dp
if ((abs(rlat1-rlat2).lt.0.0003).and. &
(abs(rlon1-rlon2).lt.0.0003)) then
distance2=0.0_dp
else
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
clat1=cos(real(rlat1,kind=dp))
slat1=sin(real(rlat1,kind=dp))
clat2=cos(real(rlat2,kind=dp))
slat2=sin(real(rlat2,kind=dp))
cdlon=cos(real(rlon1-rlon2,kind=dp))
crd=slat1*slat2+clat1*clat2*cdlon
distance2=real(rerth*acos(crd)/1000.0_dp)
endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end function distance2
|
