source: flexpart.git/src/redist.f90 @ 4c64400

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 4c64400 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

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