| 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364 |
- ! First include the set of model-wide compiler flags
- #include "tm5.inc"
- module misctools
- implicit none
- private
- public :: dist
- contains
- !---------------------
- ! subroutine dist
- ! compute distance between points on the globe
- !---------------------
- subroutine dist( x1g, y1g, x2g, y2g, ddg )
- !---------------------------------------------------------------
- !!! compute distance of two points on the globe.
- !
- !!! this formula is correct, also for points close together.
- !
- !!! x1g, y1g : longitude and latitude of first point (degrees)
- !!! x2g, y2g : longitude and latitude of second point (degrees)
- !!! ddg : distance between these points (km)
- !----------------------------------------------------------------
- implicit none
- !
- real,intent(in) :: x1g, x2g, y1g, y2g
- real,intent(out) :: ddg
- !
- real :: x1, x2, y1, y2, dd
- real :: cc, dx, dy, dx2, dy2, dd2
- real,parameter :: pi=3.141593, rearth=6378.5
- !
- cc = pi / 180.0
- !
- x1 = x1g * cc
- x2 = x2g * cc
- y1 = y1g * cc
- y2 = y2g * cc
- !
- ! print *,'x1,x2,y1,y2', x1,x2,y1,y2
- dy = 0.5 * ( y2 - y1 )
- dy = sin(dy)
- dy2 = dy * dy
- !
- dx = 0.5 * ( x2 - x1 )
- dx = sin(dx)
- dx2 = dx * dx * cos(y1) * cos(y2)
- !
- ! print *,'dx,dy,dx2,dy2',dx,dy,dx2,dy2
- dd2 = dx2 + dy2
- dd = sqrt(dd2)
- dd = 2.0 * asin(dd)
- !
- ddg = dd / cc
- ddg = (ddg * 2.0 * pi * rearth) / 360.0
- ! print *,'dd2,dd,ddg', dd2,dd,ddg
- !
- end subroutine dist
- end module misctools
|
