69 integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
70 real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
71 real(kind=dp) :: jul1,jul2,
juldate
72 character(len=50) :: line
77 open(unitreleases,file=path(1)(1:length(1))//
'RELEASES',status=
'old', &
85 call
skplin(12,unitreleases)
86 read (unitreleases,901) line
88 if (index(line,
'Total') .eq. 0)
then
99 call
skplin(11,unitreleases)
102 read(unitreleases,*,err=998) nspec
103 if (old) call
skplin(2,unitreleases)
105 read(unitreleases,*,err=998) specnum_rel
106 if (old) call
skplin(2,unitreleases)
110 100 numpoint=numpoint+1
111 read(unitreleases,*,end=25)
112 read(unitreleases,*,err=998,end=25) idum,idum
113 if (old) call
skplin(2,unitreleases)
114 read(unitreleases,*,err=998) idum,idum
115 if (old) call
skplin(2,unitreleases)
116 read(unitreleases,*,err=998) xdum
117 if (old) call
skplin(2,unitreleases)
118 read(unitreleases,*,err=998) xdum
119 if (old) call
skplin(2,unitreleases)
120 read(unitreleases,*,err=998) xdum
121 if (old) call
skplin(2,unitreleases)
122 read(unitreleases,*,err=998) xdum
123 if (old) call
skplin(2,unitreleases)
124 read(unitreleases,*,err=998) idum
125 if (old) call
skplin(2,unitreleases)
126 read(unitreleases,*,err=998) xdum
127 if (old) call
skplin(2,unitreleases)
128 read(unitreleases,*,err=998) xdum
129 if (old) call
skplin(2,unitreleases)
130 read(unitreleases,*,err=998) idum
131 if (old) call
skplin(2,unitreleases)
133 read(unitreleases,*,err=998) xdum
134 if (old) call
skplin(2,unitreleases)
137 read(unitreleases,
'(a40)',err=998) compoint(1)(1:40)
138 if (old) call
skplin(1,unitreleases)
142 25 numpoint=numpoint-1
145 allocate(ireleasestart(numpoint) &
147 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
148 allocate(ireleaseend(numpoint) &
150 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
151 allocate(xpoint1(numpoint) &
153 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
154 allocate(xpoint2(numpoint) &
156 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
157 allocate(ypoint1(numpoint) &
159 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
160 allocate(ypoint2(numpoint) &
162 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
163 allocate(zpoint1(numpoint) &
165 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
166 allocate(zpoint2(numpoint) &
168 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
169 allocate(kindz(numpoint) &
171 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
172 allocate(xmass(numpoint,maxspec) &
174 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
175 allocate(rho_rel(numpoint) &
177 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
178 allocate(npart(numpoint) &
180 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
181 allocate(xmasssave(numpoint) &
183 if (stat.ne.0)
write(*,*)
'ERROR: could not allocate RELEASPOINT'
185 write (*,*)
' Releasepoints allocated: ', numpoint
197 drydepspec(i)=.false.
206 call
skplin(11,unitreleases)
212 read(unitreleases,*,err=998) nspec
213 if (nspec.gt.maxspec) goto 994
214 if (old) call
skplin(2,unitreleases)
216 read(unitreleases,*,err=998) specnum_rel
217 if (old) call
skplin(2,unitreleases)
234 if (((iout.eq.2).or.(iout.eq.3)).and. &
235 (weightmolar(i).lt.0.))
then
236 write(*,*)
'For mixing ratio output, valid molar weight'
237 write(*,*)
'must be specified for all simulated species.'
238 write(*,*)
'Check table SPECIES or choose concentration'
239 write(*,*)
'output instead if molar weight is not known.'
247 decay(i)=0.693147/decay(i)
253 if (reldiff(i).gt.0.) &
254 rm(i)=1./(henry(i)/3000.+100.*f0(i))
261 dquer(i)=dquer(i)*1000000.
262 if (density(i).gt.0.)
then
263 call
part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
268 cunningham(i)=cunningham(i)+cun*fract(i,j)
269 vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
271 write(*,*)
'Average setting velocity: ',i,vsetaver(i)
277 dryvel(i)=dryvel(i)*0.01
281 if (weta(i).gt.0.)
then
283 write (*,*)
'Wetdeposition switched on: ',weta(i),i
285 if (ohreact(i).gt.0)
then
287 write (*,*)
'OHreaction switched on: ',ohreact(i),i
291 if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
292 (dryvel(i).gt.0.))
then
299 if (wetdep.or.drydep) dep=.true.
307 1000 numpoint=numpoint+1
308 read(unitreleases,*,end=250)
309 read(unitreleases,*,err=998,end=250) id1,it1
310 if (old) call
skplin(2,unitreleases)
311 read(unitreleases,*,err=998) id2,it2
312 if (old) call
skplin(2,unitreleases)
313 read(unitreleases,*,err=998) xpoint1(numpoint)
314 if (old) call
skplin(2,unitreleases)
315 read(unitreleases,*,err=998) ypoint1(numpoint)
316 if (old) call
skplin(2,unitreleases)
317 read(unitreleases,*,err=998) xpoint2(numpoint)
318 if (old) call
skplin(2,unitreleases)
319 read(unitreleases,*,err=998) ypoint2(numpoint)
320 if (old) call
skplin(2,unitreleases)
321 read(unitreleases,*,err=998) kindz(numpoint)
322 if (old) call
skplin(2,unitreleases)
323 read(unitreleases,*,err=998) zpoint1(numpoint)
324 if (old) call
skplin(2,unitreleases)
325 read(unitreleases,*,err=998) zpoint2(numpoint)
326 if (old) call
skplin(2,unitreleases)
327 read(unitreleases,*,err=998) npart(numpoint)
328 if (old) call
skplin(2,unitreleases)
330 read(unitreleases,*,err=998) xmass(numpoint,i)
331 if (old) call
skplin(2,unitreleases)
334 if (numpoint.le.1000)
then
335 read(unitreleases,
'(a40)',err=998) compoint(numpoint)(1:40)
337 read(unitreleases,
'(a40)',err=998) compoint(1001)(1:40)
339 if (old) call
skplin(1,unitreleases)
340 if (numpoint.le.1000)
then
341 if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
342 (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
343 (compoint(numpoint)(1:8).eq.
' ')) goto 250
345 if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
346 (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
353 if (npart(numpoint).eq.0)
then
354 write(*,*)
'FLEXPART MODEL ERROR'
355 write(*,*)
'RELEASES file is corrupt.'
356 write(*,*)
'At least for one release point, there are zero'
357 write(*,*)
'particles released. Make changes to RELEASES.'
364 if (xpoint1(numpoint).lt.xlon0) &
365 xpoint1(numpoint)=xpoint1(numpoint)+360.
366 if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
367 xpoint1(numpoint)=xpoint1(numpoint)-360.
368 if (xpoint2(numpoint).lt.xlon0) &
369 xpoint2(numpoint)=xpoint2(numpoint)+360.
370 if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
371 xpoint2(numpoint)=xpoint2(numpoint)-360.
378 if (jul1.gt.jul2)
then
379 write(*,*)
'FLEXPART MODEL ERROR'
380 write(*,*)
'Release stops before it begins.'
381 write(*,*)
'Make changes to file RELEASES.'
384 if (mdomainfill.eq.0)
then
385 if (ldirect.eq.1)
then
386 if ((jul1.lt.bdate).or.(jul2.gt.edate))
then
387 write(*,*)
'FLEXPART MODEL ERROR'
388 write(*,*)
'Release starts before simulation begins or ends'
389 write(*,*)
'after simulation stops.'
390 write(*,*)
'Make files COMMAND and RELEASES consistent.'
393 ireleasestart(numpoint)=int((jul1-bdate)*86400.)
394 ireleaseend(numpoint)=int((jul2-bdate)*86400.)
395 else if (ldirect.eq.-1)
then
396 if ((jul1.lt.edate).or.(jul2.gt.bdate))
then
397 write(*,*)
'FLEXPART MODEL ERROR'
398 write(*,*)
'Release starts before simulation begins or ends'
399 write(*,*)
'after simulation stops.'
400 write(*,*)
'Make files COMMAND and RELEASES consistent.'
403 ireleasestart(numpoint)=int((jul1-bdate)*86400.)
404 ireleaseend(numpoint)=int((jul2-bdate)*86400.)
417 if (ireleasestart(numpoint).ne.ireleaseend(numpoint))
then
418 releaserate=releaserate+
real(npart(numpoint))/ &
419 real(ireleaseend(numpoint)-ireleasestart(numpoint))
423 numpartmax=numpartmax+npart(numpoint)
427 250
close(unitreleases)
429 write (*,*)
' Particles allocated for this run: ',maxpart,
', released in simulation: ', numpartmax
432 if (ioutputforeachrelease.eq.1)
then
433 maxpointspec_act=numpoint
438 if (releaserate.gt. &
439 0.99*
real(maxpart)/
real(lage(nageclass))) then
440 if (numpartmax.gt.maxpart)
then
441 write(*,*)
'#####################################################'
442 write(*,*)
'#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
443 write(*,*)
'#### ####'
444 write(*,*)
'####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'
445 write(*,*)
'#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'
446 write(*,*)
'#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'
447 write(*,*)
'#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'
448 write(*,*)
'#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'
449 write(*,*)
'#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'
450 write(*,*)
'#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'
451 write(*,*)
'#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'
452 write(*,*)
'#####################################################'
453 write(*,*)
'Maximum release rate may be: ',releaserate, &
454 ' particles per second'
455 write(*,*)
'Maximum allowed release rate is: ', &
456 real(maxpart)/
real(lage(nageclass)),
' particles per second'
458 'Total number of particles released during the simulation is: ', &
460 write(*,*)
'Maximum allowed number of particles is: ',maxpart
466 994
write(*,*)
'#####################################################'
467 write(*,*)
'#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
468 write(*,*)
'#### ####'
469 write(*,*)
'#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
470 write(*,*)
'#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
471 write(*,*)
'#####################################################'
474 998
write(*,*)
'#####################################################'
475 write(*,*)
'#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'
476 write(*,*)
'#### ####'
477 write(*,*)
'#### FATAL ERROR - FILE "RELEASES" IS ####'
478 write(*,*)
'#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####'
479 write(*,*)
'#### MISTAKES OR GET A NEW "RELEASES"- ####'
480 write(*,*)
'#### FILE ... ####'
481 write(*,*)
'#####################################################'
485 999
write(*,*)
'#####################################################'
486 write(*,*)
' FLEXPART MODEL SUBROUTINE READRELEASES: '
488 write(*,*)
'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS'
489 write(*,*)
'POINTS IS NOT AVAILABLE OR YOU ARE NOT'
490 write(*,*)
'PERMITTED FOR ANY ACCESS'
491 write(*,*)
'#####################################################'
subroutine skplin(nlines, iunit)
subroutine readspecies(id_spec, pos_spec)
real(kind=dp) function juldate(yyyymmdd, hhmiss)
subroutine part0(dquer, dsigma, density, fract, schmi, cun, vsh)