Changeset d7935de in flexpart.git for src/readreleases.f90

Timestamp:

Sep 11, 2018, 6:06:31 PM (6 years ago)

Author:

pesei <petra seibert at univie ac at>

Branches:

univie

Children:

93786a1

Parents:

34f1452

Message:

modify most input read subroutines

changed some variable names (mostly for I-N reasons)
includes two names appearing also in timemanager, com_mod
corrected a few mistakes
simplified some parts of code
changed options/RELEASES which is in nml fmt correspondingly

File:

• src/readreleases.f90 (modified) (27 diffs)

src/readreleases.f90

@@ -33 +33 @@
   !     Update: 29 January 2001                                                *
   !     Release altitude can be either in magl or masl                         *
-  !     HSO, 12 August 2013
-  !     Added optional namelist input
+  !     HSO, 12 August 2013: Added optional namelist input                     *
+  !     PS, 2/2015: access= -> position= (ported from 9.2 review)              *
+  !     PS, 8/2018: check for uninitialised mass (ticket:204),                 *
+  !          rename int/real variables, a bit of code layout improvement       *
   !                                                                            *
   !*****************************************************************************
@@ -74 +76 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat,irel,ispc,nsettle
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,iallocstat,irel,ispc,nsettle
   integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
   real(kind=dp) :: jul1,jul2,julm,juldate
-  real,parameter :: eps2=1.e-9
   character(len=50) :: line
   logical :: old
 
   ! help variables for namelist reading
-  integer :: numpoints, parts, readerror
-  integer*2 :: zkind
+  integer :: numpoints, nparts, ios
+  integer*2 :: ikindz
   integer :: idate1, itime1, idate2, itime2
-  real :: lon1,lon2,lat1,lat2,z1,z2
+  real :: rlon1,rlon2,rlat1,rlat2,z1,z2
   character*40 :: comment
   integer,parameter :: unitreleasesout=2
-  real,allocatable, dimension (:) :: mass
-  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
+  real,allocatable, dimension (:) :: rmass
+  integer,allocatable, dimension (:) :: numspec_rel,numspec_rel2
 
   ! declare namelists
   namelist /releases_ctrl/ &
        nspec, &
-       specnum_rel
+       numspec_rel
 
   namelist /release/ &
        idate1, itime1, &
        idate2, itime2, &
-       lon1, lon2, &
-       lat1, lat2, &
+       rlon1, rlon2, &
+       rlat1, rlat2, &
        z1, z2, &
-       zkind, &
-       mass, &
-       parts, &
+       ikindz, &
+       rmass, &
+       nparts, &
        comment
+! initialiase mass
 
   numpoint=0
 
-  ! allocate with maxspec for first input loop
-  allocate(mass(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel(maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-
-  ! presetting namelist releases_ctrl
+! allocate with maxspec for first input loop
+  allocate(rmass(maxspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate rmass'
+  allocate(numspec_rel(maxspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel'
+
+! presetting namelist releases_ctrl
   nspec = -1  ! use negative value to determine failed namelist input
-  specnum_rel = 0
-
-  !sec, read release to find how many releasepoints should be allocated
-  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
-
-  ! check if namelist input provided
-  read(unitreleases,releases_ctrl,iostat=readerror)
-
-  ! prepare namelist output if requested
+  numspec_rel = 0
+
+!SEC: read release to find how many releasepoints should be allocated
+  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
+    form='formatted',err=999)
+
+! prepare namelist output if requested
   if (nmlout.and.lroot) then
-    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='replace',err=1000)
-  endif
-
-  if ((readerror.ne.0).or.(nspec.lt.0)) then
-
-  ! no namelist format, close file and allow reopening in old format
+    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist', &
+      position='append',status='replace',err=1000)
+  endif
+
+! check whether namelist input is available
+  read(unitreleases,releases_ctrl,iostat=ios)
+
+  if (ios .ne. 0 .or. nspec .lt. 0 ) then
+
+  ! not in namelist format, close file and reopen as simple text file
+
     close(unitreleases)
-    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
-
-    readerror=1 ! indicates old format
-
-  ! Check the format of the RELEASES file (either in free format,
-  ! or using a formatted mask)
+    open(unitreleases,file=path(1)(1:length(1))//'RELEASES', status='old', &
+      err=999)
+
+    ios=1 ! indicates old format
+
+  ! Check the format of the RELEASES file
+  ! (either in free format, or using a formatted mask)
   ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
   !**************************************************************************
+! Note PS: something has changed, now Total is checked. Is this compatible??
 
     call skplin(12,unitreleases)
@@ -163 +170 @@
     if (old) call skplin(2,unitreleases)
     do i=1,nspec
-      read(unitreleases,*,err=998) specnum_rel(i)
+      read(unitreleases,*,err=998) numspec_rel(i)
       if (old) call skplin(2,unitreleases)
     end do
@@ -202 +209 @@
 25  numpoint=numpoint-1
 
-  else 
-
-    readerror=0
-    do while (readerror.eq.0) 
+  else ! read namelist input
+
+    ios=0
+    do while (ios.eq.0)  ! loop over release points
       idate1=-1
-      read(unitreleases,release,iostat=readerror)
-      if ((idate1.lt.0).or.(readerror.ne.0)) then
-        readerror=1
+      read(unitreleases,release,iostat=ios)
+      if (idate1.lt.0 .or. ios.ne.0) then
+        ios=1
       else
         numpoint=numpoint+1
       endif
     end do
-    readerror=0
-  endif ! if namelist input
+    ios=0
+
+  endif 
 
   rewind(unitreleases)
@@ -221 +229 @@
   if (nspec.gt.maxspec) goto 994
 
-  ! allocate arrays of matching size for number of species (namelist output)
-  deallocate(mass)
-  allocate(mass(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
-  allocate(specnum_rel2(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
-  specnum_rel2=specnum_rel(1:nspec)
-  deallocate(specnum_rel) 
+! allocate arrays of matching size for number of species (namelist output)
+  deallocate(rmass)
+  allocate(rmass(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate rmass'
+  allocate(numspec_rel2(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel2'
+
+  numspec_rel2 = numspec_rel(1:nspec)
+
+  deallocate(numspec_rel)
   ! eso: BUG, crashes here for nspec=12 and maxspec=6,
   ! TODO: catch error and exit
-  allocate(specnum_rel(nspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
-  specnum_rel=specnum_rel2
-  deallocate(specnum_rel2)
-
-  !allocate memory for numpoint releaspoints
-  allocate(ireleasestart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
-  allocate(ireleaseend(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
-  allocate(xpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
-  allocate(xpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
-  allocate(ypoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
-  allocate(ypoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
-  allocate(zpoint1(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
-  allocate(zpoint2(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
-  allocate(kindz(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
-  allocate(xmass(numpoint,maxspec),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
-  allocate(rho_rel(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
-  allocate(npart(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
-  allocate(xmasssave(numpoint),stat=stat)
-  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
-
-  if (lroot) write (*,*) 'Releasepoints : ', numpoint
-
-  do i=1,numpoint
-    xmasssave(i)=0.
-  end do
+  allocate(numspec_rel(nspec),stat=iallocstat)
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate numspec_rel'
+
+  numspec_rel = numspec_rel2
+
+  deallocate(numspec_rel2)
+
+! allocate memory for numpoint releaspoints
+  allocate ( &
+    ireleasestart(numpoint), ireleaseend(numpoint),xpoint1(numpoint), &
+    xpoint2(numpoint), ypoint1(numpoint), ypoint2(numpoint), &
+    zpoint1(numpoint), zpoint2(numpoint), kindz(numpoint), &
+    xmass(numpoint,maxspec), rho_rel(numpoint), npart(numpoint), &
+    xmasssave(numpoint), &
+      stat=iallocstat )
+  if (iallocstat.ne.0) write(*,*)'ERROR: could not allocate variables' // &
+    ' dimensioned numpoint =', numpoint
+
+  if (lroot) write (*,*) 'Release points : ', numpoint
+
+  xmasssave(:)=0.
 
   !now save the information
-  DEP=.false.
-  DRYDEP=.false.
-  WETDEP=.false.
-  OHREA=.false.
-  do i=1,maxspec
-    DRYDEPSPEC(i)=.false.
-    WETDEPSPEC(i)=.false.
-  end do
-
-  if (readerror.ne.0) then
+  dep   =.false.
+  drydep=.false.
+  wetdep=.false.
+  ohrea =.false.
+  drydepspec(:)=.false.
+  wetdepspec(:)=.false.
+
+  if (ios.ne.0) then
+
   ! Skip first 11 lines (file header)
   !**********************************
@@ -292 +284 @@
     if (nspec.gt.maxspec) goto 994
     if (old) call skplin(2,unitreleases)
-  endif
-
-  ! namelist output
-  if (nmlout.and.lroot) then
-    write(unitreleasesout,nml=releases_ctrl)
-  endif
-
+
+  endif
+
+! namelist output
+  if (nmlout.and.lroot) write(unitreleasesout,nml=releases_ctrl)
+
+! species loop
   do i=1,nspec
-    if (readerror.ne.0) then
-      read(unitreleases,*,err=998) specnum_rel(i)
+  
+    if (ios.ne.0) then
+      read(unitreleases,*,err=998) numspec_rel(i)
       if (old) call skplin(2,unitreleases)
-      call readspecies(specnum_rel(i),i)
+      call readspecies(numspec_rel(i),i)
     else
-      call readspecies(specnum_rel(i),i)
+      call readspecies(numspec_rel(i),i)
     endif
 
@@ -322 +315 @@
   !*****************
 
-    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
+    if ( (iout.eq.2 .or. iout.eq.3) .and. weightmolar(i).lt.0. ) then
       write(*,*) 'For mixing ratio output, valid molar weight'
       write(*,*) 'must be specified for all simulated species.'
@@ -334 +327 @@
 
     decay(i)=0.693147/decay(i) !conversion half life to decay constant
-
+! Note PS: this should be replace by = alog(2.)/decay(i)
+! keep it for the moment so that results are not affect in test phase of v10
 
   ! Dry deposition of gases
   !************************
 
-    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
+    if (reldiff(i).gt.0.) rm(i)=1./( henry(i)/3000. + 100.*f0(i) )! mesophyll resistance
 
   ! Dry deposition of particles
   !****************************
 
-    vsetaver(i)=0.
-    cunningham(i)=0.
-    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
+    vsetaver(i) = 0.
+    cunningham(i) = 0.
+    dquer(i) = dquer(i)*1.e6 ! Conversion m to um
     if (density(i).gt.0.) then         ! Additional parameters
       call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
       do j=1,ni
-        fract(i,j)=fracth(j)
-        schmi(i,j)=schmih(j)
-        vset(i,j)=vsh(j)
-        cunningham(i)=cunningham(i)+cun*fract(i,j)
-        vsetaver(i)=vsetaver(i)-vset(i,j)*fract(i,j)
+        fract(i,j) = fracth(j)
+        schmi(i,j) = schmih(j)
+        vset(i,j) = vsh(j)
+        cunningham(i) = cunningham(i) + cun*fract(i,j)
+        vsetaver(i) = vsetaver(i) - vset(i,j)*fract(i,j)
       end do
       if (lroot) write(*,*) 'Average settling velocity: ',i,vsetaver(i)
@@ -368 +362 @@
 
   ! ESO 04.2016 check for below-cloud scavenging (gas or aerosol)
-    if ((dquer(i).le.0..and.(weta_gas(i).gt.0. .or. wetb_gas(i).gt.0.)) .or. &
-         &(dquer(i).gt.0. .and. (crain_aero(i) .gt. 0. .or. csnow_aero(i).gt.0.)))  then
-      WETDEP=.true.
-      WETDEPSPEC(i)=.true.
+    if ( &
+      dquer(i).le.0. .and. (weta_gas(i)  .gt.0. .or. wetb_gas(i)  .gt.0.) .or. &
+      dquer(i).gt.0. .and. (crain_aero(i).gt.0. .or. csnow_aero(i).gt.0.)  &
+       )  then
+      wetdep        = .true.
+      wetdepspec(i) = .true.
       if (lroot) then
         write (*,*) '  Below-cloud scavenging: ON'
@@ -381 +377 @@
 
   ! NIK 31.01.2013 + 10.12.2013 + 15.02.2015
-    if (dquer(i).gt.0..and.(ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
-      WETDEP=.true.
-      WETDEPSPEC(i)=.true.
+    if (dquer(i).gt.0. .and. (ccn_aero(i).gt.0. .or. in_aero(i).gt.0.))  then
+      wetdep        = .true.
+      wetdepspec(i) = .true.
       if (lroot) then
         write (*,*) '  In-cloud scavenging: ON'
@@ -394 +390 @@
 
     if (ohcconst(i).gt.0.) then
-      OHREA=.true.
+      ohrea = .true.
       if (lroot) write (*,*) '  OHreaction switched on: ',ohcconst(i),i
     endif
 
-    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
-      DRYDEP=.true.
-      DRYDEPSPEC(i)=.true.
+    if (reldiff(i).gt.0. .or. density(i).gt.0. .or. dryvel(i).gt.0.) then
+      drydep        = .true.
+      drydepspec(i) = .true.
     endif
 
@@ -415 +411 @@
 101 numpoint=numpoint+1
 
-  if (readerror.lt.1) then ! reading namelist format
+  if (ios.lt.1) then ! reading namelist format
 
     if (numpoint.gt.numpoints) goto 250
-    zkind = 1
-    mass = 0
-    parts = 0
+    ikindz = 1
+    rmass = -1. ! initialise as negative to catch if not given in input
+    nparts = 0
     comment = ' '
-    read(unitreleases,release,iostat=readerror)
+    read(unitreleases,release,iostat=ios)
     id1=idate1
     it1=itime1
     id2=idate2
     it2=itime2
-    xpoint1(numpoint)=lon1
-    xpoint2(numpoint)=lon2
-    ypoint1(numpoint)=lat1
-    ypoint2(numpoint)=lat2
+    xpoint1(numpoint)=rlon1
+    xpoint2(numpoint)=rlon2
+    ypoint1(numpoint)=rlat1
+    ypoint2(numpoint)=rlat2
     zpoint1(numpoint)=z1
     zpoint2(numpoint)=z2
-    kindz(numpoint)=zkind
+    kindz(numpoint)=ikindz
     do i=1,nspec
-      xmass(numpoint,i)=mass(i)
+      if (rmass(i) .lt. 0.) goto 995 ! write error and stop
+      xmass(numpoint,i)=rmass(i)
     end do
-    npart(numpoint)=parts
+    npart(numpoint)=nparts
     compoint(min(1001,numpoint))=comment
 
@@ -445 +442 @@
     endif
 
-  else
+  else ! simple text format
 
     read(unitreleases,*,end=250)
@@ -471 +468 @@
       read(unitreleases,*,err=998) xmass(numpoint,i)
       if (old) call skplin(2,unitreleases)
-      mass(i)=xmass(numpoint,i)
+      rmass(i)=xmass(numpoint,i)
     end do
   !save compoint only for the first 1000 release points
@@ -489 +486 @@
       idate2=id2
       itime2=it2
-      lon1=xpoint1(numpoint)
-      lon2=xpoint2(numpoint)
-      lat1=ypoint1(numpoint)
-      lat2=ypoint2(numpoint)
+      rlon1=xpoint1(numpoint)
+      rlon2=xpoint2(numpoint)
+      rlat1=ypoint1(numpoint)
+      rlat2=ypoint2(numpoint)
       z1=zpoint1(numpoint)
       z2=zpoint2(numpoint)
-      zkind=kindz(numpoint)
-      parts=npart(numpoint)
+      ikindz=kindz(numpoint)
+      nparts=npart(numpoint)
       write(unitreleasesout,nml=release)
     endif
 
     if (numpoint.le.1000) then
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
-           (compoint(numpoint)(1:8).eq.'        ')) goto 250
+      if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
+          xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0. .and. &
+          compoint(numpoint)(1:8).eq.'        ') goto 250
     else
-      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
+      if (xpoint1(numpoint).eq.0. .and. ypoint1(numpoint).eq.0. .and. &
+          xpoint2(numpoint).eq.0. .and. ypoint2(numpoint).eq.0.) goto 250
     endif
 
@@ -539 +536 @@
   !*********************************************************************
 
-  if (xpoint1(numpoint).lt.xlon0) &
-       xpoint1(numpoint)=xpoint1(numpoint)+360.
-  if (xpoint1(numpoint).gt.xlon0+(nxmin1)*dx) &
-       xpoint1(numpoint)=xpoint1(numpoint)-360.
-  if (xpoint2(numpoint).lt.xlon0) &
-       xpoint2(numpoint)=xpoint2(numpoint)+360.
-  if (xpoint2(numpoint).gt.xlon0+(nxmin1)*dx) &
-       xpoint2(numpoint)=xpoint2(numpoint)-360.
+  if (xpoint1(numpoint) .lt. xlon0) &
+      xpoint1(numpoint) = xpoint1(numpoint) + 360.
+  if (xpoint1(numpoint) .gt. xlon0+nxmin1*dx) &
+      xpoint1(numpoint) = xpoint1(numpoint) - 360.
+  if (xpoint2(numpoint) .lt. xlon0) &
+      xpoint2(numpoint) = xpoint2(numpoint) + 360.
+  if (xpoint2(numpoint) .gt. xlon0+nxmin1*dx) &
+      xpoint2(numpoint) = xpoint2(numpoint) - 360.
 
   ! Determine relative beginning and ending times of particle release
@@ -562 +559 @@
   if (mdomainfill.eq.0) then   ! no domain filling
     if (ldirect.eq.1) then
-      if ((jul1.lt.bdate).or.(jul2.gt.edate)) then
+      if (jul1.lt.bdate .or. jul2.gt.edate) then
         write(*,*) 'FLEXPART MODEL ERROR'
         write(*,*) 'Release starts before simulation begins or ends'
@@ -577 +574 @@
       endif
     else if (ldirect.eq.-1) then
-      if ((jul1.lt.edate).or.(jul2.gt.bdate)) then
+      if (jul1.lt.edate .or. jul2.gt.bdate) then
         write(*,*) 'FLEXPART MODEL ERROR'
         write(*,*) 'Release starts before simulation begins or ends'
@@ -586 +583 @@
       if (npart(numpoint).gt.num_min_discrete) then
         ireleasestart(numpoint)=int((jul1-bdate)*86400.)
-        ireleaseend(numpoint)=int((jul2-bdate)*86400.)
+        ireleaseend(numpoint) = int((jul2-bdate)*86400.)
       else
         ireleasestart(numpoint)=int((julm-bdate)*86400.)
-        ireleaseend(numpoint)=int((julm-bdate)*86400.)
+        ireleaseend(numpoint) = int((julm-bdate)*86400.)
       endif
     endif
@@ -600 +597 @@
 
   if (ireleasestart(numpoint).ne.ireleaseend(numpoint)) then
-    releaserate=releaserate+real(npart(numpoint))/ &
-         real(ireleaseend(numpoint)-ireleasestart(numpoint))
+    releaserate = releaserate + real( npart(numpoint) ) / &
+      real( ireleaseend(numpoint) - ireleasestart(numpoint) )
   else
     releaserate=99999999
@@ -634 +631 @@
       nsettle=0
       do ispc=1,nspec
-        if (xmass(irel,ispc).gt.eps2) nsettle=nsettle+1
+        if (xmass(irel,ispc).gt.tiny(1.)) nsettle=nsettle+1
       end do
       if (nsettle.gt.1) lsettling=.false.
@@ -641 +638 @@
 
   if (lroot) then
-    if (.not.lsettling) then 
-      write(*,*) 'WARNING: more than 1 species per release point, settling &
-           &disabled'
-    end if
+    if (.not.lsettling) &
+     write(*,*) 'WARNING: more than 1 species per release point,' &
+           // ' settling disabled'
   end if
 
@@ -651 +647 @@
   !************************************************************************
 
-  if (releaserate.gt. &
-       0.99*real(maxpart)/real(lage(nageclass))) then
+  if (releaserate .gt. 0.99*real(maxpart)/real(lage(nageclass))) then
     if (numpartmax.gt.maxpart.and.lroot) then
       write(*,*) '#####################################################'
@@ -670 +665 @@
       write(*,*) 'Maximum allowed release rate is: ', &
            real(maxpart)/real(lage(nageclass)),' particles per second'
-      write(*,*) &
-           'Total number of particles released during the simulation is: ', &
-           numpartmax
+      write(*,*) 'Total number of particles released during the simulation is:'&
+        ,numpartmax
       write(*,*) 'Maximum allowed number of particles is: ',maxpart
     endif
@@ -679 +673 @@
 
   if (lroot) then
-    write(*,FMT='(A,ES14.7)') ' Total mass released:', sum(xmass(1:numpoint,1:nspec))
+    write(*,FMT='(A,ES14.7)') ' Total mass released:', &
+      sum(xmass(1:numpoint,1:nspec))
   end if
 
   return
 
-994 write(*,*) '#####################################################'
+995 write(*,*) '#####################################################'
   write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
   write(*,*) '####                                             ####'
-  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS####'
-  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES. ####'
+  write(*,*) '#### ERROR - mass for species',i,' not given     ####'
   write(*,*) '#####################################################'
+  stop 'readreleases - undefined mass'
+
+994 write(*,*) '######################################################'
+  write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:      ####'
+  write(*,*) '####                                              ####'
+  write(*,*) '#### ERROR - MAXIMUM NUMBER OF EMITTED SPECIES IS ####'
+  write(*,*) '#### TOO LARGE. PLEASE REDUCE NUMBER OF SPECIES.  ####'
+  write(*,*) '######################################################'
   stop
 
@@ -702 +704 @@
   stop
 
-
 999 write(*,*) '#####################################################'
   write(*,*) '   FLEXPART MODEL SUBROUTINE READRELEASES: '
@@ -713 +714 @@
 
 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,*)    ' #### CANNOT BE OPENED IN THE DIRECTORY        #### '
   write(*,'(a)') path(2)(1:length(2))
   stop