Changeset 16b61a5 in flexpart.git for src/boundcond_domainfill_mpi.f90
- Timestamp:
- Oct 14, 2016, 3:19:00 PM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- 4c64400
- Parents:
- 861805a
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundcond_domainfill_mpi.f90
r0f7835d r16b61a5 21 21 22 22 subroutine boundcond_domainfill(itime,loutend) 23 ! i i 24 !***************************************************************************** 25 ! * 26 ! Particles are created by this subroutine continuously throughout the * 27 ! simulation at the boundaries of the domain-filling box. * 28 ! All particles carry the same amount of mass which alltogether comprises the* 29 ! mass of air within the box, which remains (more or less) constant. * 30 ! * 31 ! Author: A. Stohl * 32 ! * 33 ! 16 October 2002 * 34 ! * 35 !***************************************************************************** 36 ! * 37 ! Variables: * 38 ! * 39 ! nx_we(2) grid indices for western and eastern boundary of domain- * 40 ! filling trajectory calculations * 41 ! ny_sn(2) grid indices for southern and northern boundary of domain- * 42 ! filling trajectory calculations * 43 ! * 44 !***************************************************************************** 23 ! i i 24 !***************************************************************************** 25 ! * 26 ! Particles are created by this subroutine continuously throughout the * 27 ! simulation at the boundaries of the domain-filling box. * 28 ! All particles carry the same amount of mass which alltogether comprises the* 29 ! mass of air within the box, which remains (more or less) constant. * 30 ! * 31 ! Author: A. Stohl * 32 ! * 33 ! 16 October 2002 * 34 ! * 35 !***************************************************************************** 36 ! * 37 ! Variables: * 38 ! * 39 ! nx_we(2) grid indices for western and eastern boundary of domain- * 40 ! filling trajectory calculations * 41 ! ny_sn(2) grid indices for southern and northern boundary of domain- * 42 ! filling trajectory calculations * 43 ! * 44 !***************************************************************************** 45 ! CHANGES 46 ! 08/2016 eso: MPI version: 47 ! 48 ! -Root process release particles and distributes to other processes. 49 ! Temporary arrays are used, also for the non-root (receiving) processes. 50 ! -The scheme can be improved by having all processes report numpart 51 ! (keeping track of how many particles have left the domain), so that 52 ! a proportional amount of new particles can be distributed (however 53 ! we have a separate function called from timemanager that will 54 ! redistribute particles among processes if there are imbalances) 55 !***************************************************************************** 45 56 46 57 use point_mod … … 55 66 integer :: itime,in,indz,indzp,i,loutend 56 67 integer :: j,k,ix,jy,m,indzh,indexh,minpart,ipart,mmass 57 integer :: numactiveparticles 68 integer :: numactiveparticles, numpart_total, rel_counter 69 integer,allocatable,dimension(:) :: numrel_mpi !, numactiveparticles_mpi 58 70 59 71 real :: windl(2),rhol(2) … … 66 78 67 79 integer :: idummy = -11 80 integer :: mtag 68 81 logical :: first_call=.true. 69 70 ! Use different seed for each process 82 ! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if 83 ! needed. 84 real,parameter :: maxpartfract=0.1 85 integer :: tmp_size = int(maxpartfract*maxpart) 86 87 ! Use different seed for each process 71 88 if (first_call) then 72 89 idummy=idummy+mp_seed … … 75 92 76 93 77 78 94 ! If domain-filling is global, no boundary conditions are needed 95 !*************************************************************** 79 96 80 97 if (gdomainfill) return … … 82 99 accmasst=0. 83 100 numactiveparticles=0 84 85 ! Terminate trajectories that have left the domain, if domain-filling 86 ! trajectory calculation domain is not global 87 !******************************************************************** 101 ! Keep track of active particles on each process 102 allocate(numrel_mpi(0:mp_partgroup_np-1)) 103 ! numactiveparticles_mpi(0:mp_partgroup_np-1) 104 105 ! New particles to be released on each process 106 numrel_mpi(:)=0 107 108 ! Terminate trajectories that have left the domain, if domain-filling 109 ! trajectory calculation domain is not global. Done for all processes 110 !******************************************************************** 88 111 89 112 do i=1,numpart … … 98 121 numactiveparticles+1 99 122 end do 100 101 102 ! Determine auxiliary variables for time interpolation 103 !***************************************************** 104 105 dt1=real(itime-memtime(1)) 106 dt2=real(memtime(2)-itime) 107 dtt=1./(dt1+dt2) 108 109 ! Initialize auxiliary variable used to search for vacant storage space 110 !********************************************************************** 111 112 minpart=1 113 114 !*************************************** 115 ! Western and eastern boundary condition 116 !*************************************** 117 118 ! Loop from south to north 119 !************************* 120 121 do jy=ny_sn(1),ny_sn(2) 122 123 ! Loop over western (index 1) and eastern (index 2) boundary 124 !*********************************************************** 125 126 do k=1,2 127 128 ! Loop over all release locations in a column 129 !******************************************** 130 131 do j=1,numcolumn_we(k,jy) 132 133 ! Determine, for each release location, the area of the corresponding boundary 134 !***************************************************************************** 135 136 if (j.eq.1) then 137 deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2. 138 else if (j.eq.numcolumn_we(k,jy)) then 139 ! deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+ 140 ! + zcolumn_we(k,jy,j))/2. 141 ! In order to avoid taking a very high column for very many particles, 142 ! use the deltaz from one particle below instead 143 deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2. 144 else 145 deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2. 146 endif 147 if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then 148 boundarea=deltaz*111198.5/2.*dy 149 else 150 boundarea=deltaz*111198.5*dy 151 endif 152 153 154 ! Interpolate the wind velocity and density to the release location 155 !****************************************************************** 156 157 ! Determine the model level below the release position 158 !***************************************************** 159 160 do i=2,nz 161 if (height(i).gt.zcolumn_we(k,jy,j)) then 162 indz=i-1 163 indzp=i 164 goto 6 123 ! numactiveparticles_mpi(mp_partid) = numactiveparticles 124 125 126 ! Collect number of active particles from all processes 127 ! call MPI_Allgather(numactiveparticles, 1, MPI_INTEGER, & 128 ! &numactiveparticles_mpi, 1, MPI_INTEGER, mp_comm_used, mp_ierr) 129 130 131 ! Total number of new releases 132 numpart_total = 0 133 134 135 ! This section only done by root process 136 !*************************************** 137 138 if (lroot) then 139 140 ! Use separate arrays for newly released particles 141 !************************************************* 142 143 allocate(itra1_tmp(tmp_size),npoint_tmp(tmp_size),nclass_tmp(tmp_size),& 144 & idt_tmp(tmp_size),itramem_tmp(tmp_size),itrasplit_tmp(tmp_size),& 145 & xtra1_tmp(tmp_size),ytra1_tmp(tmp_size),ztra1_tmp(tmp_size),& 146 & xmass1_tmp(tmp_size, maxspec)) 147 148 ! Initialize all particles as non-existent 149 itra1_tmp(:)=-999999999 150 151 ! Determine auxiliary variables for time interpolation 152 !***************************************************** 153 154 dt1=real(itime-memtime(1)) 155 dt2=real(memtime(2)-itime) 156 dtt=1./(dt1+dt2) 157 158 ! Initialize auxiliary variable used to search for vacant storage space 159 !********************************************************************** 160 161 minpart=1 162 163 !*************************************** 164 ! Western and eastern boundary condition 165 !*************************************** 166 167 ! Loop from south to north 168 !************************* 169 170 do jy=ny_sn(1),ny_sn(2) 171 172 ! Loop over western (index 1) and eastern (index 2) boundary 173 !*********************************************************** 174 175 do k=1,2 176 177 ! Loop over all release locations in a column 178 !******************************************** 179 180 do j=1,numcolumn_we(k,jy) 181 182 ! Determine, for each release location, the area of the corresponding boundary 183 !***************************************************************************** 184 185 if (j.eq.1) then 186 deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2. 187 else if (j.eq.numcolumn_we(k,jy)) then 188 ! deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+ 189 ! + zcolumn_we(k,jy,j))/2. 190 ! In order to avoid taking a very high column for very many particles, 191 ! use the deltaz from one particle below instead 192 deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2. 193 else 194 deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2. 165 195 endif 166 end do 167 6 continue 168 169 ! Vertical distance to the level below and above current position 170 !**************************************************************** 171 172 dz1=zcolumn_we(k,jy,j)-height(indz) 173 dz2=height(indzp)-zcolumn_we(k,jy,j) 174 dz=1./(dz1+dz2) 175 176 ! Vertical and temporal interpolation 177 !************************************ 178 179 do m=1,2 180 indexh=memind(m) 181 do in=1,2 182 indzh=indz+in-1 183 windl(in)=uu(nx_we(k),jy,indzh,indexh) 184 rhol(in)=rho(nx_we(k),jy,indzh,indexh) 185 end do 186 187 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 188 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 189 end do 190 191 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 192 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 193 194 ! Calculate mass flux, divided by number of processes 195 !**************************************************** 196 197 fluxofmass=windx*rhox*boundarea*real(lsynctime)/mp_partgroup_np 198 199 200 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 201 ! if it is out of the domain, set accumulated mass flux to zero 202 !****************************************************************************** 203 204 if (k.eq.1) then 205 if (fluxofmass.ge.0.) then 206 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass 196 if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then 197 boundarea=deltaz*111198.5/2.*dy 207 198 else 208 acc_mass_we(k,jy,j)=0.199 boundarea=deltaz*111198.5*dy 209 200 endif 210 else 211 if (fluxofmass.le.0.) then 212 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass) 213 else 214 acc_mass_we(k,jy,j)=0. 215 endif 216 endif 217 accmasst=accmasst+acc_mass_we(k,jy,j) 218 219 ! If the accumulated mass exceeds half the mass that each particle shall carry, 220 ! one (or more) particle(s) is (are) released and the accumulated mass is 221 ! reduced by the mass of this (these) particle(s) 222 !****************************************************************************** 223 224 if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then 225 mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ & 226 xmassperparticle) 227 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- & 228 real(mmass)*xmassperparticle 229 else 230 mmass=0 231 endif 232 233 do m=1,mmass 234 do ipart=minpart,maxpart_mpi 235 236 ! If a vacant storage space is found, attribute everything to this array element 237 !***************************************************************************** 238 239 if (itra1(ipart).ne.itime) then 240 241 ! Assign particle positions 242 !************************** 243 244 xtra1(ipart)=real(nx_we(k)) 245 if (jy.eq.ny_sn(1)) then 246 ytra1(ipart)=real(jy)+0.5*ran1(idummy) 247 else if (jy.eq.ny_sn(2)) then 248 ytra1(ipart)=real(jy)-0.5*ran1(idummy) 249 else 250 ytra1(ipart)=real(jy)+(ran1(idummy)-.5) 251 endif 252 if (j.eq.1) then 253 ztra1(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- & 254 zcolumn_we(k,jy,1))/4. 255 else if (j.eq.numcolumn_we(k,jy)) then 256 ztra1(ipart)=(2.*zcolumn_we(k,jy,j)+ & 257 zcolumn_we(k,jy,j-1)+height(nz))/4. 258 else 259 ztra1(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* & 260 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1)) 261 endif 262 263 ! Interpolate PV to the particle position 264 !**************************************** 265 ixm=int(xtra1(ipart)) 266 jym=int(ytra1(ipart)) 267 ixp=ixm+1 268 jyp=jym+1 269 ddx=xtra1(ipart)-real(ixm) 270 ddy=ytra1(ipart)-real(jym) 271 rddx=1.-ddx 272 rddy=1.-ddy 273 p1=rddx*rddy 274 p2=ddx*rddy 275 p3=rddx*ddy 276 p4=ddx*ddy 277 do i=2,nz 278 if (height(i).gt.ztra1(ipart)) then 279 indzm=i-1 280 indzp=i 281 goto 26 282 endif 283 end do 284 26 continue 285 dz1=ztra1(ipart)-height(indzm) 286 dz2=height(indzp)-ztra1(ipart) 287 dz=1./(dz1+dz2) 288 do mm=1,2 289 indexh=memind(mm) 290 do in=1,2 291 indzh=indzm+in-1 292 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 293 +p2*pv(ixp,jym,indzh,indexh) & 294 +p3*pv(ixm,jyp,indzh,indexh) & 295 +p4*pv(ixp,jyp,indzh,indexh) 296 end do 297 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 298 end do 299 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 300 ylat=ylat0+ytra1(ipart)*dy 301 if (ylat.lt.0.) pvpart=-1.*pvpart 302 303 304 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 305 !***************************************************************************** 306 307 if (((ztra1(ipart).gt.3000.).and. & 308 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 309 nclass(ipart)=min(int(ran1(idummy)* & 310 real(nclassunc))+1,nclassunc) 311 numactiveparticles=numactiveparticles+1 312 numparticlecount=numparticlecount+1 313 npoint(ipart)=numparticlecount 314 idt(ipart)=mintime 315 itra1(ipart)=itime 316 itramem(ipart)=itra1(ipart) 317 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit 318 xmass1(ipart,1)=xmassperparticle 319 if (mdomainfill.eq.2) xmass1(ipart,1)= & 320 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9 321 else 322 goto 71 323 endif 324 325 326 ! Increase numpart, if necessary 327 !******************************* 328 329 numpart=max(numpart,ipart) 330 goto 73 ! Storage space has been found, stop searching 201 202 203 ! Interpolate the wind velocity and density to the release location 204 !****************************************************************** 205 206 ! Determine the model level below the release position 207 !***************************************************** 208 209 do i=2,nz 210 if (height(i).gt.zcolumn_we(k,jy,j)) then 211 indz=i-1 212 indzp=i 213 goto 6 331 214 endif 332 215 end do 333 if (ipart.gt.maxpart_mpi) & 334 stop 'boundcond_domainfill.f: too many particles required' 335 73 minpart=ipart+1 336 71 continue 216 6 continue 217 218 ! Vertical distance to the level below and above current position 219 !**************************************************************** 220 221 dz1=zcolumn_we(k,jy,j)-height(indz) 222 dz2=height(indzp)-zcolumn_we(k,jy,j) 223 dz=1./(dz1+dz2) 224 225 ! Vertical and temporal interpolation 226 !************************************ 227 228 do m=1,2 229 indexh=memind(m) 230 do in=1,2 231 indzh=indz+in-1 232 windl(in)=uu(nx_we(k),jy,indzh,indexh) 233 rhol(in)=rho(nx_we(k),jy,indzh,indexh) 234 end do 235 236 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 237 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 238 end do 239 240 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 241 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 242 243 ! Calculate mass flux 244 !******************** 245 246 fluxofmass=windx*rhox*boundarea*real(lsynctime) 247 248 249 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 250 ! if it is out of the domain, set accumulated mass flux to zero 251 !****************************************************************************** 252 253 if (k.eq.1) then 254 if (fluxofmass.ge.0.) then 255 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass 256 else 257 acc_mass_we(k,jy,j)=0. 258 endif 259 else 260 if (fluxofmass.le.0.) then 261 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass) 262 else 263 acc_mass_we(k,jy,j)=0. 264 endif 265 endif 266 accmasst=accmasst+acc_mass_we(k,jy,j) 267 268 ! If the accumulated mass exceeds half the mass that each particle shall carry, 269 ! one (or more) particle(s) is (are) released and the accumulated mass is 270 ! reduced by the mass of this (these) particle(s) 271 !****************************************************************************** 272 273 if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then 274 mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ & 275 xmassperparticle) 276 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- & 277 real(mmass)*xmassperparticle 278 else 279 mmass=0 280 endif 281 282 do m=1,mmass 283 do ipart=minpart,maxpart 284 285 ! If a vacant storage space is found, attribute everything to this array element 286 ! TODO: for the MPI version this test can be removed, as all 287 ! elements in _tmp arrays are initialized to zero 288 !***************************************************************************** 289 290 if (itra1_tmp(ipart).ne.itime) then 291 292 ! Assign particle positions 293 !************************** 294 295 xtra1_tmp(ipart)=real(nx_we(k)) 296 if (jy.eq.ny_sn(1)) then 297 ytra1_tmp(ipart)=real(jy)+0.5*ran1(idummy) 298 else if (jy.eq.ny_sn(2)) then 299 ytra1_tmp(ipart)=real(jy)-0.5*ran1(idummy) 300 else 301 ytra1_tmp(ipart)=real(jy)+(ran1(idummy)-.5) 302 endif 303 if (j.eq.1) then 304 ztra1_tmp(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- & 305 zcolumn_we(k,jy,1))/4. 306 else if (j.eq.numcolumn_we(k,jy)) then 307 ztra1_tmp(ipart)=(2.*zcolumn_we(k,jy,j)+ & 308 zcolumn_we(k,jy,j-1)+height(nz))/4. 309 else 310 ztra1_tmp(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* & 311 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1)) 312 endif 313 314 ! Interpolate PV to the particle position 315 !**************************************** 316 ixm=int(xtra1_tmp(ipart)) 317 jym=int(ytra1_tmp(ipart)) 318 ixp=ixm+1 319 jyp=jym+1 320 ddx=xtra1_tmp(ipart)-real(ixm) 321 ddy=ytra1_tmp(ipart)-real(jym) 322 rddx=1.-ddx 323 rddy=1.-ddy 324 p1=rddx*rddy 325 p2=ddx*rddy 326 p3=rddx*ddy 327 p4=ddx*ddy 328 do i=2,nz 329 if (height(i).gt.ztra1_tmp(ipart)) then 330 indzm=i-1 331 indzp=i 332 goto 26 333 endif 334 end do 335 26 continue 336 dz1=ztra1_tmp(ipart)-height(indzm) 337 dz2=height(indzp)-ztra1_tmp(ipart) 338 dz=1./(dz1+dz2) 339 do mm=1,2 340 indexh=memind(mm) 341 do in=1,2 342 indzh=indzm+in-1 343 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 344 +p2*pv(ixp,jym,indzh,indexh) & 345 +p3*pv(ixm,jyp,indzh,indexh) & 346 +p4*pv(ixp,jyp,indzh,indexh) 347 end do 348 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 349 end do 350 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 351 ylat=ylat0+ytra1_tmp(ipart)*dy 352 if (ylat.lt.0.) pvpart=-1.*pvpart 353 354 355 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 356 !***************************************************************************** 357 358 if (((ztra1_tmp(ipart).gt.3000.).and. & 359 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 360 nclass_tmp(ipart)=min(int(ran1(idummy)* & 361 real(nclassunc))+1,nclassunc) 362 numactiveparticles=numactiveparticles+1 363 numparticlecount=numparticlecount+1 364 npoint_tmp(ipart)=numparticlecount 365 idt_tmp(ipart)=mintime 366 itra1_tmp(ipart)=itime 367 itramem_tmp(ipart)=itra1_tmp(ipart) 368 itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit 369 xmass1_tmp(ipart,1)=xmassperparticle 370 if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= & 371 xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9 372 else 373 goto 71 374 endif 375 376 377 ! Increase numpart, if necessary 378 !******************************* 379 380 numpart_total=max(numpart_total,ipart) 381 goto 73 ! Storage space has been found, stop searching 382 endif 383 end do 384 if (ipart.gt.tmp_size) & 385 stop 'boundcond_domainfill_mpi.f90: too many particles required' 386 73 minpart=ipart+1 387 71 continue 388 end do 389 390 337 391 end do 338 339 340 392 end do 341 393 end do 342 end do 343 344 345 !***************************************** 346 ! Southern and northern boundary condition 347 !***************************************** 348 349 ! Loop from west to east 350 !*********************** 351 352 do ix=nx_we(1),nx_we(2) 353 354 ! Loop over southern (index 1) and northern (index 2) boundary 355 !************************************************************* 356 357 do k=1,2 358 ylat=ylat0+real(ny_sn(k))*dy 359 cosfact=cos(ylat*pi180) 360 361 ! Loop over all release locations in a column 362 !******************************************** 363 364 do j=1,numcolumn_sn(k,ix) 365 366 ! Determine, for each release location, the area of the corresponding boundary 367 !***************************************************************************** 368 369 if (j.eq.1) then 370 deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2. 371 else if (j.eq.numcolumn_sn(k,ix)) then 372 ! deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+ 373 ! + zcolumn_sn(k,ix,j))/2. 374 ! In order to avoid taking a very high column for very many particles, 375 ! use the deltaz from one particle below instead 376 deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2. 377 else 378 deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2. 379 endif 380 if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then 381 boundarea=deltaz*111198.5/2.*cosfact*dx 382 else 383 boundarea=deltaz*111198.5*cosfact*dx 384 endif 385 386 387 ! Interpolate the wind velocity and density to the release location 388 !****************************************************************** 389 390 ! Determine the model level below the release position 391 !***************************************************** 392 393 do i=2,nz 394 if (height(i).gt.zcolumn_sn(k,ix,j)) then 395 indz=i-1 396 indzp=i 397 goto 16 394 395 396 !***************************************** 397 ! Southern and northern boundary condition 398 !***************************************** 399 400 ! Loop from west to east 401 !*********************** 402 403 do ix=nx_we(1),nx_we(2) 404 405 ! Loop over southern (index 1) and northern (index 2) boundary 406 !************************************************************* 407 408 do k=1,2 409 ylat=ylat0+real(ny_sn(k))*dy 410 cosfact=cos(ylat*pi180) 411 412 ! Loop over all release locations in a column 413 !******************************************** 414 415 do j=1,numcolumn_sn(k,ix) 416 417 ! Determine, for each release location, the area of the corresponding boundary 418 !***************************************************************************** 419 420 if (j.eq.1) then 421 deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2. 422 else if (j.eq.numcolumn_sn(k,ix)) then 423 ! deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+ 424 ! + zcolumn_sn(k,ix,j))/2. 425 ! In order to avoid taking a very high column for very many particles, 426 ! use the deltaz from one particle below instead 427 deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2. 428 else 429 deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2. 398 430 endif 399 end do 400 16 continue 401 402 ! Vertical distance to the level below and above current position 403 !**************************************************************** 404 405 dz1=zcolumn_sn(k,ix,j)-height(indz) 406 dz2=height(indzp)-zcolumn_sn(k,ix,j) 407 dz=1./(dz1+dz2) 408 409 ! Vertical and temporal interpolation 410 !************************************ 411 412 do m=1,2 413 indexh=memind(m) 414 do in=1,2 415 indzh=indz+in-1 416 windl(in)=vv(ix,ny_sn(k),indzh,indexh) 417 rhol(in)=rho(ix,ny_sn(k),indzh,indexh) 418 end do 419 420 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 421 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 422 end do 423 424 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 425 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 426 427 ! Calculate mass flux, divided by number of processes 428 !**************************************************** 429 430 fluxofmass=windx*rhox*boundarea*real(lsynctime)/mp_partgroup_np 431 432 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 433 ! if it is out of the domain, set accumulated mass flux to zero 434 !****************************************************************************** 435 436 if (k.eq.1) then 437 if (fluxofmass.ge.0.) then 438 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass 431 if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then 432 boundarea=deltaz*111198.5/2.*cosfact*dx 439 433 else 440 acc_mass_sn(k,ix,j)=0.434 boundarea=deltaz*111198.5*cosfact*dx 441 435 endif 442 else 443 if (fluxofmass.le.0.) then 444 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass) 445 else 446 acc_mass_sn(k,ix,j)=0. 447 endif 448 endif 449 accmasst=accmasst+acc_mass_sn(k,ix,j) 450 451 ! If the accumulated mass exceeds half the mass that each particle shall carry, 452 ! one (or more) particle(s) is (are) released and the accumulated mass is 453 ! reduced by the mass of this (these) particle(s) 454 !****************************************************************************** 455 456 if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then 457 mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ & 458 xmassperparticle) 459 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- & 460 real(mmass)*xmassperparticle 461 else 462 mmass=0 463 endif 464 465 do m=1,mmass 466 do ipart=minpart,maxpart_mpi 467 468 ! If a vacant storage space is found, attribute everything to this array element 469 !***************************************************************************** 470 471 if (itra1(ipart).ne.itime) then 472 473 ! Assign particle positions 474 !************************** 475 476 ytra1(ipart)=real(ny_sn(k)) 477 if (ix.eq.nx_we(1)) then 478 xtra1(ipart)=real(ix)+0.5*ran1(idummy) 479 else if (ix.eq.nx_we(2)) then 480 xtra1(ipart)=real(ix)-0.5*ran1(idummy) 481 else 482 xtra1(ipart)=real(ix)+(ran1(idummy)-.5) 483 endif 484 if (j.eq.1) then 485 ztra1(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- & 486 zcolumn_sn(k,ix,1))/4. 487 else if (j.eq.numcolumn_sn(k,ix)) then 488 ztra1(ipart)=(2.*zcolumn_sn(k,ix,j)+ & 489 zcolumn_sn(k,ix,j-1)+height(nz))/4. 490 else 491 ztra1(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* & 492 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1)) 493 endif 494 495 496 ! Interpolate PV to the particle position 497 !**************************************** 498 ixm=int(xtra1(ipart)) 499 jym=int(ytra1(ipart)) 500 ixp=ixm+1 501 jyp=jym+1 502 ddx=xtra1(ipart)-real(ixm) 503 ddy=ytra1(ipart)-real(jym) 504 rddx=1.-ddx 505 rddy=1.-ddy 506 p1=rddx*rddy 507 p2=ddx*rddy 508 p3=rddx*ddy 509 p4=ddx*ddy 510 do i=2,nz 511 if (height(i).gt.ztra1(ipart)) then 512 indzm=i-1 513 indzp=i 514 goto 126 515 endif 516 end do 517 126 continue 518 dz1=ztra1(ipart)-height(indzm) 519 dz2=height(indzp)-ztra1(ipart) 520 dz=1./(dz1+dz2) 521 do mm=1,2 522 indexh=memind(mm) 523 do in=1,2 524 indzh=indzm+in-1 525 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 526 +p2*pv(ixp,jym,indzh,indexh) & 527 +p3*pv(ixm,jyp,indzh,indexh) & 528 +p4*pv(ixp,jyp,indzh,indexh) 529 end do 530 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 531 end do 532 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 533 if (ylat.lt.0.) pvpart=-1.*pvpart 534 535 536 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 537 !***************************************************************************** 538 539 if (((ztra1(ipart).gt.3000.).and. & 540 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 541 nclass(ipart)=min(int(ran1(idummy)* & 542 real(nclassunc))+1,nclassunc) 543 numactiveparticles=numactiveparticles+1 544 numparticlecount=numparticlecount+1 545 npoint(ipart)=numparticlecount 546 idt(ipart)=mintime 547 itra1(ipart)=itime 548 itramem(ipart)=itra1(ipart) 549 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit 550 xmass1(ipart,1)=xmassperparticle 551 if (mdomainfill.eq.2) xmass1(ipart,1)= & 552 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9 553 else 554 goto 171 555 endif 556 557 558 ! Increase numpart, if necessary 559 !******************************* 560 numpart=max(numpart,ipart) 561 goto 173 ! Storage space has been found, stop searching 436 437 438 ! Interpolate the wind velocity and density to the release location 439 !****************************************************************** 440 441 ! Determine the model level below the release position 442 !***************************************************** 443 444 do i=2,nz 445 if (height(i).gt.zcolumn_sn(k,ix,j)) then 446 indz=i-1 447 indzp=i 448 goto 16 562 449 endif 563 450 end do 564 if (ipart.gt.maxpart_mpi) & 565 stop 'boundcond_domainfill.f: too many particles required' 566 173 minpart=ipart+1 567 171 continue 451 16 continue 452 453 ! Vertical distance to the level below and above current position 454 !**************************************************************** 455 456 dz1=zcolumn_sn(k,ix,j)-height(indz) 457 dz2=height(indzp)-zcolumn_sn(k,ix,j) 458 dz=1./(dz1+dz2) 459 460 ! Vertical and temporal interpolation 461 !************************************ 462 463 do m=1,2 464 indexh=memind(m) 465 do in=1,2 466 indzh=indz+in-1 467 windl(in)=vv(ix,ny_sn(k),indzh,indexh) 468 rhol(in)=rho(ix,ny_sn(k),indzh,indexh) 469 end do 470 471 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz 472 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz 473 end do 474 475 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt 476 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt 477 478 ! Calculate mass flux 479 !******************** 480 481 fluxofmass=windx*rhox*boundarea*real(lsynctime) 482 483 ! If the mass flux is directed into the domain, add it to previous mass fluxes; 484 ! if it is out of the domain, set accumulated mass flux to zero 485 !****************************************************************************** 486 487 if (k.eq.1) then 488 if (fluxofmass.ge.0.) then 489 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass 490 else 491 acc_mass_sn(k,ix,j)=0. 492 endif 493 else 494 if (fluxofmass.le.0.) then 495 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass) 496 else 497 acc_mass_sn(k,ix,j)=0. 498 endif 499 endif 500 accmasst=accmasst+acc_mass_sn(k,ix,j) 501 502 ! If the accumulated mass exceeds half the mass that each particle shall carry, 503 ! one (or more) particle(s) is (are) released and the accumulated mass is 504 ! reduced by the mass of this (these) particle(s) 505 !****************************************************************************** 506 507 if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then 508 mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ & 509 xmassperparticle) 510 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- & 511 real(mmass)*xmassperparticle 512 else 513 mmass=0 514 endif 515 516 do m=1,mmass 517 do ipart=minpart,maxpart 518 519 ! If a vacant storage space is found, attribute everything to this array element 520 !***************************************************************************** 521 522 if (itra1_tmp(ipart).ne.itime) then 523 524 ! Assign particle positions 525 !************************** 526 527 ytra1_tmp(ipart)=real(ny_sn(k)) 528 if (ix.eq.nx_we(1)) then 529 xtra1_tmp(ipart)=real(ix)+0.5*ran1(idummy) 530 else if (ix.eq.nx_we(2)) then 531 xtra1_tmp(ipart)=real(ix)-0.5*ran1(idummy) 532 else 533 xtra1_tmp(ipart)=real(ix)+(ran1(idummy)-.5) 534 endif 535 if (j.eq.1) then 536 ztra1_tmp(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- & 537 zcolumn_sn(k,ix,1))/4. 538 else if (j.eq.numcolumn_sn(k,ix)) then 539 ztra1_tmp(ipart)=(2.*zcolumn_sn(k,ix,j)+ & 540 zcolumn_sn(k,ix,j-1)+height(nz))/4. 541 else 542 ztra1_tmp(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* & 543 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1)) 544 endif 545 546 547 ! Interpolate PV to the particle position 548 !**************************************** 549 ixm=int(xtra1_tmp(ipart)) 550 jym=int(ytra1_tmp(ipart)) 551 ixp=ixm+1 552 jyp=jym+1 553 ddx=xtra1_tmp(ipart)-real(ixm) 554 ddy=ytra1_tmp(ipart)-real(jym) 555 rddx=1.-ddx 556 rddy=1.-ddy 557 p1=rddx*rddy 558 p2=ddx*rddy 559 p3=rddx*ddy 560 p4=ddx*ddy 561 do i=2,nz 562 if (height(i).gt.ztra1_tmp(ipart)) then 563 indzm=i-1 564 indzp=i 565 goto 126 566 endif 567 end do 568 126 continue 569 dz1=ztra1_tmp(ipart)-height(indzm) 570 dz2=height(indzp)-ztra1_tmp(ipart) 571 dz=1./(dz1+dz2) 572 do mm=1,2 573 indexh=memind(mm) 574 do in=1,2 575 indzh=indzm+in-1 576 y1(in)=p1*pv(ixm,jym,indzh,indexh) & 577 +p2*pv(ixp,jym,indzh,indexh) & 578 +p3*pv(ixm,jyp,indzh,indexh) & 579 +p4*pv(ixp,jyp,indzh,indexh) 580 end do 581 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz 582 end do 583 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt 584 if (ylat.lt.0.) pvpart=-1.*pvpart 585 586 587 ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere 588 !***************************************************************************** 589 590 if (((ztra1_tmp(ipart).gt.3000.).and. & 591 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then 592 nclass_tmp(ipart)=min(int(ran1(idummy)* & 593 real(nclassunc))+1,nclassunc) 594 numactiveparticles=numactiveparticles+1 595 numparticlecount=numparticlecount+1 596 npoint_tmp(ipart)=numparticlecount 597 idt_tmp(ipart)=mintime 598 itra1_tmp(ipart)=itime 599 itramem_tmp(ipart)=itra1_tmp(ipart) 600 itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit 601 xmass1_tmp(ipart,1)=xmassperparticle 602 if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= & 603 xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9 604 else 605 goto 171 606 endif 607 608 609 ! Increase numpart, if necessary 610 !******************************* 611 numpart_total=max(numpart_total,ipart) 612 goto 173 ! Storage space has been found, stop searching 613 endif 614 end do 615 if (ipart.gt.tmp_size) & 616 stop 'boundcond_domainfill.f: too many particles required' 617 173 minpart=ipart+1 618 171 continue 619 end do 620 621 568 622 end do 569 570 571 623 end do 572 624 end do 625 626 627 ! xm=0. 628 ! do i=1,numpart_total 629 ! if (itra1_tmp(i).eq.itime) xm=xm+xmass1(i,1) 630 ! end do 631 632 !write(*,*) itime,numactiveparticles,numparticlecount,numpart, 633 ! +xm,accmasst,xm+accmasst 634 635 end if ! if lroot 636 637 ! Distribute the number of particles to be released 638 ! ************************************************* 639 call MPI_Bcast(numpart_total, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr) 640 641 do i=0, mp_partgroup_np-1 642 numrel_mpi(i) = numpart_total/mp_partgroup_np 643 if (i.lt.mod(numpart_total,mp_partgroup_np)) numrel_mpi(i) = numrel_mpi(i) + 1 573 644 end do 574 645 575 576 xm=0. 577 do i=1,numpart 578 if (itra1(i).eq.itime) xm=xm+xmass1(i,1) 646 ! Allocate temporary arrays for receiving processes 647 if (.not.lroot) then 648 allocate(itra1_tmp(numrel_mpi(mp_partid)),& 649 & npoint_tmp(numrel_mpi(mp_partid)),& 650 & nclass_tmp(numrel_mpi(mp_partid)),& 651 & idt_tmp(numrel_mpi(mp_partid)),& 652 & itramem_tmp(numrel_mpi(mp_partid)),& 653 & itrasplit_tmp(numrel_mpi(mp_partid)),& 654 & xtra1_tmp(numrel_mpi(mp_partid)),& 655 & ytra1_tmp(numrel_mpi(mp_partid)),& 656 & ztra1_tmp(numrel_mpi(mp_partid)),& 657 & xmass1_tmp(numrel_mpi(mp_partid),maxspec)) 658 659 ! Initialize all particles as non-existent 660 itra1_tmp(:)=-999999999 661 end if 662 663 ! Distribute particles 664 ! Keep track of released particles so far 665 rel_counter = 0 666 mtag = 1000 667 668 do i=0, mp_partgroup_np-1 669 670 ! For root process, nothing to do except update release count 671 if (i.eq.0) then 672 rel_counter = rel_counter + numrel_mpi(i) 673 cycle 674 end if 675 676 ! Send particles from root to non-root processes 677 if (lroot.and.numrel_mpi(i).gt.0) then 678 679 call MPI_SEND(nclass_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 680 &numrel_mpi(i),MPI_INTEGER,i,mtag+1*i,mp_comm_used,mp_ierr) 681 682 call MPI_SEND(npoint_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 683 &numrel_mpi(i),MPI_INTEGER,i,mtag+2*i,mp_comm_used,mp_ierr) 684 685 call MPI_SEND(itra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 686 &numrel_mpi(i),MPI_INTEGER,i,mtag+3*i,mp_comm_used,mp_ierr) 687 688 call MPI_SEND(idt_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 689 &numrel_mpi(i),MPI_INTEGER,i,mtag+4*i,mp_comm_used,mp_ierr) 690 691 call MPI_SEND(itramem_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 692 &numrel_mpi(i),MPI_INTEGER,i,mtag+5*i,mp_comm_used,mp_ierr) 693 694 call MPI_SEND(itrasplit_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 695 &numrel_mpi(i),MPI_INTEGER,i,mtag+6*i,mp_comm_used,mp_ierr) 696 697 call MPI_SEND(xtra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 698 &numrel_mpi(i),mp_dp,i,mtag+7*i,mp_comm_used,mp_ierr) 699 700 call MPI_SEND(ytra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 701 &numrel_mpi(i),mp_dp,i,mtag+8*i,mp_comm_used,mp_ierr) 702 703 call MPI_SEND(ztra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),& 704 &numrel_mpi(i),mp_sp,i,mtag+9*i,mp_comm_used,mp_ierr) 705 706 do j=1,nspec 707 call MPI_SEND(xmass1_tmp(rel_counter+1:rel_counter+numrel_mpi(i),j),& 708 &numrel_mpi(i),mp_sp,i,mtag+(9+j)*i,mp_comm_used,mp_ierr) 709 end do 710 711 ! Non-root processes issue receive requests 712 else if (i.eq.mp_partid.and.numrel_mpi(i).gt.0) then 713 call MPI_RECV(nclass_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 714 &MPI_INTEGER,id_root,mtag+1*i,mp_comm_used,mp_status,mp_ierr) 715 716 call MPI_RECV(npoint_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 717 &MPI_INTEGER,id_root,mtag+2*i,mp_comm_used,mp_status,mp_ierr) 718 719 call MPI_RECV(itra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 720 &MPI_INTEGER,id_root,mtag+3*i,mp_comm_used,mp_status,mp_ierr) 721 722 call MPI_RECV(idt_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 723 &MPI_INTEGER,id_root,mtag+4*i,mp_comm_used,mp_status,mp_ierr) 724 725 call MPI_RECV(itramem_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 726 &MPI_INTEGER,id_root,mtag+5*i,mp_comm_used,mp_status,mp_ierr) 727 728 call MPI_RECV(itrasplit_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 729 &MPI_INTEGER,id_root,mtag+6*i,mp_comm_used,mp_status,mp_ierr) 730 731 call MPI_RECV(xtra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 732 &mp_dp,id_root,mtag+7*i,mp_comm_used,mp_status,mp_ierr) 733 734 call MPI_RECV(ytra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 735 &mp_dp,id_root,mtag+8*i,mp_comm_used,mp_status,mp_ierr) 736 737 call MPI_RECV(ztra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),& 738 &mp_sp,id_root,mtag+9*i,mp_comm_used,mp_status,mp_ierr) 739 740 do j=1,nspec 741 call MPI_RECV(xmass1_tmp(1:numrel_mpi(i),j),numrel_mpi(i),& 742 &mp_sp,id_root,mtag+(9+j)*i,mp_comm_used,mp_status,mp_ierr) 743 744 end do 745 end if 746 rel_counter = rel_counter + numrel_mpi(i) 579 747 end do 580 748 581 !write(*,*) itime,numactiveparticles,numparticlecount,numpart, 582 ! +xm,accmasst,xm+accmasst 583 584 585 ! If particles shall be dumped, then accumulated masses at the domain boundaries 586 ! must be dumped, too, to be used for later runs 587 !***************************************************************************** 588 589 ! :TODO: eso parallelize 749 ! Find free storage space for the new particles. 750 ! This section is independent of the redistribution scheme used 751 ! ******************************************************************** 752 753 ! Keep track of released particles so far 754 minpart=1 755 756 ! The algorithm should be correct also for root process 757 do i=1, numrel_mpi(mp_partid) 758 do ipart=minpart, maxpart 759 if (itra1(ipart).ne.itime) then 760 itra1(ipart) = itra1_tmp(i) 761 npoint(ipart) = npoint_tmp(i) 762 nclass(ipart) = nclass_tmp(i) 763 idt(ipart) = idt_tmp(i) 764 itramem(ipart) = itramem_tmp(i) 765 itrasplit(ipart) = itrasplit_tmp(i) 766 xtra1(ipart) = xtra1_tmp(i) 767 ytra1(ipart) = ytra1_tmp(i) 768 ztra1(ipart) = ztra1_tmp(i) 769 xmass1(ipart,:) = xmass1_tmp(i,:) 770 ! Increase numpart, if necessary 771 numpart=max(numpart,ipart) 772 goto 200 ! Storage space has been found, stop searching 773 end if 774 end do 775 200 minpart=ipart+1 776 end do 777 778 ! If particles shall be dumped, then accumulated masses at the domain boundaries 779 ! must be dumped, too, to be used for later runs 780 !***************************************************************************** 781 590 782 if ((ipout.gt.0).and.(itime.eq.loutend)) then 591 783 if (lroot) then 784 call mpif_mtime('iotime',0) 592 785 open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', & 593 786 form='unformatted') … … 595 788 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn 596 789 close(unitboundcond) 790 call mpif_mtime('iotime',1) 597 791 end if 598 792 endif 599 793 794 ! Deallocate temporary arrays 795 deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,& 796 & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp,numrel_mpi) 797 ! numactiveparticles_mpi 798 799 600 800 end subroutine boundcond_domainfill
Note: See TracChangeset
for help on using the changeset viewer.