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