source: trunk/src/convmix.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 9 years ago
File size: 9.9 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 convmix(itime)
23  !                     i
24  !**************************************************************
25  !handles all the calculations related to convective mixing
26  !Petra Seibert, Bernd C. Krueger, Feb 2001
27  !nested grids included, Bernd C. Krueger, May 2001
28  !
29  !Changes by Caroline Forster, April 2004 - February 2005:
30  !  convmix called every lsynctime seconds
31  !CHANGES by A. Stohl:
32  !  various run-time optimizations - February 2005
33  !**************************************************************
34
35  use flux_mod
36  use par_mod
37  use com_mod
38  use conv_mod
39
40  implicit none
41
42  integer :: igr,igrold, ipart, itime, ix, j, inest
43  integer :: ipconv
44  integer :: jy, kpart, ktop, ngrid,kz
45  integer :: igrid(maxpart), ipoint(maxpart), igridn(maxpart,maxnests)
46  ! itime [s]                 current time
47  ! igrid(maxpart)            horizontal grid position of each particle
48  ! igridn(maxpart,maxnests)  dto. for nested grids
49  ! ipoint(maxpart)           pointer to access particles according to grid position
50
51  logical :: lconv
52  real :: x, y, xtn,ytn, ztold, delt
53  real :: dt1,dt2,dtt
54  integer :: mind1,mind2
55  ! dt1,dt2,dtt,mind1,mind2       variables used for time interpolation
56  integer :: itage,nage
57  real,parameter :: eps=nxmax/3.e5
58
59
60  !monitoring variables
61  !real sumconv,sumall
62
63
64  ! Calculate auxiliary variables for time interpolation
65  !*****************************************************
66
67  dt1=real(itime-memtime(1))
68  dt2=real(memtime(2)-itime)
69  dtt=1./(dt1+dt2)
70  mind1=memind(1)
71  mind2=memind(2)
72  delt=real(abs(lsynctime))
73
74
75  lconv = .false.
76
77  ! if no particles are present return after initialization
78  !********************************************************
79
80  if (numpart.le.0) return
81
82  ! Assign igrid and igridn, which are pseudo grid numbers indicating particles
83  ! that are outside the part of the grid under consideration
84  ! (e.g. particles near the poles or particles in other nests).
85  ! Do this for all nests but use only the innermost nest; for all others
86  ! igrid shall be -1
87  ! Also, initialize index vector ipoint
88  !************************************************************************
89
90  do ipart=1,numpart
91    igrid(ipart)=-1
92    do j=numbnests,1,-1
93      igridn(ipart,j)=-1
94    end do
95    ipoint(ipart)=ipart
96  ! do not consider particles that are (yet) not part of simulation
97    if (itra1(ipart).ne.itime) goto 20
98    x = xtra1(ipart)
99    y = ytra1(ipart)
100
101  ! Determine which nesting level to be used
102  !**********************************************************
103
104    ngrid=0
105    do j=numbnests,1,-1
106      if ( x.gt.xln(j)+eps .and. x.lt.xrn(j)-eps .and. &
107           y.gt.yln(j)+eps .and. y.lt.yrn(j)-eps ) then
108        ngrid=j
109        goto 23
110      endif
111    end do
112 23   continue
113
114  ! Determine nested grid coordinates
115  !**********************************
116
117    if (ngrid.gt.0) then
118  ! nested grids
119      xtn=(x-xln(ngrid))*xresoln(ngrid)
120      ytn=(y-yln(ngrid))*yresoln(ngrid)
121      ix=nint(xtn)
122      jy=nint(ytn)
123      igridn(ipart,ngrid) = 1 + jy*nxn(ngrid) + ix
124    else if(ngrid.eq.0) then
125  ! mother grid
126      ix=nint(x)
127      jy=nint(y)
128      igrid(ipart) = 1 + jy*nx + ix
129    endif
130
131 20 continue
132  end do
133
134  !sumall = 0.
135  !sumconv = 0.
136
137  !*****************************************************************************
138  ! 1. Now, do everything for the mother domain and, later, for all of the nested domains
139  ! While all particles have to be considered for redistribution, the Emanuel convection
140  ! scheme only needs to be called once for every grid column where particles are present.
141  ! Therefore, particles are sorted according to their grid position. Whenever a new grid
142  ! cell is encountered by looping through the sorted particles, the convection scheme is called.
143  !*****************************************************************************
144
145  ! sort particles according to horizontal position and calculate index vector IPOINT
146
147  call sort2(numpart,igrid,ipoint)
148
149  ! Now visit all grid columns where particles are present
150  ! by going through the sorted particles
151
152  igrold = -1
153  do kpart=1,numpart
154    igr = igrid(kpart)
155    if (igr .eq. -1) goto 50
156    ipart = ipoint(kpart)
157  !  sumall = sumall + 1
158    if (igr .ne. igrold) then
159  ! we are in a new grid column
160      jy = (igr-1)/nx
161      ix = igr - jy*nx - 1
162
163  ! Interpolate all meteorological data needed for the convection scheme
164      psconv=(ps(ix,jy,1,mind1)*dt2+ps(ix,jy,1,mind2)*dt1)*dtt
165      tt2conv=(tt2(ix,jy,1,mind1)*dt2+tt2(ix,jy,1,mind2)*dt1)*dtt
166      td2conv=(td2(ix,jy,1,mind1)*dt2+td2(ix,jy,1,mind2)*dt1)*dtt
167!!$      do kz=1,nconvlev+1      !old
168        do kz=1,nuvz-1           !bugfix
169        tconv(kz)=(tth(ix,jy,kz+1,mind1)*dt2+ &
170             tth(ix,jy,kz+1,mind2)*dt1)*dtt
171        qconv(kz)=(qvh(ix,jy,kz+1,mind1)*dt2+ &
172             qvh(ix,jy,kz+1,mind2)*dt1)*dtt
173      end do
174
175  ! Calculate translocation matrix
176      call calcmatrix(lconv,delt,cbaseflux(ix,jy))
177      igrold = igr
178      ktop = 0
179    endif
180
181  ! treat particle only if column has convection
182    if (lconv .eqv. .true.) then
183  ! assign new vertical position to particle
184
185      ztold=ztra1(ipart)
186      call redist(ipart,ktop,ipconv)
187  !    if (ipconv.le.0) sumconv = sumconv+1
188
189  ! Calculate the gross fluxes across layer interfaces
190  !***************************************************
191
192      if (iflux.eq.1) then
193        itage=abs(itra1(ipart)-itramem(ipart))
194        do nage=1,nageclass
195          if (itage.lt.lage(nage)) goto 37
196        end do
197 37     continue
198
199        if (nage.le.nageclass) &
200             call calcfluxes(nage,ipart,real(xtra1(ipart)), &
201             real(ytra1(ipart)),ztold)
202      endif
203
204    endif   !(lconv .eqv. .true)
205 50 continue
206  end do
207
208
209  !*****************************************************************************
210  ! 2. Nested domains
211  !*****************************************************************************
212
213  ! sort particles according to horizontal position and calculate index vector IPOINT
214
215  do inest=1,numbnests
216    do ipart=1,numpart
217      ipoint(ipart)=ipart
218      igrid(ipart) = igridn(ipart,inest)
219    enddo
220    call sort2(numpart,igrid,ipoint)
221
222  ! Now visit all grid columns where particles are present
223  ! by going through the sorted particles
224
225    igrold = -1
226    do kpart=1,numpart
227      igr = igrid(kpart)
228      if (igr .eq. -1) goto 60
229      ipart = ipoint(kpart)
230  !    sumall = sumall + 1
231      if (igr .ne. igrold) then
232  ! we are in a new grid column
233        jy = (igr-1)/nxn(inest)
234        ix = igr - jy*nxn(inest) - 1
235
236  ! Interpolate all meteorological data needed for the convection scheme
237        psconv=(psn(ix,jy,1,mind1,inest)*dt2+ &
238             psn(ix,jy,1,mind2,inest)*dt1)*dtt
239        tt2conv=(tt2n(ix,jy,1,mind1,inest)*dt2+ &
240             tt2n(ix,jy,1,mind2,inest)*dt1)*dtt
241        td2conv=(td2n(ix,jy,1,mind1,inest)*dt2+ &
242             td2n(ix,jy,1,mind2,inest)*dt1)*dtt
243!!$        do kz=1,nconvlev+1    !old
244        do kz=1,nuvz-1           !bugfix
245          tconv(kz)=(tthn(ix,jy,kz+1,mind1,inest)*dt2+ &
246               tthn(ix,jy,kz+1,mind2,inest)*dt1)*dtt
247          qconv(kz)=(qvhn(ix,jy,kz+1,mind1,inest)*dt2+ &
248               qvhn(ix,jy,kz+1,mind2,inest)*dt1)*dtt
249        end do
250
251  ! calculate translocation matrix
252  !*******************************
253        call calcmatrix(lconv,delt,cbasefluxn(ix,jy,inest))
254        igrold = igr
255        ktop = 0
256      endif
257
258  ! treat particle only if column has convection
259      if (lconv .eqv. .true.) then
260  ! assign new vertical position to particle
261        ztold=ztra1(ipart)
262        call redist(ipart,ktop,ipconv)
263  !      if (ipconv.le.0) sumconv = sumconv+1
264
265  ! Calculate the gross fluxes across layer interfaces
266  !***************************************************
267
268        if (iflux.eq.1) then
269          itage=abs(itra1(ipart)-itramem(ipart))
270          do nage=1,nageclass
271            if (itage.lt.lage(nage)) goto 47
272          end do
273 47       continue
274
275          if (nage.le.nageclass) &
276               call calcfluxes(nage,ipart,real(xtra1(ipart)), &
277               real(ytra1(ipart)),ztold)
278        endif
279
280      endif !(lconv .eqv. .true.)
281
282
28360    continue
284    end do
285  end do
286  !--------------------------------------------------------------------------
287  !write(*,*)'############################################'
288  !write(*,*)'TIME=',
289  !    &  itime
290  !write(*,*)'fraction of particles under convection',
291  !    &  sumconv/(sumall+0.001)
292  !write(*,*)'total number of particles',
293  !    &  sumall
294  !write(*,*)'number of particles under convection',
295  !    &  sumconv
296  !write(*,*)'############################################'
297
298  return
299end subroutine convmix
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG