source: branches/flexpart91_hasod/src_parallel/clustering.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

File size: 8.0 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,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  integer,allocatable,dimension(:) :: nclust
68  integer :: stat
69
70
71  ! allocate local arrays on heap
72  allocate(nclust(maxpart),stat=stat)
73  if (stat.ne.0) print*,'ERROR: could not allocate nclust'
74
75  if (n.lt.ncluster) return
76  rmsold=-5.
77
78  ! Convert longitude and latitude from degrees to radians
79  !*******************************************************
80
81  do i=1,n
82    nclust(i)=i
83    xl(i)=xl(i)*pi180
84    yl(i)=yl(i)*pi180
85  end do
86
87
88  ! Generate a seed for each cluster
89  !*********************************
90
91  do j=1,ncluster
92    zclust(j)=0.
93    xclust(j)=xl(j*n/ncluster)
94    yclust(j)=yl(j*n/ncluster)
95  end do
96
97
98  ! Iterative loop to compute the cluster means
99  !********************************************
100
101  do l=1,100
102
103  ! Assign each particle to a cluster: criterion minimum distance to the
104  ! cluster mean position
105  !*********************************************************************
106
107
108    do i=1,n
109      distancemin=10.**10.
110      do j=1,ncluster
111        distances=distance2(yl(i),xl(i),yclust(j),xclust(j))
112        if (distances.lt.distancemin) then
113          distancemin=distances
114          ncl=j
115        endif
116      end do
117      nclust(i)=ncl
118    end do
119
120
121  ! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates,
122  ! calculate mean position, and re-project this point onto the Earth's surface
123  !*****************************************************************************
124
125    do j=1,ncluster
126      xav(j)=0.
127      yav(j)=0.
128      zav(j)=0.
129      rmsclust(j)=0.
130      numb(j)=0
131    end do
132    rms=0.
133
134    do i=1,n
135      numb(nclust(i))=numb(nclust(i))+1
136      distances=distance2(yl(i),xl(i), &
137           yclust(nclust(i)),xclust(nclust(i)))
138
139  ! rms is the total rms of all particles
140  ! rmsclust is the rms for a particular cluster
141  !*********************************************
142
143      rms=rms+distances*distances
144      rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
145
146  ! Calculate Cartesian 3D coordinates from longitude and latitude
147  !***************************************************************
148
149      x = cos(yl(i))*sin(xl(i))
150      y = -1.*cos(yl(i))*cos(xl(i))
151      z = sin(yl(i))
152      xav(nclust(i))=xav(nclust(i))+x
153      yav(nclust(i))=yav(nclust(i))+y
154      zav(nclust(i))=zav(nclust(i))+z
155    end do
156
157    rms=sqrt(rms/real(n))
158
159
160  ! Find the mean location in Cartesian coordinates
161  !************************************************
162
163    do j=1,ncluster
164      if (numb(j).gt.0) then
165        rmsclust(j)=sqrt(rmsclust(j)/real(numb(j)))
166        xav(j)=xav(j)/real(numb(j))
167        yav(j)=yav(j)/real(numb(j))
168        zav(j)=zav(j)/real(numb(j))
169
170  ! Project the point back onto Earth's surface
171  !********************************************
172
173        xclust(j)=atan2(xav(j),-1.*yav(j))
174        yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
175      endif
176    end do
177
178
179  ! Leave the loop if the RMS distance decreases only slightly between 2 iterations
180  !*****************************************************************************
181
182    if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
183    rmsold=rms
184
185  end do
186
18799   continue
188
189  ! Convert longitude and latitude from radians to degrees
190  !*******************************************************
191
192  do i=1,n
193    xl(i)=xl(i)/pi180
194    yl(i)=yl(i)/pi180
195    zclust(nclust(i))=zclust(nclust(i))+zl(i)
196  end do
197
198  do j=1,ncluster
199    xclust(j)=xclust(j)/pi180
200    yclust(j)=yclust(j)/pi180
201    if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j))
202    fclust(j)=100.*real(numb(j))/real(n)
203  end do
204
205  ! Determine total vertical RMS deviation
206  !***************************************
207
208  zrms=0.
209  do i=1,n
210    zdist=zl(i)-zclust(nclust(i))
211    zrms=zrms+zdist*zdist
212  end do
213  if (zrms.gt.0.) zrms=sqrt(zrms/real(n))
214
215  deallocate(nclust)
216
217end subroutine clustering
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG