source: flexpart.git/src/clustering.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 6.5 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
5       rmsclust,zrms)
6  !                      i  i  i  i   o      o      o      o     o
7  !   o      o
8  !*****************************************************************************
9  !                                                                            *
10  !   This routine clusters the particle position into ncluster custers.       *
11  !   Input are the longitudes (xl) and latitudes (yl) of the individual       *
12  !   points, output are the cluster mean positions (xclust,yclust).           *
13  !   Vertical positions are not directly used for the clustering.             *
14  !                                                                            *
15  !   For clustering, the procedure described in Dorling et al. (1992) is used.*
16  !                                                                            *
17  !   Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992):                     *
18  !   Cluster analysis: a technique for estimating the synoptic meteorological *
19  !   controls on air and precipitation chemistry - method and applications.   *
20  !   Atmospheric Environment 26A, 2575-2581.                                  *
21  !                                                                            *
22  !                                                                            *
23  !     Author: A. Stohl                                                       *
24  !                                                                            *
25  !     1 February 2002                                                        *
26  !                                                                            *
27  ! Variables:                                                                 *
28  ! fclust          fraction of particles belonging to each cluster            *
29  ! ncluster        number of clusters to be used                              *
30  ! rms             total horizontal rms distance after clustering             *
31  ! rmsclust        horizontal rms distance for each individual cluster        *
32  ! zrms            total vertical rms distance after clustering               *
33  ! xclust,yclust,  Cluster centroid positions                                 *
34  ! zclust                                                                     *
35  ! xl,yl,zl        particle positions                                         *
36  !                                                                            *
37  !*****************************************************************************
38
39  use par_mod
40
41  implicit none
42
43  integer :: n,i,j,l,nclust(maxpart),numb(ncluster),ncl
44  real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z
45  real :: zclust(ncluster),distance2,distances,distancemin,rms,rmsold
46  real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster)
47  real :: rmsclust(ncluster)
48  real :: zdist,zrms
49
50
51
52  if (n.lt.ncluster) return
53  rmsold=-5.
54
55  ! Convert longitude and latitude from degrees to radians
56  !*******************************************************
57
58  do i=1,n
59    nclust(i)=i
60    xl(i)=xl(i)*pi180
61    yl(i)=yl(i)*pi180
62  end do
63
64
65  ! Generate a seed for each cluster
66  !*********************************
67
68  do j=1,ncluster
69    zclust(j)=0.
70    xclust(j)=xl(j*n/ncluster)
71    yclust(j)=yl(j*n/ncluster)
72  end do
73
74
75  ! Iterative loop to compute the cluster means
76  !********************************************
77
78  do l=1,100
79
80  ! Assign each particle to a cluster: criterion minimum distance to the
81  ! cluster mean position
82  !*********************************************************************
83
84
85    do i=1,n
86      distancemin=10.**10.
87      do j=1,ncluster
88        distances=distance2(yl(i),xl(i),yclust(j),xclust(j))
89        if (distances.lt.distancemin) then
90          distancemin=distances
91          ncl=j
92        endif
93      end do
94      nclust(i)=ncl
95    end do
96
97
98  ! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates,
99  ! calculate mean position, and re-project this point onto the Earth's surface
100  !*****************************************************************************
101
102    do j=1,ncluster
103      xav(j)=0.
104      yav(j)=0.
105      zav(j)=0.
106      rmsclust(j)=0.
107      numb(j)=0
108    end do
109    rms=0.
110
111    do i=1,n
112      numb(nclust(i))=numb(nclust(i))+1
113      distances=distance2(yl(i),xl(i), &
114           yclust(nclust(i)),xclust(nclust(i)))
115
116  ! rms is the total rms of all particles
117  ! rmsclust is the rms for a particular cluster
118  !*********************************************
119
120      rms=rms+distances*distances
121      rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
122
123  ! Calculate Cartesian 3D coordinates from longitude and latitude
124  !***************************************************************
125
126      x = cos(yl(i))*sin(xl(i))
127      y = -1.*cos(yl(i))*cos(xl(i))
128      z = sin(yl(i))
129      xav(nclust(i))=xav(nclust(i))+x
130      yav(nclust(i))=yav(nclust(i))+y
131      zav(nclust(i))=zav(nclust(i))+z
132    end do
133
134    rms=sqrt(rms/real(n))
135
136
137  ! Find the mean location in Cartesian coordinates
138  !************************************************
139
140    do j=1,ncluster
141      if (numb(j).gt.0) then
142        rmsclust(j)=sqrt(rmsclust(j)/real(numb(j)))
143        xav(j)=xav(j)/real(numb(j))
144        yav(j)=yav(j)/real(numb(j))
145        zav(j)=zav(j)/real(numb(j))
146
147  ! Project the point back onto Earth's surface
148  !********************************************
149
150        xclust(j)=atan2(xav(j),-1.*yav(j))
151        yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
152      endif
153    end do
154
155
156  ! Leave the loop if the RMS distance decreases only slightly between 2 iterations
157  !*****************************************************************************
158
159    if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
160    rmsold=rms
161
162  end do
163
16499   continue
165
166  ! Convert longitude and latitude from radians to degrees
167  !*******************************************************
168
169  do i=1,n
170    xl(i)=xl(i)/pi180
171    yl(i)=yl(i)/pi180
172    zclust(nclust(i))=zclust(nclust(i))+zl(i)
173  end do
174
175  do j=1,ncluster
176    xclust(j)=xclust(j)/pi180
177    yclust(j)=yclust(j)/pi180
178    if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j))
179    fclust(j)=100.*real(numb(j))/real(n)
180  end do
181
182  ! Determine total vertical RMS deviation
183  !***************************************
184
185  zrms=0.
186  do i=1,n
187    zdist=zl(i)-zclust(nclust(i))
188    zrms=zrms+zdist*zdist
189  end do
190  if (zrms.gt.0.) zrms=sqrt(zrms/real(n))
191
192end subroutine clustering
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG