source: trunk/src/plumetraj.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 10 years ago
File size: 8.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 plumetraj(itime)
23  !                       i
24  !*****************************************************************************
25  !                                                                            *
26  ! Determines a plume centroid trajectory for each release site, and manages  *
27  ! clustering of particle locations. Certain parameters (average PV,          *
28  ! tropopause height, etc., are provided along the plume trajectories.        *
29  ! At the end, output is written to file 'trajectories.txt'.                  *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     24 January 2002                                                        *
34  !                                                                            *
35  ! Variables:                                                                 *
36  ! fclust          fraction of particles belonging to each cluster            *
37  ! hmixcenter      mean mixing height for all particles                       *
38  ! ncluster        number of clusters to be used                              *
39  ! pvcenter        mean PV for all particles                                  *
40  ! pvfract         fraction of particles with PV<2pvu                         *
41  ! rms             total horizontal rms distance after clustering             *
42  ! rmsdist         total horizontal rms distance before clustering            *
43  ! rmsclust        horizontal rms distance for each individual cluster        *
44  ! topocenter      mean topography underlying all particles                   *
45  ! tropocenter     mean tropopause height at the positions of particles       *
46  ! tropofract      fraction of particles within the troposphere               *
47  ! zrms            total vertical rms distance after clustering               *
48  ! zrmsdist        total vertical rms distance before clustering              *
49  ! xclust,yclust,  Cluster centroid positions                                 *
50  ! zclust                                                                     *
51  !                                                                            *
52  !*****************************************************************************
53
54  use point_mod
55  use par_mod
56  use com_mod
57
58  implicit none
59
60  integer :: itime,ix,jy,ixp,jyp,indexh,i,j,k,m,n,il,ind,indz,indzp
61  real :: xl(maxpart),yl(maxpart),zl(maxpart)
62  real :: xcenter,ycenter,zcenter,dist,distance,rmsdist,zrmsdist
63
64  real :: xclust(ncluster),yclust(ncluster),zclust(ncluster)
65  real :: fclust(ncluster),rms,rmsclust(ncluster),zrms
66
67  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
68  real :: topo,topocenter,hm(2),hmixi,hmixfract,hmixcenter
69  real :: pv1(2),pvprof(2),pvi,pvcenter,pvfract,tr(2),tri,tropofract
70  real :: tropocenter
71
72
73  dt1=real(itime-memtime(1))
74  dt2=real(memtime(2)-itime)
75  dtt=1./(dt1+dt2)
76
77
78  ! Loop about all release points
79  !******************************
80
81  do j=1,numpoint
82    if (abs(ireleasestart(j)-itime).gt.lage(nageclass)) goto 10
83    topocenter=0.
84    hmixcenter=0.
85    hmixfract=0.
86    tropocenter=0.
87    tropofract=0.
88    pvfract=0.
89    pvcenter=0.
90    rmsdist=0.
91    zrmsdist=0.
92
93    n=0
94    do i=1,numpart
95      if (itra1(i).ne.itime) goto 20
96      if (npoint(i).ne.j) goto 20
97      n=n+1
98      xl(n)=xlon0+xtra1(i)*dx
99      yl(n)=ylat0+ytra1(i)*dy
100      zl(n)=ztra1(i)
101
102
103  ! Interpolate PBL height, PV, and tropopause height to each
104  ! particle position in order to determine fraction of particles
105  ! within the PBL, above tropopause height, and average PV.
106  ! Interpolate topography, too, and convert to altitude asl
107  !**************************************************************
108
109      ix=int(xtra1(i))
110      jy=int(ytra1(i))
111      ixp=ix+1
112      jyp=jy+1
113      ddx=xtra1(i)-real(ix)
114      ddy=ytra1(i)-real(jy)
115      rddx=1.-ddx
116      rddy=1.-ddy
117      p1=rddx*rddy
118      p2=ddx*rddy
119      p3=rddx*ddy
120      p4=ddx*ddy
121
122  ! Topography
123  !***********
124
125      topo=p1*oro(ix ,jy) &
126           + p2*oro(ixp,jy) &
127           + p3*oro(ix ,jyp) &
128           + p4*oro(ixp,jyp)
129      topocenter=topocenter+topo
130
131  ! Potential vorticity
132  !********************
133
134      do il=2,nz
135        if (height(il).gt.zl(n)) then
136          indz=il-1
137          indzp=il
138          goto 6
139        endif
140      end do
1416     continue
142
143      dz1=zl(n)-height(indz)
144      dz2=height(indzp)-zl(n)
145      dz=1./(dz1+dz2)
146
147
148      do ind=indz,indzp
149        do m=1,2
150          indexh=memind(m)
151          pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
152               +p2*pv(ixp,jy ,ind,indexh) &
153               +p3*pv(ix ,jyp,ind,indexh) &
154               +p4*pv(ixp,jyp,ind,indexh)
155        end do
156        pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
157      end do
158      pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
159      pvcenter=pvcenter+pvi
160      if (yl(n).gt.0.) then
161        if (pvi.lt.2.) pvfract=pvfract+1.
162      else
163        if (pvi.gt.-2.) pvfract=pvfract+1.
164      endif
165
166
167  ! Tropopause and PBL height
168  !**************************
169
170      do m=1,2
171        indexh=memind(m)
172
173        tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
174             + p2*tropopause(ixp,jy ,1,indexh) &
175             + p3*tropopause(ix ,jyp,1,indexh) &
176             + p4*tropopause(ixp,jyp,1,indexh)
177
178        hm(m)=p1*hmix(ix ,jy ,1,indexh) &
179             + p2*hmix(ixp,jy ,1,indexh) &
180             + p3*hmix(ix ,jyp,1,indexh) &
181             + p4*hmix(ixp,jyp,1,indexh)
182      end do
183
184      hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
185      tri=(tr(1)*dt2+tr(2)*dt1)*dtt
186      if (zl(n).lt.tri) tropofract=tropofract+1.
187      tropocenter=tropocenter+tri+topo
188      if (zl(n).lt.hmixi) hmixfract=hmixfract+1.
189      zl(n)=zl(n)+topo        ! convert to height asl
190      hmixcenter=hmixcenter+hmixi
191
192
19320    continue
194    end do
195
196
197  ! Make statistics for all plumes with n>0 particles
198  !**************************************************
199
200    if (n.gt.0) then
201      topocenter=topocenter/real(n)
202      hmixcenter=hmixcenter/real(n)
203      pvcenter=pvcenter/real(n)
204      tropocenter=tropocenter/real(n)
205      hmixfract=100.*hmixfract/real(n)
206      pvfract=100.*pvfract/real(n)
207      tropofract=100.*tropofract/real(n)
208
209  ! Cluster the particle positions
210  !*******************************
211
212      call clustering(xl,yl,zl,n,xclust,yclust,zclust,fclust,rms, &
213           rmsclust,zrms)
214
215
216  ! Determine center of mass position on earth and average height
217  !**************************************************************
218
219      call centerofmass(xl,yl,n,xcenter,ycenter)
220      call mean(zl,zcenter,zrmsdist,n)
221
222  ! Root mean square distance from center of mass
223  !**********************************************
224
225      do k=1,n
226        dist=distance(yl(k),xl(k),ycenter,xcenter)
227        rmsdist=rmsdist+dist*dist
228      end do
229      if (rmsdist.gt.0.) rmsdist=sqrt(rmsdist/real(n))
230      rmsdist=max(rmsdist,0.)
231
232  ! Write out results in trajectory data file
233  !******************************************
234
235      write(unitouttraj,'(i5,i8,2f9.4,4f8.1,f8.2,4f8.1,3f6.1,&
236           &5(2f8.3,f7.0,f6.1,f8.1))')&
237           &j,itime-(ireleasestart(j)+ireleaseend(j))/2, &
238           xcenter,ycenter,zcenter,topocenter,hmixcenter,tropocenter, &
239           pvcenter,rmsdist,rms,zrmsdist,zrms,hmixfract,pvfract, &
240           tropofract, &
241           (xclust(k),yclust(k),zclust(k),fclust(k),rmsclust(k), &
242           k=1,ncluster)
243    endif
244
245
24610  continue
247  end do
248
249
250end subroutine plumetraj
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG