Changeset 27

Timestamp:

Jul 3, 2014, 2:55:50 PM (9 years ago)

Author:

hasod

Message:

• Implemented optional namelist input for COMMAND, RELEASES, SPECIES, AGECLASSES,OUTGRID,OUTGRID_NEST,RECEPTORS
• Implemented com_mod switch nmlout to write input files as namelist to the output directory (.true. by default)
• Proposed updated startup and runtime output (may change back to previous info if desired)

Location:

trunk/src

Files:

• FLEXPART.f90 (modified) (21 diffs)
• com_mod.f90 (modified) (1 diff)
• gridcheck.f90 (modified) (1 diff)
• gridcheck_gfs.f90 (modified) (2 diffs)
• gridcheck_nests.f90 (modified) (1 diff)
• outgrid_init.f90 (modified) (1 diff)
• par_mod.f90 (modified) (2 diffs)
• readageclasses.f90 (modified) (6 diffs)
• readavailable.f90 (modified) (1 diff)
• readcommand.f90 (modified) (8 diffs)
• readoutgrid.f90 (modified) (6 diffs)
• readoutgrid_nest.f90 (modified) (2 diffs)
• readpaths.f90 (modified) (1 diff)
• readreceptors.f90 (modified) (5 diffs)
• readreleases.f90 (modified) (11 diffs)
• readspecies.f90 (modified) (6 diffs)
• timemanager.f90 (modified) (2 diffs)

trunk/src/FLEXPART.f90

@@ -60 +60 @@
   call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
 
-  ! 
-  flexversion='Version 9.2 beta (2014-05-23)'
-  !verbosity=0
+  ! FLEXPART version string
+  flexversion='Version 9.2 beta (2014-07-01)'
+  verbosity=0
+
   ! Read the pathnames where input/output files are stored
   !*******************************************************
@@ -76 +77 @@
     call getarg(1,arg1)
     pathfile=arg1
-    verbosity=0
     if (arg1(1:1).eq.'-') then
-        write(pathfile,'(a11)') './pathnames'
-        inline_options=arg1 
+      write(pathfile,'(a11)') './pathnames'
+      inline_options=arg1 
     endif
   case (0)
     write(pathfile,'(a11)') './pathnames'
-    verbosity=0
   end select
   
-    if (inline_options(1:1).eq.'-') then
-      print*, 'inline options=', inline_options       
-      if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
-         print*, 'verbose mode 1: additional information will be displayed'
-         verbosity=1
-      endif
-      if (trim(inline_options).eq.'-v2') then
-         print*, 'verbose mode 2: additional information will be displayed'
-         verbosity=2
-      endif
-      if (trim(inline_options).eq.'-i') then
-         print*, 'info mode: will provide run specific information and stop'
-         verbosity=1
-         info_flag=1
-      endif
-      if (trim(inline_options).eq.'-i2') then
-         print*, 'info mode: will provide run specific information and stop'
-         verbosity=2
-         info_flag=1
-      endif
-    endif
-
-
   ! Print the GPL License statement
   !*******************************************************
   print*,'Welcome to FLEXPART ', trim(flexversion)
-  print*,'FLEXPART is free software released under the GNU Genera'// &
-       'l Public License.'
-            
-  if (verbosity.gt.0) then
-        WRITE(*,*) 'call readpaths'
+  print*,'FLEXPART is free software released under the GNU General Public License.'
+ 
+  if (inline_options(1:1).eq.'-') then
+    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
+       print*, 'Verbose mode 1: display detailed information during run'
+       verbosity=1
+    endif
+    if (trim(inline_options).eq.'-v2') then
+       print*, 'Verbose mode 2: display more detailed information during run'
+       verbosity=2
+    endif
+    if (trim(inline_options).eq.'-i') then
+       print*, 'Info mode: provide detailed run specific information and stop'
+       verbosity=1
+       info_flag=1
+    endif
+    if (trim(inline_options).eq.'-i2') then
+       print*, 'Info mode: provide more detailed run specific information and stop'
+       verbosity=2
+       info_flag=1
+    endif
+  endif
+           
+  if (verbosity.gt.0) then
+    write(*,*) 'call readpaths'
   endif 
   call readpaths(pathfile)
-   
  
   if (verbosity.gt.1) then !show clock info 
@@ -131 +126 @@
   endif
 
-
   ! Read the user specifications for the current model run
   !*******************************************************
 
   if (verbosity.gt.0) then
-        WRITE(*,*) 'call readcommand'
+    write(*,*) 'call readcommand'
   endif
   call readcommand
   if (verbosity.gt.0) then
-        WRITE(*,*) '    ldirect=', ldirect
-        WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
-        WRITE(*,*) '    iedate,ietime=', iedate,ietime
-        if (verbosity.gt.1) then   
-                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
-        endif     
-  endif
-
+    write(*,*) '    ldirect=', ldirect
+    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
+    write(*,*) '    iedate,ietime=', iedate,ietime
+    if (verbosity.gt.1) then   
+      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+    endif     
+  endif
 
   ! Read the age classes to be used
   !********************************
   if (verbosity.gt.0) then
-        WRITE(*,*) 'call readageclasses'
+    write(*,*) 'call readageclasses'
   endif
   call readageclasses
 
   if (verbosity.gt.1) then   
-                CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-                WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
   endif     
- 
-
 
   ! Read, which wind fields are available within the modelling period
@@ -168 +159 @@
 
   if (verbosity.gt.0) then
-        WRITE(*,*) 'call readavailable'
+    write(*,*) 'call readavailable'
   endif  
   call readavailable
@@ -177 +168 @@
  
   if (verbosity.gt.0) then
-     WRITE(*,*) 'call gridcheck'
+     write(*,*) 'call gridcheck'
   endif
 
@@ -183 +174 @@
 
   if (verbosity.gt.1) then   
-     CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-     WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
+    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
   endif      
-  
-
-  if (verbosity.gt.0) then
-        WRITE(*,*) 'call gridcheck_nests'
+
+  if (verbosity.gt.0) then
+    write(*,*) 'call gridcheck_nests'
   endif  
   call gridcheck_nests
 
-
   ! Read the output grid specifications
   !************************************
 
   if (verbosity.gt.0) then
-        WRITE(*,*) 'call readoutgrid'
+    write(*,*) 'call readoutgrid'
   endif
 
@@ -204 +193 @@
 
   if (nested_output.eq.1) then
-          call readoutgrid_nest
+    call readoutgrid_nest
     if (verbosity.gt.0) then
-        WRITE(*,*) '# readoutgrid_nest'
+      write(*,*) '# readoutgrid_nest'
     endif
   endif
@@ -232 +221 @@
   call readlanduse
 
-
   ! Assign fractional cover of landuse classes to each ECMWF grid point
   !********************************************************************
@@ -241 +229 @@
   call assignland
 
-
-
   ! Read the coordinates of the release locations
   !**********************************************
@@ -251 +237 @@
   call readreleases
 
-
   ! Read and compute surface resistances to dry deposition of gases
   !****************************************************************
@@ -267 +252 @@
     print*,'call coordtrafo'
   endif
-
 
   ! Initialize all particles to non-existent
@@ -295 +279 @@
   endif
 
-
   ! Calculate volume, surface area, etc., of all output grid cells
   ! Allocate fluxes and OHfield if necessary
   !***************************************************************
 
-
   if (verbosity.gt.0) then
     print*,'call outgrid_init'
@@ -307 +289 @@
   if (nested_output.eq.1) call outgrid_init_nest
 
-
   ! Read the OH field
   !******************
@@ -313 +294 @@
   if (OHREA.eqv..TRUE.) then
     if (verbosity.gt.0) then
-       print*,'call readOHfield'
-    endif
-       call readOHfield
+      print*,'call readOHfield'
+    endif
+    call readOHfield
   endif
 
@@ -321 +302 @@
   ! and open files that are to be kept open throughout the simulation
   !******************************************************************
-
 
   if (verbosity.gt.0) then
@@ -332 +312 @@
   !if (nested_output.eq.1) call writeheader_nest
   if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
-
   if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
   if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
 
-
-
   !open(unitdates,file=path(2)(1:length(2))//'dates')
 
@@ -346 +323 @@
   if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
 
-
   ! Releases can only start and end at discrete times (multiples of lsynctime)
   !***************************************************************************
@@ -354 +330 @@
   endif
   do i=1,numpoint
-    ireleasestart(i)=nint(real(ireleasestart(i))/ &
-         real(lsynctime))*lsynctime
-    ireleaseend(i)=nint(real(ireleaseend(i))/ &
-         real(lsynctime))*lsynctime
-  end do
-
+    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
+    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
+  end do
 
   ! Initialize cloud-base mass fluxes for the convection scheme
@@ -381 +354 @@
   end do
 
-
   ! Calculate particle trajectories
   !********************************
@@ -388 +360 @@
      if (verbosity.gt.1) then   
        CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-       WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
      endif
      if (info_flag.eq.1) then
-         print*, 'info only mode (stop)'    
-         stop
+       print*, 'info only mode (stop)'    
+       stop
      endif
      print*,'call timemanager'
@@ -399 +371 @@
   call timemanager
 
-
-  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
-       &XPART MODEL RUN!'
+  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
 
 end program flexpart

trunk/src/com_mod.f90

@@ -682 +682 @@
 
   !********************
-  ! Verbosity, testing flags
+  ! Verbosity, testing flags, namelist I/O
   !********************   
   integer :: verbosity=0
   integer :: info_flag=0
-  INTEGER :: count_clock, count_clock0,  count_rate, count_max
+  integer :: count_clock, count_clock0,  count_rate, count_max
+  logical :: nmlout=.true.
    
 

trunk/src/gridcheck.f90

@@ -443 +443 @@
   !********************
 
-  write(*,*)
-  write(*,*)
-  write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
+  write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
        nuvz+1,nwz
   write(*,*)
-  write(*,'(a)') 'Mother domain:'
+  write(*,'(a)') ' Mother domain:'
   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
        xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
+  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
        ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
   write(*,*)

trunk/src/gridcheck_gfs.f90

@@ -423 +423 @@
   write(*,*)
   write(*,*)
-  write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
+  write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
        nuvz,nwz
   write(*,*)
@@ -429 +429 @@
   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
        xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
-  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
+  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
        ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
   write(*,*)

trunk/src/gridcheck_nests.f90

@@ -350 +350 @@
   !********************
 
-  write(*,'(a,i2)') 'Nested domain #: ',l
+  write(*,'(a,i2,a)') ' Nested domain ',l,':'
   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
        xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
        '   Grid distance: ',dxn(l)
-  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
+  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
        ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
        '   Grid distance: ',dyn(l)

trunk/src/outgrid_init.f90

@@ -225 +225 @@
   !     maxspec,maxpointspec_act,nclassunc,maxageclass
 
-  write (*,*) ' Allocating fields for nested and global output (x,y): ', &
-       max(numxgrid,numxgridn),max(numygrid,numygridn) 
+  write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
+  write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn 
 
   ! allocate fields for concoutput with maximum dimension of outgrid

trunk/src/par_mod.f90

@@ -122 +122 @@
   
   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL XF
-  integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 
-  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
+  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new 
+  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
   !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
@@ -198 +198 @@
   !**************************************************
 
-  integer,parameter :: maxpart=15000000
+  integer,parameter :: maxpart=150000
   integer,parameter :: maxspec=4
 

trunk/src/readageclasses.f90

@@ -28 +28 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !                                                                            *
   !     20 March 2000                                                          *
+  !     HSO, 1 July 2014                                                       *
+  !     Added optional namelist input                                          *
   !                                                                            *
   !*****************************************************************************
@@ -46 +47 @@
   integer :: i
 
+  ! namelist help variables
+  integer :: readerror
+
+  ! namelist declaration
+  namelist /ageclass/ &
+    nageclass, &
+    lage
+
+  nageclass=-1 ! preset to negative value to identify failed namelist input
 
   ! If age spectra calculation is switched off, set number of age classes
@@ -57 +67 @@
   endif
 
-
   ! If age spectra claculation is switched on,
   ! open the AGECLASSSES file and read user options
   !************************************************
 
-  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
-       status='old',err=999)
+  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
 
-  do i=1,13
-    read(unitageclasses,*)
-  end do
-  read(unitageclasses,*) nageclass
+  ! try to read in as a namelist
+  read(unitageclasses,ageclass,iostat=readerror)
+  close(unitageclasses)
 
+  if ((nageclass.lt.0).or.(readerror.ne.0)) then
+    open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',status='old',err=999)
+    do i=1,13
+      read(unitageclasses,*)
+    end do
+    read(unitageclasses,*) nageclass
+    read(unitageclasses,*) lage(1)
+    do i=2,nageclass
+      read(unitageclasses,*) lage(i)
+    end do
+    close(unitageclasses)
+  endif
+
+  ! write ageclasses file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist',err=1000)
+    write(unitageclasses,nml=ageclass)
+    close(unitageclasses)
+  endif
 
   if (nageclass.gt.maxageclass) then
@@ -80 +106 @@
   endif
 
-  read(unitageclasses,*) lage(1)
   if (lage(1).le.0) then
     write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
@@ -89 +114 @@
 
   do i=2,nageclass
-    read(unitageclasses,*) lage(i)
     if (lage(i).le.lage(i-1)) then
       write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
@@ -105 +129 @@
   stop
 
+1000  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
+
+
 end subroutine readageclasses

trunk/src/readavailable.f90

@@ -123 +123 @@
 
   do k=1,numbnests
-  print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
-  print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
+  !print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
+  !print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     open(unitavailab,file=path(numpath+2*(k-1)+2) &
          (1:length(numpath+2*(k-1)+2)),status='old',err=998)

trunk/src/readcommand.f90

@@ -29 +29 @@
   !                                                                            *
   !     18 May 1996                                                            *
+  !     HSO, 1 July 2014                                                       *
+  !     Added optional namelist input                                          *
   !                                                                            *
   !*****************************************************************************
@@ -80 +82 @@
   character(len=50) :: line
   logical :: old
-  logical :: nml_COMMAND=.true. , nmlout=.true.  !.false.
   integer :: readerror
 
@@ -111 +112 @@
 
   ! Presetting namelist command
-  ldirect=1
+  ldirect=0
   ibdate=20000101
   ibtime=0
@@ -142 +143 @@
   ! Namelist input first: try to read as namelist file
   !**************************************************************************
-  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-         form='formatted',iostat=readerror)
-   ! If fail, check if file does not exist
-   if (readerror.ne.0) then
-     print*,'***ERROR: file COMMAND not found in ' 
-     print*, path(1)(1:length(1))//'COMMAND'
-     print*, 'Check your pathnames file.'
-     stop
-   endif
-
-   ! print error code
-   !write(*,*) 'readcommand > readerror open=' , readerror
-   !probe first line  
-   read (unitcommand,901) line
-   !write(*,*) 'index(line,COMMAND) =', index(line,'COMMAND') 
-
+  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old',form='formatted',err=999)
+
+  ! try namelist input (default)
+  read(unitcommand,command,iostat=readerror)
+  close(unitcommand)
+
+  ! distinguish namelist from fixed text input
+  if ((readerror.ne.0).or.(ldirect.eq.0)) then ! parse as text file format
  
-   !default is namelist input 
-   ! distinguish namelist from fixed text input
-   if (index(line,'COMMAND') .eq. 0) then
-   nml_COMMAND = .false. 
-   !write(*,*) 'COMMAND file does not contain the string COMMAND in the first line'      
-   endif
-   !write(*,*) 'readcommand > read as namelist? ' , nml_COMMAND
-   rewind(unitcommand)
-   read(unitcommand,command,iostat=readerror)
-
-  close(unitcommand)
-
-  !write(*,*) 'readcommand > readerror read=' , readerror
-  ! If error in namelist format, try to open with old input code
-  ! if (readerror.ne.0) then
-  ! IP 21/5/2014 the previous line cause the old long format 
-  ! to be confused with namelist input
-  
-  ! use text input 
-  if (nml_COMMAND .eqv. .false.) then 
-
-    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-         err=999)
+    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', err=999)
 
     ! Check the format of the COMMAND file (either in free format,
@@ -193 +164 @@
     if (index(line,'LDIRECT') .eq. 0) then
       old = .false.
-    !write(*,*) 'readcommand old short'
+      write(*,*) 'COMMAND in old short format, please update to namelist format'
     else
       old = .true.
-    !write(*,*) 'readcommand old long'
+      write(*,*) 'COMMAND in old long format, please update to namelist format'
     endif
     rewind(unitcommand)
@@ -206 +177 @@
     call skplin(7,unitcommand)
     if (old) call skplin(1,unitcommand)
-
     read(unitcommand,*) ldirect
     if (old) call skplin(3,unitcommand)
@@ -265 +235 @@
     write(unitcommand,nml=command)
     close(unitcommand)
-     ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
-     ! write(unitheader,NML=COMMAND)
-     !close(unitheader) 
   endif
 
@@ -616 +583 @@
 
 1000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
-       write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-        write(*,'(a)') path(2)(1:length(1))
-        stop
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(1))
+  stop
 end subroutine readcommand

trunk/src/readoutgrid.f90

@@ -29 +29 @@
   !                                                                            *
   !     4 June 1996                                                            *
+  !     HSO, 1 July 2014
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -53 +55 @@
   real,parameter :: eps=1.e-4
 
-
+  ! namelist variables
+  integer, parameter :: maxoutlev=500
+  integer :: readerror
+  real,allocatable, dimension (:) :: outheights
+
+  ! declare namelist
+  namelist /outgrid/ &
+    outlon0,outlat0, &
+    numxgrid,numygrid, &
+    dxout,dyout, &
+    outheights
+
+  ! allocate large array for reading input
+  allocate(outheights(maxoutlev),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
+
+  ! helps identifying failed namelist input
+  dxout=-1.0
+  outheights=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
-       err=999)
-
-
-  call skplin(5,unitoutgrid)
-
-
-  ! 1.  Read horizontal grid specifications
-  !****************************************
-
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlon0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlat0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numxgrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numygrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dxout
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dyout
-
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
+
+  ! try namelist input
+  read(unitoutgrid,outgrid,iostat=readerror)
+  close(unitoutgrid)
+
+  if ((dxout.le.0).or.(readerror.ne.0)) then
+
+    readerror=1
+
+    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
+
+    call skplin(5,unitoutgrid)
+
+    ! 1.  Read horizontal grid specifications
+    !****************************************
+
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlon0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlat0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numxgrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numygrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dxout
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dyout
+
+  endif
 
   ! Check validity of output grid (shall be within model domain)
@@ -91 +119 @@
   if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
        .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
-    write(*,*) 'outlon0,outlat0:'
     write(*,*) outlon0,outlat0
-    write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
     write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
     write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
@@ -104 +130 @@
   ! 2. Count Vertical levels of output grid
   !****************************************
-  j=0
-100   j=j+1
+
+  if (readerror.ne.0) then
+    j=0
+100 j=j+1
     do i=1,3
       read(unitoutgrid,*,end=99)
@@ -112 +140 @@
     if (outhelp.eq.0.) goto 99
     goto 100
-99   numzgrid=j-1
-
-    allocate(outheight(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(outheighthalf(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-
-
-  rewind(unitoutgrid)
-  call skplin(29,unitoutgrid)
+99  numzgrid=j-1
+  else
+    do i=1,maxoutlev
+      if (outheights(i).lt.0) exit
+    end do
+    numzgrid=i-1
+  end if
+
+  allocate(outheight(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
+  allocate(outheighthalf(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
 
   ! 2. Vertical levels of output grid
   !**********************************
 
-  j=0
-1000   j=j+1
-    do i=1,3
-      read(unitoutgrid,*,end=990)
-    end do
-    read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
-    if (outhelp.eq.0.) goto 99
-    outheight(j)=outhelp
-    goto 1000
-990   numzgrid=j-1
-
+  if (readerror.ne.0) then
+
+    rewind(unitoutgrid)
+    call skplin(29,unitoutgrid)
+
+    do j=1,numzgrid
+      do i=1,3
+        read(unitoutgrid,*)
+      end do
+      read(unitoutgrid,'(4x,f7.1)') outhelp
+      outheight(j)=outhelp
+      outheights(j)=outhelp
+    end do
+    close(unitoutgrid)
+
+  else
+
+    do j=1,numzgrid
+      outheight(j)=outheights(j)
+    end do
+
+  endif
+
+  ! write outgrid file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    ! reallocate outheights with actually required dimension for namelist writing
+    deallocate(outheights)
+    allocate(outheights(numzgrid),stat=stat)
+    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
+
+    do j=1,numzgrid
+      outheights(j)=outheight(j)
+    end do
+
+    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
+    write(unitoutgrid,nml=outgrid)
+    close(unitoutgrid)
+  endif
 
   ! Check whether vertical levels are specified in ascending order
@@ -160 +215 @@
   end do
 
-
   xoutshift=xlon0-outlon0
   youtshift=ylat0-outlat0
-  close(unitoutgrid)
-
-    allocate(oroout(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(area(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+
+  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
+  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
+  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
+  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
+  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
   return
 
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,*) ' #### xxx/flexpart/options                    #### '
+  write(*,'(a)') path(1)(1:length(1))
   stop
 
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
+
 end subroutine readoutgrid

trunk/src/readoutgrid_nest.f90

@@ -53 +53 @@
   real,parameter :: eps=1.e-4
 
+  integer :: readerror
 
+  ! declare namelist
+  namelist /outgridn/ &
+    outlon0n,outlat0n, &
+    numxgridn,numygridn, &
+    dxoutn,dyoutn
+
+  ! helps identifying failed namelist input
+  dxoutn=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
   !**********************************************************
 
-  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST', &
-       status='old',err=999)
+  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
 
+  ! try namelist input
+  read(unitoutgrid,outgridn,iostat=readerror)
+  close(unitoutgrid)
 
-  call skplin(5,unitoutgrid)
+  if ((dxoutn.le.0).or.(readerror.ne.0)) then
 
+    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
+    call skplin(5,unitoutgrid)
 
-  ! 1.  Read horizontal grid specifications
-  !****************************************
+    ! 1.  Read horizontal grid specifications
+    !****************************************
 
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlon0n
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlat0n
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numxgridn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numygridn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dxoutn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dyoutn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlon0n
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlat0n
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numxgridn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numygridn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dxoutn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dyoutn
 
+    close(unitoutgrid)
+  endif
 
-    allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(arean(0:numxgridn-1,0:numygridn-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+  ! write outgrid_nest file in namelist format to output directory if requested
+  if (nmlout.eqv..true.) then
+    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',err=1000)
+    write(unitoutgrid,nml=outgridn)
+    close(unitoutgrid)
+  endif
+
+  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
+  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
+  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
 
   ! Check validity of output grid (shall be within model domain)
@@ -110 +129 @@
   xoutshiftn=xlon0-outlon0n
   youtshiftn=ylat0-outlat0n
-  close(unitoutgrid)
   return
 
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(1)(1:length(1))
+  stop
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
-  write(*,*) ' #### xxx/flexpart/options                    #### '
+  write(*,'(a)') path(2)(1:length(2))
   stop
 

trunk/src/readpaths.f90

@@ -88 +88 @@
     length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
   end do
-  print*,length(5),length(6)
 
 

trunk/src/readreceptors.f90

@@ -27 +27 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !                                                                            *
   !     1 August 1996                                                          *
+  !     HSO, 14 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -51 +52 @@
   character(len=16) :: receptor
 
+  integer :: readerror
+  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
+  integer,parameter :: unitreceptorout=2
+
+  ! declare namelist
+  namelist /receptors/ &
+    receptor, lon, lat
+
+  lon=-999.9
 
   ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
@@ -64 +74 @@
   !************************************************************
 
-  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS', &
-       status='old',err=999)
+  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
 
-  call skplin(5,unitreceptor)
+  ! try namelist input
+  read(unitreceptor,receptors,iostat=readerror)
 
+  ! prepare namelist output if requested
+  if (nmlout.eqv..true.) then
+    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='new',err=1000)
+  endif
 
-  ! Read the names and coordinates of the receptors
-  !************************************************
+  if ((lon.lt.-900).or.(readerror.ne.0)) then
 
-  j=0
-100   j=j+1
+    close(unitreceptor)
+    open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
+    call skplin(5,unitreceptor)
+
+    ! Read the names and coordinates of the receptors
+    !************************************************
+
+    j=0
+100 j=j+1
     read(unitreceptor,*,end=99)
     read(unitreceptor,*,end=99)
@@ -89 +109 @@
     endif
     if (j.gt.maxreceptor) then
-    write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
-    write(*,*) ' #### POINTS ARE GIVEN.                       #### '
-    write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
-    write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
+      write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
+      write(*,*) ' #### POINTS ARE GIVEN.                       #### '
+      write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
+      write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     endif
     receptorname(j)=receptor
@@ -100 +120 @@
     ym=r_earth*dy/180.*pi
     receptorarea(j)=xm*ym
+
+    ! write receptors file in namelist format to output directory if requested
+    if (nmlout.eqv..true.) then
+      lon=x
+      lat=y
+      write(unitreceptorout,nml=receptors)
+    endif
+
     goto 100
 
-99   numreceptor=j-1
-  close(unitreceptor)
+99  numreceptor=j-1
+
+  else ! continue with namelist input
+
+    j=0
+    do while (readerror.eq.0)
+      j=j+1
+      lon=-999.9
+      read(unitreceptor,receptors,iostat=readerror)
+      if ((lon.lt.-900).or.(readerror.ne.0)) then
+        readerror=1
+      else
+        if (j.gt.maxreceptor) then
+          write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
+          write(*,*) ' #### POINTS ARE GIVEN.                       #### '
+          write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
+          write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
+        endif
+        receptorname(j)=receptor
+        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
+        yreceptor(j)=(lat-ylat0)/dy
+        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
+        ym=r_earth*dy/180.*pi
+        receptorarea(j)=xm*ym
+      endif
+
+      ! write receptors file in namelist format to output directory if requested
+      if (nmlout.eqv..true.) then
+        write(unitreceptorout,nml=receptors)
+      endif
+
+    end do
+    numreceptor=j-1
+    close(unitreceptor)
+
+  endif
+
+  if (nmlout.eqv..true.) then
+    close(unitreceptorout)
+  endif
+
   return
 
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
   write(*,'(a)') path(1)(1:length(1))
   stop
 
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
+
 end subroutine readreceptors

trunk/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
   !                                                                            *
   !*****************************************************************************
@@ -71 +73 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
   integer,parameter :: num_min_discrete=100
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
@@ -78 +80 @@
   logical :: old
 
+  ! help variables for namelist reading
+  integer :: numpoints, parts, readerror
+  integer*2 :: zkind
+  integer :: idate1, itime1, idate2, itime2
+  real :: lon1,lon2,lat1,lat2,z1,z2
+  character*40 :: comment
+  integer,parameter :: unitreleasesout=2
+  real,allocatable, dimension (:) :: mass
+  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
+
+  ! declare namelists
+  namelist /releases_ctrl/ &
+    nspec, &
+    specnum_rel
+
+  namelist /release/ &
+    idate1, itime1, &
+    idate2, itime2, &
+    lon1, lon2, &
+    lat1, lat2, &
+    z1, z2, &
+    zkind, &
+    mass, &
+    parts, &
+    comment
+
+  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
+  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', &
-       err=999)
-
-  ! 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'
-  !**************************************************************************
-
-  call skplin(12,unitreleases)
-  read (unitreleases,901) line
-901   format (a)
-  if (index(line,'Total') .eq. 0) then
-    old = .false.
-  else
-    old = .true.
-  endif
+  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
+  if (nmlout.eqv..true.) then
+    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
+  endif
+
+  if ((readerror.ne.0).or.(nspec.lt.0)) then
+
+    ! no namelist format, close file and allow reopening in old format
+    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)
+    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
+    !**************************************************************************
+
+    call skplin(12,unitreleases)
+    read (unitreleases,901) line
+901 format (a)
+    if (index(line,'Total') .eq. 0) then
+      old = .false.
+    else
+      old = .true.
+    endif
+    rewind(unitreleases)
+
+
+    ! Skip first 11 lines (file header)
+    !**********************************
+
+    call skplin(11,unitreleases)
+
+    read(unitreleases,*,err=998) nspec
+    if (old) call skplin(2,unitreleases)
+    do i=1,nspec
+      read(unitreleases,*,err=998) specnum_rel(i)
+      if (old) call skplin(2,unitreleases)
+    end do
+
+    numpoint=0
+100 numpoint=numpoint+1
+    read(unitreleases,*,end=25)
+    read(unitreleases,*,err=998,end=25) idum,idum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) idum,idum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) idum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xdum
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) idum
+    if (old) call skplin(2,unitreleases)
+    do i=1,nspec
+      read(unitreleases,*,err=998) xdum
+      if (old) call skplin(2,unitreleases)
+    end do
+    !save compoint only for the first 1000 release points
+    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
+    if (old) call skplin(1,unitreleases)
+
+    goto 100
+
+25  numpoint=numpoint-1
+
+  else 
+
+    readerror=0
+    do while (readerror.eq.0) 
+      idate1=-1
+      read(unitreleases,release,iostat=readerror)
+      if ((idate1.lt.0).or.(readerror.ne.0)) then
+        readerror=1
+      else
+        numpoint=numpoint+1
+      endif
+    end do
+    readerror=0
+  endif ! if namelist input
+
   rewind(unitreleases)
 
-
-  ! Skip first 11 lines (file header)
-  !**********************************
-
-  call skplin(11,unitreleases)
-
-
-  read(unitreleases,*,err=998) nspec
-  if (old) call skplin(2,unitreleases)
-  do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
-    if (old) call skplin(2,unitreleases)
+  ! 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(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'
+
+  write (*,*) 'Releasepoints : ', numpoint
+
+  do i=1,numpoint
+    xmasssave(i)=0.
   end do
-
-  numpoint=0
-100   numpoint=numpoint+1
-  read(unitreleases,*,end=25)
-  read(unitreleases,*,err=998,end=25) idum,idum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) idum,idum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) idum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xdum
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) idum
-  if (old) call skplin(2,unitreleases)
-  do i=1,nspec
-    read(unitreleases,*,err=998) xdum
-    if (old) call skplin(2,unitreleases)
-  end do
-  !save compoint only for the first 1000 release points
-  read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
-  if (old) call skplin(1,unitreleases)
-
-  goto 100
-
-25   numpoint=numpoint-1
-
-  !allocate memory for numpoint releaspoint
-    allocate(ireleasestart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ireleaseend(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(kindz(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmass(numpoint,maxspec) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(rho_rel(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(npart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmasssave(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-
-   write (*,*) ' Releasepoints allocated: ', numpoint
-
-   do i=1,numpoint
-     xmasssave(i)=0.
-   end do
 
   !now save the information
@@ -203 +273 @@
   end do
 
-  rewind(unitreleases)
-
-
-  ! Skip first 11 lines (file header)
-  !**********************************
-
-  call skplin(11,unitreleases)
-
-
-  ! Assign species-specific parameters needed for physical processes
-  !*************************************************************************
-
-  read(unitreleases,*,err=998) nspec
-  if (nspec.gt.maxspec) goto 994
-  if (old) call skplin(2,unitreleases)
+  if (readerror.ne.0) then
+    ! Skip first 11 lines (file header)
+    !**********************************
+
+    call skplin(11,unitreleases)
+
+    ! Assign species-specific parameters needed for physical processes
+    !*************************************************************************
+
+    read(unitreleases,*,err=998) nspec
+    if (nspec.gt.maxspec) goto 994
+    if (old) call skplin(2,unitreleases)
+  endif
+
+  ! namelist output
+  if (nmlout.eqv..true.) then
+    write(unitreleasesout,nml=releases_ctrl)
+  endif
+
   do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
-    if (old) call skplin(2,unitreleases)
-    call readspecies(specnum_rel,i)
-
-  ! For backward runs, only 1 species is allowed
-  !*********************************************
-
-  !if ((ldirect.lt.0).and.(nspec.gt.1)) then
-  !write(*,*) '#####################################################'
-  !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
-  !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
-  !write(*,*) '#####################################################'
-  !  stop
-  !endif
-
-  ! Molecular weight
-  !*****************
-
-    if (((iout.eq.2).or.(iout.eq.3)).and. &
-         (weightmolar(i).lt.0.)) then
+    if (readerror.ne.0) then
+      read(unitreleases,*,err=998) specnum_rel(i)
+      if (old) call skplin(2,unitreleases)
+      call readspecies(specnum_rel(i),i)
+    else
+      call readspecies(specnum_rel(i),i)
+    endif
+
+    ! For backward runs, only 1 species is allowed
+    !*********************************************
+
+    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
+    !write(*,*) '#####################################################'
+    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
+    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
+    !write(*,*) '#####################################################'
+    !  stop
+    !endif
+
+    ! Molecular weight
+    !*****************
+
+    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.'
@@ -246 +323 @@
     endif
 
-
-  ! Radioactive decay
-  !******************
+    ! Radioactive decay
+    !******************
 
     decay(i)=0.693147/decay(i) !conversion half life to decay constant
 
 
-  ! Dry deposition of gases
-  !************************
-
-    if (reldiff(i).gt.0.) &
-         rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
-
-  ! Dry deposition of particles
-  !****************************
+    ! Dry deposition of gases
+    !************************
+
+    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
-    if (density(i).gt.0.) then                  ! Additional parameters
-     call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
+    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)
@@ -277 +352 @@
     endif
 
-  ! Dry deposition for constant deposition velocity
-  !************************************************
+    ! Dry deposition for constant deposition velocity
+    !************************************************
 
     dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
 
-  ! Check if wet deposition or OH reaction shall be calculated
-  !***********************************************************
+    ! Check if wet deposition or OH reaction shall be calculated
+    !***********************************************************
     if (weta(i).gt.0.)  then
       WETDEP=.true.
-      write (*,*) 'Below-cloud scavenging is switched on'
+      write (*,*) 'Below-cloud scavenging: ON'
       write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
     else
-      write (*,*) 'Below-cloud scavenging is switched OFF'
+      write (*,*) 'Below-cloud scavenging: OFF'
     endif
     
-! NIK 31.01.2013 + 10.12.2013
+    ! NIK 31.01.2013 + 10.12.2013
     if (weta_in(i).gt.0.)  then
       WETDEP=.true.
-      write (*,*) 'In-cloud scavenging is switched on'
+      write (*,*) 'In-cloud scavenging: ON'
       write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
     else
-      write (*,*) 'In-cloud scavenging is switched OFF' 
+      write (*,*) 'In-cloud scavenging: OFF' 
     endif
 
     if (ohreact(i).gt.0) then
       OHREA=.true.
-      write (*,*) 'OHreaction switched on: ',ohreact(i),i
-    endif
-
-
-    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
-         (dryvel(i).gt.0.)) then
+      write (*,*) 'OHreaction: ON (',ohreact(i),i,')'
+    endif
+
+    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
       DRYDEP=.true.
       DRYDEPSPEC(i)=.true.
@@ -315 +388 @@
   end do
 
-    if (WETDEP.or.DRYDEP) DEP=.true.
+  if (WETDEP.or.DRYDEP) DEP=.true.
 
   ! Read specifications for each release point
   !*******************************************
-
+  numpoints=numpoint
   numpoint=0
   numpartmax=0
   releaserate=0.
-1000   numpoint=numpoint+1
-  read(unitreleases,*,end=250)
-  read(unitreleases,*,err=998,end=250) id1,it1
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) id2,it2
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xpoint1(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) ypoint1(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) xpoint2(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) ypoint2(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) kindz(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) zpoint1(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) zpoint2(numpoint)
-  if (old) call skplin(2,unitreleases)
-  read(unitreleases,*,err=998) npart(numpoint)
-  if (old) call skplin(2,unitreleases)
-  do i=1,nspec
-    read(unitreleases,*,err=998) xmass(numpoint,i)
-    if (old) call skplin(2,unitreleases)
-  end do
-  !save compoint only for the first 1000 release points
-  if (numpoint.le.1000) then
-    read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
+101   numpoint=numpoint+1
+
+  if (readerror.lt.1) then ! reading namelist format
+
+    if (numpoint.gt.numpoints) goto 250
+    zkind = 1
+    mass = 0
+    parts = 0
+    comment = ' '
+    read(unitreleases,release,iostat=readerror)
+    id1=idate1
+    it1=itime1
+    id2=idate2
+    it2=itime2
+    xpoint1(numpoint)=lon1
+    xpoint2(numpoint)=lon2
+    ypoint1(numpoint)=lat1
+    ypoint2(numpoint)=lat2
+    zpoint1(numpoint)=z1
+    zpoint2(numpoint)=z2
+    kindz(numpoint)=zkind
+    do i=1,nspec
+      xmass(numpoint,i)=mass(i)
+    end do
+    npart(numpoint)=parts
+    compoint(min(1001,numpoint))=comment
+
+    ! namelist output
+    if (nmlout.eqv..true.) then
+      write(unitreleasesout,nml=release)
+    endif
+
   else
-    read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
-  endif
-  if (old) call skplin(1,unitreleases)
-  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
-  else
-    if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
-         (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
-  endif
-
+
+    read(unitreleases,*,end=250)
+    read(unitreleases,*,err=998,end=250) id1,it1
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) id2,it2
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xpoint1(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) ypoint1(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) xpoint2(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) ypoint2(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) kindz(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) zpoint1(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) zpoint2(numpoint)
+    if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) npart(numpoint)
+    if (old) call skplin(2,unitreleases)
+    do i=1,nspec
+      read(unitreleases,*,err=998) xmass(numpoint,i)
+      if (old) call skplin(2,unitreleases)
+      mass(i)=xmass(numpoint,i)
+    end do
+    !save compoint only for the first 1000 release points
+    if (numpoint.le.1000) then
+      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
+      comment=compoint(numpoint)(1:40)
+    else
+      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
+      comment=compoint(1001)(1:40)
+    endif
+    if (old) call skplin(1,unitreleases)
+
+    ! namelist output
+    if (nmlout.eqv..true.) then
+      idate1=id1
+      itime1=it1
+      idate2=id2
+      itime2=it2
+      lon1=xpoint1(numpoint)
+      lon2=xpoint2(numpoint)
+      lat1=ypoint1(numpoint)
+      lat2=ypoint2(numpoint)
+      z1=zpoint1(numpoint)
+      z2=zpoint2(numpoint)
+      zkind=kindz(numpoint)
+      parts=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
+    else
+      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
+           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
+    endif
+
+  endif ! if namelist format
 
   ! If a release point contains no particles, stop and issue error message
@@ -435 +563 @@
   endif
 
-
-  ! Check, whether the total number of particles may exceed totally allowed
-  ! number of particles at some time during the simulation
-  !************************************************************************
-
   ! Determine the release rate (particles per second) and total number
   ! of particles released during the simulation
@@ -451 +574 @@
   endif
   numpartmax=numpartmax+npart(numpoint)
-  goto 1000
-
+  goto 101
 
 250   close(unitreleases)
 
-  write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
+  if (nmlout.eqv..true.) then
+    close(unitreleasesout)
+  endif
+
+  write (*,*) 'Particles allocated (maxpart)  : ',maxpart
+  write (*,*) 'Particles released (numpartmax): ',numpartmax
   numpoint=numpoint-1
 
@@ -464 +591 @@
     maxpointspec_act=1
   endif
+
+  ! Check, whether the total number of particles may exceed totally allowed
+  ! number of particles at some time during the simulation
+  !************************************************************************
 
   if (releaserate.gt. &
@@ -521 +652 @@
   stop
 
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
+
 end subroutine readreleases

trunk/src/readspecies.f90

@@ -30 +30 @@
   !                                                                            *
   !   11 July 1996                                                             *
-  !                                                                            *
+  !
   !   Changes:                                                                 *
   !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
+  !                                                                            *
+  !   HSO, 13 August 2013
+  !   added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -62 +65 @@
   logical :: spec_found
 
+  character(len=16) :: pspecies
+  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
+  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
+  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
+  integer :: readerror
+
+  ! declare namelist
+  namelist /species_params/ &
+   pspecies, pdecay, pweta, pwetb, &
+   pweta_in, pwetb_in, pwetc_in, pwetd_in, &
+   preldiff, phenry, pf0, pdensity, pdquer, &
+   pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
+
+  pspecies=" "
+  pdecay=-999.9
+  pweta=-9.9E-09
+  pwetb=0.0
+  pweta_in=-9.9E-09
+  pwetb_in=-9.9E-09
+  pwetc_in=-9.9E-09
+  pwetd_in=-9.9E-09
+  preldiff=-9.9
+  phenry=0.0
+  pf0=0.0
+  pdensity=-9.9E09
+  pdquer=0.0
+  pdsigma=0.0
+  pdryvel=-9.99
+  pohreact=-9.9E-09
+  pspec_ass=-9
+  pkao=-99.99
+  pweightmolar=-789.0 ! read failure indicator value
+
   ! Open the SPECIES file and read species names and properties
   !************************************************************
   specnum(pos_spec)=id_spec
   write(aspecnumb,'(i3.3)') specnum(pos_spec)
-  open(unitspecies,file= &
-       path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
-       err=998)
+  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
   !write(*,*) 'reading SPECIES',specnum(pos_spec)
 
   ASSSPEC=.FALSE.
 
-  do i=1,6
-    read(unitspecies,*)
-  end do
+  ! try namelist input
+  read(unitspecies,species_params,iostat=readerror)
+  close(unitspecies)
+   
+  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
+
+    readerror=1
+
+    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
+
+    do i=1,6
+      read(unitspecies,*)
+    end do
 
     read(unitspecies,'(a10)',end=22) species(pos_spec)
@@ -86 +130 @@
   !  write(*,*) wetb(pos_spec)
 
-!*** NIK 31.01.2013: including in-cloud scavening parameters
+  !*** NIK 31.01.2013: including in-cloud scavening parameters
    read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
   !  write(*,*) weta_in(pos_spec)
@@ -118 +162 @@
     read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
   !       write(*,*) kao(pos_spec)
-    i=pos_spec
-
-    if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
-       if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
+
+    pspecies=species(pos_spec)
+    pdecay=decay(pos_spec)
+    pweta=weta(pos_spec)
+    pwetb=wetb(pos_spec)
+    pweta_in=weta_in(pos_spec)
+    pwetb_in=wetb_in(pos_spec)
+    pwetc_in=wetc_in(pos_spec)
+    pwetd_in=wetd_in(pos_spec)
+    preldiff=reldiff(pos_spec)
+    phenry=henry(pos_spec)
+    pf0=f0(pos_spec)
+    pdensity=density(pos_spec)
+    pdquer=dquer(pos_spec)
+    pdsigma=dsigma(pos_spec)
+    pdryvel=dryvel(pos_spec)
+    pweightmolar=weightmolar(pos_spec)
+    pohreact=ohreact(pos_spec)
+    pspec_ass=spec_ass(pos_spec)
+    pkao=kao(pos_spec)
+
+  else
+
+    species(pos_spec)=pspecies
+    decay(pos_spec)=pdecay
+    weta(pos_spec)=pweta
+    wetb(pos_spec)=pwetb
+    weta_in(pos_spec)=pweta_in
+    wetb_in(pos_spec)=pwetb_in
+    wetc_in(pos_spec)=pwetc_in
+    wetd_in(pos_spec)=pwetd_in
+    reldiff(pos_spec)=preldiff
+    henry(pos_spec)=phenry
+    f0(pos_spec)=pf0
+    density(pos_spec)=pdensity
+    dquer(pos_spec)=pdquer
+    dsigma(pos_spec)=pdsigma
+    dryvel(pos_spec)=pdryvel
+    weightmolar(pos_spec)=pweightmolar
+    ohreact(pos_spec)=pohreact
+    spec_ass(pos_spec)=pspec_ass
+    kao(pos_spec)=pkao
+
+  endif
+
+  i=pos_spec
+
+  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
+   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
+  endif
+
+  if (spec_ass(pos_spec).gt.0) then
+    spec_found=.FALSE.
+    do j=1,pos_spec-1
+      if (spec_ass(pos_spec).eq.specnum(j)) then
+        spec_ass(pos_spec)=j
+        spec_found=.TRUE.
+        ASSSPEC=.TRUE.
+      endif
+    end do
+    if (spec_found.eqv..FALSE.) then
+      goto 997
     endif
-
-    if (spec_ass(pos_spec).gt.0) then
-       spec_found=.FALSE.
-       do j=1,pos_spec-1
-          if (spec_ass(pos_spec).eq.specnum(j)) then
-             spec_ass(pos_spec)=j
-             spec_found=.TRUE.
-             ASSSPEC=.TRUE.
-          endif
-       end do
-       if (spec_found.eqv..FALSE.) then
-          goto 997
-       endif
-    endif
-
-    if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
-    if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
-
-    if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
-      write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
-      write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
-      write(*,*) '#### PARTICLE AND GAS.                       ####'
-      write(*,*) '#### SPECIES NUMBER',aspecnumb
-      stop
-    endif
+  endif
+
+  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
+  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
+
+  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
+    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
+    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
+    write(*,*) '#### PARTICLE AND GAS.                       ####'
+    write(*,*) '#### SPECIES NUMBER',aspecnumb
+    stop
+  endif
 20   continue
 
-
-  ! Read in daily and day-of-week variation of emissions, if available
-  !*******************************************************************
+  if (readerror.ne.0) then ! text format input
+
+    ! Read in daily and day-of-week variation of emissions, if available
+    !*******************************************************************
+    ! HSO: This is not yet implemented as namelist parameters
 
     do j=1,24           ! initialize everything to no variation
@@ -172 +263 @@
     end do
 
-22   close(unitspecies)
-
-   return
+  endif
+
+22 close(unitspecies)
+
+  ! namelist output if requested
+  if (nmlout.eqv..true.) then
+    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='new',err=1000)
+    write(unitspecies,nml=species_params)
+    close(unitspecies)
+  endif
+
+  return
 
 996   write(*,*) '#####################################################'
@@ -203 +303 @@
   stop
 
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
+  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+  write(*,'(a)') path(2)(1:length(2))
+  stop
 
 end subroutine readspecies

trunk/src/timemanager.f90

@@ -129 +129 @@
 
 
+  !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
+  write(*,46) float(itime)/3600,itime,numpart
   if (verbosity.gt.0) then
     write (*,*) 'timemanager> starting simulation'
@@ -371 +373 @@
         if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
         if (iflux.eq.1) call fluxoutput(itime)
-        write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
-             drygridtotalunc
-45      format(i9,' SECONDS SIMULATED: ',i8, &
-             ' PARTICLES:    Uncertainty: ',3f7.3)
+        !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
+        write(*,46) float(itime)/3600,itime,numpart
+45      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
+46      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
         if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
         loutnext=loutnext+loutstep