1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine redist (ipart,ktop,ipconv)
!**************************************************************************
! Do the redistribution of particles due to convection
! This subroutine is called for each particle which is assigned
! a new vertical position randomly, based on the convective redistribution
! matrix
!**************************************************************************
! Petra Seibert, Feb 2001, Apr 2001, May 2001, Jan 2002, Nov 2002 and
! Andreas Frank, Nov 2002
! Caroline Forster: November 2004 - February 2005
use par_mod
use com_mod
use conv_mod
implicit none
real,parameter :: const=r_air/ga
integer :: ipart, ktop,ipconv
integer :: k, kz, levnew, levold
real,save :: uvzlev(nuvzmax)
real :: wsub(nuvzmax)
real :: totlevmass, wsubpart
real :: temp_levold,temp_levold1
real :: sub_levold,sub_levold1
real :: pint, pold, rn, tv, tvold, dlevfrac
real :: ew,ran3, ztold,ffraction
real :: tv1, tv2, dlogp, dz, dz1, dz2
integer :: iseed = -88
! ipart ... number of particle to be treated
ipconv=1
! determine height of the eta half-levels (uvzlev)
! do that only once for each grid column
! i.e. when ktop.eq.1
!**************************************************************
if (ktop .le. 1) then
tvold=tt2conv*(1.+0.378*ew(td2conv)/psconv)
pold=psconv
uvzlev(1)=0.
pint = phconv(2)
! determine next virtual temperatures
tv1 = tconv(1)*(1.+0.608*qconv(1))
tv2 = tconv(2)*(1.+0.608*qconv(2))
! interpolate virtual temperature to half-level
tv = tv1 + (tv2-tv1)*(pconv(1)-phconv(2))/(pconv(1)-pconv(2))
if (abs(tv-tvold).gt.0.2) then
uvzlev(2) = uvzlev(1) + &
const*log(pold/pint)* &
(tv-tvold)/log(tv/tvold)
else
uvzlev(2) = uvzlev(1)+ &
const*log(pold/pint)*tv
endif
tvold=tv
tv1=tv2
pold=pint
! integrate profile (calculation of height agl of eta layers) as required
do kz = 3, nconvtop+1
! note that variables defined in calcmatrix.f (pconv,tconv,qconv)
! start at the first real ECMWF model level whereas kz and
! thus uvzlev(kz) starts at the surface. uvzlev is defined at the
! half-levels (between the tconv, qconv etc. values !)
! Thus, uvzlev(kz) is the lower boundary of the tconv(kz) cell.
pint = phconv(kz)
! determine next virtual temperatures
tv2 = tconv(kz)*(1.+0.608*qconv(kz))
! interpolate virtual temperature to half-level
tv = tv1 + (tv2-tv1)*(pconv(kz-1)-phconv(kz))/ &
(pconv(kz-1)-pconv(kz))
if (abs(tv-tvold).gt.0.2) then
uvzlev(kz) = uvzlev(kz-1) + &
const*log(pold/pint)* &
(tv-tvold)/log(tv/tvold)
else
uvzlev(kz) = uvzlev(kz-1)+ &
const*log(pold/pint)*tv
endif
tvold=tv
tv1=tv2
pold=pint
end do
ktop = 2
endif ! (if ktop .le. 1) then
! determine vertical grid position of particle in the eta system
!****************************************************************
ztold = ztra1(abs(ipart))
! find old particle grid position
do kz = 2, nconvtop
if (uvzlev(kz) .ge. ztold ) then
levold = kz-1
goto 30
endif
end do
! Particle is above the potentially convective domain. Skip it.
goto 90
30 continue
! now redistribute particles
!****************************
! Choose a random number and find corresponding level of destination
! Random numbers to be evenly distributed in [0,1]
rn = ran3(iseed)
! initialize levnew
levnew = levold
ffraction = 0.
totlevmass=dpr(levold)/ga
do k = 1,nconvtop
! for backward runs use the transposed matrix
if (ldirect.eq.1) then
ffraction=ffraction+fmassfrac(levold,k) &
/totlevmass
else
ffraction=ffraction+fmassfrac(k,levold) &
/totlevmass
endif
if (rn.le.ffraction) then
levnew=k
! avoid division by zero or a too small number
! if division by zero or a too small number happens the
! particle is assigned to the center of the grid cell
if (ffraction.gt.1.e-20) then
if (ldirect.eq.1) then
dlevfrac = (ffraction-rn) / fmassfrac(levold,k) * totlevmass
else
dlevfrac = (ffraction-rn) / fmassfrac(k,levold) * totlevmass
endif
else
dlevfrac = 0.5
endif
goto 40
endif
end do
40 continue
! now assign new position to particle
if (levnew.le.nconvtop) then
if (levnew.eq.levold) then
ztra1(abs(ipart)) = ztold
else
dlogp = (1.-dlevfrac)* &
(log(phconv(levnew+1))-log(phconv(levnew)))
pint = log(phconv(levnew))+dlogp
dz1 = pint - log(phconv(levnew))
dz2 = log(phconv(levnew+1)) - pint
dz = dz1 + dz2
ztra1(abs(ipart)) = (uvzlev(levnew)*dz2+uvzlev(levnew+1)*dz1)/dz
if (ztra1(abs(ipart)).lt.0.) &
ztra1(abs(ipart))=-1.*ztra1(abs(ipart))
if (ipconv.gt.0) ipconv=-1
endif
endif
! displace particle according to compensating subsidence
! this is done to those particles, that were not redistributed
! by the matrix
!**************************************************************
if (levnew.le.nconvtop.and.levnew.eq.levold) then
ztold = ztra1(abs(ipart))
! determine compensating vertical velocity at the levels
! above and below the particel position
! increase compensating subsidence by the fraction that
! is displaced by convection to this level
if (levold.gt.1) then
temp_levold = tconv(levold-1) + &
(tconv(levold)-tconv(levold-1)) &
*(pconv(levold-1)-phconv(levold))/ &
(pconv(levold-1)-pconv(levold))
sub_levold = sub(levold)/(1.-sub(levold)/dpr(levold)*ga)
wsub(levold)=-1.*sub_levold*r_air*temp_levold/(phconv(levold))
else
wsub(levold)=0.
endif
temp_levold1 = tconv(levold) + &
(tconv(levold+1)-tconv(levold)) &
*(pconv(levold)-phconv(levold+1))/ &
(pconv(levold)-pconv(levold+1))
sub_levold1 = sub(levold+1)/(1.-sub(levold+1)/dpr(levold+1)*ga)
wsub(levold+1)=-1.*sub_levold1*r_air*temp_levold1/ &
(phconv(levold+1))
! interpolate wsub to the vertical particle position
dz1 = ztold - uvzlev(levold)
dz2 = uvzlev(levold+1) - ztold
dz = dz1 + dz2
wsubpart = (dz2*wsub(levold)+dz1*wsub(levold+1))/dz
ztra1(abs(ipart)) = ztold+wsubpart*real(lsynctime)
if (ztra1(abs(ipart)).lt.0.) then
ztra1(abs(ipart))=-1.*ztra1(abs(ipart))
endif
endif !(levnew.le.nconvtop.and.levnew.eq.levold)
! Maximum altitude .5 meter below uppermost model level
!*******************************************************
90 continue
if (ztra1(abs(ipart)) .gt. height(nz)-0.5) &
ztra1(abs(ipart)) = height(nz)-0.5
end subroutine redist
|
