source: flexpart.git/src/convmix.f90 @ 027e844

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 027e844 was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Merged changes from CTBTO project

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