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