Ticket #277: init_domainfill.2.f90

File init_domainfill.2.f90, 14.5 KB (added by nehstand, 4 years ago)
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine init_domainfill
5!
6!*****************************************************************************
7!                                                                            *
8! Initializes particles equally distributed over the first release location  *
9! specified in file RELEASES. This box is assumed to be the domain for doing *
10! domain-filling trajectory calculations.                                    *
11! All particles carry the same amount of mass which alltogether comprises the*
12! mass of air within the box.                                                *
13!                                                                            *
14!     Author: A. Stohl                                                       *
15!                                                                            *
16!     15 October 2002                                                        *
17!                                                                            *
18!*****************************************************************************
19!                                                                            *
20! Variables:                                                                 *
21!                                                                            *
22! numparticlecount    consecutively counts the number of particles released  *
23! nx_we(2)       grid indices for western and eastern boundary of domain-    *
24!                filling trajectory calculations                             *
25! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
26!                filling trajectory calculations                             *
27!                                                                            *
28!*****************************************************************************
29
30  use point_mod
31  use par_mod
32  use com_mod
33  use random_mod
34
35  implicit none
36
37  integer :: j,ix,jy,kz,ncolumn,numparttot
38  real :: gridarea(0:nymax-1),pp(nzmax),ylat,ylatp,ylatm,hzone
39  real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,fractus
40  real,parameter :: pih=pi/180.
41  real :: colmass(0:nxmax-1,0:nymax-1),colmasstotal,zposition
42
43  integer :: ixm,ixp,jym,jyp,indzm,indzp,in,indzh,i,jj
44  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2)
45
46  integer :: idummy = -11
47
48
49! Determine the release region (only full grid cells), over which particles
50! shall be initialized
51! Use 2 fields for west/east and south/north boundary
52!**************************************************************************
53
54  nx_we(1)=max(int(xpoint1(1)),0)
55  nx_we(2)=min((int(xpoint2(1))+1),nxmin1)
56  ny_sn(1)=max(int(ypoint1(1)),0)
57  ny_sn(2)=min((int(ypoint2(1))+1),nymin1)
58
59! For global simulations (both global wind data and global domain-filling),
60! set a switch, such that no boundary conditions are used
61!**************************************************************************
62  if (xglobal.and.sglobal.and.nglobal) then
63    if ((nx_we(1).eq.0).and.(nx_we(2).eq.nxmin1).and. &
64         (ny_sn(1).eq.0).and.(ny_sn(2).eq.nymin1)) then
65      gdomainfill=.true.
66    else
67      gdomainfill=.false.
68    endif
69  endif
70
71! Exit here if resuming a run from particle dump
72!***********************************************
73  if (gdomainfill.and.ipin.ne.0) return
74
75! Do not release particles twice (i.e., not at both in the leftmost and rightmost
76! grid cell) for a global domain
77!*****************************************************************************
78  if (xglobal) nx_we(2)=min(nx_we(2),nx-2)
79
80
81! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
82! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
83!************************************************************
84
85  do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
86    ylat=ylat0+real(jy)*dy
87    ylatp=ylat+0.5*dy
88    ylatm=ylat-0.5*dy
89    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
90      hzone=1./dyconst
91    else
92      cosfactp=cos(ylatp*pih)*r_earth
93      cosfactm=cos(ylatm*pih)*r_earth
94      if (cosfactp.lt.cosfactm) then
95        hzone=sqrt(r_earth**2-cosfactp**2)- &
96             sqrt(r_earth**2-cosfactm**2)
97      else
98        hzone=sqrt(r_earth**2-cosfactm**2)- &
99             sqrt(r_earth**2-cosfactp**2)
100      endif
101    endif
102    gridarea(jy)=2.*pi*r_earth*hzone*dx/360.
103  end do
104
105! Do the same for the south pole
106
107  if (sglobal) then
108    ylat=ylat0
109    ylatp=ylat+0.5*dy
110    ylatm=ylat
111    cosfactm=0.
112    cosfactp=cos(ylatp*pih)*r_earth
113    hzone=sqrt(r_earth**2-cosfactm**2)- &
114         sqrt(r_earth**2-cosfactp**2)
115    gridarea(0)=2.*pi*r_earth*hzone*dx/360.
116  endif
117
118! Do the same for the north pole
119
120  if (nglobal) then
121    ylat=ylat0+real(nymin1)*dy
122    ylatp=ylat
123    ylatm=ylat-0.5*dy
124    cosfactp=0.
125    cosfactm=cos(ylatm*pih)*r_earth
126    hzone=sqrt(r_earth**2-cosfactp**2)- &
127         sqrt(r_earth**2-cosfactm**2)
128    gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360.
129  endif
130
131
132! Calculate total mass of each grid column and of the whole atmosphere
133!*********************************************************************
134
135  colmasstotal=0.
136  do jy=ny_sn(1),ny_sn(2)          ! loop about latitudes
137    do ix=nx_we(1),nx_we(2)      ! loop about longitudes
138      pp(1)=rho(ix,jy,1,1)*r_air*tt(ix,jy,1,1)
139      pp(nz)=rho(ix,jy,nz,1)*r_air*tt(ix,jy,nz,1)
140      colmass(ix,jy)=(pp(1)-pp(nz))/ga*gridarea(jy)
141      colmasstotal=colmasstotal+colmass(ix,jy)
142    end do
143  end do
144
145  write(*,*) 'Atm. mass: ',colmasstotal
146
147
148  if (ipin.eq.0) numpart=0
149
150! Determine the particle positions
151!*********************************
152
153  numparttot=0
154  numcolumn=0
155  do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
156    ylat=ylat0+real(jy)*dy
157    do ix=nx_we(1),nx_we(2)      ! loop about longitudes
158      ncolumn=nint(0.999*real(npart(1))*colmass(ix,jy)/ &
159           colmasstotal)
160      if (ncolumn.eq.0) goto 30
161      if (ncolumn.gt.numcolumn) numcolumn=ncolumn
162
163! Calculate pressure at the altitudes of model surfaces, using the air density
164! information, which is stored as a 3-d field
165!*****************************************************************************
166
167      do kz=1,nz
168        pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
169      end do
170
171
172      deltacol=(pp(1)-pp(nz))/real(ncolumn)
173      pnew=pp(1)+deltacol/2.
174      jj=0
175      do j=1,ncolumn
176        jj=jj+1
177
178
179! For columns with many particles (i.e. around the equator), distribute
180! the particles equally, for columns with few particles (i.e. around the
181! poles), distribute the particles randomly
182!***********************************************************************
183
184
185        if (ncolumn.gt.20) then
186          pnew=pnew-deltacol
187        else
188          pnew=pp(1)-ran1(idummy)*(pp(1)-pp(nz))
189        endif
190
191        do kz=1,nz-1
192          if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
193            dz1=pp(kz)-pnew
194            dz2=pnew-pp(kz+1)
195            dz=1./(dz1+dz2)
196
197! Assign particle position
198!*************************
199! Do the following steps only if particles are not read in from previous model run
200!*****************************************************************************
201            if (ipin.eq.0) then
202              xtra1(numpart+jj)=real(ix)-0.5+ran1(idummy)
203              if (ix.eq.0) xtra1(numpart+jj)=ran1(idummy)
204              if (ix.eq.nxmin1) xtra1(numpart+jj)= &
205                   real(nxmin1)-ran1(idummy)
206              ytra1(numpart+jj)=real(jy)-0.5+ran1(idummy)
207              ztra1(numpart+jj)=(height(kz)*dz2+height(kz+1)*dz1)*dz
208              if (ztra1(numpart+jj).gt.height(nz)-0.5) &
209                   ztra1(numpart+jj)=height(nz)-0.5
210
211
212! Interpolate PV to the particle position
213!****************************************
214              ixm=int(xtra1(numpart+jj))
215              jym=int(ytra1(numpart+jj))
216              ixp=ixm+1
217              jyp=jym+1
218              ddx=xtra1(numpart+jj)-real(ixm)
219              ddy=ytra1(numpart+jj)-real(jym)
220              rddx=1.-ddx
221              rddy=1.-ddy
222              p1=rddx*rddy
223              p2=ddx*rddy
224              p3=rddx*ddy
225              p4=ddx*ddy
226              do i=2,nz
227                if (height(i).gt.ztra1(numpart+jj)) then
228                  indzm=i-1
229                  indzp=i
230                  goto 6
231                endif
232              end do
2336             continue
234              dz1=ztra1(numpart+jj)-height(indzm)
235              dz2=height(indzp)-ztra1(numpart+jj)
236              dz=1./(dz1+dz2)
237              do in=1,2
238                indzh=indzm+in-1
239                y1(in)=p1*pv(ixm,jym,indzh,1) &
240                     +p2*pv(ixp,jym,indzh,1) &
241                     +p3*pv(ixm,jyp,indzh,1) &
242                     +p4*pv(ixp,jyp,indzh,1)
243              end do
244              pvpart=(dz2*y1(1)+dz1*y1(2))*dz
245              if (ylat.lt.0.) pvpart=-1.*pvpart
246
247
248! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
249!*****************************************************************************
250
251              if (((ztra1(numpart+jj).gt.3000.).and. &
252                   (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
253
254! Assign certain properties to the particle
255!******************************************
256                nclass(numpart+jj)=min(int(ran1(idummy)* &
257                     real(nclassunc))+1,nclassunc)
258                numparticlecount=numparticlecount+1
259                npoint(numpart+jj)=numparticlecount
260                idt(numpart+jj)=mintime
261                itra1(numpart+jj)=0
262                itramem(numpart+jj)=0
263                itrasplit(numpart+jj)=itra1(numpart+jj)+ldirect* &
264                     itsplit
265                xmass1(numpart+jj,1)=colmass(ix,jy)/real(ncolumn)
266                if (mdomainfill.eq.2) xmass1(numpart+jj,1)= &
267                     xmass1(numpart+jj,1)*pvpart*48./29.*ozonescale/10.**9
268              else
269                jj=jj-1
270              endif
271            endif
272          endif
273        end do
274      end do
275      numparttot=numparttot+ncolumn
276      if (ipin.eq.0) numpart=numpart+jj
27730    continue
278    end do
279  end do
280
281
282! Check whether numpart is really smaller than maxpart
283!*****************************************************
284
285! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier
286  if (numpart.gt.maxpart) then
287    write(*,*) 'numpart too large: change source in init_atm_mass.f'
288    write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart
289  endif
290
291
292  xmassperparticle=colmasstotal/real(numparttot)
293
294
295! Make sure that all particles are within domain
296!***********************************************
297
298  do j=1,numpart
299    if ((xtra1(j).lt.0.).or.(xtra1(j).ge.real(nxmin1)).or. &
300         (ytra1(j).lt.0.).or.(ytra1(j).ge.real(nymin1))) then
301      itra1(j)=-999999999
302    endif
303  end do
304
305
306
307
308! For boundary conditions, we need fewer particle release heights per column,
309! because otherwise it takes too long until enough mass has accumulated to
310! release a particle at the boundary (would take dx/u seconds), leading to
311! relatively large position errors of the order of one grid distance.
312! It's better to release fewer particles per column, but to do so more often.
313! Thus, use on the order of nz starting heights per column.
314! We thus repeat the above to determine fewer starting heights, that are
315! used furtheron in subroutine boundcond_domainfill.f.
316!****************************************************************************
317
318  fractus=real(numcolumn)/real(nz)
319  write(*,*) 'Total number of particles at model start: ',numpart
320  write(*,*) 'Maximum number of particles per column: ',numcolumn
321  write(*,*) 'If ',fractus,' <1, better use more particles'
322  fractus=sqrt(max(fractus,1.))/2.
323
324  do jy=ny_sn(1),ny_sn(2)      ! loop about latitudes
325    do ix=nx_we(1),nx_we(2)      ! loop about longitudes
326      ncolumn=nint(0.999/fractus*real(npart(1))*colmass(ix,jy) &
327           /colmasstotal)
328      if (ncolumn.gt.maxcolumn) stop 'maxcolumn too small'
329      if (ncolumn.eq.0) goto 80
330
331
332! Memorize how many particles per column shall be used for all boundaries
333! This is further used in subroutine boundcond_domainfill.f
334! Use 2 fields for west/east and south/north boundary
335!************************************************************************
336
337      if (ix.eq.nx_we(1)) numcolumn_we(1,jy)=ncolumn
338      if (ix.eq.nx_we(2)) numcolumn_we(2,jy)=ncolumn
339      if (jy.eq.ny_sn(1)) numcolumn_sn(1,ix)=ncolumn
340      if (jy.eq.ny_sn(2)) numcolumn_sn(2,ix)=ncolumn
341
342! Calculate pressure at the altitudes of model surfaces, using the air density
343! information, which is stored as a 3-d field
344!*****************************************************************************
345
346      do kz=1,nz
347        pp(kz)=rho(ix,jy,kz,1)*r_air*tt(ix,jy,kz,1)
348      end do
349
350! Determine the reference starting altitudes
351!*******************************************
352
353      deltacol=(pp(1)-pp(nz))/real(ncolumn)
354      pnew=pp(1)+deltacol/2.
355      do j=1,ncolumn
356        pnew=pnew-deltacol
357        do kz=1,nz-1
358          if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
359            dz1=pp(kz)-pnew
360            dz2=pnew-pp(kz+1)
361            dz=1./(dz1+dz2)
362            zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
363            if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5
364
365! Memorize vertical positions where particles are introduced
366! This is further used in subroutine boundcond_domainfill.f
367!***********************************************************
368
369            if (ix.eq.nx_we(1)) zcolumn_we(1,jy,j)=zposition
370            if (ix.eq.nx_we(2)) zcolumn_we(2,jy,j)=zposition
371            if (jy.eq.ny_sn(1)) zcolumn_sn(1,ix,j)=zposition
372            if (jy.eq.ny_sn(2)) zcolumn_sn(2,ix,j)=zposition
373
374! Initialize mass that has accumulated at boundary to zero
375!*********************************************************
376
377            acc_mass_we(1,jy,j)=0.
378            acc_mass_we(2,jy,j)=0.
379            acc_mass_sn(1,jy,j)=0.
380            acc_mass_sn(2,jy,j)=0.
381          endif
382        end do
383      end do
38480    continue
385    end do
386  end do
387
388! Reduce numpart if invalid particles at end of arrays
389  do i=numpart, 1, -1
390      if (itra1(i).eq.-999999999) then
391        numpart=numpart-1
392      else
393        exit
394      end if
395    end do
396
397! If particles shall be read in to continue an existing run,
398! then the accumulated masses at the domain boundaries must be read in, too.
399! This overrides any previous calculations.
400!***************************************************************************
401
402  if ((ipin.eq.1).and.(.not.gdomainfill)) then
403    open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
404         form='unformatted')
405    read(unitboundcond) numcolumn_we,numcolumn_sn, &
406         zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
407    close(unitboundcond)
408  endif
409
410
411end subroutine init_domainfill
hosted by ZAMG