source: flexpart.git/src/redist_mpi.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago

Updated wet depo scheme

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