source: branches/flexpart91_hasod/src_parallel/init_domainfill.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago

Added parallel version of Flexpart91

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