Changeset 27


Ignore:
Timestamp:
Jul 3, 2014, 2:55:50 PM (5 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:
17 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/FLEXPART.f90

    r24 r27  
    6060  call gasdev1(idummy,rannumb(maxrand),rannumb(maxrand-1))
    6161
    62   !
    63   flexversion='Version 9.2 beta (2014-05-23)'
    64   !verbosity=0
     62  ! FLEXPART version string
     63  flexversion='Version 9.2 beta (2014-07-01)'
     64  verbosity=0
     65
    6566  ! Read the pathnames where input/output files are stored
    6667  !*******************************************************
     
    7677    call getarg(1,arg1)
    7778    pathfile=arg1
    78     verbosity=0
    7979    if (arg1(1:1).eq.'-') then
    80         write(pathfile,'(a11)') './pathnames'
    81         inline_options=arg1
     80      write(pathfile,'(a11)') './pathnames'
     81      inline_options=arg1
    8282    endif
    8383  case (0)
    8484    write(pathfile,'(a11)') './pathnames'
    85     verbosity=0
    8685  end select
    8786 
    88     if (inline_options(1:1).eq.'-') then
    89       print*, 'inline options=', inline_options       
    90       if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
    91          print*, 'verbose mode 1: additional information will be displayed'
    92          verbosity=1
    93       endif
    94       if (trim(inline_options).eq.'-v2') then
    95          print*, 'verbose mode 2: additional information will be displayed'
    96          verbosity=2
    97       endif
    98       if (trim(inline_options).eq.'-i') then
    99          print*, 'info mode: will provide run specific information and stop'
    100          verbosity=1
    101          info_flag=1
    102       endif
    103       if (trim(inline_options).eq.'-i2') then
    104          print*, 'info mode: will provide run specific information and stop'
    105          verbosity=2
    106          info_flag=1
    107       endif
    108     endif
    109 
    110 
    11187  ! Print the GPL License statement
    11288  !*******************************************************
    11389  print*,'Welcome to FLEXPART ', trim(flexversion)
    114   print*,'FLEXPART is free software released under the GNU Genera'// &
    115        'l Public License.'
    116            
    117   if (verbosity.gt.0) then
    118         WRITE(*,*) 'call readpaths'
     90  print*,'FLEXPART is free software released under the GNU General Public License.'
     91 
     92  if (inline_options(1:1).eq.'-') then
     93    if (trim(inline_options).eq.'-v'.or.trim(inline_options).eq.'-v1') then
     94       print*, 'Verbose mode 1: display detailed information during run'
     95       verbosity=1
     96    endif
     97    if (trim(inline_options).eq.'-v2') then
     98       print*, 'Verbose mode 2: display more detailed information during run'
     99       verbosity=2
     100    endif
     101    if (trim(inline_options).eq.'-i') then
     102       print*, 'Info mode: provide detailed run specific information and stop'
     103       verbosity=1
     104       info_flag=1
     105    endif
     106    if (trim(inline_options).eq.'-i2') then
     107       print*, 'Info mode: provide more detailed run specific information and stop'
     108       verbosity=2
     109       info_flag=1
     110    endif
     111  endif
     112           
     113  if (verbosity.gt.0) then
     114    write(*,*) 'call readpaths'
    119115  endif
    120116  call readpaths(pathfile)
    121    
    122117 
    123118  if (verbosity.gt.1) then !show clock info
     
    131126  endif
    132127
    133 
    134128  ! Read the user specifications for the current model run
    135129  !*******************************************************
    136130
    137131  if (verbosity.gt.0) then
    138         WRITE(*,*) 'call readcommand'
     132    write(*,*) 'call readcommand'
    139133  endif
    140134  call readcommand
    141135  if (verbosity.gt.0) then
    142         WRITE(*,*) '    ldirect=', ldirect
    143         WRITE(*,*) '    ibdate,ibtime=',ibdate,ibtime
    144         WRITE(*,*) '    iedate,ietime=', iedate,ietime
    145         if (verbosity.gt.1) then   
    146                 CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    147                 WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    148         endif     
    149   endif
    150 
     136    write(*,*) '    ldirect=', ldirect
     137    write(*,*) '    ibdate,ibtime=',ibdate,ibtime
     138    write(*,*) '    iedate,ietime=', iedate,ietime
     139    if (verbosity.gt.1) then   
     140      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     141      write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     142    endif     
     143  endif
    151144
    152145  ! Read the age classes to be used
    153146  !********************************
    154147  if (verbosity.gt.0) then
    155         WRITE(*,*) 'call readageclasses'
     148    write(*,*) 'call readageclasses'
    156149  endif
    157150  call readageclasses
    158151
    159152  if (verbosity.gt.1) then   
    160                 CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    161                 WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     153    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     154    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    162155  endif     
    163  
    164 
    165156
    166157  ! Read, which wind fields are available within the modelling period
     
    168159
    169160  if (verbosity.gt.0) then
    170         WRITE(*,*) 'call readavailable'
     161    write(*,*) 'call readavailable'
    171162  endif 
    172163  call readavailable
     
    177168 
    178169  if (verbosity.gt.0) then
    179      WRITE(*,*) 'call gridcheck'
     170     write(*,*) 'call gridcheck'
    180171  endif
    181172
     
    183174
    184175  if (verbosity.gt.1) then   
    185      CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    186      WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     176    CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
     177    write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    187178  endif     
    188  
    189 
    190   if (verbosity.gt.0) then
    191         WRITE(*,*) 'call gridcheck_nests'
     179
     180  if (verbosity.gt.0) then
     181    write(*,*) 'call gridcheck_nests'
    192182  endif 
    193183  call gridcheck_nests
    194184
    195 
    196185  ! Read the output grid specifications
    197186  !************************************
    198187
    199188  if (verbosity.gt.0) then
    200         WRITE(*,*) 'call readoutgrid'
     189    write(*,*) 'call readoutgrid'
    201190  endif
    202191
     
    204193
    205194  if (nested_output.eq.1) then
    206           call readoutgrid_nest
     195    call readoutgrid_nest
    207196    if (verbosity.gt.0) then
    208         WRITE(*,*) '# readoutgrid_nest'
     197      write(*,*) '# readoutgrid_nest'
    209198    endif
    210199  endif
     
    232221  call readlanduse
    233222
    234 
    235223  ! Assign fractional cover of landuse classes to each ECMWF grid point
    236224  !********************************************************************
     
    241229  call assignland
    242230
    243 
    244 
    245231  ! Read the coordinates of the release locations
    246232  !**********************************************
     
    251237  call readreleases
    252238
    253 
    254239  ! Read and compute surface resistances to dry deposition of gases
    255240  !****************************************************************
     
    267252    print*,'call coordtrafo'
    268253  endif
    269 
    270254
    271255  ! Initialize all particles to non-existent
     
    295279  endif
    296280
    297 
    298281  ! Calculate volume, surface area, etc., of all output grid cells
    299282  ! Allocate fluxes and OHfield if necessary
    300283  !***************************************************************
    301284
    302 
    303285  if (verbosity.gt.0) then
    304286    print*,'call outgrid_init'
     
    307289  if (nested_output.eq.1) call outgrid_init_nest
    308290
    309 
    310291  ! Read the OH field
    311292  !******************
     
    313294  if (OHREA.eqv..TRUE.) then
    314295    if (verbosity.gt.0) then
    315        print*,'call readOHfield'
    316     endif
    317        call readOHfield
     296      print*,'call readOHfield'
     297    endif
     298    call readOHfield
    318299  endif
    319300
     
    321302  ! and open files that are to be kept open throughout the simulation
    322303  !******************************************************************
    323 
    324304
    325305  if (verbosity.gt.0) then
     
    332312  !if (nested_output.eq.1) call writeheader_nest
    333313  if (nested_output.eq.1.and.surf_only.ne.1) call writeheader_nest
    334 
    335314  if (nested_output.eq.1.and.surf_only.eq.1) call writeheader_nest_surf
    336315  if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
    337316
    338 
    339 
    340317  !open(unitdates,file=path(2)(1:length(2))//'dates')
    341318
     
    346323  if ((iout.eq.4).or.(iout.eq.5)) call openouttraj
    347324
    348 
    349325  ! Releases can only start and end at discrete times (multiples of lsynctime)
    350326  !***************************************************************************
     
    354330  endif
    355331  do i=1,numpoint
    356     ireleasestart(i)=nint(real(ireleasestart(i))/ &
    357          real(lsynctime))*lsynctime
    358     ireleaseend(i)=nint(real(ireleaseend(i))/ &
    359          real(lsynctime))*lsynctime
    360   end do
    361 
     332    ireleasestart(i)=nint(real(ireleasestart(i))/real(lsynctime))*lsynctime
     333    ireleaseend(i)=nint(real(ireleaseend(i))/real(lsynctime))*lsynctime
     334  end do
    362335
    363336  ! Initialize cloud-base mass fluxes for the convection scheme
     
    381354  end do
    382355
    383 
    384356  ! Calculate particle trajectories
    385357  !********************************
     
    388360     if (verbosity.gt.1) then   
    389361       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
    390        WRITE(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
     362       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
    391363     endif
    392364     if (info_flag.eq.1) then
    393          print*, 'info only mode (stop)'   
    394          stop
     365       print*, 'info only mode (stop)'   
     366       stop
    395367     endif
    396368     print*,'call timemanager'
     
    399371  call timemanager
    400372
    401 
    402   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLE&
    403        &XPART MODEL RUN!'
     373  write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
    404374
    405375end program flexpart
  • trunk/src/com_mod.f90

    r24 r27  
    682682
    683683  !********************
    684   ! Verbosity, testing flags
     684  ! Verbosity, testing flags, namelist I/O
    685685  !********************   
    686686  integer :: verbosity=0
    687687  integer :: info_flag=0
    688   INTEGER :: count_clock, count_clock0,  count_rate, count_max
     688  integer :: count_clock, count_clock0,  count_rate, count_max
     689  logical :: nmlout=.true.
    689690   
    690691
  • trunk/src/gridcheck.f90

    r24 r27  
    443443  !********************
    444444
    445   write(*,*)
    446   write(*,*)
    447   write(*,'(a,2i7)') '# of vertical levels in ECMWF data: ', &
     445  write(*,'(a,2i7)') ' Vertical levels in ECMWF data: ', &
    448446       nuvz+1,nwz
    449447  write(*,*)
    450   write(*,'(a)') 'Mother domain:'
     448  write(*,'(a)') ' Mother domain:'
    451449  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    452450       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    453   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
     451  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
    454452       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    455453  write(*,*)
  • trunk/src/gridcheck_gfs.f90

    r4 r27  
    423423  write(*,*)
    424424  write(*,*)
    425   write(*,'(a,2i7)') '# of vertical levels in NCEP data: ', &
     425  write(*,'(a,2i7)') 'Vertical levels in NCEP data: ', &
    426426       nuvz,nwz
    427427  write(*,*)
     
    429429  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Longitude range: ', &
    430430       xlon0,' to ',xlon0+(nx-1)*dx,'   Grid distance: ',dx
    431   write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range: ', &
     431  write(*,'(a,f10.2,a1,f10.2,a,f10.2)') '  Latitude range : ', &
    432432       ylat0,' to ',ylat0+(ny-1)*dy,'   Grid distance: ',dy
    433433  write(*,*)
  • trunk/src/gridcheck_nests.f90

    r24 r27  
    350350  !********************
    351351
    352   write(*,'(a,i2)') 'Nested domain #: ',l
     352  write(*,'(a,i2,a)') ' Nested domain ',l,':'
    353353  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Longitude range: ', &
    354354       xlon0n(l),' to ',xlon0n(l)+(nxn(l)-1)*dxn(l), &
    355355       '   Grid distance: ',dxn(l)
    356   write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range: ', &
     356  write(*,'(a,f10.5,a,f10.5,a,f10.5)') '  Latitude range : ', &
    357357       ylat0n(l),' to ',ylat0n(l)+(nyn(l)-1)*dyn(l), &
    358358       '   Grid distance: ',dyn(l)
  • trunk/src/outgrid_init.f90

    r20 r27  
    225225  !     maxspec,maxpointspec_act,nclassunc,maxageclass
    226226
    227   write (*,*) ' Allocating fields for nested and global output (x,y): ', &
    228        max(numxgrid,numxgridn),max(numygrid,numygridn)
     227  write (*,*) 'Allocating fields for global output (x,y): ', numxgrid,numygrid
     228  write (*,*) 'Allocating fields for nested output (x,y): ', numxgridn,numygridn
    229229
    230230  ! allocate fields for concoutput with maximum dimension of outgrid
  • trunk/src/par_mod.f90

    r24 r27  
    122122 
    123123  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !FNL XF
    124   integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new
    125   !integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
     124  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=152,nwzmax=152,nzmax=152 !ECMWF new
     125  integer,parameter :: nxmax=361,nymax=181,nuvzmax=92,nwzmax=92,nzmax=92 !ECMWF
    126126  !integer,parameter :: nxmax=361,nymax=181,nuvzmax=26,nwzmax=26,nzmax=26
    127127  !integer,parameter :: nxmax=721,nymax=361,nuvzmax=64,nwzmax=64,nzmax=64
     
    198198  !**************************************************
    199199
    200   integer,parameter :: maxpart=15000000
     200  integer,parameter :: maxpart=150000
    201201  integer,parameter :: maxspec=4
    202202
  • trunk/src/readageclasses.f90

    r4 r27  
    2828  !                                                                            *
    2929  !     Author: A. Stohl                                                       *
    30   !                                                                            *
    3130  !     20 March 2000                                                          *
     31  !     HSO, 1 July 2014                                                       *
     32  !     Added optional namelist input                                          *
    3233  !                                                                            *
    3334  !*****************************************************************************
     
    4647  integer :: i
    4748
     49  ! namelist help variables
     50  integer :: readerror
     51
     52  ! namelist declaration
     53  namelist /ageclass/ &
     54    nageclass, &
     55    lage
     56
     57  nageclass=-1 ! preset to negative value to identify failed namelist input
    4858
    4959  ! If age spectra calculation is switched off, set number of age classes
     
    5767  endif
    5868
    59 
    6069  ! If age spectra claculation is switched on,
    6170  ! open the AGECLASSSES file and read user options
    6271  !************************************************
    6372
    64   open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES', &
    65        status='old',err=999)
     73  open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',form='formatted',status='old',err=999)
    6674
    67   do i=1,13
    68     read(unitageclasses,*)
    69   end do
    70   read(unitageclasses,*) nageclass
     75  ! try to read in as a namelist
     76  read(unitageclasses,ageclass,iostat=readerror)
     77  close(unitageclasses)
    7178
     79  if ((nageclass.lt.0).or.(readerror.ne.0)) then
     80    open(unitageclasses,file=path(1)(1:length(1))//'AGECLASSES',status='old',err=999)
     81    do i=1,13
     82      read(unitageclasses,*)
     83    end do
     84    read(unitageclasses,*) nageclass
     85    read(unitageclasses,*) lage(1)
     86    do i=2,nageclass
     87      read(unitageclasses,*) lage(i)
     88    end do
     89    close(unitageclasses)
     90  endif
     91
     92  ! write ageclasses file in namelist format to output directory if requested
     93  if (nmlout.eqv..true.) then
     94    open(unitageclasses,file=path(2)(1:length(2))//'AGECLASSES.namelist',err=1000)
     95    write(unitageclasses,nml=ageclass)
     96    close(unitageclasses)
     97  endif
    7298
    7399  if (nageclass.gt.maxageclass) then
     
    80106  endif
    81107
    82   read(unitageclasses,*) lage(1)
    83108  if (lage(1).le.0) then
    84109    write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
     
    89114
    90115  do i=2,nageclass
    91     read(unitageclasses,*) lage(i)
    92116    if (lage(i).le.lage(i-1)) then
    93117      write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
     
    105129  stop
    106130
     1311000  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "AGECLASSES" #### '
     132  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     133  write(*,'(a)') path(2)(1:length(2))
     134  stop
     135
     136
    107137end subroutine readageclasses
  • trunk/src/readavailable.f90

    r20 r27  
    123123
    124124  do k=1,numbnests
    125   print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
    126   print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
     125  !print*,length(numpath+2*(k-1)+1),length(numpath+2*(k-1)+2),length(4),length(3)
     126  !print*,path(numpath+2*(k-1)+2)(1:length(numpath+2*(k-1)+2))
    127127    open(unitavailab,file=path(numpath+2*(k-1)+2) &
    128128         (1:length(numpath+2*(k-1)+2)),status='old',err=998)
  • trunk/src/readcommand.f90

    r24 r27  
    2929  !                                                                            *
    3030  !     18 May 1996                                                            *
     31  !     HSO, 1 July 2014                                                       *
     32  !     Added optional namelist input                                          *
    3133  !                                                                            *
    3234  !*****************************************************************************
     
    8082  character(len=50) :: line
    8183  logical :: old
    82   logical :: nml_COMMAND=.true. , nmlout=.true.  !.false.
    8384  integer :: readerror
    8485
     
    111112
    112113  ! Presetting namelist command
    113   ldirect=1
     114  ldirect=0
    114115  ibdate=20000101
    115116  ibtime=0
     
    142143  ! Namelist input first: try to read as namelist file
    143144  !**************************************************************************
    144   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    145          form='formatted',iostat=readerror)
    146    ! If fail, check if file does not exist
    147    if (readerror.ne.0) then
    148      print*,'***ERROR: file COMMAND not found in '
    149      print*, path(1)(1:length(1))//'COMMAND'
    150      print*, 'Check your pathnames file.'
    151      stop
    152    endif
    153 
    154    ! print error code
    155    !write(*,*) 'readcommand > readerror open=' , readerror
    156    !probe first line 
    157    read (unitcommand,901) line
    158    !write(*,*) 'index(line,COMMAND) =', index(line,'COMMAND')
    159 
     145  open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old',form='formatted',err=999)
     146
     147  ! try namelist input (default)
     148  read(unitcommand,command,iostat=readerror)
     149  close(unitcommand)
     150
     151  ! distinguish namelist from fixed text input
     152  if ((readerror.ne.0).or.(ldirect.eq.0)) then ! parse as text file format
    160153 
    161    !default is namelist input
    162    ! distinguish namelist from fixed text input
    163    if (index(line,'COMMAND') .eq. 0) then
    164    nml_COMMAND = .false.
    165    !write(*,*) 'COMMAND file does not contain the string COMMAND in the first line'     
    166    endif
    167    !write(*,*) 'readcommand > read as namelist? ' , nml_COMMAND
    168    rewind(unitcommand)
    169    read(unitcommand,command,iostat=readerror)
    170 
    171   close(unitcommand)
    172 
    173   !write(*,*) 'readcommand > readerror read=' , readerror
    174   ! If error in namelist format, try to open with old input code
    175   ! if (readerror.ne.0) then
    176   ! IP 21/5/2014 the previous line cause the old long format
    177   ! to be confused with namelist input
    178  
    179   ! use text input
    180   if (nml_COMMAND .eqv. .false.) then
    181 
    182     open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
    183          err=999)
     154    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', err=999)
    184155
    185156    ! Check the format of the COMMAND file (either in free format,
     
    193164    if (index(line,'LDIRECT') .eq. 0) then
    194165      old = .false.
    195     !write(*,*) 'readcommand old short'
     166      write(*,*) 'COMMAND in old short format, please update to namelist format'
    196167    else
    197168      old = .true.
    198     !write(*,*) 'readcommand old long'
     169      write(*,*) 'COMMAND in old long format, please update to namelist format'
    199170    endif
    200171    rewind(unitcommand)
     
    206177    call skplin(7,unitcommand)
    207178    if (old) call skplin(1,unitcommand)
    208 
    209179    read(unitcommand,*) ldirect
    210180    if (old) call skplin(3,unitcommand)
     
    265235    write(unitcommand,nml=command)
    266236    close(unitcommand)
    267      ! open(unitheader,file=path(2)(1:length(2))//'header_nml',status='new',err=999)
    268      ! write(unitheader,NML=COMMAND)
    269      !close(unitheader)
    270237  endif
    271238
     
    616583
    6175841000   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "COMMAND"    #### '
    618        write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    619         write(*,'(a)') path(2)(1:length(1))
    620         stop
     585  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     586  write(*,'(a)') path(2)(1:length(1))
     587  stop
    621588end subroutine readcommand
  • trunk/src/readoutgrid.f90

    r20 r27  
    2929  !                                                                            *
    3030  !     4 June 1996                                                            *
     31  !     HSO, 1 July 2014
     32  !     Added optional namelist input
    3133  !                                                                            *
    3234  !*****************************************************************************
     
    5355  real,parameter :: eps=1.e-4
    5456
    55 
     57  ! namelist variables
     58  integer, parameter :: maxoutlev=500
     59  integer :: readerror
     60  real,allocatable, dimension (:) :: outheights
     61
     62  ! declare namelist
     63  namelist /outgrid/ &
     64    outlon0,outlat0, &
     65    numxgrid,numygrid, &
     66    dxout,dyout, &
     67    outheights
     68
     69  ! allocate large array for reading input
     70  allocate(outheights(maxoutlev),stat=stat)
     71  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
     72
     73  ! helps identifying failed namelist input
     74  dxout=-1.0
     75  outheights=-1.0
    5676
    5777  ! Open the OUTGRID file and read output grid specifications
    5878  !**********************************************************
    5979
    60   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old', &
    61        err=999)
    62 
    63 
    64   call skplin(5,unitoutgrid)
    65 
    66 
    67   ! 1.  Read horizontal grid specifications
    68   !****************************************
    69 
    70   call skplin(3,unitoutgrid)
    71   read(unitoutgrid,'(4x,f11.4)') outlon0
    72   call skplin(3,unitoutgrid)
    73   read(unitoutgrid,'(4x,f11.4)') outlat0
    74   call skplin(3,unitoutgrid)
    75   read(unitoutgrid,'(4x,i5)') numxgrid
    76   call skplin(3,unitoutgrid)
    77   read(unitoutgrid,'(4x,i5)') numygrid
    78   call skplin(3,unitoutgrid)
    79   read(unitoutgrid,'(4x,f12.5)') dxout
    80   call skplin(3,unitoutgrid)
    81   read(unitoutgrid,'(4x,f12.5)') dyout
    82 
     80  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',form='formatted',err=999)
     81
     82  ! try namelist input
     83  read(unitoutgrid,outgrid,iostat=readerror)
     84  close(unitoutgrid)
     85
     86  if ((dxout.le.0).or.(readerror.ne.0)) then
     87
     88    readerror=1
     89
     90    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID',status='old',err=999)
     91
     92    call skplin(5,unitoutgrid)
     93
     94    ! 1.  Read horizontal grid specifications
     95    !****************************************
     96
     97    call skplin(3,unitoutgrid)
     98    read(unitoutgrid,'(4x,f11.4)') outlon0
     99    call skplin(3,unitoutgrid)
     100    read(unitoutgrid,'(4x,f11.4)') outlat0
     101    call skplin(3,unitoutgrid)
     102    read(unitoutgrid,'(4x,i5)') numxgrid
     103    call skplin(3,unitoutgrid)
     104    read(unitoutgrid,'(4x,i5)') numygrid
     105    call skplin(3,unitoutgrid)
     106    read(unitoutgrid,'(4x,f12.5)') dxout
     107    call skplin(3,unitoutgrid)
     108    read(unitoutgrid,'(4x,f12.5)') dyout
     109
     110  endif
    83111
    84112  ! Check validity of output grid (shall be within model domain)
     
    91119  if ((outlon0+eps.lt.xlon0).or.(outlat0+eps.lt.ylat0) &
    92120       .or.(xr.gt.xr1+eps).or.(yr.gt.yr1+eps)) then
    93     write(*,*) 'outlon0,outlat0:'
    94121    write(*,*) outlon0,outlat0
    95     write(*,*) 'xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout:'
    96122    write(*,*) xr1,yr1,xlon0,ylat0,xr,yr,dxout,dyout
    97123    write(*,*) ' #### FLEXPART MODEL ERROR! PART OF OUTPUT    ####'
     
    104130  ! 2. Count Vertical levels of output grid
    105131  !****************************************
    106   j=0
    107 100   j=j+1
     132
     133  if (readerror.ne.0) then
     134    j=0
     135100 j=j+1
    108136    do i=1,3
    109137      read(unitoutgrid,*,end=99)
     
    112140    if (outhelp.eq.0.) goto 99
    113141    goto 100
    114 99   numzgrid=j-1
    115 
    116     allocate(outheight(numzgrid) &
    117          ,stat=stat)
    118     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    119     allocate(outheighthalf(numzgrid) &
    120          ,stat=stat)
    121     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    122 
    123 
    124   rewind(unitoutgrid)
    125   call skplin(29,unitoutgrid)
     14299  numzgrid=j-1
     143  else
     144    do i=1,maxoutlev
     145      if (outheights(i).lt.0) exit
     146    end do
     147    numzgrid=i-1
     148  end if
     149
     150  allocate(outheight(numzgrid),stat=stat)
     151  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
     152  allocate(outheighthalf(numzgrid),stat=stat)
     153  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
    126154
    127155  ! 2. Vertical levels of output grid
    128156  !**********************************
    129157
    130   j=0
    131 1000   j=j+1
    132     do i=1,3
    133       read(unitoutgrid,*,end=990)
    134     end do
    135     read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
    136     if (outhelp.eq.0.) goto 99
    137     outheight(j)=outhelp
    138     goto 1000
    139 990   numzgrid=j-1
    140 
     158  if (readerror.ne.0) then
     159
     160    rewind(unitoutgrid)
     161    call skplin(29,unitoutgrid)
     162
     163    do j=1,numzgrid
     164      do i=1,3
     165        read(unitoutgrid,*)
     166      end do
     167      read(unitoutgrid,'(4x,f7.1)') outhelp
     168      outheight(j)=outhelp
     169      outheights(j)=outhelp
     170    end do
     171    close(unitoutgrid)
     172
     173  else
     174
     175    do j=1,numzgrid
     176      outheight(j)=outheights(j)
     177    end do
     178
     179  endif
     180
     181  ! write outgrid file in namelist format to output directory if requested
     182  if (nmlout.eqv..true.) then
     183    ! reallocate outheights with actually required dimension for namelist writing
     184    deallocate(outheights)
     185    allocate(outheights(numzgrid),stat=stat)
     186    if (stat.ne.0) write(*,*)'ERROR: could not allocate outheights'
     187
     188    do j=1,numzgrid
     189      outheights(j)=outheight(j)
     190    end do
     191
     192    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID.namelist',err=1000)
     193    write(unitoutgrid,nml=outgrid)
     194    close(unitoutgrid)
     195  endif
    141196
    142197  ! Check whether vertical levels are specified in ascending order
     
    160215  end do
    161216
    162 
    163217  xoutshift=xlon0-outlon0
    164218  youtshift=ylat0-outlat0
    165   close(unitoutgrid)
    166 
    167     allocate(oroout(0:numxgrid-1,0:numygrid-1) &
    168          ,stat=stat)
    169     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    170     allocate(area(0:numxgrid-1,0:numygrid-1) &
    171          ,stat=stat)
    172     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    173     allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
    174          ,stat=stat)
    175     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    176     allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
    177          ,stat=stat)
    178     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    179     allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
    180          ,stat=stat)
    181     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
     219
     220  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
     221  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
     222  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
     223  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
     224  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     225  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
     226  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     227  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
     228  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
     229  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
    182230  return
    183231
    184232
    185 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     233999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    186234  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    187   write(*,*) ' #### xxx/flexpart/options                    #### '
     235  write(*,'(a)') path(1)(1:length(1))
    188236  stop
    189237
     2381000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     239  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     240  write(*,'(a)') path(2)(1:length(2))
     241  stop
     242
    190243end subroutine readoutgrid
  • trunk/src/readoutgrid_nest.f90

    r4 r27  
    5353  real,parameter :: eps=1.e-4
    5454
     55  integer :: readerror
    5556
     57  ! declare namelist
     58  namelist /outgridn/ &
     59    outlon0n,outlat0n, &
     60    numxgridn,numygridn, &
     61    dxoutn,dyoutn
     62
     63  ! helps identifying failed namelist input
     64  dxoutn=-1.0
    5665
    5766  ! Open the OUTGRID file and read output grid specifications
    5867  !**********************************************************
    5968
    60   open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST', &
    61        status='old',err=999)
     69  open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',form='formatted',status='old',err=999)
    6270
     71  ! try namelist input
     72  read(unitoutgrid,outgridn,iostat=readerror)
     73  close(unitoutgrid)
    6374
    64   call skplin(5,unitoutgrid)
     75  if ((dxoutn.le.0).or.(readerror.ne.0)) then
    6576
     77    open(unitoutgrid,file=path(1)(1:length(1))//'OUTGRID_NEST',status='old',err=999)
     78    call skplin(5,unitoutgrid)
    6679
    67   ! 1.  Read horizontal grid specifications
    68   !****************************************
     80    ! 1.  Read horizontal grid specifications
     81    !****************************************
    6982
    70   call skplin(3,unitoutgrid)
    71   read(unitoutgrid,'(4x,f11.4)') outlon0n
    72   call skplin(3,unitoutgrid)
    73   read(unitoutgrid,'(4x,f11.4)') outlat0n
    74   call skplin(3,unitoutgrid)
    75   read(unitoutgrid,'(4x,i5)') numxgridn
    76   call skplin(3,unitoutgrid)
    77   read(unitoutgrid,'(4x,i5)') numygridn
    78   call skplin(3,unitoutgrid)
    79   read(unitoutgrid,'(4x,f12.5)') dxoutn
    80   call skplin(3,unitoutgrid)
    81   read(unitoutgrid,'(4x,f12.5)') dyoutn
     83    call skplin(3,unitoutgrid)
     84    read(unitoutgrid,'(4x,f11.4)') outlon0n
     85    call skplin(3,unitoutgrid)
     86    read(unitoutgrid,'(4x,f11.4)') outlat0n
     87    call skplin(3,unitoutgrid)
     88    read(unitoutgrid,'(4x,i5)') numxgridn
     89    call skplin(3,unitoutgrid)
     90    read(unitoutgrid,'(4x,i5)') numygridn
     91    call skplin(3,unitoutgrid)
     92    read(unitoutgrid,'(4x,f12.5)') dxoutn
     93    call skplin(3,unitoutgrid)
     94    read(unitoutgrid,'(4x,f12.5)') dyoutn
    8295
     96    close(unitoutgrid)
     97  endif
    8398
    84     allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
    85          ,stat=stat)
    86     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    87     allocate(arean(0:numxgridn-1,0:numygridn-1) &
    88          ,stat=stat)
    89     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
    90     allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
    91          ,stat=stat)
    92     if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
     99  ! write outgrid_nest file in namelist format to output directory if requested
     100  if (nmlout.eqv..true.) then
     101    open(unitoutgrid,file=path(2)(1:length(2))//'OUTGRID_NEST.namelist',err=1000)
     102    write(unitoutgrid,nml=outgridn)
     103    close(unitoutgrid)
     104  endif
     105
     106  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
     107  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
     108  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
     109  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
     110  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
     111  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
    93112
    94113  ! Check validity of output grid (shall be within model domain)
     
    110129  xoutshiftn=xlon0-outlon0n
    111130  youtshiftn=ylat0-outlat0n
    112   close(unitoutgrid)
    113131  return
    114132
     133999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
     134  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     135  write(*,'(a)') path(1)(1:length(1))
     136  stop
    115137
    116 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
     1381000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
    117139  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    118   write(*,*) ' #### xxx/flexpart/options                    #### '
     140  write(*,'(a)') path(2)(1:length(2))
    119141  stop
    120142
  • trunk/src/readpaths.f90

    r20 r27  
    8888    length(numpath+2*(i-1)+2)=index(path(numpath+2*(i-1)+2),' ')-1
    8989  end do
    90   print*,length(5),length(6)
    9190
    9291
  • trunk/src/readreceptors.f90

    r4 r27  
    2727  !                                                                            *
    2828  !     Author: A. Stohl                                                       *
    29   !                                                                            *
    3029  !     1 August 1996                                                          *
     30  !     HSO, 14 August 2013
     31  !     Added optional namelist input
    3132  !                                                                            *
    3233  !*****************************************************************************
     
    5152  character(len=16) :: receptor
    5253
     54  integer :: readerror
     55  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
     56  integer,parameter :: unitreceptorout=2
     57
     58  ! declare namelist
     59  namelist /receptors/ &
     60    receptor, lon, lat
     61
     62  lon=-999.9
    5363
    5464  ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
     
    6474  !************************************************************
    6575
    66   open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS', &
    67        status='old',err=999)
     76  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
    6877
    69   call skplin(5,unitreceptor)
     78  ! try namelist input
     79  read(unitreceptor,receptors,iostat=readerror)
    7080
     81  ! prepare namelist output if requested
     82  if (nmlout.eqv..true.) then
     83    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='new',err=1000)
     84  endif
    7185
    72   ! Read the names and coordinates of the receptors
    73   !************************************************
     86  if ((lon.lt.-900).or.(readerror.ne.0)) then
    7487
    75   j=0
    76 100   j=j+1
     88    close(unitreceptor)
     89    open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
     90    call skplin(5,unitreceptor)
     91
     92    ! Read the names and coordinates of the receptors
     93    !************************************************
     94
     95    j=0
     96100 j=j+1
    7797    read(unitreceptor,*,end=99)
    7898    read(unitreceptor,*,end=99)
     
    89109    endif
    90110    if (j.gt.maxreceptor) then
    91     write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
    92     write(*,*) ' #### POINTS ARE GIVEN.                       #### '
    93     write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
    94     write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     111      write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
     112      write(*,*) ' #### POINTS ARE GIVEN.                       #### '
     113      write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
     114      write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
    95115    endif
    96116    receptorname(j)=receptor
     
    100120    ym=r_earth*dy/180.*pi
    101121    receptorarea(j)=xm*ym
     122
     123    ! write receptors file in namelist format to output directory if requested
     124    if (nmlout.eqv..true.) then
     125      lon=x
     126      lat=y
     127      write(unitreceptorout,nml=receptors)
     128    endif
     129
    102130    goto 100
    103131
    104 99   numreceptor=j-1
    105   close(unitreceptor)
     13299  numreceptor=j-1
     133
     134  else ! continue with namelist input
     135
     136    j=0
     137    do while (readerror.eq.0)
     138      j=j+1
     139      lon=-999.9
     140      read(unitreceptor,receptors,iostat=readerror)
     141      if ((lon.lt.-900).or.(readerror.ne.0)) then
     142        readerror=1
     143      else
     144        if (j.gt.maxreceptor) then
     145          write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
     146          write(*,*) ' #### POINTS ARE GIVEN.                       #### '
     147          write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
     148          write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     149        endif
     150        receptorname(j)=receptor
     151        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
     152        yreceptor(j)=(lat-ylat0)/dy
     153        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
     154        ym=r_earth*dy/180.*pi
     155        receptorarea(j)=xm*ym
     156      endif
     157
     158      ! write receptors file in namelist format to output directory if requested
     159      if (nmlout.eqv..true.) then
     160        write(unitreceptorout,nml=receptors)
     161      endif
     162
     163    end do
     164    numreceptor=j-1
     165    close(unitreceptor)
     166
     167  endif
     168
     169  if (nmlout.eqv..true.) then
     170    close(unitreceptorout)
     171  endif
     172
    106173  return
    107174
    108175
    109 999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
     176999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
    110177  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    111178  write(*,'(a)') path(1)(1:length(1))
    112179  stop
    113180
     1811000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
     182  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     183  write(*,'(a)') path(2)(1:length(2))
     184  stop
     185
    114186end subroutine readreceptors
  • trunk/src/readreleases.f90

    r24 r27  
    3333  !     Update: 29 January 2001                                                *
    3434  !     Release altitude can be either in magl or masl                         *
     35  !     HSO, 12 August 2013
     36  !     Added optional namelist input
    3537  !                                                                            *
    3638  !*****************************************************************************
     
    7173  implicit none
    7274
    73   integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
     75  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
    7476  integer,parameter :: num_min_discrete=100
    7577  real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
     
    7880  logical :: old
    7981
     82  ! help variables for namelist reading
     83  integer :: numpoints, parts, readerror
     84  integer*2 :: zkind
     85  integer :: idate1, itime1, idate2, itime2
     86  real :: lon1,lon2,lat1,lat2,z1,z2
     87  character*40 :: comment
     88  integer,parameter :: unitreleasesout=2
     89  real,allocatable, dimension (:) :: mass
     90  integer,allocatable, dimension (:) :: specnum_rel,specnum_rel2
     91
     92  ! declare namelists
     93  namelist /releases_ctrl/ &
     94    nspec, &
     95    specnum_rel
     96
     97  namelist /release/ &
     98    idate1, itime1, &
     99    idate2, itime2, &
     100    lon1, lon2, &
     101    lat1, lat2, &
     102    z1, z2, &
     103    zkind, &
     104    mass, &
     105    parts, &
     106    comment
     107
     108  numpoint=0
     109
     110  ! allocate with maxspec for first input loop
     111  allocate(mass(maxspec),stat=stat)
     112  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
     113  allocate(specnum_rel(maxspec),stat=stat)
     114  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
     115
     116  ! presetting namelist releases_ctrl
     117  nspec = -1  ! use negative value to determine failed namelist input
     118  specnum_rel = 0
     119
    80120  !sec, read release to find how many releasepoints should be allocated
    81 
    82   open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
    83        err=999)
    84 
    85   ! Check the format of the RELEASES file (either in free format,
    86   ! or using a formatted mask)
    87   ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
    88   !**************************************************************************
    89 
    90   call skplin(12,unitreleases)
    91   read (unitreleases,901) line
    92 901   format (a)
    93   if (index(line,'Total') .eq. 0) then
    94     old = .false.
    95   else
    96     old = .true.
    97   endif
     121  open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',form='formatted',err=999)
     122
     123  ! check if namelist input provided
     124  read(unitreleases,releases_ctrl,iostat=readerror)
     125
     126  ! prepare namelist output if requested
     127  if (nmlout.eqv..true.) then
     128    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
     129  endif
     130
     131  if ((readerror.ne.0).or.(nspec.lt.0)) then
     132
     133    ! no namelist format, close file and allow reopening in old format
     134    close(unitreleases)
     135    open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old',err=999)
     136
     137    readerror=1 ! indicates old format
     138
     139    ! Check the format of the RELEASES file (either in free format,
     140    ! or using a formatted mask)
     141    ! Use of formatted mask is assumed if line 10 contains the word 'DIRECTION'
     142    !**************************************************************************
     143
     144    call skplin(12,unitreleases)
     145    read (unitreleases,901) line
     146901 format (a)
     147    if (index(line,'Total') .eq. 0) then
     148      old = .false.
     149    else
     150      old = .true.
     151    endif
     152    rewind(unitreleases)
     153
     154
     155    ! Skip first 11 lines (file header)
     156    !**********************************
     157
     158    call skplin(11,unitreleases)
     159
     160    read(unitreleases,*,err=998) nspec
     161    if (old) call skplin(2,unitreleases)
     162    do i=1,nspec
     163      read(unitreleases,*,err=998) specnum_rel(i)
     164      if (old) call skplin(2,unitreleases)
     165    end do
     166
     167    numpoint=0
     168100 numpoint=numpoint+1
     169    read(unitreleases,*,end=25)
     170    read(unitreleases,*,err=998,end=25) idum,idum
     171    if (old) call skplin(2,unitreleases)
     172    read(unitreleases,*,err=998) idum,idum
     173    if (old) call skplin(2,unitreleases)
     174    read(unitreleases,*,err=998) xdum
     175    if (old) call skplin(2,unitreleases)
     176    read(unitreleases,*,err=998) xdum
     177    if (old) call skplin(2,unitreleases)
     178    read(unitreleases,*,err=998) xdum
     179    if (old) call skplin(2,unitreleases)
     180    read(unitreleases,*,err=998) xdum
     181    if (old) call skplin(2,unitreleases)
     182    read(unitreleases,*,err=998) idum
     183    if (old) call skplin(2,unitreleases)
     184    read(unitreleases,*,err=998) xdum
     185    if (old) call skplin(2,unitreleases)
     186    read(unitreleases,*,err=998) xdum
     187    if (old) call skplin(2,unitreleases)
     188    read(unitreleases,*,err=998) idum
     189    if (old) call skplin(2,unitreleases)
     190    do i=1,nspec
     191      read(unitreleases,*,err=998) xdum
     192      if (old) call skplin(2,unitreleases)
     193    end do
     194    !save compoint only for the first 1000 release points
     195    read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
     196    if (old) call skplin(1,unitreleases)
     197
     198    goto 100
     199
     20025  numpoint=numpoint-1
     201
     202  else
     203
     204    readerror=0
     205    do while (readerror.eq.0)
     206      idate1=-1
     207      read(unitreleases,release,iostat=readerror)
     208      if ((idate1.lt.0).or.(readerror.ne.0)) then
     209        readerror=1
     210      else
     211        numpoint=numpoint+1
     212      endif
     213    end do
     214    readerror=0
     215  endif ! if namelist input
     216
    98217  rewind(unitreleases)
    99218
    100 
    101   ! Skip first 11 lines (file header)
    102   !**********************************
    103 
    104   call skplin(11,unitreleases)
    105 
    106 
    107   read(unitreleases,*,err=998) nspec
    108   if (old) call skplin(2,unitreleases)
    109   do i=1,nspec
    110     read(unitreleases,*,err=998) specnum_rel
    111     if (old) call skplin(2,unitreleases)
     219  ! allocate arrays of matching size for number of species (namelist output)
     220  deallocate(mass)
     221  allocate(mass(nspec),stat=stat)
     222  if (stat.ne.0) write(*,*)'ERROR: could not allocate mass'
     223  allocate(specnum_rel2(nspec),stat=stat)
     224  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel2'
     225  specnum_rel2=specnum_rel(1:nspec)
     226  deallocate(specnum_rel)
     227  allocate(specnum_rel(nspec),stat=stat)
     228  if (stat.ne.0) write(*,*)'ERROR: could not allocate specnum_rel'
     229  specnum_rel=specnum_rel2
     230  deallocate(specnum_rel2)
     231
     232  !allocate memory for numpoint releaspoints
     233  allocate(ireleasestart(numpoint),stat=stat)
     234  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleasestart'
     235  allocate(ireleaseend(numpoint),stat=stat)
     236  if (stat.ne.0) write(*,*)'ERROR: could not allocate ireleaseend'
     237  allocate(xpoint1(numpoint),stat=stat)
     238  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint1'
     239  allocate(xpoint2(numpoint),stat=stat)
     240  if (stat.ne.0) write(*,*)'ERROR: could not allocate xpoint2'
     241  allocate(ypoint1(numpoint),stat=stat)
     242  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint1'
     243  allocate(ypoint2(numpoint),stat=stat)
     244  if (stat.ne.0) write(*,*)'ERROR: could not allocate ypoint2'
     245  allocate(zpoint1(numpoint),stat=stat)
     246  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint1'
     247  allocate(zpoint2(numpoint),stat=stat)
     248  if (stat.ne.0) write(*,*)'ERROR: could not allocate zpoint2'
     249  allocate(kindz(numpoint),stat=stat)
     250  if (stat.ne.0) write(*,*)'ERROR: could not allocate kindz'
     251  allocate(xmass(numpoint,maxspec),stat=stat)
     252  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmass'
     253  allocate(rho_rel(numpoint),stat=stat)
     254  if (stat.ne.0) write(*,*)'ERROR: could not allocate rho_rel'
     255  allocate(npart(numpoint),stat=stat)
     256  if (stat.ne.0) write(*,*)'ERROR: could not allocate npart'
     257  allocate(xmasssave(numpoint),stat=stat)
     258  if (stat.ne.0) write(*,*)'ERROR: could not allocate xmasssave'
     259
     260  write (*,*) 'Releasepoints : ', numpoint
     261
     262  do i=1,numpoint
     263    xmasssave(i)=0.
    112264  end do
    113 
    114   numpoint=0
    115 100   numpoint=numpoint+1
    116   read(unitreleases,*,end=25)
    117   read(unitreleases,*,err=998,end=25) idum,idum
    118   if (old) call skplin(2,unitreleases)
    119   read(unitreleases,*,err=998) idum,idum
    120   if (old) call skplin(2,unitreleases)
    121   read(unitreleases,*,err=998) xdum
    122   if (old) call skplin(2,unitreleases)
    123   read(unitreleases,*,err=998) xdum
    124   if (old) call skplin(2,unitreleases)
    125   read(unitreleases,*,err=998) xdum
    126   if (old) call skplin(2,unitreleases)
    127   read(unitreleases,*,err=998) xdum
    128   if (old) call skplin(2,unitreleases)
    129   read(unitreleases,*,err=998) idum
    130   if (old) call skplin(2,unitreleases)
    131   read(unitreleases,*,err=998) xdum
    132   if (old) call skplin(2,unitreleases)
    133   read(unitreleases,*,err=998) xdum
    134   if (old) call skplin(2,unitreleases)
    135   read(unitreleases,*,err=998) idum
    136   if (old) call skplin(2,unitreleases)
    137   do i=1,nspec
    138     read(unitreleases,*,err=998) xdum
    139     if (old) call skplin(2,unitreleases)
    140   end do
    141   !save compoint only for the first 1000 release points
    142   read(unitreleases,'(a40)',err=998) compoint(1)(1:40)
    143   if (old) call skplin(1,unitreleases)
    144 
    145   goto 100
    146 
    147 25   numpoint=numpoint-1
    148 
    149   !allocate memory for numpoint releaspoint
    150     allocate(ireleasestart(numpoint) &
    151          ,stat=stat)
    152     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    153     allocate(ireleaseend(numpoint) &
    154          ,stat=stat)
    155     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    156     allocate(xpoint1(numpoint) &
    157          ,stat=stat)
    158     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    159     allocate(xpoint2(numpoint) &
    160          ,stat=stat)
    161     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    162     allocate(ypoint1(numpoint) &
    163          ,stat=stat)
    164     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    165     allocate(ypoint2(numpoint) &
    166          ,stat=stat)
    167     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    168     allocate(zpoint1(numpoint) &
    169          ,stat=stat)
    170     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    171     allocate(zpoint2(numpoint) &
    172          ,stat=stat)
    173     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    174     allocate(kindz(numpoint) &
    175          ,stat=stat)
    176     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    177     allocate(xmass(numpoint,maxspec) &
    178          ,stat=stat)
    179     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    180     allocate(rho_rel(numpoint) &
    181          ,stat=stat)
    182     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    183     allocate(npart(numpoint) &
    184          ,stat=stat)
    185     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    186     allocate(xmasssave(numpoint) &
    187          ,stat=stat)
    188     if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
    189 
    190    write (*,*) ' Releasepoints allocated: ', numpoint
    191 
    192    do i=1,numpoint
    193      xmasssave(i)=0.
    194    end do
    195265
    196266  !now save the information
     
    203273  end do
    204274
    205   rewind(unitreleases)
    206 
    207 
    208   ! Skip first 11 lines (file header)
    209   !**********************************
    210 
    211   call skplin(11,unitreleases)
    212 
    213 
    214   ! Assign species-specific parameters needed for physical processes
    215   !*************************************************************************
    216 
    217   read(unitreleases,*,err=998) nspec
    218   if (nspec.gt.maxspec) goto 994
    219   if (old) call skplin(2,unitreleases)
     275  if (readerror.ne.0) then
     276    ! Skip first 11 lines (file header)
     277    !**********************************
     278
     279    call skplin(11,unitreleases)
     280
     281    ! Assign species-specific parameters needed for physical processes
     282    !*************************************************************************
     283
     284    read(unitreleases,*,err=998) nspec
     285    if (nspec.gt.maxspec) goto 994
     286    if (old) call skplin(2,unitreleases)
     287  endif
     288
     289  ! namelist output
     290  if (nmlout.eqv..true.) then
     291    write(unitreleasesout,nml=releases_ctrl)
     292  endif
     293
    220294  do i=1,nspec
    221     read(unitreleases,*,err=998) specnum_rel
    222     if (old) call skplin(2,unitreleases)
    223     call readspecies(specnum_rel,i)
    224 
    225   ! For backward runs, only 1 species is allowed
    226   !*********************************************
    227 
    228   !if ((ldirect.lt.0).and.(nspec.gt.1)) then
    229   !write(*,*) '#####################################################'
    230   !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
    231   !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
    232   !write(*,*) '#####################################################'
    233   !  stop
    234   !endif
    235 
    236   ! Molecular weight
    237   !*****************
    238 
    239     if (((iout.eq.2).or.(iout.eq.3)).and. &
    240          (weightmolar(i).lt.0.)) then
     295    if (readerror.ne.0) then
     296      read(unitreleases,*,err=998) specnum_rel(i)
     297      if (old) call skplin(2,unitreleases)
     298      call readspecies(specnum_rel(i),i)
     299    else
     300      call readspecies(specnum_rel(i),i)
     301    endif
     302
     303    ! For backward runs, only 1 species is allowed
     304    !*********************************************
     305
     306    !if ((ldirect.lt.0).and.(nspec.gt.1)) then
     307    !write(*,*) '#####################################################'
     308    !write(*,*) '#### FLEXPART MODEL SUBROUTINE READRELEASES:     ####'
     309    !write(*,*) '#### FOR BACKWARD RUNS, ONLY 1 SPECIES IS ALLOWED####'
     310    !write(*,*) '#####################################################'
     311    !  stop
     312    !endif
     313
     314    ! Molecular weight
     315    !*****************
     316
     317    if (((iout.eq.2).or.(iout.eq.3)).and.(weightmolar(i).lt.0.)) then
    241318      write(*,*) 'For mixing ratio output, valid molar weight'
    242319      write(*,*) 'must be specified for all simulated species.'
     
    246323    endif
    247324
    248 
    249   ! Radioactive decay
    250   !******************
     325    ! Radioactive decay
     326    !******************
    251327
    252328    decay(i)=0.693147/decay(i) !conversion half life to decay constant
    253329
    254330
    255   ! Dry deposition of gases
    256   !************************
    257 
    258     if (reldiff(i).gt.0.) &
    259          rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
    260 
    261   ! Dry deposition of particles
    262   !****************************
     331    ! Dry deposition of gases
     332    !************************
     333
     334    if (reldiff(i).gt.0.) rm(i)=1./(henry(i)/3000.+100.*f0(i))    ! mesophyll resistance
     335
     336    ! Dry deposition of particles
     337    !****************************
    263338
    264339    vsetaver(i)=0.
    265340    cunningham(i)=0.
    266341    dquer(i)=dquer(i)*1000000.         ! Conversion m to um
    267     if (density(i).gt.0.) then                  ! Additional parameters
    268      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
     342    if (density(i).gt.0.) then         ! Additional parameters
     343      call part0(dquer(i),dsigma(i),density(i),fracth,schmih,cun,vsh)
    269344      do j=1,ni
    270345        fract(i,j)=fracth(j)
     
    277352    endif
    278353
    279   ! Dry deposition for constant deposition velocity
    280   !************************************************
     354    ! Dry deposition for constant deposition velocity
     355    !************************************************
    281356
    282357    dryvel(i)=dryvel(i)*0.01         ! conversion to m/s
    283358
    284   ! Check if wet deposition or OH reaction shall be calculated
    285   !***********************************************************
     359    ! Check if wet deposition or OH reaction shall be calculated
     360    !***********************************************************
    286361    if (weta(i).gt.0.)  then
    287362      WETDEP=.true.
    288       write (*,*) 'Below-cloud scavenging is switched on'
     363      write (*,*) 'Below-cloud scavenging: ON'
    289364      write (*,*) 'Below-cloud scavenging coefficients: ',weta(i),i
    290365    else
    291       write (*,*) 'Below-cloud scavenging is switched OFF'
     366      write (*,*) 'Below-cloud scavenging: OFF'
    292367    endif
    293368   
    294 ! NIK 31.01.2013 + 10.12.2013
     369    ! NIK 31.01.2013 + 10.12.2013
    295370    if (weta_in(i).gt.0.)  then
    296371      WETDEP=.true.
    297       write (*,*) 'In-cloud scavenging is switched on'
     372      write (*,*) 'In-cloud scavenging: ON'
    298373      write (*,*) 'In-cloud scavenging coefficients: ',weta_in(i),wetb_in(i), wetc_in(i), wetd_in(i),i
    299374    else
    300       write (*,*) 'In-cloud scavenging is switched OFF'
     375      write (*,*) 'In-cloud scavenging: OFF'
    301376    endif
    302377
    303378    if (ohreact(i).gt.0) then
    304379      OHREA=.true.
    305       write (*,*) 'OHreaction switched on: ',ohreact(i),i
    306     endif
    307 
    308 
    309     if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or. &
    310          (dryvel(i).gt.0.)) then
     380      write (*,*) 'OHreaction: ON (',ohreact(i),i,')'
     381    endif
     382
     383    if ((reldiff(i).gt.0.).or.(density(i).gt.0.).or.(dryvel(i).gt.0.)) then
    311384      DRYDEP=.true.
    312385      DRYDEPSPEC(i)=.true.
     
    315388  end do
    316389
    317     if (WETDEP.or.DRYDEP) DEP=.true.
     390  if (WETDEP.or.DRYDEP) DEP=.true.
    318391
    319392  ! Read specifications for each release point
    320393  !*******************************************
    321 
     394  numpoints=numpoint
    322395  numpoint=0
    323396  numpartmax=0
    324397  releaserate=0.
    325 1000   numpoint=numpoint+1
    326   read(unitreleases,*,end=250)
    327   read(unitreleases,*,err=998,end=250) id1,it1
    328   if (old) call skplin(2,unitreleases)
    329   read(unitreleases,*,err=998) id2,it2
    330   if (old) call skplin(2,unitreleases)
    331   read(unitreleases,*,err=998) xpoint1(numpoint)
    332   if (old) call skplin(2,unitreleases)
    333   read(unitreleases,*,err=998) ypoint1(numpoint)
    334   if (old) call skplin(2,unitreleases)
    335   read(unitreleases,*,err=998) xpoint2(numpoint)
    336   if (old) call skplin(2,unitreleases)
    337   read(unitreleases,*,err=998) ypoint2(numpoint)
    338   if (old) call skplin(2,unitreleases)
    339   read(unitreleases,*,err=998) kindz(numpoint)
    340   if (old) call skplin(2,unitreleases)
    341   read(unitreleases,*,err=998) zpoint1(numpoint)
    342   if (old) call skplin(2,unitreleases)
    343   read(unitreleases,*,err=998) zpoint2(numpoint)
    344   if (old) call skplin(2,unitreleases)
    345   read(unitreleases,*,err=998) npart(numpoint)
    346   if (old) call skplin(2,unitreleases)
    347   do i=1,nspec
    348     read(unitreleases,*,err=998) xmass(numpoint,i)
    349     if (old) call skplin(2,unitreleases)
    350   end do
    351   !save compoint only for the first 1000 release points
    352   if (numpoint.le.1000) then
    353     read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
     398101   numpoint=numpoint+1
     399
     400  if (readerror.lt.1) then ! reading namelist format
     401
     402    if (numpoint.gt.numpoints) goto 250
     403    zkind = 1
     404    mass = 0
     405    parts = 0
     406    comment = ' '
     407    read(unitreleases,release,iostat=readerror)
     408    id1=idate1
     409    it1=itime1
     410    id2=idate2
     411    it2=itime2
     412    xpoint1(numpoint)=lon1
     413    xpoint2(numpoint)=lon2
     414    ypoint1(numpoint)=lat1
     415    ypoint2(numpoint)=lat2
     416    zpoint1(numpoint)=z1
     417    zpoint2(numpoint)=z2
     418    kindz(numpoint)=zkind
     419    do i=1,nspec
     420      xmass(numpoint,i)=mass(i)
     421    end do
     422    npart(numpoint)=parts
     423    compoint(min(1001,numpoint))=comment
     424
     425    ! namelist output
     426    if (nmlout.eqv..true.) then
     427      write(unitreleasesout,nml=release)
     428    endif
     429
    354430  else
    355     read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
    356   endif
    357   if (old) call skplin(1,unitreleases)
    358   if (numpoint.le.1000) then
    359     if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
    360          (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
    361          (compoint(numpoint)(1:8).eq.'        ')) goto 250
    362   else
    363     if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
    364          (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
    365   endif
    366 
     431
     432    read(unitreleases,*,end=250)
     433    read(unitreleases,*,err=998,end=250) id1,it1
     434    if (old) call skplin(2,unitreleases)
     435    read(unitreleases,*,err=998) id2,it2
     436    if (old) call skplin(2,unitreleases)
     437    read(unitreleases,*,err=998) xpoint1(numpoint)
     438    if (old) call skplin(2,unitreleases)
     439    read(unitreleases,*,err=998) ypoint1(numpoint)
     440    if (old) call skplin(2,unitreleases)
     441    read(unitreleases,*,err=998) xpoint2(numpoint)
     442    if (old) call skplin(2,unitreleases)
     443    read(unitreleases,*,err=998) ypoint2(numpoint)
     444    if (old) call skplin(2,unitreleases)
     445    read(unitreleases,*,err=998) kindz(numpoint)
     446    if (old) call skplin(2,unitreleases)
     447    read(unitreleases,*,err=998) zpoint1(numpoint)
     448    if (old) call skplin(2,unitreleases)
     449    read(unitreleases,*,err=998) zpoint2(numpoint)
     450    if (old) call skplin(2,unitreleases)
     451    read(unitreleases,*,err=998) npart(numpoint)
     452    if (old) call skplin(2,unitreleases)
     453    do i=1,nspec
     454      read(unitreleases,*,err=998) xmass(numpoint,i)
     455      if (old) call skplin(2,unitreleases)
     456      mass(i)=xmass(numpoint,i)
     457    end do
     458    !save compoint only for the first 1000 release points
     459    if (numpoint.le.1000) then
     460      read(unitreleases,'(a40)',err=998) compoint(numpoint)(1:40)
     461      comment=compoint(numpoint)(1:40)
     462    else
     463      read(unitreleases,'(a40)',err=998) compoint(1001)(1:40)
     464      comment=compoint(1001)(1:40)
     465    endif
     466    if (old) call skplin(1,unitreleases)
     467
     468    ! namelist output
     469    if (nmlout.eqv..true.) then
     470      idate1=id1
     471      itime1=it1
     472      idate2=id2
     473      itime2=it2
     474      lon1=xpoint1(numpoint)
     475      lon2=xpoint2(numpoint)
     476      lat1=ypoint1(numpoint)
     477      lat2=ypoint2(numpoint)
     478      z1=zpoint1(numpoint)
     479      z2=zpoint2(numpoint)
     480      zkind=kindz(numpoint)
     481      parts=npart(numpoint)
     482      write(unitreleasesout,nml=release)
     483    endif
     484
     485    if (numpoint.le.1000) then
     486      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
     487           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.).and. &
     488           (compoint(numpoint)(1:8).eq.'        ')) goto 250
     489    else
     490      if((xpoint1(numpoint).eq.0.).and.(ypoint1(numpoint).eq.0.).and. &
     491           (xpoint2(numpoint).eq.0.).and.(ypoint2(numpoint).eq.0.)) goto 250
     492    endif
     493
     494  endif ! if namelist format
    367495
    368496  ! If a release point contains no particles, stop and issue error message
     
    435563  endif
    436564
    437 
    438   ! Check, whether the total number of particles may exceed totally allowed
    439   ! number of particles at some time during the simulation
    440   !************************************************************************
    441 
    442565  ! Determine the release rate (particles per second) and total number
    443566  ! of particles released during the simulation
     
    451574  endif
    452575  numpartmax=numpartmax+npart(numpoint)
    453   goto 1000
    454 
     576  goto 101
    455577
    456578250   close(unitreleases)
    457579
    458   write (*,*) ' Particles allocated for this run: ',maxpart, ', released in simulation: ',  numpartmax
     580  if (nmlout.eqv..true.) then
     581    close(unitreleasesout)
     582  endif
     583
     584  write (*,*) 'Particles allocated (maxpart)  : ',maxpart
     585  write (*,*) 'Particles released (numpartmax): ',numpartmax
    459586  numpoint=numpoint-1
    460587
     
    464591    maxpointspec_act=1
    465592  endif
     593
     594  ! Check, whether the total number of particles may exceed totally allowed
     595  ! number of particles at some time during the simulation
     596  !************************************************************************
    466597
    467598  if (releaserate.gt. &
     
    521652  stop
    522653
     6541000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RELEASES"    #### '
     655  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     656  write(*,'(a)') path(2)(1:length(2))
     657  stop
     658
    523659end subroutine readreleases
  • trunk/src/readspecies.f90

    r24 r27  
    3030  !                                                                            *
    3131  !   11 July 1996                                                             *
    32   !                                                                            *
     32  !
    3333  !   Changes:                                                                 *
    3434  !   N. Kristiansen, 31.01.2013: Including parameters for in-cloud scavenging *
     35  !                                                                            *
     36  !   HSO, 13 August 2013
     37  !   added optional namelist input
    3538  !                                                                            *
    3639  !*****************************************************************************
     
    6265  logical :: spec_found
    6366
     67  character(len=16) :: pspecies
     68  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
     69  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
     70  real :: pweta_in, pwetb_in, pwetc_in, pwetd_in
     71  integer :: readerror
     72
     73  ! declare namelist
     74  namelist /species_params/ &
     75   pspecies, pdecay, pweta, pwetb, &
     76   pweta_in, pwetb_in, pwetc_in, pwetd_in, &
     77   preldiff, phenry, pf0, pdensity, pdquer, &
     78   pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
     79
     80  pspecies=" "
     81  pdecay=-999.9
     82  pweta=-9.9E-09
     83  pwetb=0.0
     84  pweta_in=-9.9E-09
     85  pwetb_in=-9.9E-09
     86  pwetc_in=-9.9E-09
     87  pwetd_in=-9.9E-09
     88  preldiff=-9.9
     89  phenry=0.0
     90  pf0=0.0
     91  pdensity=-9.9E09
     92  pdquer=0.0
     93  pdsigma=0.0
     94  pdryvel=-9.99
     95  pohreact=-9.9E-09
     96  pspec_ass=-9
     97  pkao=-99.99
     98  pweightmolar=-789.0 ! read failure indicator value
     99
    64100  ! Open the SPECIES file and read species names and properties
    65101  !************************************************************
    66102  specnum(pos_spec)=id_spec
    67103  write(aspecnumb,'(i3.3)') specnum(pos_spec)
    68   open(unitspecies,file= &
    69        path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
    70        err=998)
     104  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',form='formatted',err=998)
    71105  !write(*,*) 'reading SPECIES',specnum(pos_spec)
    72106
    73107  ASSSPEC=.FALSE.
    74108
    75   do i=1,6
    76     read(unitspecies,*)
    77   end do
     109  ! try namelist input
     110  read(unitspecies,species_params,iostat=readerror)
     111  close(unitspecies)
     112   
     113  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
     114
     115    readerror=1
     116
     117    open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
     118
     119    do i=1,6
     120      read(unitspecies,*)
     121    end do
    78122
    79123    read(unitspecies,'(a10)',end=22) species(pos_spec)
     
    86130  !  write(*,*) wetb(pos_spec)
    87131
    88 !*** NIK 31.01.2013: including in-cloud scavening parameters
     132  !*** NIK 31.01.2013: including in-cloud scavening parameters
    89133   read(unitspecies,'(e18.1)',end=22) weta_in(pos_spec)
    90134  !  write(*,*) weta_in(pos_spec)
     
    118162    read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
    119163  !       write(*,*) kao(pos_spec)
    120     i=pos_spec
    121 
    122     if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
    123        if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
     164
     165    pspecies=species(pos_spec)
     166    pdecay=decay(pos_spec)
     167    pweta=weta(pos_spec)
     168    pwetb=wetb(pos_spec)
     169    pweta_in=weta_in(pos_spec)
     170    pwetb_in=wetb_in(pos_spec)
     171    pwetc_in=wetc_in(pos_spec)
     172    pwetd_in=wetd_in(pos_spec)
     173    preldiff=reldiff(pos_spec)
     174    phenry=henry(pos_spec)
     175    pf0=f0(pos_spec)
     176    pdensity=density(pos_spec)
     177    pdquer=dquer(pos_spec)
     178    pdsigma=dsigma(pos_spec)
     179    pdryvel=dryvel(pos_spec)
     180    pweightmolar=weightmolar(pos_spec)
     181    pohreact=ohreact(pos_spec)
     182    pspec_ass=spec_ass(pos_spec)
     183    pkao=kao(pos_spec)
     184
     185  else
     186
     187    species(pos_spec)=pspecies
     188    decay(pos_spec)=pdecay
     189    weta(pos_spec)=pweta
     190    wetb(pos_spec)=pwetb
     191    weta_in(pos_spec)=pweta_in
     192    wetb_in(pos_spec)=pwetb_in
     193    wetc_in(pos_spec)=pwetc_in
     194    wetd_in(pos_spec)=pwetd_in
     195    reldiff(pos_spec)=preldiff
     196    henry(pos_spec)=phenry
     197    f0(pos_spec)=pf0
     198    density(pos_spec)=pdensity
     199    dquer(pos_spec)=pdquer
     200    dsigma(pos_spec)=pdsigma
     201    dryvel(pos_spec)=pdryvel
     202    weightmolar(pos_spec)=pweightmolar
     203    ohreact(pos_spec)=pohreact
     204    spec_ass(pos_spec)=pspec_ass
     205    kao(pos_spec)=pkao
     206
     207  endif
     208
     209  i=pos_spec
     210
     211  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
     212   if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
     213  endif
     214
     215  if (spec_ass(pos_spec).gt.0) then
     216    spec_found=.FALSE.
     217    do j=1,pos_spec-1
     218      if (spec_ass(pos_spec).eq.specnum(j)) then
     219        spec_ass(pos_spec)=j
     220        spec_found=.TRUE.
     221        ASSSPEC=.TRUE.
     222      endif
     223    end do
     224    if (spec_found.eqv..FALSE.) then
     225      goto 997
    124226    endif
    125 
    126     if (spec_ass(pos_spec).gt.0) then
    127        spec_found=.FALSE.
    128        do j=1,pos_spec-1
    129           if (spec_ass(pos_spec).eq.specnum(j)) then
    130              spec_ass(pos_spec)=j
    131              spec_found=.TRUE.
    132              ASSSPEC=.TRUE.
    133           endif
    134        end do
    135        if (spec_found.eqv..FALSE.) then
    136           goto 997
    137        endif
    138     endif
    139 
    140     if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
    141     if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
    142 
    143     if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
    144       write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
    145       write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
    146       write(*,*) '#### PARTICLE AND GAS.                       ####'
    147       write(*,*) '#### SPECIES NUMBER',aspecnumb
    148       stop
    149     endif
     227  endif
     228
     229  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
     230  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
     231
     232  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
     233    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
     234    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
     235    write(*,*) '#### PARTICLE AND GAS.                       ####'
     236    write(*,*) '#### SPECIES NUMBER',aspecnumb
     237    stop
     238  endif
    15023920   continue
    151240
    152 
    153   ! Read in daily and day-of-week variation of emissions, if available
    154   !*******************************************************************
     241  if (readerror.ne.0) then ! text format input
     242
     243    ! Read in daily and day-of-week variation of emissions, if available
     244    !*******************************************************************
     245    ! HSO: This is not yet implemented as namelist parameters
    155246
    156247    do j=1,24           ! initialize everything to no variation
     
    172263    end do
    173264
    174 22   close(unitspecies)
    175 
    176    return
     265  endif
     266
     26722 close(unitspecies)
     268
     269  ! namelist output if requested
     270  if (nmlout.eqv..true.) then
     271    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='new',err=1000)
     272    write(unitspecies,nml=species_params)
     273    close(unitspecies)
     274  endif
     275
     276  return
    177277
    178278996   write(*,*) '#####################################################'
     
    203303  stop
    204304
     3051000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "SPECIES_',aspecnumb,'.namelist'
     306  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
     307  write(*,'(a)') path(2)(1:length(2))
     308  stop
    205309
    206310end subroutine readspecies
  • trunk/src/timemanager.f90

    r20 r27  
    129129
    130130
     131  !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
     132  write(*,46) float(itime)/3600,itime,numpart
    131133  if (verbosity.gt.0) then
    132134    write (*,*) 'timemanager> starting simulation'
     
    371373        if ((iout.eq.4).or.(iout.eq.5)) call plumetraj(itime)
    372374        if (iflux.eq.1) call fluxoutput(itime)
    373         write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc, &
    374              drygridtotalunc
    375 45      format(i9,' SECONDS SIMULATED: ',i8, &
    376              ' PARTICLES:    Uncertainty: ',3f7.3)
     375        !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
     376        write(*,46) float(itime)/3600,itime,numpart
     37745      format(i9,' SECONDS SIMULATED: ',i8, ' PARTICLES:    Uncertainty: ',3f7.3)
     37846      format(' Simulated ',f7.1,' hours (',i9,' s), ',i8, ' particles')
    377379        if (ipout.ge.1) call partoutput(itime)    ! dump particle positions
    378380        loutnext=loutnext+loutstep
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG