centerofmass.f90 3.37 KB
 Ignacio Pisso committed Nov 14, 2019 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ``````!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** `````` Matthias Langer committed Nov 28, 2012 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 ``````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``````