Changeset 861805a in flexpart.git for src/releaseparticles_mpi.f90
- Timestamp:
- Sep 6, 2016, 9:59:11 AM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 16b61a5, 93786a1
- Parents:
- 0f7835d
- File:
-
- 1 edited
-
src/releaseparticles_mpi.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/releaseparticles_mpi.f90
r7e52e2e r861805a 21 21 22 22 subroutine releaseparticles(itime) 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 ! addpartextra particle assigned to MPI process if *49 50 23 ! o 24 !***************************************************************************** 25 ! * 26 ! This subroutine releases particles from the release locations. * 27 ! * 28 ! It searches for a "vacant" storage space and assigns all particle * 29 ! information to that space. A space is vacant either when no particle * 30 ! is yet assigned to it, or when it's particle is expired and, thus, * 31 ! the storage space is made available to a new particle. * 32 ! * 33 ! Author: A. Stohl * 34 ! * 35 ! 29 June 2002 * 36 ! * 37 ! CHANGES * 38 ! 12/2014 eso: MPI version * 39 ! * 40 !***************************************************************************** 41 ! * 42 ! Variables: * 43 ! itime [s] current time * 44 ! ireleasestart, ireleaseend start and end times of all releases * 45 ! npart(maxpoint) number of particles to be released in total * 46 ! numrel number of particles to be released during this time * 47 ! step * 48 ! addone extra particle assigned to MPI process if * 49 ! mod(npart(i),mp_partgroup_np) .ne. 0) * 50 !***************************************************************************** 51 51 52 52 use point_mod … … 59 59 implicit none 60 60 61 61 !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint) 62 62 real :: xaux,yaux,zaux,rfraction 63 63 real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4 … … 73 73 74 74 integer :: idummy = -7 75 76 75 !save idummy,xmasssave 76 !data idummy/-7/,xmasssave/maxpoint*0./ 77 77 78 78 logical :: first_call=.true. 79 79 80 81 80 ! Use different seed for each process. 81 !**************************************************************************** 82 82 if (first_call) then 83 83 idummy=idummy+mp_seed … … 87 87 mind2=memind(2) 88 88 89 90 89 ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time) 90 !***************************************************************************** 91 91 92 92 julmonday=juldate(19000101,0) ! this is a Monday … … 97 97 98 98 99 100 99 ! For every release point, check whether we are in the release time interval 100 !*************************************************************************** 101 101 102 102 minpart=1 … … 105 105 (itime.le.ireleaseend(i))) then 106 106 107 108 107 ! Determine the local day and time 108 !********************************* 109 109 110 110 xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx ! longitude needed to determine local time … … 124 124 endif 125 125 126 127 128 129 126 ! Calculate a species- and time-dependent correction factor, distinguishing between 127 ! area (those with release starting at surface) and point (release starting above surface) sources 128 ! Also, calculate an average time correction factor (species independent) 129 !***************************************************************************** 130 130 average_timecorrect=0. 131 131 do k=1,nspec … … 139 139 average_timecorrect=average_timecorrect/real(nspec) 140 140 141 142 143 141 ! Determine number of particles to be released this time; at start and at end of release, 142 ! only half the particles are released 143 !***************************************************************************** 144 144 145 145 if (ireleasestart(i).ne.ireleaseend(i)) then … … 149 149 (itime.eq.ireleaseend(i))) rfraction=rfraction/2. 150 150 151 152 153 154 151 ! Take the species-average time correction factor in order to scale the 152 ! number of particles released this time 153 ! Also scale by number of MPI processes 154 !********************************************************************** 155 155 156 156 rfraction=rfraction*average_timecorrect … … 158 158 rfraction=rfraction+xmasssave(i) ! number to be released at this time 159 159 160 160 ! number to be released for one process 161 161 if (mp_partid.lt.mod(int(rfraction),mp_partgroup_np)) then 162 162 addone=1 … … 180 180 numrel=npart(i)/mp_partgroup_np+addone 181 181 endif 182 182 183 183 xaux=xpoint2(i)-xpoint1(i) 184 184 yaux=ypoint2(i)-ypoint1(i) … … 187 187 do ipart=minpart,maxpart_mpi ! search for free storage space 188 188 189 190 189 ! If a free storage space is found, attribute everything to this array element 190 !***************************************************************************** 191 191 192 192 if (itra1(ipart).ne.itime) then 193 193 194 195 196 197 198 194 ! Particle coordinates are determined by using a random position within the release volume 195 !***************************************************************************** 196 197 ! Determine horizontal particle position 198 !*************************************** 199 199 200 200 xtra1(ipart)=xpoint1(i)+ran1(idummy)*xaux … … 207 207 ytra1(ipart)=ypoint1(i)+ran1(idummy)*yaux 208 208 209 210 211 212 213 214 209 ! Assign mass to particle: Total mass divided by total number of particles. 210 ! Time variation has partly been taken into account already by a species-average 211 ! correction factor, by which the number of particles released this time has been 212 ! scaled. Adjust the mass per particle by the species-dependent time correction factor 213 ! divided by the species-average one 214 !***************************************************************************** 215 215 do k=1,nspec 216 217 218 219 220 216 xmass1(ipart,k)=xmass(i,k)/real(npart(i)) & 217 *timecorrect(k)/average_timecorrect 218 ! write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k 219 ! Assign certain properties to particle 220 !************************************** 221 221 end do 222 222 nclass(ipart)=min(int(ran1(idummy)*real(nclassunc))+1, & … … 234 234 235 235 236 237 236 ! Determine vertical particle position 237 !************************************* 238 238 239 239 ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux 240 240 241 242 243 244 245 241 ! Interpolation of topography and density 242 !**************************************** 243 244 ! Determine the nest we are in 245 !***************************** 246 246 247 247 ngrid=0 … … 257 257 43 continue 258 258 259 260 259 ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation 260 !***************************************************************************** 261 261 262 262 if (ngrid.gt.0) then … … 294 294 endif 295 295 296 297 296 ! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters 297 !***************************************************************************** 298 298 if (kindz(i).eq.3) then 299 299 presspart=ztra1(ipart) … … 337 337 endif 338 338 339 340 341 339 ! If release positions are given in meters above sea level, subtract the 340 ! topography from the starting height 341 !*********************************************************************** 342 342 343 343 if (kindz(i).eq.2) ztra1(ipart)=ztra1(ipart)-topo … … 348 348 349 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 350 ! For special simulations, multiply particle concentration air density; 351 ! Simply take the 2nd field in memory to do this (accurate enough) 352 !*********************************************************************** 353 !AF IND_SOURCE switches between different units for concentrations at the source 354 !Af NOTE that in backward simulations the release of particles takes place at the 355 !Af receptor and the sampling at the source. 356 !Af 1="mass" 357 !Af 2="mass mixing ratio" 358 !Af IND_RECEPTOR switches between different units for concentrations at the receptor 359 !Af 1="mass" 360 !Af 2="mass mixing ratio" 361 362 !Af switches for the releasefile: 363 !Af IND_REL = 1 : xmass * rho 364 !Af IND_REL = 0 : xmass * 1 365 366 !Af ind_rel is defined in readcommand.f 367 367 368 368 if (ind_rel .eq. 1) then 369 369 370 371 370 ! Interpolate the air density 371 !**************************** 372 372 373 373 do ii=2,nz … … 403 403 404 404 405 406 405 ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density 406 !******************************************************************** 407 407 408 408 do k=1,nspec … … 416 416 endif 417 417 end do 418 if (ipart.gt.maxpart) goto 996 418 ! ESO TODO: Here we could use dynamic allocation and increase array sizes 419 if (ipart.gt.maxpart_mpi) goto 996 419 420 420 421 34 minpart=ipart+1 421 422 end do 422 423 endif 423 424 end do 424 425 … … 426 427 return 427 428 428 996 429 996 continue 429 430 write(*,*) '#####################################################' 430 431 write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
Note: See TracChangeset
for help on using the changeset viewer.
