Changeset 5184a7c in flexpart.git
- Timestamp:
- Jun 20, 2017, 9:09:35 AM (7 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- b5127f9
- Parents:
- a76d954
- Location:
- src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
src/FLEXPART_MPI.f90
rc8fc724 r5184a7c 80 80 ! FLEXPART version string 81 81 flexversion_major = '10' ! Major version number, also used for species file names 82 ! flexversion='Ver. 10 Beta MPI (2015-05-01)'83 82 flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)' 84 83 verbosity=0 -
src/advance.f90
re52967c r5184a7c 535 535 536 536 if (mdomainfill.eq.0) then 537 ! ESO 05.2015 Changed this to fix MQUASILAG option, where nrelpoint is 538 ! particle number and thus xmass array goes out of bounds 539 ! do nsp=1,nspec 540 ! if (xmass(nrelpoint,nsp).gt.eps2) goto 887 541 ! end do 542 ! 887 nsp=min(nsp,nspec) 543 if (nspec.eq.1.and.density(1).gt.0.) then 544 call get_settling(itime,real(xt),real(yt),zt,nspec,settling) !bugfix 537 if (lsettling) then 538 do nsp=1,nspec 539 if (xmass(nrelpoint,nsp).gt.eps2) exit 540 end do 541 if (nsp.gt.nspec) then 542 ! This should never happen 543 write(*,*) 'advance.f90: ERROR: could not find releasepoint' 544 stop 545 end if 546 if (density(nsp).gt.0.) then 547 call get_settling(itime,real(xt),real(yt),zt,nsp,settling) !bugfix 548 w=w+settling 549 end if 545 550 end if 546 w=w+settling547 551 endif 548 552 … … 700 704 !************************************************************************* 701 705 702 703 704 if (mdomainfill.eq.0) then 705 ! ESO 05.2015 Changed this to fix MQUASILAG option, where nrelpoint is 706 ! particle number and thus xmass array goes out of bounds 707 708 ! do nsp=1,nspec 709 ! if (xmass(nrelpoint,nsp).gt.eps2) goto 888 710 ! end do 711 ! 888 nsp=min(nsp,nspec) 712 ! if (density(nsp).gt.0.) then 713 if (nspec.eq.1.and.density(1).gt.0.) then 714 call get_settling(itime,real(xt),real(yt),zt,nspec,settling) !bugfix 706 if (mdomainfill.eq.0) then 707 if (lsettling) then 708 do nsp=1,nspec 709 if (xmass(nrelpoint,nsp).gt.eps2) exit 710 end do 711 if (nsp.gt.nspec) then 712 ! This should never happen 713 write(*,*) 'advance.f90: ERROR: could not find releasepoint' 714 stop 715 715 end if 716 w=w+settling 716 if (density(nsp).gt.0.) then 717 call get_settling(itime,real(xt),real(yt),zt,nsp,settling) !bugfix 718 w=w+settling 719 end if 717 720 endif 721 end if 722 718 723 719 724 ! Calculate position at time step itime+lsynctime … … 910 915 911 916 if (mdomainfill.eq.0) then 912 ! ESO 05.2015 Changed this to fix MQUASILAG option, where nrelpoint is 913 ! particle number and thus xmass array goes out of bounds 914 ! do nsp=1,nspec 915 ! if (xmass(nrelpoint,nsp).gt.eps2) goto 889 916 ! end do 917 ! 889 nsp=min(nsp,nspec) 918 ! if (density(nsp).gt.0.) then 919 if (nspec.eq.1.and.density(1).gt.0.) then 920 call get_settling(itime+ldt,real(xt),real(yt),zt,nspec,settling) !bugfix 921 end if 922 w=w+settling 923 endif 917 if (lsettling) then 918 do nsp=1,nspec 919 if (xmass(nrelpoint,nsp).gt.eps2) exit 920 end do 921 if (nsp.gt.nspec) then 922 ! This should never happen 923 write(*,*) 'advance.f90: ERROR: could not find releasepoint' 924 stop 925 end if 926 if (density(nsp).gt.0.) then 927 call get_settling(itime+ldt,real(xt),real(yt),zt,nsp,settling) !bugfix 928 w=w+settling 929 end if 930 endif 931 end if 924 932 925 933 -
src/com_mod.f90
r2bec33e r5184a7c 135 135 !ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed) 136 136 logical :: sumclouds=.false. 137 138 !ESO: Disable settling if more than 1 species per release point 139 logical :: lsettling=.true. 137 140 138 141 logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest -
src/readreleases.f90
r6985a98 r5184a7c 74 74 implicit none 75 75 76 integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat 76 integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle 77 77 integer,parameter :: num_min_discrete=100 78 78 real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun 79 79 real(kind=dp) :: jul1,jul2,julm,juldate 80 real,parameter :: eps2=1.e-9 80 81 character(len=50) :: line 81 82 logical :: old … … 93 94 ! declare namelists 94 95 namelist /releases_ctrl/ & 95 nspec, &96 specnum_rel96 nspec, & 97 specnum_rel 97 98 98 99 namelist /release/ & 99 idate1, itime1, &100 idate2, itime2, &101 lon1, lon2, &102 lat1, lat2, &103 z1, z2, &104 zkind, &105 mass, &106 parts, &107 comment100 idate1, itime1, & 101 idate2, itime2, & 102 lon1, lon2, & 103 lat1, lat2, & 104 z1, z2, & 105 zkind, & 106 mass, & 107 parts, & 108 comment 108 109 109 110 numpoint=0 … … 132 133 if ((readerror.ne.0).or.(nspec.lt.0)) then 133 134 134 135 ! no namelist format, close file and allow reopening in old format 135 136 close(unitreleases) 136 137 open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999) … … 138 139 readerror=1 ! indicates old format 139 140 140 141 142 143 141 ! Check the format of the RELEASES file (either in free format, 142 ! or using a formatted mask) 143 ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION' 144 !************************************************************************** 144 145 145 146 call skplin(12,unitreleases) … … 154 155 155 156 156 157 157 ! Skip first 11 lines (file header) 158 !********************************** 158 159 159 160 call skplin(11,unitreleases) … … 193 194 if (old) call skplin(2,unitreleases) 194 195 end do 195 196 !save compoint only for the first 1000 release points 196 197 read(unitreleases,'(a40)',err=998) compoint(1)(1:40) 197 198 if (old) call skplin(1,unitreleases) … … 228 229 specnum_rel2=specnum_rel(1:nspec) 229 230 deallocate(specnum_rel) 230 ! eso: BUG, crashes here for nspec=12 and maxspec=6,231 ! TODO: catch error and exit231 ! eso: BUG, crashes here for nspec=12 and maxspec=6, 232 ! TODO: catch error and exit 232 233 allocate(specnum_rel(nspec),stat=stat) 233 234 if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel' … … 280 281 281 282 if (readerror.ne.0) then 282 283 283 ! Skip first 11 lines (file header) 284 !********************************** 284 285 285 286 call skplin(11,unitreleases) 286 287 287 288 288 ! Assign species-specific parameters needed for physical processes 289 !************************************************************************* 289 290 290 291 read(unitreleases,*,err=998) nspec … … 307 308 endif 308 309 309 310 311 312 313 314 315 316 317 318 319 320 321 310 ! For backward runs, only 1 species is allowed 311 !********************************************* 312 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 320 321 ! Molecular weight 322 !***************** 322 323 323 324 if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then … … 329 330 endif 330 331 331 332 332 ! Radioactive decay 333 !****************** 333 334 334 335 decay(i)=0.693147/decay(i) !conversion half life to decay constant 335 336 336 337 337 338 338 ! Dry deposition of gases 339 !************************ 339 340 340 341 if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i)) ! mesophyll resistance 341 342 342 343 343 ! Dry deposition of particles 344 !**************************** 344 345 345 346 vsetaver(i)=0. … … 358 359 endif 359 360 360 361 361 ! Dry deposition for constant deposition velocity 362 !************************************************ 362 363 363 364 dryvel(i)=dryvel(i)*0.01 ! conversion to m/s … … 373 374 if (lroot) then 374 375 write (*,*) ' Below-cloud scavenging: ON' 375 ! write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i376 ! write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i 376 377 end if 377 378 else 378 379 if (lroot) write (*,*) ' Below-cloud scavenging: OFF' 379 380 endif 380 381 ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015381 382 ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015 382 383 if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.)) then 383 384 WETDEP=.true. … … 385 386 if (lroot) then 386 387 write (*,*) ' In-cloud scavenging: ON' 387 ! write (*,*) 'In-cloud scavenging coefficients: ',&388 ! &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i388 ! write (*,*) 'In-cloud scavenging coefficients: ',& 389 ! &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i 389 390 end if 390 391 else … … 412 413 numpartmax=0 413 414 releaserate=0. 414 101 415 101 numpoint=numpoint+1 415 416 416 417 if (readerror.lt.1) then ! reading namelist format … … 439 440 compoint(min(1001,numpoint))=comment 440 441 441 442 ! namelist output 442 443 if (nmlout.and.lroot) then 443 444 write(unitreleasesout,nml=release) … … 472 473 mass(i)=xmass(numpoint,i) 473 474 end do 474 475 !save compoint only for the first 1000 release points 475 476 if (numpoint.le.1000) then 476 477 read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40) … … 482 483 if (old) call skplin(1,unitreleases) 483 484 484 485 ! namelist output 485 486 if (nmlout.and.lroot) then 486 487 idate1=id1 … … 524 525 !********************************************************************* 525 526 if (WETBKDEP) then 526 527 528 527 zpoint1(numpoint)=0. 528 zpoint2(numpoint)=20000. 529 kindz(numpoint)=1 529 530 endif 530 531 if (DRYBKDEP) then 531 532 533 532 zpoint1(numpoint)=0. 533 zpoint2(numpoint)=2.*href 534 kindz(numpoint)=1 534 535 endif 535 536 … … 538 539 !********************************************************************* 539 540 540 541 542 543 544 545 546 547 541 if (xpoint1(numpoint).lt.xlon0) & 542 xpoint1(numpoint)=xpoint1(numpoint)+360. 543 if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) & 544 xpoint1(numpoint)=xpoint1(numpoint)-360. 545 if (xpoint2(numpoint).lt.xlon0) & 546 xpoint2(numpoint)=xpoint2(numpoint)+360. 547 if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) & 548 xpoint2(numpoint)=xpoint2(numpoint)-360. 548 549 549 550 ! Determine relative beginning and ending times of particle release … … 607 608 goto 101 608 609 609 250 610 250 close(unitreleases) 610 611 611 612 if (nmlout.and.lroot) then … … 623 624 endif 624 625 626 ! Disable settling if more than 1 species at any release point 627 ! or if MQUASILAG and more than one species 628 !************************************************************* 629 630 if (mquasilag.ne.0) then 631 if (nspec.gt.1) lsettling=.false. 632 else 633 do irel=1,numpoint 634 nsettle=0 635 do ispc=1,nspec 636 if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1 637 end do 638 if (nsettle.gt.1) lsettling=.false. 639 end do 640 end if 641 642 if (lroot) then 643 if (.not.lsettling) then 644 write(*,*) 'WARNING: more than 1 species per release point, settling & 645 &disabled' 646 end if 647 end if 648 625 649 ! Check, whether the total number of particles may exceed totally allowed 626 650 ! number of particles at some time during the simulation … … 630 654 0.99*real(maxpart)/real(lage(nageclass))) then 631 655 if (numpartmax.gt.maxpart.and.lroot) then 632 write(*,*) '#####################################################'633 write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####'634 write(*,*) '#### ####'635 write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####'636 write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####'637 write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####'638 write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####'639 write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####'640 write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####'641 write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####'642 write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####'643 write(*,*) '#####################################################'656 write(*,*) '#####################################################' 657 write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' 658 write(*,*) '#### ####' 659 write(*,*) '####WARNING - TOTAL NUMBER OF PARTICLES SPECIFIED####' 660 write(*,*) '#### IN FILE "RELEASES" MAY AT SOME POINT DURING ####' 661 write(*,*) '#### THE SIMULATION EXCEED THE MAXIMUM ALLOWED ####' 662 write(*,*) '#### NUMBER (MAXPART).IF RELEASES DO NOT OVERLAP,####' 663 write(*,*) '#### FLEXPART CAN POSSIBLY COMPLETE SUCCESSFULLY.####' 664 write(*,*) '#### HOWEVER, FLEXPART MAY HAVE TO STOP ####' 665 write(*,*) '#### AT SOME TIME DURING THE SIMULATION. PLEASE ####' 666 write(*,*) '#### MAKE SURE THAT YOUR SETTINGS ARE CORRECT. ####' 667 write(*,*) '#####################################################' 644 668 write(*,*) 'Maximum release rate may be: ',releaserate, & 645 669 ' particles per second' … … 653 677 endif 654 678 679 655 680 if (lroot) then 656 681 write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec)) … … 659 684 return 660 685 661 994 686 994 write(*,*) '#####################################################' 662 687 write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' 663 688 write(*,*) '#### ####' … … 667 692 stop 668 693 669 998 694 998 write(*,*) '#####################################################' 670 695 write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES: ####' 671 696 write(*,*) '#### ####' … … 678 703 679 704 680 999 705 999 write(*,*) '#####################################################' 681 706 write(*,*) ' FLEXPART MODEL SUBROUTINE READRELEASES: ' 682 707 write(*,*)
Note: See TracChangeset
for help on using the changeset viewer.