| 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 | |
|---|
| 22 | subroutine readreleases |
|---|
| 23 | |
|---|
| 24 | !***************************************************************************** |
|---|
| 25 | ! * |
|---|
| 26 | ! This routine reads the release point specifications for the current * |
|---|
| 27 | ! model run. Several release points can be used at the same time. * |
|---|
| 28 | ! * |
|---|
| 29 | ! Author: A. Stohl * |
|---|
| 30 | ! * |
|---|
| 31 | ! 18 May 1996 * |
|---|
| 32 | ! * |
|---|
| 33 | ! Update: 29 January 2001 * |
|---|
| 34 | ! Release altitude can be either in magl or masl * |
|---|
| 35 | ! * |
|---|
| 36 | !***************************************************************************** |
|---|
| 37 | ! * |
|---|
| 38 | ! Variables: * |
|---|
| 39 | ! decay decay constant of species * |
|---|
| 40 | ! dquer [um] mean particle diameters * |
|---|
| 41 | ! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass* |
|---|
| 42 | ! are between 0.1*dquer and 10*dquer * |
|---|
| 43 | ! ireleasestart, ireleaseend [s] starting time and ending time of each * |
|---|
| 44 | ! release * |
|---|
| 45 | ! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint* |
|---|
| 46 | ! is in hPa * |
|---|
| 47 | ! npart number of particles to be released * |
|---|
| 48 | ! nspec number of species to be released * |
|---|
| 49 | ! density [kg/m3] density of the particles * |
|---|
| 50 | ! rm [s/m] Mesophyll resistance * |
|---|
| 51 | ! species name of species * |
|---|
| 52 | ! xmass total mass of each species * |
|---|
| 53 | ! xpoint1,ypoint1 geograf. coordinates of lower left corner of release * |
|---|
| 54 | ! area * |
|---|
| 55 | ! xpoint2,ypoint2 geograf. coordinates of upper right corner of release * |
|---|
| 56 | ! area * |
|---|
| 57 | ! weta, wetb parameters to determine the below-cloud scavenging * |
|---|
| 58 | ! weta_in, wetb_in parameters to determine the in-cloud scavenging * |
|---|
| 59 | ! wetc_in, wetd_in parameters to determine the in-cloud scavenging * |
|---|
| 60 | ! zpoint1,zpoint2 height range, over which release takes place * |
|---|
| 61 | ! num_min_discrete if less, release cannot be randomized and happens at * |
|---|
| 62 | ! time mid-point of release interval * |
|---|
| 63 | ! * |
|---|
| 64 | !***************************************************************************** |
|---|
| 65 | |
|---|
| 66 | use point_mod |
|---|
| 67 | use xmass_mod |
|---|
| 68 | use par_mod |
|---|
| 69 | use com_mod |
|---|
| 70 | |
|---|
| 71 | implicit none |
|---|
| 72 | |
|---|
| 73 | integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat |
|---|
| 74 | integer,parameter :: num_min_discrete=100 |
|---|
| 75 | real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun |
|---|
| 76 | real(kind=dp) :: jul1,jul2,julm,juldate |
|---|
| 77 | character(len=50) :: line |
|---|
| 78 | logical :: old |
|---|
| 79 | |
|---|
| 80 | !sec, read release to find how many releasepoints should be allocated |
|---|
| 81 | |
|---|
| 82 | open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', & |
|---|
| 83 | err=999) |
|---|
| 84 | |
|---|
| 85 | ! Check the format of the RELEASES file (either in free format, |
|---|
| 86 | ! or using a formatted mask) |
|---|
| 87 | ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION' |
|---|
| 88 | !************************************************************************** |
|---|
| 89 | |
|---|
| 90 | call skplin(12,unitreleases) |
|---|
| 91 | read (unitreleases,901) line |
|---|
| 92 | 901 format (a) |
|---|
| 93 | if (index(line,'Total') .eq. 0) then |
|---|
| 94 | old = .false. |
|---|
| 95 | else |
|---|
| 96 | old = .true. |
|---|
| 97 | endif |
|---|
| 98 | rewind(unitreleases) |
|---|
| 99 | |
|---|
| 100 | |
|---|
| 101 | ! Skip first 11 lines (file header) |
|---|
| 102 | !********************************** |
|---|
| 103 | |
|---|
| 104 | call skplin(11,unitreleases) |
|---|
| 105 | |
|---|
| 106 | |
|---|
| 107 | read(unitreleases,*,err=998) nspec |
|---|
| 108 | if (old) call skplin(2,unitreleases) |
|---|
| 109 | do i=1,nspec |
|---|
| 110 | read(unitreleases,*,err=998) specnum_rel |
|---|
| 111 | if (old) call skplin(2,unitreleases) |
|---|
| 112 | end do |
|---|
| 113 | |
|---|
| 114 | numpoint=0 |
|---|
| 115 | 100 numpoint=numpoint+1 |
|---|
| 116 | read(unitreleases,*,end=25) |
|---|
| 117 | read(unitreleases,*,err=998,end=25) idum,idum |
|---|
| 118 | if (old) call skplin(2,unitreleases) |
|---|
| 119 | read(unitreleases,*,err=998) idum,idum |
|---|
| 120 | if (old) call skplin(2,unitreleases) |
|---|
| 121 | read(unitreleases,*,err=998) xdum |
|---|
| 122 | if (old) call skplin(2,unitreleases) |
|---|
| 123 | read(unitreleases,*,err=998) xdum |
|---|
| 124 | if (old) call skplin(2,unitreleases) |
|---|
| 125 | read(unitreleases,*,err=998) xdum |
|---|
| 126 | if (old) call skplin(2,unitreleases) |
|---|
| 127 | read(unitreleases,*,err=998) xdum |
|---|
| 128 | if (old) call skplin(2,unitreleases) |
|---|
| 129 | read(unitreleases,*,err=998) idum |
|---|
| 130 | if (old) call skplin(2,unitreleases) |
|---|
| 131 | read(unitreleases,*,err=998) xdum |
|---|
| 132 | if (old) call skplin(2,unitreleases) |
|---|
| 133 | read(unitreleases,*,err=998) xdum |
|---|
| 134 | if (old) call skplin(2,unitreleases) |
|---|
| 135 | read(unitreleases,*,err=998) idum |
|---|
| 136 | if (old) call skplin(2,unitreleases) |
|---|
| 137 | do i=1,nspec |
|---|
| 138 | read(unitreleases,*,err=998) xdum |
|---|
| 139 | if (old) call skplin(2,unitreleases) |
|---|
| 140 | end do |
|---|
| 141 | !save compoint only for the first 1000 release points |
|---|
| 142 | read(unitreleases,'(a40)',err=998) compoint(1)(1:40) |
|---|
| 143 | if (old) call skplin(1,unitreleases) |
|---|
| 144 | |
|---|
| 145 | goto 100 |
|---|
| 146 | |
|---|
| 147 | 25 numpoint=numpoint-1 |
|---|
| 148 | |
|---|
| 149 | !allocate memory for numpoint releaspoint |
|---|
| 150 | allocate(ireleasestart(numpoint) & |
|---|
| 151 | ,stat=stat) |
|---|
| 152 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 153 | allocate(ireleaseend(numpoint) & |
|---|
| 154 | ,stat=stat) |
|---|
| 155 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 156 | allocate(xpoint1(numpoint) & |
|---|
| 157 | ,stat=stat) |
|---|
| 158 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 159 | allocate(xpoint2(numpoint) & |
|---|
| 160 | ,stat=stat) |
|---|
| 161 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 162 | allocate(ypoint1(numpoint) & |
|---|
| 163 | ,stat=stat) |
|---|
| 164 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 165 | allocate(ypoint2(numpoint) & |
|---|
| 166 | ,stat=stat) |
|---|
| 167 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 168 | allocate(zpoint1(numpoint) & |
|---|
| 169 | ,stat=stat) |
|---|
| 170 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 171 | allocate(zpoint2(numpoint) & |
|---|
| 172 | ,stat=stat) |
|---|
| 173 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 174 | allocate(kindz(numpoint) & |
|---|
| 175 | ,stat=stat) |
|---|
| 176 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 177 | allocate(xmass(numpoint,maxspec) & |
|---|
| 178 | ,stat=stat) |
|---|
| 179 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 180 | allocate(rho_rel(numpoint) & |
|---|
| 181 | ,stat=stat) |
|---|
| 182 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 183 | allocate(npart(numpoint) & |
|---|
| 184 | ,stat=stat) |
|---|
| 185 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 186 | allocate(xmasssave(numpoint) & |
|---|
| 187 | ,stat=stat) |
|---|
| 188 | if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT' |
|---|
| 189 | |
|---|
| 190 | write (*,*) ' Releasepoints allocated: ', numpoint |
|---|
| 191 | |
|---|
| 192 | do i=1,numpoint |
|---|
| 193 | xmasssave(i)=0. |
|---|
| 194 | end do |
|---|
| 195 | |
|---|
| 196 | !now save the information |
|---|
| 197 | DEP=.false. |
|---|
| 198 | DRYDEP=.false. |
|---|
| 199 | WETDEP=.false. |
|---|
| 200 | OHREA=.false. |
|---|
| 201 | do i=1,maxspec |
|---|
| 202 | DRYDEPSPEC(i)=.false. |
|---|
| 203 | end do |
|---|
| 204 | |
|---|
| 205 | rewind(unitreleases) |
|---|
| 206 | |
|---|
| 207 | |
|---|
| 208 | ! Skip first 11 lines (file header) |
|---|
| 209 | !********************************** |
|---|
| 210 | |
|---|
| 211 | call skplin(11,unitreleases) |
|---|
| 212 | |
|---|
| 213 | |
|---|
| 214 | ! Assign species-specific parameters needed for physical processes |
|---|
| 215 | !************************************************************************* |
|---|
| 216 | |
|---|
| 217 | read(unitreleases,*,err=998) nspec |
|---|
| 218 | if (nspec.gt.maxspec) goto 994 |
|---|
| 219 | if (old) call skplin(2,unitreleases) |
|---|
| 220 | do i=1,nspec |
|---|
| 221 | read(unitreleases,*,err=998) specnum_rel |
|---|
| 222 | if (old) call skplin(2,unitreleases) |
|---|
| 223 | call readspecies(specnum_rel,i) |
|---|
| 224 | |
|---|
| 225 | ! For backward runs, only 1 species is allowed |
|---|
| 226 | !********************************************* |
|---|
| 227 | |
|---|
| 228 | !if ((ldirect.lt.0).and.(nspec.gt.1)) then |
|---|
| 229 | !write(*,*) '#####################################################' |
|---|
| 230 | !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
|---|
| 231 | !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####' |
|---|
| 232 | !write(*,*) '#####################################################' |
|---|
| 233 | ! stop |
|---|
| 234 | !endif |
|---|
| 235 | |
|---|
| 236 | ! Molecular weight |
|---|
| 237 | !***************** |
|---|
| 238 | |
|---|
| 239 | if (((iout.eq.2).or.(iout.eq.3)).and. & |
|---|
| 240 | (weightmolar(i).lt.0.)) then |
|---|
| 241 | write(*,*) 'For mixing ratio output, valid molar weight' |
|---|
| 242 | write(*,*) 'must be specified for all simulated species.' |
|---|
| 243 | write(*,*) 'Check table SPECIES or choose concentration' |
|---|
| 244 | write(*,*) 'output instead if molar weight is not known.' |
|---|
| 245 | stop |
|---|
| 246 | endif |
|---|
| 247 | |
|---|
| 248 | |
|---|
| 249 | ! Radioactive decay |
|---|
| 250 | !****************** |
|---|
| 251 | |
|---|
| 252 | decay(i)=0.693147/decay(i) !conversion half life to decay constant |
|---|
| 253 | |
|---|
| 254 | |
|---|
| 255 | ! Dry deposition of gases |
|---|
| 256 | !************************ |
|---|
| 257 | |
|---|
| 258 | if (reldiff(i).gt.0.) & |
|---|
| 259 | rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance |
|---|
| 260 | |
|---|
| 261 | ! Dry deposition of particles |
|---|
| 262 | !**************************** |
|---|
| 263 | |
|---|
| 264 | vsetaver(i)=0. |
|---|
| 265 | cunningham(i)=0. |
|---|
| 266 | dquer(i)=dquer(i)*1000000. ! Conversion m to um |
|---|
| 267 | if (density(i).gt.0.) then ! Additional parameters |
|---|
| 268 | call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh) |
|---|
| 269 | do j=1,ni |
|---|
| 270 | fract(i,j)=fracth(j) |
|---|
| 271 | schmi(i,j)=schmih(j) |
|---|
| 272 | vset(i,j)=vsh(j) |
|---|
| 273 | cunningham(i)=cunningham(i)+cun*fract(i,j) |
|---|
| 274 | vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j) |
|---|
| 275 | end do |
|---|
| 276 | write(*,*) 'Average settling velocity: ',i,vsetaver(i) |
|---|
| 277 | endif |
|---|
| 278 | |
|---|
| 279 | ! Dry deposition for constant deposition velocity |
|---|
| 280 | !************************************************ |
|---|
| 281 | |
|---|
| 282 | dryvel(i)=dryvel(i)*0.01 ! conversion to m/s |
|---|
| 283 | |
|---|
| 284 | ! Check if wet deposition or OH reaction shall be calculated |
|---|
| 285 | !*********************************************************** |
|---|
| 286 | if (weta(i).gt.0.) then |
|---|
| 287 | WETDEP=.true. |
|---|
| 288 | write (*,*) 'Below-cloud scavenging is switched on' |
|---|
| 289 | write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i |
|---|
| 290 | else |
|---|
| 291 | write (*,*) 'Below-cloud scavenging is switched OFF' |
|---|
| 292 | endif |
|---|
| 293 | |
|---|
| 294 | ! NIK 31.01.2013 + 10.12.2013 |
|---|
| 295 | if (weta_in(i).gt.0.) then |
|---|
| 296 | WETDEP=.true. |
|---|
| 297 | write (*,*) 'In-cloud scavenging is switched on' |
|---|
| 298 | write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i |
|---|
| 299 | else |
|---|
| 300 | write (*,*) 'In-cloud scavenging is switched OFF' |
|---|
| 301 | endif |
|---|
| 302 | |
|---|
| 303 | if (ohreact(i).gt.0) then |
|---|
| 304 | OHREA=.true. |
|---|
| 305 | write (*,*) 'OHreaction switched on: ',ohreact(i),i |
|---|
| 306 | endif |
|---|
| 307 | |
|---|
| 308 | |
|---|
| 309 | if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. & |
|---|
| 310 | (dryvel(i).gt.0.)) then |
|---|
| 311 | DRYDEP=.true. |
|---|
| 312 | DRYDEPSPEC(i)=.true. |
|---|
| 313 | endif |
|---|
| 314 | |
|---|
| 315 | end do |
|---|
| 316 | |
|---|
| 317 | if (WETDEP.or.DRYDEP) DEP=.true. |
|---|
| 318 | |
|---|
| 319 | ! Read specifications for each release point |
|---|
| 320 | !******************************************* |
|---|
| 321 | |
|---|
| 322 | numpoint=0 |
|---|
| 323 | numpartmax=0 |
|---|
| 324 | releaserate=0. |
|---|
| 325 | 1000 numpoint=numpoint+1 |
|---|
| 326 | read(unitreleases,*,end=250) |
|---|
| 327 | read(unitreleases,*,err=998,end=250) id1,it1 |
|---|
| 328 | if (old) call skplin(2,unitreleases) |
|---|
| 329 | read(unitreleases,*,err=998) id2,it2 |
|---|
| 330 | if (old) call skplin(2,unitreleases) |
|---|
| 331 | read(unitreleases,*,err=998) xpoint1(numpoint) |
|---|
| 332 | if (old) call skplin(2,unitreleases) |
|---|
| 333 | read(unitreleases,*,err=998) ypoint1(numpoint) |
|---|
| 334 | if (old) call skplin(2,unitreleases) |
|---|
| 335 | read(unitreleases,*,err=998) xpoint2(numpoint) |
|---|
| 336 | if (old) call skplin(2,unitreleases) |
|---|
| 337 | read(unitreleases,*,err=998) ypoint2(numpoint) |
|---|
| 338 | if (old) call skplin(2,unitreleases) |
|---|
| 339 | read(unitreleases,*,err=998) kindz(numpoint) |
|---|
| 340 | if (old) call skplin(2,unitreleases) |
|---|
| 341 | read(unitreleases,*,err=998) zpoint1(numpoint) |
|---|
| 342 | if (old) call skplin(2,unitreleases) |
|---|
| 343 | read(unitreleases,*,err=998) zpoint2(numpoint) |
|---|
| 344 | if (old) call skplin(2,unitreleases) |
|---|
| 345 | read(unitreleases,*,err=998) npart(numpoint) |
|---|
| 346 | if (old) call skplin(2,unitreleases) |
|---|
| 347 | do i=1,nspec |
|---|
| 348 | read(unitreleases,*,err=998) xmass(numpoint,i) |
|---|
| 349 | if (old) call skplin(2,unitreleases) |
|---|
| 350 | end do |
|---|
| 351 | !save compoint only for the first 1000 release points |
|---|
| 352 | if (numpoint.le.1000) then |
|---|
| 353 | read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40) |
|---|
| 354 | else |
|---|
| 355 | read(unitreleases,'(a40)',err=998) compoint(1001)(1:40) |
|---|
| 356 | endif |
|---|
| 357 | if (old) call skplin(1,unitreleases) |
|---|
| 358 | if (numpoint.le.1000) then |
|---|
| 359 | if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & |
|---|
| 360 | (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. & |
|---|
| 361 | (compoint(numpoint)(1:8).eq.' ')) goto 250 |
|---|
| 362 | else |
|---|
| 363 | if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & |
|---|
| 364 | (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250 |
|---|
| 365 | endif |
|---|
| 366 | |
|---|
| 367 | |
|---|
| 368 | ! If a release point contains no particles, stop and issue error message |
|---|
| 369 | !*********************************************************************** |
|---|
| 370 | |
|---|
| 371 | if (npart(numpoint).eq.0) then |
|---|
| 372 | write(*,*) 'FLEXPART MODEL ERROR' |
|---|
| 373 | write(*,*) 'RELEASES file is corrupt.' |
|---|
| 374 | write(*,*) 'At least for one release point, there are zero' |
|---|
| 375 | write(*,*) 'particles released. Make changes to RELEASES.' |
|---|
| 376 | stop |
|---|
| 377 | endif |
|---|
| 378 | |
|---|
| 379 | ! Check whether x coordinates of release point are within model domain |
|---|
| 380 | !********************************************************************* |
|---|
| 381 | |
|---|
| 382 | if (xpoint1(numpoint).lt.xlon0) & |
|---|
| 383 | xpoint1(numpoint)=xpoint1(numpoint)+360. |
|---|
| 384 | if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) & |
|---|
| 385 | xpoint1(numpoint)=xpoint1(numpoint)-360. |
|---|
| 386 | if (xpoint2(numpoint).lt.xlon0) & |
|---|
| 387 | xpoint2(numpoint)=xpoint2(numpoint)+360. |
|---|
| 388 | if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) & |
|---|
| 389 | xpoint2(numpoint)=xpoint2(numpoint)-360. |
|---|
| 390 | |
|---|
| 391 | ! Determine relative beginning and ending times of particle release |
|---|
| 392 | !****************************************************************** |
|---|
| 393 | |
|---|
| 394 | jul1=juldate(id1,it1) |
|---|
| 395 | jul2=juldate(id2,it2) |
|---|
| 396 | julm=(jul1+jul2)/2. |
|---|
| 397 | if (jul1.gt.jul2) then |
|---|
| 398 | write(*,*) 'FLEXPART MODEL ERROR' |
|---|
| 399 | write(*,*) 'Release stops before it begins.' |
|---|
| 400 | write(*,*) 'Make changes to file RELEASES.' |
|---|
| 401 | stop |
|---|
| 402 | endif |
|---|
| 403 | if (mdomainfill.eq.0) then ! no domain filling |
|---|
| 404 | if (ldirect.eq.1) then |
|---|
| 405 | if ((jul1.lt.bdate).or.(jul2.gt.edate)) then |
|---|
| 406 | write(*,*) 'FLEXPART MODEL ERROR' |
|---|
| 407 | write(*,*) 'Release starts before simulation begins or ends' |
|---|
| 408 | write(*,*) 'after simulation stops.' |
|---|
| 409 | write(*,*) 'Make files COMMAND and RELEASES consistent.' |
|---|
| 410 | stop |
|---|
| 411 | endif |
|---|
| 412 | if (npart(numpoint).gt.num_min_discrete) then |
|---|
| 413 | ireleasestart(numpoint)=int((jul1-bdate)*86400.) |
|---|
| 414 | ireleaseend(numpoint)=int((jul2-bdate)*86400.) |
|---|
| 415 | else |
|---|
| 416 | ireleasestart(numpoint)=int((julm-bdate)*86400.) |
|---|
| 417 | ireleaseend(numpoint)=int((julm-bdate)*86400.) |
|---|
| 418 | endif |
|---|
| 419 | else if (ldirect.eq.-1) then |
|---|
| 420 | if ((jul1.lt.edate).or.(jul2.gt.bdate)) then |
|---|
| 421 | write(*,*) 'FLEXPART MODEL ERROR' |
|---|
| 422 | write(*,*) 'Release starts before simulation begins or ends' |
|---|
| 423 | write(*,*) 'after simulation stops.' |
|---|
| 424 | write(*,*) 'Make files COMMAND and RELEASES consistent.' |
|---|
| 425 | stop |
|---|
| 426 | endif |
|---|
| 427 | if (npart(numpoint).gt.num_min_discrete) then |
|---|
| 428 | ireleasestart(numpoint)=int((jul1-bdate)*86400.) |
|---|
| 429 | ireleaseend(numpoint)=int((jul2-bdate)*86400.) |
|---|
| 430 | else |
|---|
| 431 | ireleasestart(numpoint)=int((julm-bdate)*86400.) |
|---|
| 432 | ireleaseend(numpoint)=int((julm-bdate)*86400.) |
|---|
| 433 | endif |
|---|
| 434 | endif |
|---|
| 435 | endif |
|---|
| 436 | |
|---|
| 437 | |
|---|
| 438 | ! Check, whether the total number of particles may exceed totally allowed |
|---|
| 439 | ! number of particles at some time during the simulation |
|---|
| 440 | !************************************************************************ |
|---|
| 441 | |
|---|
| 442 | ! Determine the release rate (particles per second) and total number |
|---|
| 443 | ! of particles released during the simulation |
|---|
| 444 | !******************************************************************* |
|---|
| 445 | |
|---|
| 446 | if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then |
|---|
| 447 | releaserate=releaserate+real(npart(numpoint))/ & |
|---|
| 448 | real(ireleaseend(numpoint)-ireleasestart(numpoint)) |
|---|
| 449 | else |
|---|
| 450 | releaserate=99999999 |
|---|
| 451 | endif |
|---|
| 452 | numpartmax=numpartmax+npart(numpoint) |
|---|
| 453 | goto 1000 |
|---|
| 454 | |
|---|
| 455 | |
|---|
| 456 | 250 close(unitreleases) |
|---|
| 457 | |
|---|
| 458 | write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ', numpartmax |
|---|
| 459 | numpoint=numpoint-1 |
|---|
| 460 | |
|---|
| 461 | if (ioutputforeachrelease.eq.1) then |
|---|
| 462 | maxpointspec_act=numpoint |
|---|
| 463 | else |
|---|
| 464 | maxpointspec_act=1 |
|---|
| 465 | endif |
|---|
| 466 | |
|---|
| 467 | if (releaserate.gt. & |
|---|
| 468 | 0.99*real(maxpart)/real(lage(nageclass))) then |
|---|
| 469 | if (numpartmax.gt.maxpart) then |
|---|
| 470 | write(*,*) '#####################################################' |
|---|
| 471 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
|---|
| 472 | write(*,*) '#### ####' |
|---|
| 473 | write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####' |
|---|
| 474 | write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####' |
|---|
| 475 | write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####' |
|---|
| 476 | write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####' |
|---|
| 477 | write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####' |
|---|
| 478 | write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####' |
|---|
| 479 | write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####' |
|---|
| 480 | write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####' |
|---|
| 481 | write(*,*) '#####################################################' |
|---|
| 482 | write(*,*) 'Maximum release rate may be: ',releaserate, & |
|---|
| 483 | ' particles per second' |
|---|
| 484 | write(*,*) 'Maximum allowed release rate is: ', & |
|---|
| 485 | real(maxpart)/real(lage(nageclass)),' particles per second' |
|---|
| 486 | write(*,*) & |
|---|
| 487 | 'Total number of particles released during the simulation is: ', & |
|---|
| 488 | numpartmax |
|---|
| 489 | write(*,*) 'Maximum allowed number of particles is: ',maxpart |
|---|
| 490 | endif |
|---|
| 491 | endif |
|---|
| 492 | |
|---|
| 493 | return |
|---|
| 494 | |
|---|
| 495 | 994 write(*,*) '#####################################################' |
|---|
| 496 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
|---|
| 497 | write(*,*) '#### ####' |
|---|
| 498 | write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####' |
|---|
| 499 | write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####' |
|---|
| 500 | write(*,*) '#####################################################' |
|---|
| 501 | stop |
|---|
| 502 | |
|---|
| 503 | 998 write(*,*) '#####################################################' |
|---|
| 504 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
|---|
| 505 | write(*,*) '#### ####' |
|---|
| 506 | write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS ####' |
|---|
| 507 | write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####' |
|---|
| 508 | write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"- ####' |
|---|
| 509 | write(*,*) '#### FILE ... ####' |
|---|
| 510 | write(*,*) '#####################################################' |
|---|
| 511 | stop |
|---|
| 512 | |
|---|
| 513 | |
|---|
| 514 | 999 write(*,*) '#####################################################' |
|---|
| 515 | write(*,*) ' FLEXPART MODEL SUBROUTINE READRELEASES: ' |
|---|
| 516 | write(*,*) |
|---|
| 517 | write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS' |
|---|
| 518 | write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT' |
|---|
| 519 | write(*,*) 'PERMITTED FOR ANY ACCESS' |
|---|
| 520 | write(*,*) '#####################################################' |
|---|
| 521 | stop |
|---|
| 522 | |
|---|
| 523 | end subroutine readreleases |
|---|
