Changeset 5184a7c in flexpart.git


Ignore:
Timestamp:
Jun 20, 2017, 9:09:35 AM (7 years ago)
Author:
Espen Sollum ATMOS <eso@…>
Branches:
master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
Children:
b5127f9
Parents:
a76d954
Message:

Enable settling with multiple species if from separate releases

Location:
src
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • src/FLEXPART_MPI.f90

    rc8fc724 r5184a7c  
    8080  ! FLEXPART version string
    8181  flexversion_major = '10' ! Major version number, also used for species file names
    82 !  flexversion='Ver. 10 Beta MPI (2015-05-01)'
    8382  flexversion='Ver. '//trim(flexversion_major)//'.1beta MPI (2016-11-02)'
    8483  verbosity=0
  • src/advance.f90

    re52967c r5184a7c  
    535535
    536536      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
    545550        end if
    546         w=w+settling
    547551      endif
    548552
     
    700704  !*************************************************************************
    701705
    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
    715715      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
    717720    endif
     721  end if
     722
    718723
    719724  ! Calculate position at time step itime+lsynctime
     
    910915
    911916  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
    924932
    925933
  • src/com_mod.f90

    r2bec33e r5184a7c  
    135135!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
    136136  logical :: sumclouds=.false.
     137
     138!ESO: Disable settling if more than 1 species per release point
     139  logical :: lsettling=.true.
    137140
    138141  logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
  • src/readreleases.f90

    r6985a98 r5184a7c  
    7474  implicit none
    7575
    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
    7777  integer,parameter :: num_min_discrete=100
    7878  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
    7979  real(kind=dp) :: jul1,jul2,julm,juldate
     80  real,parameter :: eps2=1.e-9
    8081  character(len=50) :: line
    8182  logical :: old
     
    9394  ! declare namelists
    9495  namelist /releases_ctrl/ &
    95     nspec, &
    96     specnum_rel
     96       nspec, &
     97       specnum_rel
    9798
    9899  namelist /release/ &
    99     idate1, itime1, &
    100     idate2, itime2, &
    101     lon1, lon2, &
    102     lat1, lat2, &
    103     z1, z2, &
    104     zkind, &
    105     mass, &
    106     parts, &
    107     comment
     100       idate1, itime1, &
     101       idate2, itime2, &
     102       lon1, lon2, &
     103       lat1, lat2, &
     104       z1, z2, &
     105       zkind, &
     106       mass, &
     107       parts, &
     108       comment
    108109
    109110  numpoint=0
     
    132133  if ((readerror.ne.0).or.(nspec.lt.0)) then
    133134
    134     ! no namelist format, close file and allow reopening in old format
     135  ! no namelist format, close file and allow reopening in old format
    135136    close(unitreleases)
    136137    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
     
    138139    readerror=1 ! indicates old format
    139140
    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     !**************************************************************************
     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  !**************************************************************************
    144145
    145146    call skplin(12,unitreleases)
     
    154155
    155156
    156     ! Skip first 11 lines (file header)
    157     !**********************************
     157  ! Skip first 11 lines (file header)
     158  !**********************************
    158159
    159160    call skplin(11,unitreleases)
     
    193194      if (old) call skplin(2,unitreleases)
    194195    end do
    195     !save compoint only for the first 1000 release points
     196  !save compoint only for the first 1000 release points
    196197    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
    197198    if (old) call skplin(1,unitreleases)
     
    228229  specnum_rel2=specnum_rel(1:nspec)
    229230  deallocate(specnum_rel)
    230 ! eso: BUG, crashes here for nspec=12 and maxspec=6,
    231 ! TODO: catch error and exit
     231  ! eso: BUG, crashes here for nspec=12 and maxspec=6,
     232  ! TODO: catch error and exit
    232233  allocate(specnum_rel(nspec),stat=stat)
    233234  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
     
    280281
    281282  if (readerror.ne.0) then
    282     ! Skip first 11 lines (file header)
    283     !**********************************
     283  ! Skip first 11 lines (file header)
     284  !**********************************
    284285
    285286    call skplin(11,unitreleases)
    286287
    287     ! Assign species-specific parameters needed for physical processes
    288     !*************************************************************************
     288  ! Assign species-specific parameters needed for physical processes
     289  !*************************************************************************
    289290
    290291    read(unitreleases,*,err=998) nspec
     
    307308    endif
    308309
    309     ! For backward runs, only 1 species is allowed
    310     !*********************************************
    311 
    312     !if ((ldirect.lt.0).and.(nspec.gt.1)) then
    313     !write(*,*) '#####################################################'
    314     !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    315     !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
    316     !write(*,*) '#####################################################'
    317     !  stop
    318     !endif
    319 
    320     ! Molecular weight
    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  !*****************
    322323
    323324    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
     
    329330    endif
    330331
    331     ! Radioactive decay
    332     !******************
     332  ! Radioactive decay
     333  !******************
    333334
    334335    decay(i)=0.693147/decay(i) !conversion half life to decay constant
    335336
    336337
    337     ! Dry deposition of gases
    338     !************************
     338  ! Dry deposition of gases
     339  !************************
    339340
    340341    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
    341342
    342     ! Dry deposition of particles
    343     !****************************
     343  ! Dry deposition of particles
     344  !****************************
    344345
    345346    vsetaver(i)=0.
     
    358359    endif
    359360
    360     ! Dry deposition for constant deposition velocity
    361     !************************************************
     361  ! Dry deposition for constant deposition velocity
     362  !************************************************
    362363
    363364    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
     
    373374      if (lroot) then
    374375        write (*,*) '  Below-cloud scavenging: ON'
    375 !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
     376  !  write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
    376377      end if
    377378    else
    378379      if (lroot) write (*,*) '  Below-cloud scavenging: OFF'
    379380    endif
    380    
    381 ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
     381
     382  ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
    382383    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
    383384      WETDEP=.true.
     
    385386      if (lroot) then
    386387        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),i
     388  !        write (*,*) 'In-cloud scavenging coefficients: ',&
     389  !           &ccn_aero(i),in_aero(i),i !,wetc_in(i), wetd_in(i),i
    389390      end if
    390391    else
     
    412413  numpartmax=0
    413414  releaserate=0.
    414 101   numpoint=numpoint+1
     415101 numpoint=numpoint+1
    415416
    416417  if (readerror.lt.1) then ! reading namelist format
     
    439440    compoint(min(1001,numpoint))=comment
    440441
    441     ! namelist output
     442  ! namelist output
    442443    if (nmlout.and.lroot) then
    443444      write(unitreleasesout,nml=release)
     
    472473      mass(i)=xmass(numpoint,i)
    473474    end do
    474     !save compoint only for the first 1000 release points
     475  !save compoint only for the first 1000 release points
    475476    if (numpoint.le.1000) then
    476477      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
     
    482483    if (old) call skplin(1,unitreleases)
    483484
    484     ! namelist output
     485  ! namelist output
    485486    if (nmlout.and.lroot) then
    486487      idate1=id1
     
    524525  !*********************************************************************
    525526  if (WETBKDEP) then
    526       zpoint1(numpoint)=0.
    527       zpoint2(numpoint)=20000.
    528       kindz(numpoint)=1
     527    zpoint1(numpoint)=0.
     528    zpoint2(numpoint)=20000.
     529    kindz(numpoint)=1
    529530  endif
    530531  if (DRYBKDEP) then
    531       zpoint1(numpoint)=0.
    532       zpoint2(numpoint)=2.*href
    533       kindz(numpoint)=1
     532    zpoint1(numpoint)=0.
     533    zpoint2(numpoint)=2.*href
     534    kindz(numpoint)=1
    534535  endif
    535536
     
    538539  !*********************************************************************
    539540
    540    if (xpoint1(numpoint).lt.xlon0) &
    541         xpoint1(numpoint)=xpoint1(numpoint)+360.
    542    if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
    543         xpoint1(numpoint)=xpoint1(numpoint)-360.
    544    if (xpoint2(numpoint).lt.xlon0) &
    545         xpoint2(numpoint)=xpoint2(numpoint)+360.
    546    if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
    547         xpoint2(numpoint)=xpoint2(numpoint)-360.
     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.
    548549
    549550  ! Determine relative beginning and ending times of particle release
     
    607608  goto 101
    608609
    609 250   close(unitreleases)
     610250 close(unitreleases)
    610611
    611612  if (nmlout.and.lroot) then
     
    623624  endif
    624625
     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
    625649  ! Check, whether the total number of particles may exceed totally allowed
    626650  ! number of particles at some time during the simulation
     
    630654       0.99*real(maxpart)/real(lage(nageclass))) then
    631655    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(*,*) '#####################################################'
    644668      write(*,*) 'Maximum release rate may be: ',releaserate, &
    645669           ' particles per second'
     
    653677  endif
    654678
     679
    655680  if (lroot) then
    656681    write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
     
    659684  return
    660685
    661 994   write(*,*) '#####################################################'
     686994 write(*,*) '#####################################################'
    662687  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    663688  write(*,*) '####                                             ####'
     
    667692  stop
    668693
    669 998   write(*,*) '#####################################################'
     694998 write(*,*) '#####################################################'
    670695  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    671696  write(*,*) '####                                             ####'
     
    678703
    679704
    680 999   write(*,*) '#####################################################'
     705999 write(*,*) '#####################################################'
    681706  write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
    682707  write(*,*)
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG