! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later

subroutine centerofmass(xl,yl,n,xcenter,ycenter)
  !                        i  i  i    o       o
  !*****************************************************************************
  !                                                                            *
  !   This routine calculates the center of mass of n points on the Earth.     *
  !   Input are the longitudes (xl) and latitudes (yl) of the individual       *
  !   points, output is the longitude and latitude of the centre of mass.      *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     24 January 2002                                                        *
  !                                                                            *
  !*****************************************************************************

  use par_mod

  implicit none

  integer :: n,l
  real :: xl(n),yl(n),xll,yll,xav,yav,zav,x,y,z,xcenter,ycenter


  xav=0.
  yav=0.
  zav=0.

  do l=1,n

  ! Convert longitude and latitude from degrees to radians
  !*******************************************************

    xll=xl(l)*pi180
    yll=yl(l)*pi180

  ! Calculate 3D coordinates from longitude and latitude
  !*****************************************************

    x = cos(yll)*sin(xll)
    y = -1.*cos(yll)*cos(xll)
    z = sin(yll)


  ! Find the mean location in Cartesian coordinates
  !************************************************

    xav=xav+x
    yav=yav+y
    zav=zav+z
  end do

  xav=xav/real(n)
  yav=yav/real(n)
  zav=zav/real(n)


  ! Project the point back onto Earth's surface
  !********************************************

  xcenter=atan2(xav,-1.*yav)
  ycenter=atan2(zav,sqrt(xav*xav+yav*yav))

  ! Convert back to degrees
  !************************

  xcenter=xcenter/pi180
  ycenter=ycenter/pi180

end subroutine centerofmass