source: trunk/src/clustering.f90

Last change on this file was 4, checked in by mlanger, 10 years ago
File size: 7.8 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
23       rmsclust,zrms)
24  !                      i  i  i  i   o      o      o      o     o
25  !   o      o
26  !*****************************************************************************
27  !                                                                            *
28  !   This routine clusters the particle position into ncluster custers.       *
29  !   Input are the longitudes (xl) and latitudes (yl) of the individual       *
30  !   points, output are the cluster mean positions (xclust,yclust).           *
31  !   Vertical positions are not directly used for the clustering.             *
32  !                                                                            *
33  !   For clustering, the procedure described in Dorling et al. (1992) is used.*
34  !                                                                            *
35  !   Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992):                     *
36  !   Cluster analysis: a technique for estimating the synoptic meteorological *
37  !   controls on air and precipitation chemistry - method and applications.   *
38  !   Atmospheric Environment 26A, 2575-2581.                                  *
39  !                                                                            *
40  !                                                                            *
41  !     Author: A. Stohl                                                       *
42  !                                                                            *
43  !     1 February 2002                                                        *
44  !                                                                            *
45  ! Variables:                                                                 *
46  ! fclust          fraction of particles belonging to each cluster            *
47  ! ncluster        number of clusters to be used                              *
48  ! rms             total horizontal rms distance after clustering             *
49  ! rmsclust        horizontal rms distance for each individual cluster        *
50  ! zrms            total vertical rms distance after clustering               *
51  ! xclust,yclust,  Cluster centroid positions                                 *
52  ! zclust                                                                     *
53  ! xl,yl,zl        particle positions                                         *
54  !                                                                            *
55  !*****************************************************************************
56
57  use par_mod
58
59  implicit none
60
61  integer :: n,i,j,l,nclust(maxpart),numb(ncluster),ncl
62  real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z
63  real :: zclust(ncluster),distance2,distances,distancemin,rms,rmsold
64  real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster)
65  real :: rmsclust(ncluster)
66  real :: zdist,zrms
67
68
69
70  if (n.lt.ncluster) return
71  rmsold=-5.
72
73  ! Convert longitude and latitude from degrees to radians
74  !*******************************************************
75
76  do i=1,n
77    nclust(i)=i
78    xl(i)=xl(i)*pi180
79    yl(i)=yl(i)*pi180
80  end do
81
82
83  ! Generate a seed for each cluster
84  !*********************************
85
86  do j=1,ncluster
87    zclust(j)=0.
88    xclust(j)=xl(j*n/ncluster)
89    yclust(j)=yl(j*n/ncluster)
90  end do
91
92
93  ! Iterative loop to compute the cluster means
94  !********************************************
95
96  do l=1,100
97
98  ! Assign each particle to a cluster: criterion minimum distance to the
99  ! cluster mean position
100  !*********************************************************************
101
102
103    do i=1,n
104      distancemin=10.**10.
105      do j=1,ncluster
106        distances=distance2(yl(i),xl(i),yclust(j),xclust(j))
107        if (distances.lt.distancemin) then
108          distancemin=distances
109          ncl=j
110        endif
111      end do
112      nclust(i)=ncl
113    end do
114
115
116  ! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates,
117  ! calculate mean position, and re-project this point onto the Earth's surface
118  !*****************************************************************************
119
120    do j=1,ncluster
121      xav(j)=0.
122      yav(j)=0.
123      zav(j)=0.
124      rmsclust(j)=0.
125      numb(j)=0
126    end do
127    rms=0.
128
129    do i=1,n
130      numb(nclust(i))=numb(nclust(i))+1
131      distances=distance2(yl(i),xl(i), &
132           yclust(nclust(i)),xclust(nclust(i)))
133
134  ! rms is the total rms of all particles
135  ! rmsclust is the rms for a particular cluster
136  !*********************************************
137
138      rms=rms+distances*distances
139      rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
140
141  ! Calculate Cartesian 3D coordinates from longitude and latitude
142  !***************************************************************
143
144      x = cos(yl(i))*sin(xl(i))
145      y = -1.*cos(yl(i))*cos(xl(i))
146      z = sin(yl(i))
147      xav(nclust(i))=xav(nclust(i))+x
148      yav(nclust(i))=yav(nclust(i))+y
149      zav(nclust(i))=zav(nclust(i))+z
150    end do
151
152    rms=sqrt(rms/real(n))
153
154
155  ! Find the mean location in Cartesian coordinates
156  !************************************************
157
158    do j=1,ncluster
159      if (numb(j).gt.0) then
160        rmsclust(j)=sqrt(rmsclust(j)/real(numb(j)))
161        xav(j)=xav(j)/real(numb(j))
162        yav(j)=yav(j)/real(numb(j))
163        zav(j)=zav(j)/real(numb(j))
164
165  ! Project the point back onto Earth's surface
166  !********************************************
167
168        xclust(j)=atan2(xav(j),-1.*yav(j))
169        yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
170      endif
171    end do
172
173
174  ! Leave the loop if the RMS distance decreases only slightly between 2 iterations
175  !*****************************************************************************
176
177    if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
178    rmsold=rms
179
180  end do
181
18299   continue
183
184  ! Convert longitude and latitude from radians to degrees
185  !*******************************************************
186
187  do i=1,n
188    xl(i)=xl(i)/pi180
189    yl(i)=yl(i)/pi180
190    zclust(nclust(i))=zclust(nclust(i))+zl(i)
191  end do
192
193  do j=1,ncluster
194    xclust(j)=xclust(j)/pi180
195    yclust(j)=yclust(j)/pi180
196    if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j))
197    fclust(j)=100.*real(numb(j))/real(n)
198  end do
199
200  ! Determine total vertical RMS deviation
201  !***************************************
202
203  zrms=0.
204  do i=1,n
205    zdist=zl(i)-zclust(nclust(i))
206    zrms=zrms+zdist*zdist
207  end do
208  if (zrms.gt.0.) zrms=sqrt(zrms/real(n))
209
210end subroutine clustering
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG