[e200b7a] | 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 * |
---|
[b4d29ce] | 35 | ! HSO, 12 August 2013 |
---|
| 36 | ! Added optional namelist input |
---|
[e200b7a] | 37 | ! * |
---|
| 38 | !***************************************************************************** |
---|
| 39 | ! * |
---|
| 40 | ! Variables: * |
---|
| 41 | ! decay decay constant of species * |
---|
| 42 | ! dquer [um] mean particle diameters * |
---|
| 43 | ! dsigma e.g. dsigma=10 or dsigma=0.1 means that 68% of the mass* |
---|
| 44 | ! are between 0.1*dquer and 10*dquer * |
---|
| 45 | ! ireleasestart, ireleaseend [s] starting time and ending time of each * |
---|
| 46 | ! release * |
---|
| 47 | ! kindz 1: zpoint is in m agl, 2: zpoint is in m asl, 3: zpoint* |
---|
| 48 | ! is in hPa * |
---|
| 49 | ! npart number of particles to be released * |
---|
| 50 | ! nspec number of species to be released * |
---|
| 51 | ! density [kg/m3] density of the particles * |
---|
| 52 | ! rm [s/m] Mesophyll resistance * |
---|
| 53 | ! species name of species * |
---|
| 54 | ! xmass total mass of each species * |
---|
| 55 | ! xpoint1,ypoint1 geograf. coordinates of lower left corner of release * |
---|
| 56 | ! area * |
---|
| 57 | ! xpoint2,ypoint2 geograf. coordinates of upper right corner of release * |
---|
| 58 | ! area * |
---|
[4fbe7a5] | 59 | ! weta, wetb parameters to determine the below-cloud scavenging * |
---|
| 60 | ! weta_in, wetb_in parameters to determine the in-cloud scavenging * |
---|
| 61 | ! wetc_in, wetd_in parameters to determine the in-cloud scavenging * |
---|
[e200b7a] | 62 | ! zpoint1,zpoint2 height range, over which release takes place * |
---|
[f13406c] | 63 | ! num_min_discrete if less, release cannot be randomized and happens at * |
---|
| 64 | ! time mid-point of release interval * |
---|
[e200b7a] | 65 | ! * |
---|
| 66 | !***************************************************************************** |
---|
| 67 | |
---|
| 68 | use point_mod |
---|
| 69 | use xmass_mod |
---|
| 70 | use par_mod |
---|
| 71 | use com_mod |
---|
| 72 | |
---|
| 73 | implicit none |
---|
| 74 | |
---|
[b4d29ce] | 75 | integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat |
---|
[f13406c] | 76 | integer,parameter :: num_min_discrete=100 |
---|
[e200b7a] | 77 | real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun |
---|
[f13406c] | 78 | real(kind=dp) :: jul1,jul2,julm,juldate |
---|
[e200b7a] | 79 | character(len=50) :: line |
---|
| 80 | logical :: old |
---|
| 81 | |
---|
[b4d29ce] | 82 | ! help variables for namelist reading |
---|
| 83 | integer :: numpoints, parts, readerror |
---|
| 84 | integer*2 :: zkind |
---|
| 85 | integer :: idate1, itime1, idate2, itime2 |
---|
| 86 | real :: lon1,lon2,lat1,lat2,z1,z2 |
---|
| 87 | character*40 :: comment |
---|
| 88 | integer,parameter :: unitreleasesout=2 |
---|
| 89 | real,allocatable, dimension (:) :: mass |
---|
| 90 | integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2 |
---|
| 91 | |
---|
| 92 | ! declare namelists |
---|
| 93 | namelist /releases_ctrl/ & |
---|
| 94 | nspec, & |
---|
| 95 | specnum_rel |
---|
| 96 | |
---|
| 97 | namelist /release/ & |
---|
| 98 | idate1, itime1, & |
---|
| 99 | idate2, itime2, & |
---|
| 100 | lon1, lon2, & |
---|
| 101 | lat1, lat2, & |
---|
| 102 | z1, z2, & |
---|
| 103 | zkind, & |
---|
| 104 | mass, & |
---|
| 105 | parts, & |
---|
| 106 | comment |
---|
[e200b7a] | 107 | |
---|
[b4d29ce] | 108 | numpoint=0 |
---|
[e200b7a] | 109 | |
---|
[b4d29ce] | 110 | ! allocate with maxspec for first input loop |
---|
| 111 | allocate(mass(maxspec),stat=stat) |
---|
| 112 | if (stat.ne.0) write(*,*)'ERROR: could not allocate mass' |
---|
| 113 | allocate(specnum_rel(maxspec),stat=stat) |
---|
| 114 | if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel' |
---|
[e200b7a] | 115 | |
---|
[b4d29ce] | 116 | ! presetting namelist releases_ctrl |
---|
| 117 | nspec = -1 ! use negative value to determine failed namelist input |
---|
| 118 | specnum_rel = 0 |
---|
| 119 | |
---|
| 120 | !sec, read release to find how many releasepoints should be allocated |
---|
| 121 | open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999) |
---|
| 122 | |
---|
| 123 | ! check if namelist input provided |
---|
| 124 | read(unitreleases,releases_ctrl,iostat=readerror) |
---|
| 125 | |
---|
| 126 | ! prepare namelist output if requested |
---|
| 127 | if (nmlout.eqv..true.) then |
---|
[242571d] | 128 | !open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000) |
---|
| 129 | open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',err=1000) |
---|
[e200b7a] | 130 | endif |
---|
| 131 | |
---|
[b4d29ce] | 132 | if ((readerror.ne.0).or.(nspec.lt.0)) then |
---|
[e200b7a] | 133 | |
---|
[b4d29ce] | 134 | ! no namelist format, close file and allow reopening in old format |
---|
| 135 | close(unitreleases) |
---|
| 136 | open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999) |
---|
[e200b7a] | 137 | |
---|
[b4d29ce] | 138 | readerror=1 ! indicates old format |
---|
[e200b7a] | 139 | |
---|
[b4d29ce] | 140 | ! Check the format of the RELEASES file (either in free format, |
---|
| 141 | ! or using a formatted mask) |
---|
| 142 | ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION' |
---|
| 143 | !************************************************************************** |
---|
[e200b7a] | 144 | |
---|
[b4d29ce] | 145 | call skplin(12,unitreleases) |
---|
| 146 | read (unitreleases,901) line |
---|
| 147 | 901 format (a) |
---|
| 148 | if (index(line,'Total') .eq. 0) then |
---|
| 149 | old = .false. |
---|
| 150 | else |
---|
| 151 | old = .true. |
---|
| 152 | endif |
---|
| 153 | rewind(unitreleases) |
---|
[e200b7a] | 154 | |
---|
[b4d29ce] | 155 | |
---|
| 156 | ! Skip first 11 lines (file header) |
---|
| 157 | !********************************** |
---|
| 158 | |
---|
| 159 | call skplin(11,unitreleases) |
---|
| 160 | |
---|
| 161 | read(unitreleases,*,err=998) nspec |
---|
| 162 | if (old) call skplin(2,unitreleases) |
---|
| 163 | do i=1,nspec |
---|
| 164 | read(unitreleases,*,err=998) specnum_rel(i) |
---|
| 165 | if (old) call skplin(2,unitreleases) |
---|
| 166 | end do |
---|
| 167 | |
---|
| 168 | numpoint=0 |
---|
| 169 | 100 numpoint=numpoint+1 |
---|
| 170 | read(unitreleases,*,end=25) |
---|
| 171 | read(unitreleases,*,err=998,end=25) idum,idum |
---|
| 172 | if (old) call skplin(2,unitreleases) |
---|
| 173 | read(unitreleases,*,err=998) idum,idum |
---|
| 174 | if (old) call skplin(2,unitreleases) |
---|
| 175 | read(unitreleases,*,err=998) xdum |
---|
| 176 | if (old) call skplin(2,unitreleases) |
---|
| 177 | read(unitreleases,*,err=998) xdum |
---|
| 178 | if (old) call skplin(2,unitreleases) |
---|
[e200b7a] | 179 | read(unitreleases,*,err=998) xdum |
---|
| 180 | if (old) call skplin(2,unitreleases) |
---|
[b4d29ce] | 181 | read(unitreleases,*,err=998) xdum |
---|
| 182 | if (old) call skplin(2,unitreleases) |
---|
| 183 | read(unitreleases,*,err=998) idum |
---|
| 184 | if (old) call skplin(2,unitreleases) |
---|
| 185 | read(unitreleases,*,err=998) xdum |
---|
| 186 | if (old) call skplin(2,unitreleases) |
---|
| 187 | read(unitreleases,*,err=998) xdum |
---|
| 188 | if (old) call skplin(2,unitreleases) |
---|
| 189 | read(unitreleases,*,err=998) idum |
---|
| 190 | if (old) call skplin(2,unitreleases) |
---|
| 191 | do i=1,nspec |
---|
| 192 | read(unitreleases,*,err=998) xdum |
---|
| 193 | if (old) call skplin(2,unitreleases) |
---|
| 194 | end do |
---|
| 195 | !save compoint only for the first 1000 release points |
---|
| 196 | read(unitreleases,'(a40)',err=998) compoint(1)(1:40) |
---|
| 197 | if (old) call skplin(1,unitreleases) |
---|
| 198 | |
---|
| 199 | goto 100 |
---|
| 200 | |
---|
| 201 | 25 numpoint=numpoint-1 |
---|
| 202 | |
---|
| 203 | else |
---|
| 204 | |
---|
| 205 | readerror=0 |
---|
| 206 | do while (readerror.eq.0) |
---|
| 207 | idate1=-1 |
---|
| 208 | read(unitreleases,release,iostat=readerror) |
---|
| 209 | if ((idate1.lt.0).or.(readerror.ne.0)) then |
---|
| 210 | readerror=1 |
---|
| 211 | else |
---|
| 212 | numpoint=numpoint+1 |
---|
| 213 | endif |
---|
| 214 | end do |
---|
| 215 | readerror=0 |
---|
| 216 | endif ! if namelist input |
---|
| 217 | |
---|
| 218 | rewind(unitreleases) |
---|
| 219 | |
---|
| 220 | ! allocate arrays of matching size for number of species (namelist output) |
---|
| 221 | deallocate(mass) |
---|
| 222 | allocate(mass(nspec),stat=stat) |
---|
| 223 | if (stat.ne.0) write(*,*)'ERROR: could not allocate mass' |
---|
| 224 | allocate(specnum_rel2(nspec),stat=stat) |
---|
| 225 | if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2' |
---|
| 226 | specnum_rel2=specnum_rel(1:nspec) |
---|
| 227 | deallocate(specnum_rel) |
---|
| 228 | allocate(specnum_rel(nspec),stat=stat) |
---|
| 229 | if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel' |
---|
| 230 | specnum_rel=specnum_rel2 |
---|
| 231 | deallocate(specnum_rel2) |
---|
| 232 | |
---|
| 233 | !allocate memory for numpoint releaspoints |
---|
| 234 | allocate(ireleasestart(numpoint),stat=stat) |
---|
| 235 | if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart' |
---|
| 236 | allocate(ireleaseend(numpoint),stat=stat) |
---|
| 237 | if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend' |
---|
| 238 | allocate(xpoint1(numpoint),stat=stat) |
---|
| 239 | if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1' |
---|
| 240 | allocate(xpoint2(numpoint),stat=stat) |
---|
| 241 | if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2' |
---|
| 242 | allocate(ypoint1(numpoint),stat=stat) |
---|
| 243 | if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1' |
---|
| 244 | allocate(ypoint2(numpoint),stat=stat) |
---|
| 245 | if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2' |
---|
| 246 | allocate(zpoint1(numpoint),stat=stat) |
---|
| 247 | if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1' |
---|
| 248 | allocate(zpoint2(numpoint),stat=stat) |
---|
| 249 | if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2' |
---|
| 250 | allocate(kindz(numpoint),stat=stat) |
---|
| 251 | if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz' |
---|
| 252 | allocate(xmass(numpoint,maxspec),stat=stat) |
---|
| 253 | if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass' |
---|
| 254 | allocate(rho_rel(numpoint),stat=stat) |
---|
| 255 | if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel' |
---|
| 256 | allocate(npart(numpoint),stat=stat) |
---|
| 257 | if (stat.ne.0) write(*,*)'ERROR: could not allocate npart' |
---|
| 258 | allocate(xmasssave(numpoint),stat=stat) |
---|
| 259 | if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave' |
---|
| 260 | |
---|
[b7ae015] | 261 | if (verbosity.gt.0) then |
---|
| 262 | write (*,*) 'readreleases> Releasepoints : ', numpoint |
---|
| 263 | endif |
---|
[b4d29ce] | 264 | |
---|
| 265 | do i=1,numpoint |
---|
| 266 | xmasssave(i)=0. |
---|
[e200b7a] | 267 | end do |
---|
| 268 | |
---|
| 269 | !now save the information |
---|
| 270 | DEP=.false. |
---|
| 271 | DRYDEP=.false. |
---|
| 272 | WETDEP=.false. |
---|
| 273 | OHREA=.false. |
---|
| 274 | do i=1,maxspec |
---|
| 275 | DRYDEPSPEC(i)=.false. |
---|
| 276 | end do |
---|
| 277 | |
---|
[b4d29ce] | 278 | if (readerror.ne.0) then |
---|
| 279 | ! Skip first 11 lines (file header) |
---|
| 280 | !********************************** |
---|
[e200b7a] | 281 | |
---|
[b4d29ce] | 282 | call skplin(11,unitreleases) |
---|
[e200b7a] | 283 | |
---|
[b4d29ce] | 284 | ! Assign species-specific parameters needed for physical processes |
---|
| 285 | !************************************************************************* |
---|
[e200b7a] | 286 | |
---|
[b4d29ce] | 287 | read(unitreleases,*,err=998) nspec |
---|
| 288 | if (nspec.gt.maxspec) goto 994 |
---|
| 289 | if (old) call skplin(2,unitreleases) |
---|
| 290 | endif |
---|
[e200b7a] | 291 | |
---|
[b4d29ce] | 292 | ! namelist output |
---|
| 293 | if (nmlout.eqv..true.) then |
---|
| 294 | write(unitreleasesout,nml=releases_ctrl) |
---|
| 295 | endif |
---|
[e200b7a] | 296 | |
---|
| 297 | do i=1,nspec |
---|
[242571d] | 298 | if (verbosity.gt.0) then |
---|
| 299 | print*, 'readreleases> call readspecies', i |
---|
| 300 | endif |
---|
| 301 | |
---|
[b4d29ce] | 302 | if (readerror.ne.0) then |
---|
| 303 | read(unitreleases,*,err=998) specnum_rel(i) |
---|
| 304 | if (old) call skplin(2,unitreleases) |
---|
| 305 | call readspecies(specnum_rel(i),i) |
---|
| 306 | else |
---|
| 307 | call readspecies(specnum_rel(i),i) |
---|
| 308 | endif |
---|
[e200b7a] | 309 | |
---|
[b4d29ce] | 310 | ! For backward runs, only 1 species is allowed |
---|
| 311 | !********************************************* |
---|
[e200b7a] | 312 | |
---|
[b4d29ce] | 313 | !if ((ldirect.lt.0).and.(nspec.gt.1)) then |
---|
| 314 | !write(*,*) '#####################################################' |
---|
| 315 | !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
---|
| 316 | !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####' |
---|
| 317 | !write(*,*) '#####################################################' |
---|
| 318 | ! stop |
---|
| 319 | !endif |
---|
[e200b7a] | 320 | |
---|
[b4d29ce] | 321 | ! Molecular weight |
---|
| 322 | !***************** |
---|
[e200b7a] | 323 | |
---|
[b4d29ce] | 324 | if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then |
---|
[e200b7a] | 325 | write(*,*) 'For mixing ratio output, valid molar weight' |
---|
| 326 | write(*,*) 'must be specified for all simulated species.' |
---|
| 327 | write(*,*) 'Check table SPECIES or choose concentration' |
---|
| 328 | write(*,*) 'output instead if molar weight is not known.' |
---|
| 329 | stop |
---|
| 330 | endif |
---|
| 331 | |
---|
[b4d29ce] | 332 | ! Radioactive decay |
---|
| 333 | !****************** |
---|
[e200b7a] | 334 | |
---|
| 335 | decay(i)=0.693147/decay(i) !conversion half life to decay constant |
---|
| 336 | |
---|
| 337 | |
---|
[b4d29ce] | 338 | ! Dry deposition of gases |
---|
| 339 | !************************ |
---|
[e200b7a] | 340 | |
---|
[b4d29ce] | 341 | if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance |
---|
[e200b7a] | 342 | |
---|
[b4d29ce] | 343 | ! Dry deposition of particles |
---|
| 344 | !**************************** |
---|
[e200b7a] | 345 | |
---|
| 346 | vsetaver(i)=0. |
---|
| 347 | cunningham(i)=0. |
---|
| 348 | dquer(i)=dquer(i)*1000000. ! Conversion m to um |
---|
[b4d29ce] | 349 | if (density(i).gt.0.) then ! Additional parameters |
---|
| 350 | call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh) |
---|
[e200b7a] | 351 | do j=1,ni |
---|
| 352 | fract(i,j)=fracth(j) |
---|
| 353 | schmi(i,j)=schmih(j) |
---|
| 354 | vset(i,j)=vsh(j) |
---|
| 355 | cunningham(i)=cunningham(i)+cun*fract(i,j) |
---|
| 356 | vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j) |
---|
| 357 | end do |
---|
[4fbe7a5] | 358 | write(*,*) 'Average settling velocity: ',i,vsetaver(i) |
---|
[e200b7a] | 359 | endif |
---|
| 360 | |
---|
[b4d29ce] | 361 | ! Dry deposition for constant deposition velocity |
---|
| 362 | !************************************************ |
---|
[e200b7a] | 363 | |
---|
| 364 | dryvel(i)=dryvel(i)*0.01 ! conversion to m/s |
---|
| 365 | |
---|
[b4d29ce] | 366 | ! Check if wet deposition or OH reaction shall be calculated |
---|
| 367 | !*********************************************************** |
---|
[e200b7a] | 368 | if (weta(i).gt.0.) then |
---|
| 369 | WETDEP=.true. |
---|
[b4d29ce] | 370 | write (*,*) 'Below-cloud scavenging: ON' |
---|
[b7ae015] | 371 | if (verbosity.gt.0) then |
---|
[4fbe7a5] | 372 | write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i |
---|
[b7ae015] | 373 | endif |
---|
[4fbe7a5] | 374 | else |
---|
[b7ae015] | 375 | if (verbosity.gt.0) then |
---|
[b4d29ce] | 376 | write (*,*) 'Below-cloud scavenging: OFF' |
---|
[b7ae015] | 377 | endif |
---|
[e200b7a] | 378 | endif |
---|
[4fbe7a5] | 379 | |
---|
[b4d29ce] | 380 | ! NIK 31.01.2013 + 10.12.2013 |
---|
[4fbe7a5] | 381 | if (weta_in(i).gt.0.) then |
---|
| 382 | WETDEP=.true. |
---|
[b4d29ce] | 383 | write (*,*) 'In-cloud scavenging: ON' |
---|
[b7ae015] | 384 | if (verbosity.gt.0) then |
---|
[4fbe7a5] | 385 | write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i |
---|
[b7ae015] | 386 | endif |
---|
| 387 | else |
---|
| 388 | if (verbosity.gt.0) then |
---|
| 389 | write (*,*) 'In-cloud scavenging: OFF' |
---|
| 390 | endif |
---|
[4fbe7a5] | 391 | endif |
---|
| 392 | |
---|
[e200b7a] | 393 | if (ohreact(i).gt.0) then |
---|
| 394 | OHREA=.true. |
---|
[b4d29ce] | 395 | write (*,*) 'OHreaction: ON (',ohreact(i),i,')' |
---|
[e200b7a] | 396 | endif |
---|
| 397 | |
---|
[b4d29ce] | 398 | if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then |
---|
[e200b7a] | 399 | DRYDEP=.true. |
---|
| 400 | DRYDEPSPEC(i)=.true. |
---|
| 401 | endif |
---|
| 402 | |
---|
| 403 | end do |
---|
| 404 | |
---|
[b4d29ce] | 405 | if (WETDEP.or.DRYDEP) DEP=.true. |
---|
[e200b7a] | 406 | |
---|
| 407 | ! Read specifications for each release point |
---|
| 408 | !******************************************* |
---|
[b4d29ce] | 409 | numpoints=numpoint |
---|
[e200b7a] | 410 | numpoint=0 |
---|
| 411 | numpartmax=0 |
---|
| 412 | releaserate=0. |
---|
[b4d29ce] | 413 | 101 numpoint=numpoint+1 |
---|
| 414 | |
---|
| 415 | if (readerror.lt.1) then ! reading namelist format |
---|
| 416 | |
---|
| 417 | if (numpoint.gt.numpoints) goto 250 |
---|
| 418 | zkind = 1 |
---|
| 419 | mass = 0 |
---|
| 420 | parts = 0 |
---|
| 421 | comment = ' ' |
---|
| 422 | read(unitreleases,release,iostat=readerror) |
---|
| 423 | id1=idate1 |
---|
| 424 | it1=itime1 |
---|
| 425 | id2=idate2 |
---|
| 426 | it2=itime2 |
---|
| 427 | xpoint1(numpoint)=lon1 |
---|
| 428 | xpoint2(numpoint)=lon2 |
---|
| 429 | ypoint1(numpoint)=lat1 |
---|
| 430 | ypoint2(numpoint)=lat2 |
---|
| 431 | zpoint1(numpoint)=z1 |
---|
| 432 | zpoint2(numpoint)=z2 |
---|
| 433 | kindz(numpoint)=zkind |
---|
| 434 | do i=1,nspec |
---|
| 435 | xmass(numpoint,i)=mass(i) |
---|
| 436 | end do |
---|
| 437 | npart(numpoint)=parts |
---|
| 438 | compoint(min(1001,numpoint))=comment |
---|
| 439 | |
---|
| 440 | ! namelist output |
---|
| 441 | if (nmlout.eqv..true.) then |
---|
| 442 | write(unitreleasesout,nml=release) |
---|
| 443 | endif |
---|
| 444 | |
---|
[e200b7a] | 445 | else |
---|
| 446 | |
---|
[b4d29ce] | 447 | read(unitreleases,*,end=250) |
---|
| 448 | read(unitreleases,*,err=998,end=250) id1,it1 |
---|
| 449 | if (old) call skplin(2,unitreleases) |
---|
| 450 | read(unitreleases,*,err=998) id2,it2 |
---|
| 451 | if (old) call skplin(2,unitreleases) |
---|
| 452 | read(unitreleases,*,err=998) xpoint1(numpoint) |
---|
| 453 | if (old) call skplin(2,unitreleases) |
---|
| 454 | read(unitreleases,*,err=998) ypoint1(numpoint) |
---|
| 455 | if (old) call skplin(2,unitreleases) |
---|
| 456 | read(unitreleases,*,err=998) xpoint2(numpoint) |
---|
| 457 | if (old) call skplin(2,unitreleases) |
---|
| 458 | read(unitreleases,*,err=998) ypoint2(numpoint) |
---|
| 459 | if (old) call skplin(2,unitreleases) |
---|
| 460 | read(unitreleases,*,err=998) kindz(numpoint) |
---|
| 461 | if (old) call skplin(2,unitreleases) |
---|
| 462 | read(unitreleases,*,err=998) zpoint1(numpoint) |
---|
| 463 | if (old) call skplin(2,unitreleases) |
---|
| 464 | read(unitreleases,*,err=998) zpoint2(numpoint) |
---|
| 465 | if (old) call skplin(2,unitreleases) |
---|
| 466 | read(unitreleases,*,err=998) npart(numpoint) |
---|
| 467 | if (old) call skplin(2,unitreleases) |
---|
| 468 | do i=1,nspec |
---|
| 469 | read(unitreleases,*,err=998) xmass(numpoint,i) |
---|
| 470 | if (old) call skplin(2,unitreleases) |
---|
| 471 | mass(i)=xmass(numpoint,i) |
---|
| 472 | end do |
---|
| 473 | !save compoint only for the first 1000 release points |
---|
| 474 | if (numpoint.le.1000) then |
---|
| 475 | read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40) |
---|
| 476 | comment=compoint(numpoint)(1:40) |
---|
| 477 | else |
---|
| 478 | read(unitreleases,'(a40)',err=998) compoint(1001)(1:40) |
---|
| 479 | comment=compoint(1001)(1:40) |
---|
| 480 | endif |
---|
| 481 | if (old) call skplin(1,unitreleases) |
---|
| 482 | |
---|
| 483 | ! namelist output |
---|
| 484 | if (nmlout.eqv..true.) then |
---|
| 485 | idate1=id1 |
---|
| 486 | itime1=it1 |
---|
| 487 | idate2=id2 |
---|
| 488 | itime2=it2 |
---|
| 489 | lon1=xpoint1(numpoint) |
---|
| 490 | lon2=xpoint2(numpoint) |
---|
| 491 | lat1=ypoint1(numpoint) |
---|
| 492 | lat2=ypoint2(numpoint) |
---|
| 493 | z1=zpoint1(numpoint) |
---|
| 494 | z2=zpoint2(numpoint) |
---|
| 495 | zkind=kindz(numpoint) |
---|
| 496 | parts=npart(numpoint) |
---|
| 497 | write(unitreleasesout,nml=release) |
---|
| 498 | endif |
---|
| 499 | |
---|
| 500 | if (numpoint.le.1000) then |
---|
| 501 | if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & |
---|
| 502 | (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. & |
---|
| 503 | (compoint(numpoint)(1:8).eq.' ')) goto 250 |
---|
| 504 | else |
---|
| 505 | if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & |
---|
| 506 | (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250 |
---|
| 507 | endif |
---|
| 508 | |
---|
| 509 | endif ! if namelist format |
---|
[e200b7a] | 510 | |
---|
[242571d] | 511 | |
---|
| 512 | if (verbosity.gt.1 .and. numpoint.eq.1) then ! verbosity 2 or larger |
---|
| 513 | write(*,*) 'numpoint=', numpoint |
---|
| 514 | print*, id1,it1 |
---|
| 515 | print*, id2,it2 |
---|
| 516 | print*, xpoint1(numpoint) |
---|
| 517 | print*, ypoint1(numpoint) |
---|
| 518 | print*, xpoint2(numpoint) |
---|
| 519 | print*, ypoint2(numpoint) |
---|
| 520 | print*, 'kindz=' , kindz(numpoint) |
---|
| 521 | print*, zpoint1(numpoint) |
---|
| 522 | print*, zpoint2(numpoint) |
---|
| 523 | print*, npart(numpoint) |
---|
| 524 | do i=1,nspec |
---|
| 525 | !mass(i)= |
---|
| 526 | print*, 'xmass=', xmass(numpoint,i) |
---|
| 527 | end do |
---|
| 528 | print*, compoint(numpoint) |
---|
| 529 | endif |
---|
| 530 | |
---|
| 531 | |
---|
[e200b7a] | 532 | ! If a release point contains no particles, stop and issue error message |
---|
| 533 | !*********************************************************************** |
---|
| 534 | |
---|
| 535 | if (npart(numpoint).eq.0) then |
---|
| 536 | write(*,*) 'FLEXPART MODEL ERROR' |
---|
| 537 | write(*,*) 'RELEASES file is corrupt.' |
---|
| 538 | write(*,*) 'At least for one release point, there are zero' |
---|
| 539 | write(*,*) 'particles released. Make changes to RELEASES.' |
---|
| 540 | stop |
---|
| 541 | endif |
---|
| 542 | |
---|
| 543 | ! Check whether x coordinates of release point are within model domain |
---|
| 544 | !********************************************************************* |
---|
| 545 | |
---|
| 546 | if (xpoint1(numpoint).lt.xlon0) & |
---|
| 547 | xpoint1(numpoint)=xpoint1(numpoint)+360. |
---|
| 548 | if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) & |
---|
| 549 | xpoint1(numpoint)=xpoint1(numpoint)-360. |
---|
| 550 | if (xpoint2(numpoint).lt.xlon0) & |
---|
| 551 | xpoint2(numpoint)=xpoint2(numpoint)+360. |
---|
| 552 | if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) & |
---|
| 553 | xpoint2(numpoint)=xpoint2(numpoint)-360. |
---|
| 554 | |
---|
| 555 | ! Determine relative beginning and ending times of particle release |
---|
| 556 | !****************************************************************** |
---|
| 557 | |
---|
| 558 | jul1=juldate(id1,it1) |
---|
| 559 | jul2=juldate(id2,it2) |
---|
[f13406c] | 560 | julm=(jul1+jul2)/2. |
---|
[e200b7a] | 561 | if (jul1.gt.jul2) then |
---|
| 562 | write(*,*) 'FLEXPART MODEL ERROR' |
---|
| 563 | write(*,*) 'Release stops before it begins.' |
---|
| 564 | write(*,*) 'Make changes to file RELEASES.' |
---|
| 565 | stop |
---|
| 566 | endif |
---|
| 567 | if (mdomainfill.eq.0) then ! no domain filling |
---|
| 568 | if (ldirect.eq.1) then |
---|
| 569 | if ((jul1.lt.bdate).or.(jul2.gt.edate)) then |
---|
| 570 | write(*,*) 'FLEXPART MODEL ERROR' |
---|
| 571 | write(*,*) 'Release starts before simulation begins or ends' |
---|
| 572 | write(*,*) 'after simulation stops.' |
---|
| 573 | write(*,*) 'Make files COMMAND and RELEASES consistent.' |
---|
[242571d] | 574 | write(*,*) jul1, ' < ' , bdate |
---|
| 575 | write(*,*) ' .or. ' |
---|
| 576 | write(*,*) jul2 , ' > ', edate |
---|
| 577 | |
---|
[e200b7a] | 578 | stop |
---|
| 579 | endif |
---|
[f13406c] | 580 | if (npart(numpoint).gt.num_min_discrete) then |
---|
| 581 | ireleasestart(numpoint)=int((jul1-bdate)*86400.) |
---|
| 582 | ireleaseend(numpoint)=int((jul2-bdate)*86400.) |
---|
| 583 | else |
---|
| 584 | ireleasestart(numpoint)=int((julm-bdate)*86400.) |
---|
| 585 | ireleaseend(numpoint)=int((julm-bdate)*86400.) |
---|
| 586 | endif |
---|
[e200b7a] | 587 | else if (ldirect.eq.-1) then |
---|
| 588 | if ((jul1.lt.edate).or.(jul2.gt.bdate)) then |
---|
| 589 | write(*,*) 'FLEXPART MODEL ERROR' |
---|
| 590 | write(*,*) 'Release starts before simulation begins or ends' |
---|
| 591 | write(*,*) 'after simulation stops.' |
---|
| 592 | write(*,*) 'Make files COMMAND and RELEASES consistent.' |
---|
| 593 | stop |
---|
| 594 | endif |
---|
[f13406c] | 595 | if (npart(numpoint).gt.num_min_discrete) then |
---|
| 596 | ireleasestart(numpoint)=int((jul1-bdate)*86400.) |
---|
| 597 | ireleaseend(numpoint)=int((jul2-bdate)*86400.) |
---|
| 598 | else |
---|
| 599 | ireleasestart(numpoint)=int((julm-bdate)*86400.) |
---|
| 600 | ireleaseend(numpoint)=int((julm-bdate)*86400.) |
---|
| 601 | endif |
---|
[e200b7a] | 602 | endif |
---|
| 603 | endif |
---|
| 604 | |
---|
[242571d] | 605 | if (verbosity.gt.1 .and. numpoint.eq.1) then ! verbosity 2 or larger |
---|
| 606 | print*, 'ireleasestart(',numpoint,')', ireleasestart(numpoint) |
---|
| 607 | print*, 'ireleaseend(',numpoint,')', ireleaseend(numpoint) |
---|
| 608 | endif |
---|
| 609 | |
---|
[e200b7a] | 610 | ! Determine the release rate (particles per second) and total number |
---|
| 611 | ! of particles released during the simulation |
---|
| 612 | !******************************************************************* |
---|
| 613 | |
---|
| 614 | if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then |
---|
| 615 | releaserate=releaserate+real(npart(numpoint))/ & |
---|
| 616 | real(ireleaseend(numpoint)-ireleasestart(numpoint)) |
---|
| 617 | else |
---|
| 618 | releaserate=99999999 |
---|
| 619 | endif |
---|
| 620 | numpartmax=numpartmax+npart(numpoint) |
---|
[b4d29ce] | 621 | goto 101 |
---|
[e200b7a] | 622 | |
---|
| 623 | 250 close(unitreleases) |
---|
| 624 | |
---|
[b4d29ce] | 625 | if (nmlout.eqv..true.) then |
---|
| 626 | close(unitreleasesout) |
---|
| 627 | endif |
---|
| 628 | |
---|
| 629 | write (*,*) 'Particles allocated (maxpart) : ',maxpart |
---|
| 630 | write (*,*) 'Particles released (numpartmax): ',numpartmax |
---|
[e200b7a] | 631 | numpoint=numpoint-1 |
---|
| 632 | |
---|
| 633 | if (ioutputforeachrelease.eq.1) then |
---|
| 634 | maxpointspec_act=numpoint |
---|
| 635 | else |
---|
| 636 | maxpointspec_act=1 |
---|
| 637 | endif |
---|
| 638 | |
---|
[b4d29ce] | 639 | ! Check, whether the total number of particles may exceed totally allowed |
---|
| 640 | ! number of particles at some time during the simulation |
---|
| 641 | !************************************************************************ |
---|
| 642 | |
---|
[e200b7a] | 643 | if (releaserate.gt. & |
---|
| 644 | 0.99*real(maxpart)/real(lage(nageclass))) then |
---|
| 645 | if (numpartmax.gt.maxpart) then |
---|
| 646 | write(*,*) '#####################################################' |
---|
| 647 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
---|
| 648 | write(*,*) '#### ####' |
---|
| 649 | write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####' |
---|
| 650 | write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####' |
---|
| 651 | write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####' |
---|
| 652 | write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####' |
---|
| 653 | write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####' |
---|
| 654 | write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####' |
---|
| 655 | write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####' |
---|
| 656 | write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####' |
---|
| 657 | write(*,*) '#####################################################' |
---|
| 658 | write(*,*) 'Maximum release rate may be: ',releaserate, & |
---|
| 659 | ' particles per second' |
---|
| 660 | write(*,*) 'Maximum allowed release rate is: ', & |
---|
| 661 | real(maxpart)/real(lage(nageclass)),' particles per second' |
---|
| 662 | write(*,*) & |
---|
| 663 | 'Total number of particles released during the simulation is: ', & |
---|
| 664 | numpartmax |
---|
| 665 | write(*,*) 'Maximum allowed number of particles is: ',maxpart |
---|
| 666 | endif |
---|
| 667 | endif |
---|
| 668 | |
---|
| 669 | return |
---|
| 670 | |
---|
| 671 | 994 write(*,*) '#####################################################' |
---|
| 672 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
---|
| 673 | write(*,*) '#### ####' |
---|
| 674 | write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####' |
---|
| 675 | write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####' |
---|
| 676 | write(*,*) '#####################################################' |
---|
| 677 | stop |
---|
| 678 | |
---|
| 679 | 998 write(*,*) '#####################################################' |
---|
| 680 | write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' |
---|
| 681 | write(*,*) '#### ####' |
---|
| 682 | write(*,*) '#### FATAL ERROR - FILE "RELEASES" IS ####' |
---|
| 683 | write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR ####' |
---|
| 684 | write(*,*) '#### MISTAKES OR GET A NEW "RELEASES"- ####' |
---|
| 685 | write(*,*) '#### FILE ... ####' |
---|
| 686 | write(*,*) '#####################################################' |
---|
| 687 | stop |
---|
| 688 | |
---|
| 689 | |
---|
| 690 | 999 write(*,*) '#####################################################' |
---|
| 691 | write(*,*) ' FLEXPART MODEL SUBROUTINE READRELEASES: ' |
---|
| 692 | write(*,*) |
---|
| 693 | write(*,*) 'FATAL ERROR - FILE CONTAINING PARTICLE RELEASE POINTS' |
---|
| 694 | write(*,*) 'POINTS IS NOT AVAILABLE OR YOU ARE NOT' |
---|
| 695 | write(*,*) 'PERMITTED FOR ANY ACCESS' |
---|
| 696 | write(*,*) '#####################################################' |
---|
| 697 | stop |
---|
| 698 | |
---|
[b4d29ce] | 699 | 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES" #### ' |
---|
| 700 | write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' |
---|
| 701 | write(*,'(a)') path(2)(1:length(2)) |
---|
| 702 | stop |
---|
| 703 | |
---|
[e200b7a] | 704 | end subroutine readreleases |
---|