Changeset b4d29ce in flexpart.git for src/readreleases.f90
- Timestamp:
- Jul 3, 2014, 2:55:50 PM (10 years ago)
- Branches:
- master, 10.4.1_pesei, FPv9.3.1, FPv9.3.1b_testing, FPv9.3.2, GFS_025, NetCDF, bugfixes+enhancements, deposition, dev, fp9.3.1-20161214-nc4, grib2nc4_repair, inputlist, laptop, release-10, release-10.4.1, scaling-bug, svn-petra, svn-trunk, univie
- Children:
- 4b093cc, 75dfd42, 326b31b
- Parents:
- 87910af
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/readreleases.f90
r4fbe7a5 rb4d29ce 33 33 ! Update: 29 January 2001 * 34 34 ! Release altitude can be either in magl or masl * 35 ! HSO, 12 August 2013 36 ! Added optional namelist input 35 37 ! * 36 38 !***************************************************************************** … … 71 73 implicit none 72 74 73 integer :: numpartmax,i,j,id1,it1,id2,it2, specnum_rel,idum,stat75 integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat 74 76 integer,parameter :: num_min_discrete=100 75 77 real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun … … 78 80 logical :: old 79 81 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 107 108 numpoint=0 109 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' 115 116 ! presetting namelist releases_ctrl 117 nspec = -1 ! use negative value to determine failed namelist input 118 specnum_rel = 0 119 80 120 !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 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 128 open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000) 129 endif 130 131 if ((readerror.ne.0).or.(nspec.lt.0)) then 132 133 ! no namelist format, close file and allow reopening in old format 134 close(unitreleases) 135 open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999) 136 137 readerror=1 ! indicates old format 138 139 ! Check the format of the RELEASES file (either in free format, 140 ! or using a formatted mask) 141 ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION' 142 !************************************************************************** 143 144 call skplin(12,unitreleases) 145 read (unitreleases,901) line 146 901 format (a) 147 if (index(line,'Total') .eq. 0) then 148 old = .false. 149 else 150 old = .true. 151 endif 152 rewind(unitreleases) 153 154 155 ! Skip first 11 lines (file header) 156 !********************************** 157 158 call skplin(11,unitreleases) 159 160 read(unitreleases,*,err=998) nspec 161 if (old) call skplin(2,unitreleases) 162 do i=1,nspec 163 read(unitreleases,*,err=998) specnum_rel(i) 164 if (old) call skplin(2,unitreleases) 165 end do 166 167 numpoint=0 168 100 numpoint=numpoint+1 169 read(unitreleases,*,end=25) 170 read(unitreleases,*,err=998,end=25) idum,idum 171 if (old) call skplin(2,unitreleases) 172 read(unitreleases,*,err=998) idum,idum 173 if (old) call skplin(2,unitreleases) 174 read(unitreleases,*,err=998) xdum 175 if (old) call skplin(2,unitreleases) 176 read(unitreleases,*,err=998) xdum 177 if (old) call skplin(2,unitreleases) 178 read(unitreleases,*,err=998) xdum 179 if (old) call skplin(2,unitreleases) 180 read(unitreleases,*,err=998) xdum 181 if (old) call skplin(2,unitreleases) 182 read(unitreleases,*,err=998) idum 183 if (old) call skplin(2,unitreleases) 184 read(unitreleases,*,err=998) xdum 185 if (old) call skplin(2,unitreleases) 186 read(unitreleases,*,err=998) xdum 187 if (old) call skplin(2,unitreleases) 188 read(unitreleases,*,err=998) idum 189 if (old) call skplin(2,unitreleases) 190 do i=1,nspec 191 read(unitreleases,*,err=998) xdum 192 if (old) call skplin(2,unitreleases) 193 end do 194 !save compoint only for the first 1000 release points 195 read(unitreleases,'(a40)',err=998) compoint(1)(1:40) 196 if (old) call skplin(1,unitreleases) 197 198 goto 100 199 200 25 numpoint=numpoint-1 201 202 else 203 204 readerror=0 205 do while (readerror.eq.0) 206 idate1=-1 207 read(unitreleases,release,iostat=readerror) 208 if ((idate1.lt.0).or.(readerror.ne.0)) then 209 readerror=1 210 else 211 numpoint=numpoint+1 212 endif 213 end do 214 readerror=0 215 endif ! if namelist input 216 98 217 rewind(unitreleases) 99 218 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) 219 ! allocate arrays of matching size for number of species (namelist output) 220 deallocate(mass) 221 allocate(mass(nspec),stat=stat) 222 if (stat.ne.0) write(*,*)'ERROR: could not allocate mass' 223 allocate(specnum_rel2(nspec),stat=stat) 224 if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2' 225 specnum_rel2=specnum_rel(1:nspec) 226 deallocate(specnum_rel) 227 allocate(specnum_rel(nspec),stat=stat) 228 if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel' 229 specnum_rel=specnum_rel2 230 deallocate(specnum_rel2) 231 232 !allocate memory for numpoint releaspoints 233 allocate(ireleasestart(numpoint),stat=stat) 234 if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart' 235 allocate(ireleaseend(numpoint),stat=stat) 236 if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend' 237 allocate(xpoint1(numpoint),stat=stat) 238 if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1' 239 allocate(xpoint2(numpoint),stat=stat) 240 if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2' 241 allocate(ypoint1(numpoint),stat=stat) 242 if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1' 243 allocate(ypoint2(numpoint),stat=stat) 244 if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2' 245 allocate(zpoint1(numpoint),stat=stat) 246 if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1' 247 allocate(zpoint2(numpoint),stat=stat) 248 if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2' 249 allocate(kindz(numpoint),stat=stat) 250 if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz' 251 allocate(xmass(numpoint,maxspec),stat=stat) 252 if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass' 253 allocate(rho_rel(numpoint),stat=stat) 254 if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel' 255 allocate(npart(numpoint),stat=stat) 256 if (stat.ne.0) write(*,*)'ERROR: could not allocate npart' 257 allocate(xmasssave(numpoint),stat=stat) 258 if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave' 259 260 write (*,*) 'Releasepoints : ', numpoint 261 262 do i=1,numpoint 263 xmasssave(i)=0. 112 264 end do 113 114 numpoint=0115 100 numpoint=numpoint+1116 read(unitreleases,*,end=25)117 read(unitreleases,*,err=998,end=25) idum,idum118 if (old) call skplin(2,unitreleases)119 read(unitreleases,*,err=998) idum,idum120 if (old) call skplin(2,unitreleases)121 read(unitreleases,*,err=998) xdum122 if (old) call skplin(2,unitreleases)123 read(unitreleases,*,err=998) xdum124 if (old) call skplin(2,unitreleases)125 read(unitreleases,*,err=998) xdum126 if (old) call skplin(2,unitreleases)127 read(unitreleases,*,err=998) xdum128 if (old) call skplin(2,unitreleases)129 read(unitreleases,*,err=998) idum130 if (old) call skplin(2,unitreleases)131 read(unitreleases,*,err=998) xdum132 if (old) call skplin(2,unitreleases)133 read(unitreleases,*,err=998) xdum134 if (old) call skplin(2,unitreleases)135 read(unitreleases,*,err=998) idum136 if (old) call skplin(2,unitreleases)137 do i=1,nspec138 read(unitreleases,*,err=998) xdum139 if (old) call skplin(2,unitreleases)140 end do141 !save compoint only for the first 1000 release points142 read(unitreleases,'(a40)',err=998) compoint(1)(1:40)143 if (old) call skplin(1,unitreleases)144 145 goto 100146 147 25 numpoint=numpoint-1148 149 !allocate memory for numpoint releaspoint150 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: ', numpoint191 192 do i=1,numpoint193 xmasssave(i)=0.194 end do195 265 196 266 !now save the information … … 203 273 end do 204 274 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) 275 if (readerror.ne.0) then 276 ! Skip first 11 lines (file header) 277 !********************************** 278 279 call skplin(11,unitreleases) 280 281 ! Assign species-specific parameters needed for physical processes 282 !************************************************************************* 283 284 read(unitreleases,*,err=998) nspec 285 if (nspec.gt.maxspec) goto 994 286 if (old) call skplin(2,unitreleases) 287 endif 288 289 ! namelist output 290 if (nmlout.eqv..true.) then 291 write(unitreleasesout,nml=releases_ctrl) 292 endif 293 220 294 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 295 if (readerror.ne.0) then 296 read(unitreleases,*,err=998) specnum_rel(i) 297 if (old) call skplin(2,unitreleases) 298 call readspecies(specnum_rel(i),i) 299 else 300 call readspecies(specnum_rel(i),i) 301 endif 302 303 ! For backward runs, only 1 species is allowed 304 !********************************************* 305 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 313 314 ! Molecular weight 315 !***************** 316 317 if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then 241 318 write(*,*) 'For mixing ratio output, valid molar weight' 242 319 write(*,*) 'must be specified for all simulated species.' … … 246 323 endif 247 324 248 249 ! Radioactive decay 250 !****************** 325 ! Radioactive decay 326 !****************** 251 327 252 328 decay(i)=0.693147/decay(i) !conversion half life to decay constant 253 329 254 330 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 !**************************** 331 ! Dry deposition of gases 332 !************************ 333 334 if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance 335 336 ! Dry deposition of particles 337 !**************************** 263 338 264 339 vsetaver(i)=0. 265 340 cunningham(i)=0. 266 341 dquer(i)=dquer(i)*1000000. ! Conversion m to um 267 if (density(i).gt.0.) then 268 call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)342 if (density(i).gt.0.) then ! Additional parameters 343 call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh) 269 344 do j=1,ni 270 345 fract(i,j)=fracth(j) … … 277 352 endif 278 353 279 ! Dry deposition for constant deposition velocity280 !************************************************354 ! Dry deposition for constant deposition velocity 355 !************************************************ 281 356 282 357 dryvel(i)=dryvel(i)*0.01 ! conversion to m/s 283 358 284 ! Check if wet deposition or OH reaction shall be calculated285 !***********************************************************359 ! Check if wet deposition or OH reaction shall be calculated 360 !*********************************************************** 286 361 if (weta(i).gt.0.) then 287 362 WETDEP=.true. 288 write (*,*) 'Below-cloud scavenging is switched on'363 write (*,*) 'Below-cloud scavenging: ON' 289 364 write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i 290 365 else 291 write (*,*) 'Below-cloud scavenging is switchedOFF'366 write (*,*) 'Below-cloud scavenging: OFF' 292 367 endif 293 368 294 ! NIK 31.01.2013 + 10.12.2013369 ! NIK 31.01.2013 + 10.12.2013 295 370 if (weta_in(i).gt.0.) then 296 371 WETDEP=.true. 297 write (*,*) 'In-cloud scavenging is switched on'372 write (*,*) 'In-cloud scavenging: ON' 298 373 write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i 299 374 else 300 write (*,*) 'In-cloud scavenging is switchedOFF'375 write (*,*) 'In-cloud scavenging: OFF' 301 376 endif 302 377 303 378 if (ohreact(i).gt.0) then 304 379 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 380 write (*,*) 'OHreaction: ON (',ohreact(i),i,')' 381 endif 382 383 if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then 311 384 DRYDEP=.true. 312 385 DRYDEPSPEC(i)=.true. … … 315 388 end do 316 389 317 390 if (WETDEP.or.DRYDEP) DEP=.true. 318 391 319 392 ! Read specifications for each release point 320 393 !******************************************* 321 394 numpoints=numpoint 322 395 numpoint=0 323 396 numpartmax=0 324 397 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) 398 101 numpoint=numpoint+1 399 400 if (readerror.lt.1) then ! reading namelist format 401 402 if (numpoint.gt.numpoints) goto 250 403 zkind = 1 404 mass = 0 405 parts = 0 406 comment = ' ' 407 read(unitreleases,release,iostat=readerror) 408 id1=idate1 409 it1=itime1 410 id2=idate2 411 it2=itime2 412 xpoint1(numpoint)=lon1 413 xpoint2(numpoint)=lon2 414 ypoint1(numpoint)=lat1 415 ypoint2(numpoint)=lat2 416 zpoint1(numpoint)=z1 417 zpoint2(numpoint)=z2 418 kindz(numpoint)=zkind 419 do i=1,nspec 420 xmass(numpoint,i)=mass(i) 421 end do 422 npart(numpoint)=parts 423 compoint(min(1001,numpoint))=comment 424 425 ! namelist output 426 if (nmlout.eqv..true.) then 427 write(unitreleasesout,nml=release) 428 endif 429 354 430 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 431 432 read(unitreleases,*,end=250) 433 read(unitreleases,*,err=998,end=250) id1,it1 434 if (old) call skplin(2,unitreleases) 435 read(unitreleases,*,err=998) id2,it2 436 if (old) call skplin(2,unitreleases) 437 read(unitreleases,*,err=998) xpoint1(numpoint) 438 if (old) call skplin(2,unitreleases) 439 read(unitreleases,*,err=998) ypoint1(numpoint) 440 if (old) call skplin(2,unitreleases) 441 read(unitreleases,*,err=998) xpoint2(numpoint) 442 if (old) call skplin(2,unitreleases) 443 read(unitreleases,*,err=998) ypoint2(numpoint) 444 if (old) call skplin(2,unitreleases) 445 read(unitreleases,*,err=998) kindz(numpoint) 446 if (old) call skplin(2,unitreleases) 447 read(unitreleases,*,err=998) zpoint1(numpoint) 448 if (old) call skplin(2,unitreleases) 449 read(unitreleases,*,err=998) zpoint2(numpoint) 450 if (old) call skplin(2,unitreleases) 451 read(unitreleases,*,err=998) npart(numpoint) 452 if (old) call skplin(2,unitreleases) 453 do i=1,nspec 454 read(unitreleases,*,err=998) xmass(numpoint,i) 455 if (old) call skplin(2,unitreleases) 456 mass(i)=xmass(numpoint,i) 457 end do 458 !save compoint only for the first 1000 release points 459 if (numpoint.le.1000) then 460 read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40) 461 comment=compoint(numpoint)(1:40) 462 else 463 read(unitreleases,'(a40)',err=998) compoint(1001)(1:40) 464 comment=compoint(1001)(1:40) 465 endif 466 if (old) call skplin(1,unitreleases) 467 468 ! namelist output 469 if (nmlout.eqv..true.) then 470 idate1=id1 471 itime1=it1 472 idate2=id2 473 itime2=it2 474 lon1=xpoint1(numpoint) 475 lon2=xpoint2(numpoint) 476 lat1=ypoint1(numpoint) 477 lat2=ypoint2(numpoint) 478 z1=zpoint1(numpoint) 479 z2=zpoint2(numpoint) 480 zkind=kindz(numpoint) 481 parts=npart(numpoint) 482 write(unitreleasesout,nml=release) 483 endif 484 485 if (numpoint.le.1000) then 486 if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & 487 (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. & 488 (compoint(numpoint)(1:8).eq.' ')) goto 250 489 else 490 if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. & 491 (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250 492 endif 493 494 endif ! if namelist format 367 495 368 496 ! If a release point contains no particles, stop and issue error message … … 435 563 endif 436 564 437 438 ! Check, whether the total number of particles may exceed totally allowed439 ! number of particles at some time during the simulation440 !************************************************************************441 442 565 ! Determine the release rate (particles per second) and total number 443 566 ! of particles released during the simulation … … 451 574 endif 452 575 numpartmax=numpartmax+npart(numpoint) 453 goto 1000 454 576 goto 101 455 577 456 578 250 close(unitreleases) 457 579 458 write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ', numpartmax 580 if (nmlout.eqv..true.) then 581 close(unitreleasesout) 582 endif 583 584 write (*,*) 'Particles allocated (maxpart) : ',maxpart 585 write (*,*) 'Particles released (numpartmax): ',numpartmax 459 586 numpoint=numpoint-1 460 587 … … 464 591 maxpointspec_act=1 465 592 endif 593 594 ! Check, whether the total number of particles may exceed totally allowed 595 ! number of particles at some time during the simulation 596 !************************************************************************ 466 597 467 598 if (releaserate.gt. & … … 521 652 stop 522 653 654 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES" #### ' 655 write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY #### ' 656 write(*,'(a)') path(2)(1:length(2)) 657 stop 658 523 659 end subroutine readreleases
Note: See TracChangeset
for help on using the changeset viewer.