1 | ! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt |
---|
2 | ! SPDX-License-Identifier: GPL-3.0-or-later |
---|
3 | |
---|
4 | !----------------------------------------------------------------------- |
---|
5 | function distance(rlat1,rlon1,rlat2,rlon2) |
---|
6 | |
---|
7 | !$$$ SUBPROGRAM DOCUMENTATION BLOCK |
---|
8 | ! |
---|
9 | ! SUBPROGRAM: GCDIST COMPUTE GREAT CIRCLE DISTANCE |
---|
10 | ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 |
---|
11 | ! |
---|
12 | ! ABSTRACT: THIS SUBPROGRAM COMPUTES GREAT CIRCLE DISTANCE |
---|
13 | ! BETWEEN TWO POINTS ON THE EARTH. |
---|
14 | ! |
---|
15 | ! PROGRAM HISTORY LOG: |
---|
16 | ! 96-04-10 IREDELL |
---|
17 | ! |
---|
18 | ! USAGE: ...GCDIST(RLAT1,RLON1,RLAT2,RLON2) |
---|
19 | ! |
---|
20 | ! INPUT ARGUMENT LIST: |
---|
21 | !rlat1 - REAL LATITUDE OF POINT 1 IN DEGREES |
---|
22 | !rlon1 - REAL LONGITUDE OF POINT 1 IN DEGREES |
---|
23 | !rlat2 - REAL LATITUDE OF POINT 2 IN DEGREES |
---|
24 | !rlon2 - REAL LONGITUDE OF POINT 2 IN DEGREES |
---|
25 | ! |
---|
26 | ! OUTPUT ARGUMENT LIST: |
---|
27 | !distance - REAL GREAT CIRCLE DISTANCE IN KILOMETERS |
---|
28 | ! |
---|
29 | ! ATTRIBUTES: |
---|
30 | ! LANGUAGE: Fortran 90 |
---|
31 | ! |
---|
32 | !$$$ |
---|
33 | |
---|
34 | use par_mod, only: dp |
---|
35 | |
---|
36 | implicit none |
---|
37 | |
---|
38 | real :: rlat1,rlon1,rlat2,rlon2,distance |
---|
39 | real(kind=dp) :: clat1,clat2,slat1,slat2,cdlon,crd |
---|
40 | real(kind=dp),parameter :: rerth=6.3712e6_dp |
---|
41 | real(kind=dp),parameter :: pi=3.14159265358979_dp, dpr=180.0_dp/pi |
---|
42 | ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
---|
43 | if ((abs(rlat1-rlat2).lt.0.03).and. & |
---|
44 | (abs(rlon1-rlon2).lt.0.03)) then |
---|
45 | distance=0. |
---|
46 | else |
---|
47 | clat1=cos(real(rlat1,kind=dp)/dpr) |
---|
48 | slat1=sin(real(rlat1,kind=dp)/dpr) |
---|
49 | clat2=cos(real(rlat2,kind=dp)/dpr) |
---|
50 | slat2=sin(real(rlat2,kind=dp)/dpr) |
---|
51 | cdlon=cos(real((rlon1-rlon2),kind=dp)/dpr) |
---|
52 | crd=slat1*slat2+clat1*clat2*cdlon |
---|
53 | distance=real(rerth*acos(crd)/1000.0_dp) |
---|
54 | endif |
---|
55 | ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
---|
56 | end function distance |
---|