Changeset 3405994 in flexpart.git


Ignore:
Timestamp:
May 28, 2020, 6:10:55 PM (4 years ago)
Author:
pesei <petra.seibert at univie.ac.at>
Branches:
10.4.1_pesei, bugfixes+enhancements, release-10.4.1, scaling-bug
Children:
477a094
Parents:
5cbd51b
Message:

ticket:270 fix reading of receptors in namelist format

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/readreceptors.f90

    r92fab65 r3405994  
    1010  !     Author: A. Stohl                                                       *
    1111  !     1 August 1996                                                          *
    12   !     HSO, 14 August 2013
    13   !     Added optional namelist input
     12  !     HSO, 14 August 2013: Added optional namelist input                     *
     13  !     PS, 2020-05-28: correct bug in nml input, code cosmetics               *
    1414  !                                                                            *
    1515  !*****************************************************************************
    … …  
    3434  character(len=16) :: receptor
    3535
    36   integer :: readerror
     36  integer :: ierr
    3737  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
    38   integer,parameter :: unitreceptorout=2
     38  integer,parameter :: iunitreceptorout=2
    3939
    4040  ! declare namelist
    41   namelist /receptors/ &
    42     receptor, lon, lat
    43 
    44   lon=-999.9
     41  namelist /receptors/ receptor, lon, lat
    4542
    4643  ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
    … …  
    5249  endif
    5350
     51! prepare namelist output if requested
     52  if (nmlout .and. lroot) &
     53    open(iunitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist', &
     54      status='replace',err=1000)
    5455
    55   ! Open the RECEPTORS file and read output grid specifications
    56   !************************************************************
     56! Open the RECEPTORS file and read output grid specifications
     57!************************************************************
     58! try namelist input
     59  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
     60  read(unitreceptor,receptors,iostat=ierr)
     61  close(unitreceptor)
    5762
    58   open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',form='formatted',status='old',err=999)
     63  open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS')
    5964
    60   ! try namelist input
    61   read(unitreceptor,receptors,iostat=readerror)
    62 
    63   ! prepare namelist output if requested
    64   if (nmlout.and.lroot) then
    65     open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',&
    66          &access='append',status='replace',err=1000)
    67   endif
    68 
    69   if ((lon.lt.-900).or.(readerror.ne.0)) then
    70 
    71     close(unitreceptor)
    72     open(unitreceptor,file=path(1)(1:length(1))//'RECEPTORS',status='old',err=999)
     65  if (ierr.ne.0) then ! not namelist
     66 
    7367    call skplin(5,unitreceptor)
    7468
    … …  
    7872    j=0
    7973100 j=j+1
    80     read(unitreceptor,*,end=99)
    81     read(unitreceptor,*,end=99)
    82     read(unitreceptor,*,end=99)
    83     read(unitreceptor,'(4x,a16)',end=99) receptor
    84     call skplin(3,unitreceptor)
    85     read(unitreceptor,'(4x,f11.4)',end=99) x
    86     call skplin(3,unitreceptor)
    87     read(unitreceptor,'(4x,f11.4)',end=99) y
    88     if ((x.eq.0.).and.(y.eq.0.).and. &
    89          (receptor.eq.'                ')) then
    90       j=j-1
    91       goto 100
    92     endif
    93     if (j.gt.maxreceptor) then
    94       write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
    95       write(*,*) ' #### POINTS ARE GIVEN.                       #### '
    96       write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
    97       write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
    98     endif
    99     receptorname(j)=receptor
    100     xreceptor(j)=(x-xlon0)/dx       ! transform to grid coordinates
    101     yreceptor(j)=(y-ylat0)/dy
    102     xm=r_earth*cos(y*pi/180.)*dx/180.*pi
    103     ym=r_earth*dy/180.*pi
    104     receptorarea(j)=xm*ym
     74      read(unitreceptor,*,end=99)
     75      read(unitreceptor,*,end=99)
     76      read(unitreceptor,*,end=99)
     77      read(unitreceptor,'(4x,a16)',end=99) receptor
     78      call skplin(3,unitreceptor)
     79      read(unitreceptor,'(4x,f11.4)',end=99) x
     80      call skplin(3,unitreceptor)
     81      read(unitreceptor,'(4x,f11.4)',end=99) y
     82      if (x.eq.0. .and. y.eq.0. .and. receptor.eq.'                ') then
     83        j=j-1
     84        goto 100
     85      endif
     86      if (j.gt.maxreceptor) goto 998 ! ERROR - STOP
     87      receptorname(j)=receptor
     88      xreceptor(j)=(x-xlon0)/dx       ! transform to grid coordinates
     89      yreceptor(j)=(y-ylat0)/dy
     90      xm=r_earth*cos(y*pi/180.)*dx/180.*pi
     91      ym=r_earth*dy/180.*pi
     92      receptorarea(j)=xm*ym
    10593
    106     ! write receptors file in namelist format to output directory if requested
    107     if (nmlout.and.lroot) then
    108       lon=x
    109       lat=y
    110       write(unitreceptorout,nml=receptors)
    111     endif
     94      ! write receptors file in namelist format to output directory if requested
     95      if (nmlout .and. lroot) then
     96        lon=x
     97        lat=y
     98        write(iunitreceptorout,nml=receptors)
     99      endif
    112100
    113     goto 100
     101    goto 100 ! read next
    114102
    115 99  numreceptor=j-1
     10399  numreceptor=j-1 ! read all
    116104
    117105  else ! continue with namelist input
    118106
    119107    j=0
    120     do while (readerror.eq.0)
     108    do while (ierr.eq.0)
    121109      j=j+1
    122       lon=-999.9
    123       read(unitreceptor,receptors,iostat=readerror)
    124       if ((lon.lt.-900).or.(readerror.ne.0)) then
    125         readerror=1
    126       else
    127         if (j.gt.maxreceptor) then
    128           write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
    129           write(*,*) ' #### POINTS ARE GIVEN.                       #### '
    130           write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
    131           write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
    132         endif
     110      read(unitreceptor,receptors,iostat=ierr)
     111      if (ierr.eq.0) then
     112        if (j.gt.maxreceptor) goto 998 ! ERROR - STOP
    133113        receptorname(j)=receptor
    134114        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
    … …  
    137117        ym=r_earth*dy/180.*pi
    138118        receptorarea(j)=xm*ym
     119      ! write receptors file in namelist format to output directory if requested
     120        if (nmlout.and.lroot) &
     121        write(iunitreceptorout,nml=receptors)
    139122      endif
    140 
    141       ! write receptors file in namelist format to output directory if requested
    142       if (nmlout.and.lroot) then
    143         write(unitreceptorout,nml=receptors)
    144       endif
    145 
    146123    end do
    147124    numreceptor=j-1
    148125    close(unitreceptor)
    149126
    150   endif
     127  endif ! end reading nml input
    151128
    152   if (nmlout.and.lroot) then
    153     close(unitreceptorout)
    154   endif
     129  if (nmlout.and.lroot) &
     130    close(iunitreceptorout)
    155131
    156132  return
    157133
     134998 continue
     135  write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
     136  write(*,*) ' #### POINTS ARE GIVEN.                       #### '
     137  write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
     138  write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     139  stop
    158140
    159 999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
     141999 continue
     142  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
    160143  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    161144  write(*,'(a)') path(1)(1:length(1))
    162145  stop
    163146
    164 1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
     1471000 continue
     148  write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
    165149  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
    166150  write(*,'(a)') path(2)(1:length(2))
Note: See TracChangeset for help on using the changeset viewer.
hosted by ZAMG