source: flexpart.git/src/centerofmass.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 2.1 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine centerofmass(xl,yl,n,xcenter,ycenter)
5  !                        i  i  i    o       o
6  !*****************************************************************************
7  !                                                                            *
8  !   This routine calculates the center of mass of n points on the Earth.     *
9  !   Input are the longitudes (xl) and latitudes (yl) of the individual       *
10  !   points, output is the longitude and latitude of the centre of mass.      *
11  !                                                                            *
12  !     Author: A. Stohl                                                       *
13  !                                                                            *
14  !     24 January 2002                                                        *
15  !                                                                            *
16  !*****************************************************************************
17
18  use par_mod
19
20  implicit none
21
22  integer :: n,l
23  real :: xl(n),yl(n),xll,yll,xav,yav,zav,x,y,z,xcenter,ycenter
24
25
26  xav=0.
27  yav=0.
28  zav=0.
29
30  do l=1,n
31
32  ! Convert longitude and latitude from degrees to radians
33  !*******************************************************
34
35    xll=xl(l)*pi180
36    yll=yl(l)*pi180
37
38  ! Calculate 3D coordinates from longitude and latitude
39  !*****************************************************
40
41    x = cos(yll)*sin(xll)
42    y = -1.*cos(yll)*cos(xll)
43    z = sin(yll)
44
45
46  ! Find the mean location in Cartesian coordinates
47  !************************************************
48
49    xav=xav+x
50    yav=yav+y
51    zav=zav+z
52  end do
53
54  xav=xav/real(n)
55  yav=yav/real(n)
56  zav=zav/real(n)
57
58
59  ! Project the point back onto Earth's surface
60  !********************************************
61
62  xcenter=atan2(xav,-1.*yav)
63  ycenter=atan2(zav,sqrt(xav*xav+yav*yav))
64
65  ! Convert back to degrees
66  !************************
67
68  xcenter=xcenter/pi180
69  ycenter=ycenter/pi180
70
71end subroutine centerofmass
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG