source: flexpart.git/src/centerofmass.f90 @ 3481cc1

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

move license from headers to a different file

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