!********************************************************************** ! 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 . * !********************************************************************** subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, & rmsclust,zrms) ! i i i i o o o o o ! o o !***************************************************************************** ! * ! This routine clusters the particle position into ncluster custers. * ! Input are the longitudes (xl) and latitudes (yl) of the individual * ! points, output are the cluster mean positions (xclust,yclust). * ! Vertical positions are not directly used for the clustering. * ! * ! For clustering, the procedure described in Dorling et al. (1992) is used.* ! * ! Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992): * ! Cluster analysis: a technique for estimating the synoptic meteorological * ! controls on air and precipitation chemistry - method and applications. * ! Atmospheric Environment 26A, 2575-2581. * ! * ! * ! Author: A. Stohl * ! * ! 1 February 2002 * ! * ! Variables: * ! fclust fraction of particles belonging to each cluster * ! ncluster number of clusters to be used * ! rms total horizontal rms distance after clustering * ! rmsclust horizontal rms distance for each individual cluster * ! zrms total vertical rms distance after clustering * ! xclust,yclust, Cluster centroid positions * ! zclust * ! xl,yl,zl particle positions * ! * !***************************************************************************** use par_mod implicit none integer :: n,i,j,l,nclust(maxpart),numb(ncluster),ncl real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z real :: zclust(ncluster),distance2,distances,distancemin,rms,rmsold real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster) real :: rmsclust(ncluster) real :: zdist,zrms if (n.lt.ncluster) return rmsold=-5. ! Convert longitude and latitude from degrees to radians !******************************************************* do i=1,n nclust(i)=i xl(i)=xl(i)*pi180 yl(i)=yl(i)*pi180 end do ! Generate a seed for each cluster !********************************* do j=1,ncluster zclust(j)=0. xclust(j)=xl(j*n/ncluster) yclust(j)=yl(j*n/ncluster) end do ! Iterative loop to compute the cluster means !******************************************** do l=1,100 ! Assign each particle to a cluster: criterion minimum distance to the ! cluster mean position !********************************************************************* do i=1,n distancemin=10.**10. do j=1,ncluster distances=distance2(yl(i),xl(i),yclust(j),xclust(j)) if (distances.lt.distancemin) then distancemin=distances ncl=j endif end do nclust(i)=ncl end do ! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates, ! calculate mean position, and re-project this point onto the Earth's surface !***************************************************************************** do j=1,ncluster xav(j)=0. yav(j)=0. zav(j)=0. rmsclust(j)=0. numb(j)=0 end do rms=0. do i=1,n numb(nclust(i))=numb(nclust(i))+1 distances=distance2(yl(i),xl(i), & yclust(nclust(i)),xclust(nclust(i))) ! rms is the total rms of all particles ! rmsclust is the rms for a particular cluster !********************************************* rms=rms+distances*distances rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances ! Calculate Cartesian 3D coordinates from longitude and latitude !*************************************************************** x = cos(yl(i))*sin(xl(i)) y = -1.*cos(yl(i))*cos(xl(i)) z = sin(yl(i)) xav(nclust(i))=xav(nclust(i))+x yav(nclust(i))=yav(nclust(i))+y zav(nclust(i))=zav(nclust(i))+z end do rms=sqrt(rms/real(n)) ! Find the mean location in Cartesian coordinates !************************************************ do j=1,ncluster if (numb(j).gt.0) then rmsclust(j)=sqrt(rmsclust(j)/real(numb(j))) xav(j)=xav(j)/real(numb(j)) yav(j)=yav(j)/real(numb(j)) zav(j)=zav(j)/real(numb(j)) ! Project the point back onto Earth's surface !******************************************** xclust(j)=atan2(xav(j),-1.*yav(j)) yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j))) endif end do ! Leave the loop if the RMS distance decreases only slightly between 2 iterations !***************************************************************************** if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99 rmsold=rms end do 99 continue ! Convert longitude and latitude from radians to degrees !******************************************************* do i=1,n xl(i)=xl(i)/pi180 yl(i)=yl(i)/pi180 zclust(nclust(i))=zclust(nclust(i))+zl(i) end do do j=1,ncluster xclust(j)=xclust(j)/pi180 yclust(j)=yclust(j)/pi180 if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j)) fclust(j)=100.*real(numb(j))/real(n) end do ! Determine total vertical RMS deviation !*************************************** zrms=0. do i=1,n zdist=zl(i)-zclust(nclust(i)) zrms=zrms+zdist*zdist end do if (zrms.gt.0.) zrms=sqrt(zrms/real(n)) end subroutine clustering