# source:trunk/src/centerofmass.f90@28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 3.4 KB
Line
1!**********************************************************************
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine centerofmass(xl,yl,n,xcenter,ycenter)
23  !                        i  i  i    o       o
24  !*****************************************************************************
25  !                                                                            *
26  !   This routine calculates the center of mass of n points on the Earth.     *
27  !   Input are the longitudes (xl) and latitudes (yl) of the individual       *
28  !   points, output is the longitude and latitude of the centre of mass.      *
29  !                                                                            *
30  !     Author: A. Stohl                                                       *
31  !                                                                            *
32  !     24 January 2002                                                        *
33  !                                                                            *
34  !*****************************************************************************
35
36  use par_mod
37
38  implicit none
39
40  integer :: n,l
41  real :: xl(n),yl(n),xll,yll,xav,yav,zav,x,y,z,xcenter,ycenter
42
43
44  xav=0.
45  yav=0.
46  zav=0.
47
48  do l=1,n
49
50  ! Convert longitude and latitude from degrees to radians
51  !*******************************************************
52
53    xll=xl(l)*pi180
54    yll=yl(l)*pi180
55
56  ! Calculate 3D coordinates from longitude and latitude
57  !*****************************************************
58
59    x = cos(yll)*sin(xll)
60    y = -1.*cos(yll)*cos(xll)
61    z = sin(yll)
62
63
64  ! Find the mean location in Cartesian coordinates
65  !************************************************
66
67    xav=xav+x
68    yav=yav+y
69    zav=zav+z
70  end do
71
72  xav=xav/real(n)
73  yav=yav/real(n)
74  zav=zav/real(n)
75
76
77  ! Project the point back onto Earth's surface
78  !********************************************
79
80  xcenter=atan2(xav,-1.*yav)
81  ycenter=atan2(zav,sqrt(xav*xav+yav*yav))
82
83  ! Convert back to degrees
84  !************************
85
86  xcenter=xcenter/pi180
87  ycenter=ycenter/pi180
88
89end subroutine centerofmass
Note: See TracBrowser for help on using the repository browser.