source: flexpart.git/src/convmix.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

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