source: trunk/src/redist.f90 @ 28

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