source: branches/jerome/src_flexwrf_v3.1/clustering.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1
File size: 10.2 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23      subroutine clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
24      rmsclust,zrms)
25!                           i  i  i  i   o      o      o      o     o
26!        o      o
27!*******************************************************************************
28!                                                                              *
29!     Note:  This is the FLEXPART_WRF version of subroutine clustering.        *
30!            The computational grid is the WRF x-y grid rather than lat-lon.   *
31!                                                                              *
32!   This routine clusters the particle position into ncluster custers.         *
33!   Input are the longitudes (xl) and latitudes (yl) of the individual         *
34!   points, output are the cluster mean positions (xclust,yclust).             *
35!   Vertical positions are not directly used for the clustering.               *
36!                                                                              *
37!   For clustering, the procedure described in Dorling et al. (1992) is used.  *
38!                                                                              *
39!   Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992):                       *
40!   Cluster analysis: a technique for estimating the synoptic meteorological   *
41!   controls on air and precipitation chemistry - method and applications.     *
42!   Atmospheric Environment 26A, 2575-2581.                                    *
43!                                                                              *
44!                                                                              *
45!     Author: A. Stohl                                                         *
46!                                                                              *
47!     1 February 2002                                                          *
48!                                                                              *
49!    26 Oct 2005, R. Easter - changes associated with WRF horizontal grid.     *
50!                 x and y coordinates are in m, so the clustering              *
51!                 calculations are simpler, with no coordinate conversions.    *
52!    10 Mar 2006, R. Easter - bug fix at (new) lines 131-2                     *
53!                 change "yclust(j)" to "yclust(nclust(i))", same for xclust   *
54!                                                                              *
55!*******************************************************************************
56!                                                                              *
57! Variables:                                                                   *
58! fclust          fraction of particles belonging to each cluster              *
59! ncluster        number of clusters to be used                                *
60! rms             total horizontal rms distance after clustering               *
61! rmsclust        horizontal rms distance for each individual cluster          *
62! zrms            total vertical rms distance after clustering                 *
63! xclust,yclust,  Cluster centroid positions                                   *
64! zclust                                                                       *
65! xl,yl,zl        particle positions                                           *
66!                                                                              *
67!                                                                              *
68!*******************************************************************************
69
70
71  use par_mod
72
73  implicit none
74
75  integer :: n,i,j,l,numb(ncluster),ncl,stat
76  real :: xl(n),yl(n),zl(n),xclust(ncluster),yclust(ncluster),x,y,z
77  real :: zclust(ncluster),distance2,distances,distancemin,rms,rmsold
78  real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster)
79  real :: rmsclust(ncluster)
80  real :: zdist,zrms
81  integer,allocatable, dimension (:) :: nclust
82
83    allocate(nclust(maxpart) ,stat=stat)
84
85      if (n.lt.ncluster) return
86      rmsold=-5.
87
88! Convert longitude and latitude from degrees to radians
89!*******************************************************
90
91      do i=1,n
92        nclust(i)=i
93
94! for FLEXPART_WRF, x & y coords are in meters
95!        xl(i)=xl(i)*pi180
96!5       yl(i)=yl(i)*pi180
97  end do
98
99! Generate a seed for each cluster
100!*********************************
101
102      do j=1,ncluster
103        zclust(j)=0.
104        xclust(j)=xl(j*n/ncluster)
105      yclust(j)=yl(j*n/ncluster)
106  end do
107
108
109! Iterative loop to compute the cluster means
110!********************************************
111
112      do l=1,100
113
114! Assign each particle to a cluster: criterion minimum distance to the
115! cluster mean position
116!*********************************************************************
117
118     
119        do i=1,n
120          distancemin=10.**10.
121          do j=1,ncluster
122
123! for FLEXPART_WRF, x & y coords are in meters, so calc distance directly
124!           distances=distance2(yl(i),xl(i),yclust(j),xclust(j))
125            distances=sqrt( (yl(i)-yclust(j))**2 +  &
126                            (xl(i)-xclust(j))**2 )
127
128            if (distances.lt.distancemin) then
129              distancemin=distances
130              ncl=j
131            endif
132      end do
133          nclust(i)=ncl
134      end do
135
136
137! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates,
138! calculate mean position, and re-project this point onto the Earth's surface
139!********************************************************************************
140
141        do j=1,ncluster
142          xav(j)=0.
143          yav(j)=0.
144          zav(j)=0.
145          rmsclust(j)=0.
146          numb(j)=0
147    end do
148
149        rms=0.
150
151        do i=1,n
152          numb(nclust(i))=numb(nclust(i))+1
153
154! for FLEXPART_WRF, x & y coords are in meters, so calc distance directly
155!          distances=distance2(yl(i),xl(i),
156!     +    yclust(nclust(i)),xclust(nclust(i)))
157! 10-mar-2006 rce - bug fix - change "yclust(j)" to
158!    "yclust(nclust(i))", same for xclust
159          distances=sqrt( (yl(i)-yclust(nclust(i)))**2 +  &
160                          (xl(i)-xclust(nclust(i)))**2 )
161
162! rms is the total rms of all particles
163! rmsclust is the rms for a particular cluster
164!*********************************************
165
166          rms=rms+distances*distances
167          rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances
168
169! Calculate Cartesian 3D coordinates from longitude and latitude
170!***************************************************************
171
172! for FLEXPART_WRF, x & y coords are in meters,
173! so no conversion is needed
174!          x = cos(yl(i))*sin(xl(i))
175!          y = -1.*cos(yl(i))*cos(xl(i))
176!          z = sin(yl(i))
177!          xav(nclust(i))=xav(nclust(i))+x
178!          yav(nclust(i))=yav(nclust(i))+y
179!50        zav(nclust(i))=zav(nclust(i))+z
180          xav(nclust(i))=xav(nclust(i))+xl(i)
181          yav(nclust(i))=yav(nclust(i))+yl(i)
182          zav(nclust(i))=0.0
183    end do
184
185
186        rms=sqrt(rms/real(n))
187
188
189! Find the mean location in Cartesian coordinates
190!************************************************
191
192        do j=1,ncluster
193          if (numb(j).gt.0) then
194            rmsclust(j)=sqrt(rmsclust(j)/real(numb(j)))
195            xav(j)=xav(j)/real(numb(j))
196            yav(j)=yav(j)/real(numb(j))
197            zav(j)=zav(j)/real(numb(j))
198
199! Project the point back onto Earth's surface
200!********************************************
201
202! for FLEXPART_WRF, x & y coords are in meters,
203! so no conversion is needed
204!            xclust(j)=atan2(xav(j),-1.*yav(j))
205!            yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j)))
206            xclust(j)=xav(j)
207            yclust(j)=yav(j)
208          endif
209    end do
210
211! Leave the loop if the RMS distance decreases only slightly between 2 iterations
212!********************************************************************************
213
214        if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) goto 99
215        rmsold=rms
216
217  end do
218
21999    continue
220
221! Convert longitude and latitude from radians to degrees
222!*******************************************************
223
224      do i=1,n
225! for FLEXPART_WRF, x & y coords are in meters
226!        xl(i)=xl(i)/pi180
227!        yl(i)=yl(i)/pi180
228       zclust(nclust(i))=zclust(nclust(i))+zl(i)
229  end do
230
231      do j=1,ncluster
232! for FLEXPART_WRF, x & y coords are in meters
233!        xclust(j)=xclust(j)/pi180
234!        yclust(j)=yclust(j)/pi180
235        if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j))
236        fclust(j)=100.*real(numb(j))/real(n)
237  end do
238
239
240! Determine total vertical RMS deviation
241!***************************************
242
243      zrms=0.
244      do i=1,n
245        zdist=zl(i)-zclust(nclust(i))
246       zrms=zrms+zdist*zdist
247  end do
248
249      if (zrms.gt.0.) zrms=sqrt(zrms/real(n))
250
251end subroutine clustering
252
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG