source: trunk/src/convmix_gfs.f90

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