centerofmass.f90 1.96 KB
Newer Older
 Matthias Langer committed Nov 28, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 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 ``````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``````