WebSVN – lagranto.ecmwf – Diff – /trunk/select/select.f


      PROGRAM select
 
c     **************************************************************
c     * Select trajectories from LSL file                          *
c     * Michael Sprenger / January, February 2008                  *
c     **************************************************************
 
      implicit none
 
c     -------------------------------------------------------------- 
c     Declaration of parameters
c     --------------------------------------------------------------
 
c     Maximum number of columns per trajectory
      integer          maxcol
      parameter        (maxcol=100)
      
c     Maximum number of commands
      integer          maxcmd
      parameter        (maxcmd=10000)
 
c     -------------------------------------------------------------- 
c     Declaration of variables
c     --------------------------------------------------------------
 
c     Input and output files
      character*120     inp_lslfile                  ! Input lsl file
      character*120     out_lslfile                  ! Output lsl file
      character*120     inp_criteria                 ! Input file with criteria
      integer           inpmode                      ! Format of input file
      integer           outmode                      ! Format of output file
      character*80      outformat                    ! Trajectory/Boolean/Index of output   
      character*80      regionf                      ! Name of the regionfile
 
c     Trajectories
      integer          ntim                          ! Number of times 
      integer          ncol                          ! Number of columns (including time...)
      integer          ntrainp,ntraout               ! Number of trajectories
      character*80     vars(maxcol)                  ! Names of trajectory columns
      integer          basetime(6)                   ! Base time of trajectory (first line in lsl)
      real,allocatable, dimension (:,:,:) :: trainp  ! Input trajectories
      real,allocatable, dimension (:,:,:) :: traout  ! Output trajectories
      integer,allocatable,dimension (:)   :: selflag ! Flag for selection
      real,allocatable, dimension (:)     :: time    ! Times of the trajectory
      integer,allocatable,dimension (:,:) :: trigger ! Trigger column 
      integer          itrigger                      ! Column index for trigger
      character*80     addtrigger                    ! Flag whether to add trigger column
 
c     Command stack
      real             cmd(maxcmd)                   ! Decoded slection criterion
      integer          ncmd                          ! Number of commands
 
c     Common block for initialisation of polygon check
      real    tlonv(2000),vlat_c(2000),vlon_c(2000)
      real    xlat_c,xlon_c
      integer ibndry,nv_c
      data    ibndry/0/
      common  /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
c     Auxiliary variables
      integer          stat                          ! Logical (state) variable
      integer          fid,fod,fcr                   ! File identifier for input and output
      integer          i,j,k                         ! Index counter
      integer          isok                          ! Flag for selected trajectory
      real             param(1000)                   ! List of parameters
      integer          nparam                        ! Number of parameters
      character*80     specialstr                    ! Name of special command
      integer          len
      integer,allocatable,dimension (:) :: trigger1  ! Trigger column 
      real,allocatable,dimension (:,:) :: trainp1    ! A single trajectory
      character        ch
 
c     -------------------------------------------------------------- 
c     Preparations
c     --------------------------------------------------------------
 
c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM SELECT ***'
      print*
 
c     Read parameter file
      open(10,file='select.param')
       read(10,*) inp_lslfile
       read(10,*) out_lslfile
       read(10,*) outformat
       read(10,*) inp_criteria
       read(10,*) ntrainp
       read(10,*) ntim
       read(10,*) ncol
       read(10,*) regionf
       read(10,*) addtrigger
      close(10)
      
c     Set the formats of the input and output files
      call mode_tra(inpmode,inp_lslfile)
      if (inpmode.eq.-1) inpmode=1
      call mode_tra(outmode,out_lslfile)
      if ( (outmode.eq.-1).and.(outformat.ne.'startf') ) then
         outmode=1
      endif
 
c     Allocate memory for a single trajectory
      allocate(trainp(ntrainp,ntim,ncol),stat=stat)
      if (stat.ne.0) stop '*** error allocating array trainp   ***'
      allocate(time(ntim),stat=stat)
      if (stat.ne.0) stop '*** error allocating array time     ***'
      allocate(selflag(ntrainp),stat=stat)
      if (stat.ne.0) stop '*** error allocating array selflag  ***'
      allocate(trigger(ntrainp,ntim),stat=stat)
      if (stat.ne.0) stop '*** error allocating array trigger  ***'
      allocate(trigger1(ntim),stat=stat)
      if (stat.ne.0) stop '*** error allocating array trigger1 ***'
      allocate(trainp1(ntim,ncol),stat=stat)
      if (stat.ne.0) stop '*** error allocating array trainp1  ***'
 
c     Read the input trajectory file
      call ropen_tra(fid,inp_lslfile,ntrainp,ntim,ncol,
     >               basetime,vars,inpmode)
      call read_tra (fid,trainp,ntrainp,ntim,ncol,inpmode)
      call close_tra(fid,inpmode)
 
c     Check that first four columns correspond to time,lon,lat,p
      if ( (vars(1).ne.'time' ).or.
     >     (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.
     >     (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.
     >     (vars(4).ne.'ppos' ).and.(vars(4).ne.'p'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
      vars(1) = 'time'
      vars(2) = 'lon'
      vars(3) = 'lat'
      vars(4) = 'p'
 
c     Get the trajectory times from first trajectory
      do i=1,ntim
         time(i)=trainp(1,i,1)
      enddo
 
c     Init the trigger field - first check whether it is already available
      itrigger = 0
      do i=1,ncol
         if ( vars(i).eq.'TRIGGER' ) itrigger = i
      enddo
 
      if ( itrigger.ne.0 ) then
         do i=1,ntrainp
            do j=1,ntim
               trigger(i,j) = nint( trainp(i,j,itrigger) )
            enddo
         enddo
      else
         do i=1,ntrainp
            do j=1,ntim
               trigger(i,j) = 0
            enddo
         enddo
      endif
 
c     Write some info about the trajectory
      print*,'---- INPUT PARAMETERS -----------------------------------'
      write(*,*) 
      write(*,*) 'Input file    : ',trim(inp_lslfile)
      write(*,*) 'Output file   : ',trim(out_lslfile)
      write(*,*) 'Output format : ',trim(outformat)
      write(*,*) 'Criteria file : ',trim(inp_criteria)
      write(*,*) '# tra         : ',ntrainp
      write(*,*) '# time        : ',ntim
      write(*,*) '# col         : ',ncol
      write(*,*) 'Region file   : ',trim(regionf)
      write(*,*) 'Add trigger   : ',trim(addtrigger)
      
      print*
      print*,'---- INPUT TRAJECTORY FILE ------------------------------'
      print*
      write(*,'(1x,a12,i4,a10)') 'Vars       : ',1,trim(vars(1))
      do i=2,ncol
         write(*,'(1x,a12,i4,a10)') '             ',i,trim(vars(i))
      enddo
      print*
      write(*,'(1x,a12,i4,f10.2)') 'Time       : ',1,time(1)
      do i=2,ntim
         write(*,'(1x,a12,i4,f10.2)') '             ',i,time(i)
      enddo
      print*
      write(*,'(1x,a12,i4,i10)') 'Base date  : ',1,basetime(1)
      write(*,'(1x,a12,i4,i10)') '             ',2,basetime(2)
      write(*,'(1x,a12,i4,i10)') '             ',3,basetime(3)
      write(*,'(1x,a12,i4,i10)') '             ',4,basetime(4)
      write(*,'(1x,a12,i4,i10)') '             ',5,basetime(5)      
      print*
      write(*,'(1x,a12,i4,i10)') 'Time range : ',6,basetime(6)
      print*
      if ( itrigger.ne.0 ) then
         print*,'TRIGGER FIELD FOUND IN COLUMN ',itrigger
         print*
      endif
 
c     Read and decode the selection criterion
      fcr = 10
      open(fcr,file=inp_criteria)
      ncmd=maxcmd
      call decode(fcr,cmd,ncmd,vars,ncol,time,ntim,regionf)
      close(fcr)
 
      print*
      print*,'---- PSEUDO CODE FOR SELECTION --------------------------'
      print*
      call dumpcode(cmd,ncmd,vars,ncol,time,ntim)
 
c     -------------------------------------------------------------- 
c     Loop over all trajectories - selection
c     --------------------------------------------------------------
 
c     Prepare string and parameters for SPECIAL commands
      if ( cmd(1).eq.0 ) then
 
c        Get command string
         j   = 2
         len = nint(cmd(j))
         specialstr = ''
         do k=1,len
            j = j + 1
            specialstr = trim(specialstr)//char(nint(cmd(j)))
         enddo
         
c        Get paramters
         j      = j + 1
         nparam = nint(cmd(j))
         do k=1,nparam
            j        = j + 1
            param(k) = cmd(j)
         enddo
         
      endif
        
c     Init the counter for selected trajectories
      ntraout = 0
 
c     Loop over all trajectories
      do i=1,ntrainp
 
c       Copy a single trajectory to <trainp1> and <trigger1>
        do j=1,ntim
           do k=1,ncol
              trainp1(j,k) = trainp(i,j,k)
           enddo
           trigger1(j) = trigger(i,j)
        enddo
 
C	    Skip the trajectory if missing data are found for positions
        isok = 1
        do j=1,ntim
           if ( trainp1(j,4).lt.0. ) isok = 0
        enddo
 
c       Decide whether the trajectory is selected (handle SPECIAL commands)
        if ( isok.eq.1 ) then
           if (cmd(1).ne.0 ) then
              call select_tra (isok,cmd,ncmd,trainp1,trigger1,ntim,ncol)
 
           else
              call special    (isok,specialstr,trainp1,ntim,ncol,
     >                         vars,time,param,nparam)
           endif
	    endif
 
c       The trigger might be changed in the selection - copy it
        do j=1,ntim
            trigger(i,j) = trigger1(j)
        enddo
 
c       Set flag for selected trajectories
        if (isok.eq.1) then
           selflag(i)          = 1
           ntraout             = ntraout + 1
        else
           selflag(i)          = 0
        endif
 
      enddo
 
c     -------------------------------------------------------------- 
c     Write output trajectories
c     --------------------------------------------------------------
 
c     ------ Write output trajectories -----------------------------
      if ( outformat.eq.'trajectory' ) then
 
 
c        Define the trigger field if it is not yet defined
         if ( ( addtrigger.eq.'-trigger' ).and.(itrigger.eq.0) ) then
            ncol       = ncol + 1
            vars(ncol) = 'TRIGGER'
            itrigger   = ncol
         endif
 
c        Allocate memory for output trajectory
         allocate(traout(ntraout,ntim,ncol),stat=stat)
         if (stat.ne.0) stop '*** error allocating array apply   ***'
 
c        Set output trajectories
         j = 0
         do i=1,ntrainp
            if (selflag(i).eq.1) then
               j             = j + 1 
               traout(j,1:ntim,1:ncol) = trainp(i,1:ntim,1:ncol)
               if ( itrigger.ne.0 ) then
                  traout(j,1:ntim,itrigger) = real(trigger(i,1:ntim))
               endif
            endif
         enddo
 
c        Write trajectories
         call wopen_tra(fod,out_lslfile,ntraout,ntim,ncol,
     >        basetime,vars,outmode)
         call write_tra(fod,traout,ntraout,ntim,ncol,outmode)
         call close_tra(fod,outmode)
         
c     ------ Write index list -------------------------------------
      elseif ( outformat.eq.'index' ) then
         
         open(10,file=out_lslfile)
          do i=1,ntrainp
             if ( selflag(i).eq.1) write(10,*) i
          enddo
         close(10)
 
c     ------ Write boolean list -----------------------------------
      elseif ( outformat.eq.'boolean' ) then
         
         open(10,file=out_lslfile)
          do i=1,ntrainp
             write(10,'(i1)') selflag(i)
          enddo
         close(10)
 
c     ------ Write count -------------------------------------------
      elseif ( outformat.eq.'count' ) then
 
         open(10,file=out_lslfile)
          write(10,'(i7)') ntraout
         close(10)
 
c     ------ Write starting positions -----------------------------
      elseif ( outformat.eq.'startf' ) then
 
c        Allocate memory for output trajectory
         allocate(traout(ntraout,1,ncol),stat=stat)
         if (stat.ne.0) stop '*** error allocating array apply   ***'
 
c        Set output trajectories
         j = 0
         do i=1,ntrainp
            if (selflag(i).eq.1) then
               j             = j + 1 
               traout(j,1,:) = trainp(i,1,:)
            endif
         enddo
 
c        Write trajectories
         if (outmode.ne.-1) then
            
           call wopen_tra(fod,out_lslfile,ntraout,1,ncol,
     >                    basetime,vars,outmode)
           call write_tra(fod,traout,ntraout,1,ncol,outmode)
           call close_tra(fod,outmode)       
 
c        Output as a triple list (corresponding to <startf> file)
         else
         
           fid = 10
           open(fid,file=out_lslfile)
            do i=1,ntraout
              write(fid,'(3f10.3)') traout(i,1,2),   ! longitude
     >                            traout(i,1,3),   ! latitude
     >                            traout(i,1,4)    ! pressure
            enddo
           close(fid)
 
        endif
  
      endif
 
c     Write some status information, and end of program message
      print*  
      print*,'---- STATUS INFORMATION --------------------------------'
      print*
      print*,' # input trajectories  : ',ntrainp
      print*,' # output trajectories : ',ntraout
      print*
      print*,'              *** END OF PROGRAM SELECT ***'
      print*,'========================================================='
 
      stop
      
c     Exception handling
 100  stop 'select: First column in input trajectory must be <time>'
 101  stop 'select: Input trajectory file is empty'
 
      end
 
c     -------------------------------------------------------------- 
c     Dump the command list
c     --------------------------------------------------------------
      
      subroutine dumpcode(out,n,vars,nvars,times,ntimes)
 
c     Write the command list to screen. The command list is decoded
c     by call to <decode>
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      n
      real         out(n)
      integer      nvars
      character*80 vars(nvars)
      integer      ntimes
      real         times(ntimes)
 
c     A single command
      character*80 cmd
      character*80 var,mode,strtim
      integer      nval
      integer      ntim
      integer      ivar,imode,icmd,itime
 
c     Auxiliary variables
      integer      i,j
 
c     Loop through the complete list
      i=0
 100  if (i.lt.n) then 
 
         write(*,*) '---------------------------------------'
 
c        Get command
         i=i+1
         icmd=nint(out(i))
 
c        Special handling of SPECIAL commands
         if ( icmd.eq.0 ) then
 
c           Write 'header' for SPECIAL command
            write(*,'(i5,f15.4,10x,a10)') i,out(i),'SPECIAL'
 
c           Write command string
            i = i + 1
            icmd = nint(out(i))
            write(*,'(i5,f15.4,10x,a10)') i,out(i),'LEN(CMD)'
            do j=1,icmd
               i    = i + 1
               ivar = nint(out(i))
               write(*,'(i5,f15.4,10x,a10)') i,out(i),char(ivar)
            enddo
 
c           Write parameters
            i    = i + 1
            nval = nint(out(i))
            write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'
            do j=1,nval
               i=i+1
               if ( var.ne.'INPOLYGON' ) then
                  write(*,'(i5,f15.4)') i,out(i)
               else
                  write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))
               endif
            enddo
 
c           Nothing else to do - exit
            goto 200
            
         endif
 
c        Set the command
         if (icmd.eq.  1) cmd='GT'
         if (icmd.eq.  2) cmd='LT'
         if (icmd.eq.  3) cmd='IN'
         if (icmd.eq.  4) cmd='OUT'
         if (icmd.eq.  5) cmd='EQ'
         if (icmd.eq.  6) cmd='TRUE'
         if (icmd.eq.  7) cmd='FALSE'
         if (icmd.eq.  8) cmd='ALL'
         if (icmd.eq.  9) cmd='ANY'
         if (icmd.eq. 10) cmd='NONE'
         if (icmd.eq. -1) cmd='BEGIN'
         if (icmd.eq. -2) cmd='END'
         if (icmd.eq. -3) cmd='AND'
         if (icmd.eq. -4) cmd='OR'
         write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(cmd)
         if (icmd.lt.0) goto 100
 
c        Get variable
         i=i+1
         ivar=nint(out(i))
         if ( ivar.eq. -1 ) then
            var = 'DIST'
         elseif ( ivar.eq. -2 ) then
            var = 'DIST0'
         elseif ( ivar.eq. -3 ) then
            var = 'INPOLYGON'
         elseif ( ivar.eq. -4 ) then
            var = 'INBOX'
         elseif ( ivar.eq. -5 ) then
            var = 'INCIRCLE'
         elseif ( ivar.eq. -6 ) then
            var = 'INREGION'
         elseif ( ivar.eq. -7 ) then
            var = 'TRIGGER'
         elseif ( ivar.eq. -8 ) then
            var = 'VERT0'
         else
            var=vars(ivar)
         endif
         write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(var)
 
c        Get variable mode
         i=i+1
         imode=nint(out(i))
         if (imode.eq.1) mode='VALUE'
         if (imode.eq.2) mode='MEAN'         
         if (imode.eq.3) mode='MAX'
         if (imode.eq.4) mode='MIN'
         if (imode.eq.5) mode='VAR'
         if (imode.eq.6) mode='SUM'
         if (imode.eq.7) mode='CHANGE'
         if (imode.eq.8) mode='DIFF'
         if (imode.eq.9) mode='RANGE'
         write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(mode)
 
c        Get values
         i=i+1
         nval=nint(out(i))
         write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'
         do j=1,nval
            i=i+1
 
            if ( var.ne.'INPOLYGON' ) then
               write(*,'(i5,f15.4)') i,out(i)
            else
               write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))
            endif
         enddo
 
c        Get the number of times
         i=i+1
         ntim=nint(out(i))
 
c        the number of times is variable - depending on TRIGGER
         if ( ntim .eq. -993 ) then
            write(*,'(i5,f15.4,7x,7x,a15)') i,out(i),'TIMES @ TRIGGER'
 
c        Get the detailed list of times
         else
            write(*,'(i5,f15.4,7x,7x,a6)') i,out(i),'#TIMES'
            do j=1,ntim
               i=i+1
               write(*,'(i5,f15.4,f7.0)') i,out(i),times(nint(out(i)))
            enddo
         endif
 
c        Get time mode
         i=i+1
         itime=nint(out(i))
         if (itime.eq.1) strtim='ALL'
         if (itime.eq.2) strtim='ANY'         
         if (itime.eq.3) strtim='NONE'
         if (itime.lt.0) strtim='TRIGGER'
 
         if ( strtim.ne.'TRIGGER' ) then
            write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(strtim)
         else
            write(*,'(i5,f15.4,10x,a10,a3,i3)') i,out(i),trim(strtim),
     >                                        ' ->',abs(itime)
         endif
 
         goto 100
         
      endif
 
c     Exit point
 200  continue
 
      write(*,*) '---------------------------------------'
 
      end
         
 
c     -------------------------------------------------------------- 
c     Read and decode a selection set
c     --------------------------------------------------------------
 
      subroutine decode(fid,out,n,vars,nvars,times,ntimes,regionf)
 
c     A selection file is opened with file id <fid> and transformed
c     into a set of commands applied to the trajectories. On input
c     <n> sets the maximum dimension of <out>, on output it gives the
c     total length of the command string. The output is a list of
c     commands with the following format:
c     
c        out(i)             = Command 
c        out(i+1)           = Column index of variable
c        out(i+2)           = Mode for variable
c        out(i+3)           = Number of parameters (n)
c        out(i+4:i+4+n)     = Parameters
c        out(i+5+n)         = Number of times 
c        out(i+6+n:i+6+n+m) = List of time indices (m)
c        out(i+7+n+m)       = Time mode
c     
c     For SPECIAL commands (to be coded in <special.f>) the format is:
c     
c        out(i)             = Length of command string (n)
c        out(i+1:i+n)       = Command string      
c        out(i+n+1)         = Number of parameters (m)
c        out(i+n+2:i+n+1+m) = List of parameters
c     
c     The following coding is used
c      
c        Command        Variable mode    Time mode
c        ----------     -------------    ---------
c        GT       1     VALUE    1       ALL      1
c        LT       2     MEAN     2       ANY      2
c        IN       3     MAX      3       NONE     3
c        OUT      4     MIN      4       TRIGGER -i (i the trigger index)
c        EQ       5     VAR      5
c        BEGIN   -1     SUM      6
c        END     -2     CHANGE   7 
c        AND     -3     DIFF     8
c        OR      -4     RANGE    9
c        TRUE     6
c        FALSE    7
c        SPECIAL  0
c        ALL      8 (TRIGGER)
c        ANY      9 (TRIGGER)
c        NONE    10 (TRIGGER)
 
 
c     Several "implicit variables" are supported - out(i+1):
c     
c        DIST      -1 : Path length of the trajectory
c        DIST0     -2 : Distance from starting position
c        INPOLYGON -3 : Specified polygon region
c        INBOX     -4 : Longitude/latitude rectangle
c        INCIRCLE  -5 : Within a specified radius 
c        INREGION  -6 : Within a specified rehion on the region file
c        TRIGGER   -7 : Trigger field
c        VERT0     -8 : Vertical distance from starting position
c         
c     For the special commands BEGIN, END, AND and OR, only one field
c     in <out> is used.
 
      implicit none
 
c     Declaration of subroutine parameters
      integer      fid
      integer      n
      real         out(n)
      integer      nvars
      character*80 vars(nvars)
      integer      ntimes
      real         times(ntimes)
      character*80 regionf
 
c     Numerical epsilon
      real         eps
      parameter    (eps=0.001)
 
c     A single command term
      character*80 cmd
      character*80 var,mode
      integer      nval
      real         val(1000)
      integer      ntim
      real         tim(1000)
      character*80 tmode
 
c     Specification of a polygon
      character*80 filename
      integer      pn                             ! Number of entries in lat/lon poly
      real         latpoly(1000)                  ! List of polygon latitudes
      real         lonpoly(1000)                  ! List of polygon longitudes
      real         loninpoly,latinpoly            ! Lon/lat inside polygon
 
c     Specification of a region
      real         xcorner(4)
      real         ycorner(4)
      integer      iregion
      character*80 string
 
c     Transformation to UPN (handling of logical operators)
      integer      nlogical
      integer      ilogical(n)
      real         tmp(n)
      integer      mlogical
      integer      isor
 
c     Auxiliary variables
      integer      i,j
      integer      j1,j2
      integer      flag(ntimes)
      integer      count
      integer      ok
      integer      itrigger,ttrigger
 
c     Common block for initialisation of polygon check
      real    tlonv(1000),vlat_c(1000),vlon_c(1000)
      real    xlat_c,xlon_c
      integer ibndry,nv_c
      common  /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
c     ------ Decode single commands ---------------------
 
c     Reset the filename for polygons
      filename='nil'
 
c     Reset the counter for logical commands
      nlogical=0
 
c     Loop through all commands
 100  continue
 
c       Read next command (handle special cases)
        read(fid,*) cmd
 
c       Special handling of SPECIAL commands
        if ( cmd.eq.'SPECIAL' ) then
 
c          Set the flag for SPECIAL command
           n      = 1
           out(n) = 0.
 
c          Read the command string
           read(fid,*) var
 
c          Add the command string
           n      = n + 1
           out(n) = len_trim(var)
           do j=1,len_trim(var)
              n      = n + 1
              out(n) = ichar(var(j:j))
           enddo
 
c          Read the parameters
           read(fid,*) nval
           read(fid,*) (val(i),i=1,nval)
 
c          Add the parameters
           n = n + 1
           out(n)=real(nval)
           do i=1,nval
              n=n+1
              out(n)=val(i)
           enddo
 
c          Goto exit point - nothing more top do
           goto 350
 
        endif
 
c       Handle structure commands
        if ( cmd.eq.'BEGIN') then
           out(1)=-1.
           n=1
           nlogical=1
           ilogical(1)=n
           goto 200
        elseif ( cmd.eq.'AND'  ) then
           n=n+1
           out(n)=-3.
           nlogical=nlogical+1
           ilogical(nlogical)=n
           goto 200
        elseif ( cmd.eq.'OR'   ) then
           n=n+1
           out(n)=-4.
           nlogical=nlogical+1
           ilogical(nlogical)=n
           goto 200
        elseif ( cmd.eq.'END'  ) then
           n=n+1
           out(n)=-2.
           nlogical=nlogical+1
           ilogical(nlogical)=n
           goto 300
        endif
 
c       Read other fields associated with the command
        read(fid,*) var,mode
 
c       Read parameter
        if ( var.eq.'INPOLYGON' ) then 
           read(fid,*) nval
           read(fid,*) filename
           filename = trim(filename)
        else
           read(fid,*) nval
           read(fid,*) (val(i),i=1,nval)
        endif
 
c       Read times (on request, change to special trigger times)
        read(fid,*) ntim
        if ( ntim.eq.-993 ) then
           ttrigger = 1
           ntim     = 1
           tim(1)   = -999.
        else
           ttrigger = 0
           read(fid,*) (tim(i),i=1,ntim)
        endif
        read(fid,*) tmode
 
c       Bring CAPITAL "TIME,LAT,LON,P" into "time,lat,lon,p"
        if (var.eq.'TIME') var='time'
        if (var.eq.'LAT' ) var='lat'
        if (var.eq.'LON' ) var='lon'
        if (var.eq.'P'   ) var='p'
 
c       If the time mode is 'TRIGGER', all times of a trajectory
c       must be considered
        if ( tmode.eq.'TRIGGER' ) then
           itrigger = nint(tim(1))
           ntim     = 1
           tim(1)   = -999.
        endif
 
c       Special times: transform into real time
        do i=1,ntim
           if ( abs(tim(i)+996.).lt.eps ) tim(i)=times(1)
           if ( abs(tim(i)+995.).lt.eps ) tim(i)=times(ntimes)
        enddo
 
c       Check whether times are valid 
        ok=0
        do i=1,ntim
           if ( (abs(tim(i)+994.).gt.eps).and.
     >          (abs(tim(i)+999.).gt.eps) )
     >     then
              do j=1,ntimes
                 if ( abs(tim(i)-times(j)).lt.eps ) then
                    ok=ok+1
                 endif
              enddo
           else
              ok=ok+1
           endif
        enddo
        if (ok.ne.ntim) goto 400
        
c       Select all times which are included in the criterion
        do i=1,ntimes
           flag(i)=0
        enddo
        i=1
 150    if (i.le.ntim) then
           
c          A list of times
           if ( (abs(tim(i)+994.).lt.eps) ) then
              j1=0
              do j=1,ntimes
                 if ( abs(tim(i-1)-times(j)).lt.eps ) then
                    j1=j
                 endif
              enddo
              j2=0
              do j=1,ntimes
                 if ( abs(tim(i+1)-times(j)).lt.eps ) then
                    j2=j
                 endif
              enddo   
              if ( (j1.eq.0).or.(j2.eq.0) ) goto 400
              do j=j1,j2
                 flag(j)=1
              enddo
              i=i+1
 
c          Explicitly given time value
           else
              do j=1,ntimes
                 if ( abs(tim(i)-times(j)).lt.eps ) then
                    flag(j)=1
                 endif
              enddo
              
           endif
 
           i=i+1
           goto 150
 
        endif
 
c       Write command identifier
        n=n+1
        if (cmd.eq.'GT'    ) out(n)= 1.
        if (cmd.eq.'LT'    ) out(n)= 2.
        if (cmd.eq.'IN'    ) out(n)= 3.
        if (cmd.eq.'OUT'   ) out(n)= 4.
        if (cmd.eq.'EQ'    ) out(n)= 5.
        if (cmd.eq.'TRUE'  ) out(n)= 6.
        if (cmd.eq.'FALSE' ) out(n)= 7.
        if (cmd.eq.'ALL  ' ) out(n)= 8.
        if (cmd.eq.'ANY  ' ) out(n)= 9.
        if (cmd.eq.'NONE ' ) out(n)=10.
 
c       Write index for variable - force implicit trigger
        ok=0
        do j=1,nvars
           if (vars(j).eq.var) ok=j
        enddo
 
        if (var.eq.'TRIGGER') ok = 0
        
        if (ok.eq.0) then
           if (var.eq.'DIST') then
              ok = -1
           elseif (var.eq.'DIST0') then
              ok = -2
           elseif (var.eq.'INPOLYGON') then
              ok = -3
           elseif (var.eq.'INBOX') then
              ok = -4
           elseif (var.eq.'INCIRCLE') then
              ok = -5
           elseif (var.eq.'INREGION') then
              ok = -6
           elseif (var.eq.'TRIGGER') then
              ok = -7
           elseif (var.eq.'VERT0') then
              ok = -8
           else
             goto 400
          endif
        endif
        n=n+1
        out(n)=real(ok)
 
c       Write mode for variable
        ok=0
        if (mode.eq.'VALUE'  ) ok=1
        if (mode.eq.'MEAN'   ) ok=2
        if (mode.eq.'MAX'    ) ok=3
        if (mode.eq.'MIN'    ) ok=4
        if (mode.eq.'VAR'    ) ok=5
        if (mode.eq.'SUM'    ) ok=6
        if (mode.eq.'CHANGE' ) ok=7
        if (mode.eq.'DIFF'   ) ok=8
        if (mode.eq.'RANGE'  ) ok=9
        if (ok.eq.0) goto 400
        n=n+1
        out(n)=real(ok)
 
c       Write the parameter values: INPOLYGON 
        if ( var.eq.'INPOLYGON' ) then
 
           n      = n+1
           out(n) = len_trim(filename)
           do j=1,len_trim(filename)
              n      = n + 1
              out(n) = ichar(filename(j:j))
           enddo
 
c       Write parameter value: INREGION 
        elseif ( var.eq.'INREGION' ) then
 
           iregion = nint(val(1))
 
           open(fid+1,file=regionf)          
 
 50        read(fid+1,*,end=51) string
 
           if ( string(1:1).ne.'#' ) then
              call regionsplit(string,i,xcorner,ycorner)
              if ( i.eq.iregion ) goto 52
           endif
 
           goto 50
           
 51        close(fid+1)        
 
           print*,' ERROR: region ',iregion,' not found on ',
     >                                         trim(regionf)
           stop
 
 52        continue
           
           n      = n + 1
           out(n) = 8                   ! Number of parameters
           do i=1,4
              n      = n + 1
              out(n) = xcorner(i)
           enddo
           do i=1,4
              n      = n + 1
              out(n) = ycorner(i)
           enddo
 
c       Write parameter values: all other cases
        else
           n=n+1
           out(n)=real(nval)
           do i=1,nval
              n=n+1
              out(n)=val(i)
           enddo
        endif
 
c       Special time handling: only trigger times are cosidered
        if ( ttrigger.eq.1 ) then
           n = n+1
           out(n)=-993.
 
c       All times are selected
        elseif ( abs(tim(1)+999.).lt.eps ) then
           n=n+1
           out(n)=real(ntimes)
           do i=1,ntimes
              n=n+1
              out(n)=real(i)
           enddo
 
c       A selection of times is given 
        else
           count=0
           do i=1,ntimes
              count=count+flag(i)
           enddo
           n=n+1
           out(n)=real(count)
           do i=1,ntimes
              if (flag(i).eq.1) then
                 n=n+1
                 out(n)=real(i)
              endif
           enddo
        endif
 
c       Write the time mode
        if ( tmode.eq.'ALL') then
           n=n+1
           out(n)=1.
        elseif ( tmode.eq.'ANY') then
           n=n+1
           out(n)=2.
        elseif ( tmode.eq.'NONE') then
           n=n+1
           out(n)=3.
        elseif ( tmode.eq.'TRIGGER') then
           n=n+1
           out(n)=-real(itrigger)
        endif
           
c     End loop: handle single command
 200  continue
      goto 100
 
c     End loop: loop over all commands 
 300  continue
 
c     ------ Read polygon file, if requested -----------
      if ( filename.ne.'nil' ) then
 
        print*
        print*,
     >     '---- POLYGON --------------------------------------------'
 
           print*
           print*,'Filename = ',trim(filename)
           print*
         
c        Read list of polygon coordinates from file
         pn = 0
         open(fid+1,file=filename)
           read(fid+1,*) loninpoly,latinpoly
           print*,'Inside (lon/lat) =',loninpoly,latinpoly
           print*
 510       continue
              pn = pn + 1
              read(fid+1,*,end=511) lonpoly(pn),
     >                              latpoly(pn)
 
              print*,pn,lonpoly(pn),latpoly(pn)
              
              goto 510
 511       continue
           pn = pn - 1
         close(fid+1)
 
c        Define the polygon boundaries
         call DefSPolyBndry(latpoly,lonpoly,pn,latinpoly,loninpoly)
 
      endif
 
c     ------ Transform to UPN --------------------------
 
c     Check whether logical commands are ok
      mlogical=nint(out(ilogical(1)))
      if ( mlogical.ne.-1) goto 400
      mlogical=nint(out(ilogical(nlogical)))
      if ( mlogical.ne.-2) goto 400
 
c     No transformation necessary if only one command
      if (nlogical.eq.2) goto 350
      
c     Copy the output to temporary list
      do i=1,n
         tmp(i)=out(i)
      enddo
 
c     Set BEGIN statement
      n=1
      out(n)=-1.
 
c     Reorder commands and transform to UPN
      isor=0
      do i=1,nlogical-1
 
c        Get the logical command
         mlogical=nint(out(ilogical(i)))
 
c        Connecting OR
         if (mlogical.eq.-4) then
            if (isor.eq.1) then
               n=n+1
               out(n)=-4.
            else
               isor=1
            endif
         endif
 
c        Put the command onto the stack
         do j=ilogical(i)+1,ilogical(i+1)-1
            n=n+1
            out(n)=tmp(j)
         enddo
 
c        Connecting AND operator
         if ( mlogical.eq.-3 ) then
            n=n+1
            out(n)=-3.
         endif
         
      enddo
      
c     Set final connecting OR
      if (isor.eq.1) then
         n=n+1
         out(n)=-4.
      endif
 
c     Set END statement
      n=n+1
      out(n)=-2.
 
c     ------ Exit point ---------------------------------
 
 350  continue
      return
 
c     ----- Exception handling --------------------------
 400  print*,'Invalid selection criterion... Stop'
      stop
 
      end
 
 
c     -------------------------------------------------------------- 
c     Decide whether a trajectory is selected or not
c     --------------------------------------------------------------
 
      subroutine select_tra (select,cmd,n,tra,trigger,ntim,ncol)
 
c     Decide whether a single trajectory is selected (<select=1>) or
c     is not selected <select=0> according to the selection criterion
c     given in <cmd(ncmd)>. The selection criterion <cmd(ncmd)> is 
c     returned from the call to the subroutine <decode>. The trajectory
c     is given in <tra(ntim,ncol)> where <ntim> is the number of times
c     and <ncol> is the number of columns.
c     
c     Important note: the structure of <tra(ntim,ncol)> must match to the
c     call parameter <vars,nvars,times,ntimes> in subroutine <decode>.
 
      implicit none
 
c     Declaration of subroutine parameters
      integer   select
      integer   n
      real      cmd(n)
      integer   ntim,ncol
      real      tra(ntim,ncol)
      integer   trigger(ntim)
 
c     Numerical epsilon (for test of equality)
      real      eps
      parameter (eps=0.000001)
 
c     A single command and the associated field
      integer   icmd,ivar,imode,itime,nsel,nval
      integer   time(ntim)
      real      param(100)
      real      var   (ntim)
      integer   intvar(ntim)
 
c     Boolean values for a single time, a single command and build-up
      integer   stack(100)
      integer   nstack
      integer   istrue(ntim)
      integer   decision
 
c     Auxiliary variables
      integer   i,j,k
      real      tmp,mea
      integer   istack1,istack2
      real      lat0,lon0,lat1,lon1
      real      length(ntim)
      integer   flag
      real      dist
      real      xcorner(4),ycorner(4)
      integer   iparam
      character ch
      real      varmin,varmax
      real      lev0,lev1
 
c     Common block for initialisation of polygon check
      real    tlonv(1000),vlat_c(1000),vlon_c(1000)
      real    xlat_c,xlon_c
      integer ibndry,nv_c
      common  /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
c     Externals
      real      sdis           ! Spherical distance
      external  sdis
      integer   inregion       ! Boolean flag for regions
      external  inregion
 
c     Reset the decision stack (with locical values)
      nstack=0
 
c     Loop through the complete command list
      i=0
 100  if (i.lt.n) then  
 
c        --- Get the command -------------------------------
         i=i+1
         icmd=nint(cmd(i))
 
c        --- Handle structural commands (BEGIN, END, AND, OR)
 
c        Handle BEGIN statement
         if ( icmd.eq.-1) then
            nstack=0
            goto 200
         endif
 
c        Handle END statement
         if (icmd.eq.-2) then
            goto 300
         endif
 
c        Handle AND statement
         if (icmd.eq.-3) then
            istack1=stack(nstack)
            nstack=nstack-1
            istack2=stack(nstack)
            if ((istack1.eq.1).and.(istack2.eq.1)) then
               stack(nstack)=1
            else
               stack(nstack)=0
            endif
            goto 200
         endif
            
c        Handle OR statement
         if (icmd.eq.-4) then
            istack1=stack(nstack)
            nstack=nstack-1
            istack2=stack(nstack)
            if ((istack1.eq.1).or.(istack2.eq.1)) then
               stack(nstack)=1
            else
               stack(nstack)=0
            endif
            goto 200
         endif
 
c        --- Get all command details (parameters, modes, times)
 
c        Get variable (<ivar> gets the column index in <tra>)
         i=i+1
         ivar=nint(cmd(i))
         
c        Get variable mode 
         i=i+1
         imode=nint(cmd(i))
 
c        Get parameter values
         i=i+1
         nval=nint(cmd(i))
         do j=1,nval
            i=i+1
            param(j)=cmd(i)
         enddo
         
c        Get times (<time(j)> gets the row indices of <tra>)
         i=i+1
         nsel=nint(cmd(i))
         if ( nsel .eq. -993 ) then
            nsel = 0
            do k=1,ntim
              if ( trigger(k).ne.0 ) then
                nsel       = nsel + 1
                time(nsel) = k
              endif
           enddo
         else
           do j=1,nsel
              i=i+1
              time(j)=nint(cmd(i))
           enddo
         endif
 
c        If no times are selected, exit with non-select status
         if ( nsel.eq.0 ) then
             stack(1) = 0
             goto 300
         endif
 
c        Get time mode
         i=i+1
         itime=nint(cmd(i))
 
c        --- Prepare field values for analysis -----------
 
c        Implicit variable: DIST
         if ( ivar.eq. -1 ) then
            length(1) = 0.
            do j=2,ntim
               lon0      = tra(j-1,2)
               lat0      = tra(j-1,3)
               lon1      = tra(j  ,2)
               lat1      = tra(j  ,3)
               length(j) = length(j-1) + sdis(lon0,lat0,lon1,lat1)
            enddo
            do j=1,nsel
               var(j) = length( time(j) )
            enddo
 
               
c        Implict variable: DIST0
         elseif ( ivar.eq. -2 ) then
            do j=1,nsel
               lon0   = tra(1      ,2)
               lat0   = tra(1      ,3)
               lon1   = tra(time(j),2)
               lat1   = tra(time(j),3)
               var(j) = sdis(lon0,lat0,lon1,lat1)
            enddo
 
c        Implict variable: INPOLYGON
         elseif ( ivar.eq. -3 ) then
            do j=1,nsel
               lon1   = tra(time(j),2)
               lat1   = tra(time(j),3)               
               call LctPtRelBndry(lat1,lon1,flag)
               if ( (flag.eq.1).or.(flag.eq.2) ) then
                  var(j) = 1.
               else
                  var(j) = 0.
               endif
            enddo
 
c        Implict variable: INBOX 
         elseif ( ivar.eq. -4 ) then
            do j=1,nsel
               lon1   = tra(time(j),2)
               lat1   = tra(time(j),3)               
               if ( ( lon1.ge.param(1) ).and.       ! lonmin
     >              ( lon1.le.param(2) ).and.       ! lonmax
     >              ( lat1.ge.param(3) ).and.       ! latmin
     >              ( lat1.le.param(4) ) )          ! latmax
     >         then
                  var(j) = 1
               else
                  var(j) = 0
               endif
            enddo
 
c        Implict variable: INCIRCLE (lonc=param(1),latc=param(2),radius=param(3))
         elseif ( ivar.eq. -5 ) then
            do j=1,nsel
               lon1   = tra(time(j),2)
               lat1   = tra(time(j),3)               
               
               dist = sdis( lon1,lat1,param(1),param(2) ) 
 
               if ( dist.le.param(3) ) then
                  var(j) = 1
               else
                  var(j) = 0
               endif
            enddo
 
c        Implict variable: INREGION (xcorner=param(1..4),ycorner=param(5..8) )
         elseif ( ivar.eq.-6 ) then
            
            do j=1,4
               xcorner(j) = param(j  )
               ycorner(j) = param(j+4)
            enddo
 
            do j=1,nsel
               lon1   = tra(time(j),2)
               lat1   = tra(time(j),3)  
               var(j) = inregion (lon1,lat1,xcorner,ycorner)               
            enddo
 
c        Implict variable: TRIGGER
         elseif ( ivar.eq. -7 ) then
            do j=1,nsel
               intvar(j) = trigger( time(j) ) 
            enddo
 
c        Implicit variable: VERT0
         elseif ( ivar.eq. -8 ) then
            do j=1,nsel
               lev0   = tra(1      ,4)
               lev1   = tra(time(j),4)
               var(j) = lev0 - lev1
            enddo
 
c        Explicit variable (column index <ivar>)
         else
            do j=1,nsel
               var(j) = tra(time(j),ivar)
            enddo
 
         endif
 
c        Take MEAN of the variable (mean of selected times)
         if (imode.eq.2) then
            tmp=0.
            do j=1,nsel
               tmp=tmp+var(j)
            enddo
            var(1)=tmp/real(nsel)
            nsel=1
 
c        Take MAX of the variable (maximum of selected times)
         elseif (imode.eq.3) then
            tmp=var(1)
            do j=2,nsel
               if (var(j).gt.tmp) tmp=var(j)
            enddo
            var(1)=tmp
            nsel=1
            
c        Take MIN of the variable (minimum of selected times)
         elseif (imode.eq.4) then
            tmp=var(1)
            do j=2,nsel
               if (var(j).lt.tmp) tmp=var(j)
            enddo
            var(1)=tmp
            nsel=1
 
c        Take VAR of the variable (variance over all selected times)
         elseif (imode.eq.5) then
            tmp=0.
            do j=1,nsel
               tmp=tmp+var(j)
            enddo
            mea=tmp/real(nsel)
            do j=1,nsel
              tmp=tmp+(var(j)-mea)**2
            enddo
            var(1)=1./real(nsel-1)*tmp
            nsel=1
 
c        Take SUM of the variable (sum over all selected times)
         elseif (imode.eq.6) then
            tmp=0.
            do j=1,nsel
               tmp=tmp+var(j)
            enddo
            var(1)=tmp
            nsel=1
            
c        Take CHANGE of the variable (absolute difference between first and last time)
         elseif (imode.eq.7) then
            var(1)=abs(var(1)-var(nsel))
            nsel=1
 
c        Take DIFF of the variable (first minus last time)
         elseif (imode.eq.8) then
            var(1)=var(1)-var(nsel)
            nsel=1
 
c        Take RANGE of the variable
         elseif (imode.eq.9) then
            varmax=var(1)
            varmin=var(1)
            do j=2,nsel
               if (var(j).gt.varmax) varmax=var(j)
               if (var(j).lt.varmin) varmin=var(j)
            enddo
            var(1) = varmax - varmin
            nsel=1
         endif
 
c        --- Apply the operators to the single values ---
 
         do j=1,nsel
 
c           GT
            if (icmd.eq.1) then
               if (var(j).gt.param(1)) then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
               
c           LT
            elseif (icmd.eq.2) then
               if (var(j).lt.param(1)) then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
 
c           IN
            elseif (icmd.eq.3) then
               if ( (var(j).gt.param(1)).and.
     >              (var(j).lt.param(2)) ) 
     >         then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
 
c           OUT
            elseif (icmd.eq.4) then
               if ( (var(j).lt.param(1)).or.
     >              (var(j).gt.param(2)) ) 
     >         then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
 
c           EQ
            elseif (icmd.eq.5) then
               if (abs(var(j)-param(1)).lt.eps) then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
 
c           TRUE
            elseif (icmd.eq.6) then
               if (abs(var(j)).lt.eps) then
                  istrue(j)=0
               else
                  istrue(j)=1
               endif
 
c           FALSE
            elseif (icmd.eq.7) then
               if (abs(var(j)).lt.eps) then
                  istrue(j)=1
               else
                  istrue(j)=0
               endif
 
c           ALL
            elseif (icmd.eq.8) then
               istrue(j) = 1
               do k=1,nval
                  iparam = nint(param(k))-1
                  if (btest(intvar(j),iparam).eqv..false.) then
                     istrue(j) = 0
                  endif
               enddo
               
c           ANY
            elseif (icmd.eq.9) then
               istrue(j) = 0
               do k=1,nval
                  iparam = nint(param(k))-1
                  if (btest(intvar(j),iparam).eqv..true.) then
                     istrue(j) = 1
                  endif
               enddo
 
c           NONE
            elseif (icmd.eq.10) then
               istrue(j) = 1
               do k=1,nval
                  iparam = nint(param(k))-1
                  if (btest(intvar(j),iparam).eqv..true.) then
                     istrue(j) = 0
                  endif
               enddo
   
            endif
 
         enddo
 
c        --- Determine the overall boolean value ----------
            
c        ALL
         if (itime.eq.1) then
            decision=1
            do j=1,nsel
               if (istrue(j).eq.0) then
                  decision=0
                  goto 110
               endif
            enddo
 110        continue
 
c        ANY
         elseif (itime.eq.2) then
            decision=0
            do j=1,nsel
               if (istrue(j).eq.1) then
                  decision=1
                  goto 120
               endif
            enddo
 120        continue
           
c        NONE
         elseif (itime.eq.3) then
            decision=1
            do j=1,nsel
               if (istrue(j).eq.1) then
                  decision=0
                  goto 130
               endif
            enddo
 130        continue
 
c        TRIGGER
         elseif (itime.lt.0) then
            decision=1
            do j=1,nsel
               if (istrue(j).eq.1) then
                  trigger(j) = ior( trigger(j), 2**(abs(itime)-1) )
               endif
            enddo
            
         endif
 
c        --- Put the new boolean value onto the stack
 
         nstack=nstack+1
         stack(nstack)=decision
 
c        Exit point for loop
 200     continue
         goto 100
 
      endif
 
c     Return the decision (selected or non-selected)
 300  continue
 
      select=stack(1)
 
      end
      
 
c     --------------------------------------------------------------------------
c     Split a region string and get corners of the domain
c     --------------------------------------------------------------------------
 
      subroutine regionsplit(string,iregion,xcorner,ycorner)
 
c     The region string comes either as <lonw,lone,lats,latn> or as <lon1,lat1,
c     lon2,lat2,lon3,lat3,lon4,lat4>: split it into ints components and get the
c     four coordinates for the region
      
      implicit none
 
c     Declaration of subroutine parameters
      character*80    string
      real            xcorner(4),ycorner(4)
      integer         iregion
 
c     Local variables
      integer         i,n
      integer         il,ir
      real            subfloat (80)
      integer         stat
      integer         len
 
c     ------- Split the string
      i    = 1
      n    = 0
      stat = 0
      il   = 1
      len  = len_trim(string)
 
 100  continue
 
c     Find start of a substring
      do while ( stat.eq.0 )
         if ( string(i:i).ne.' ' ) then
            stat = 1
            il   = i
         else
            i = i + 1
         endif
      enddo
 
c     Find end of substring
      do while ( stat.eq.1 )         
         if ( ( string(i:i).eq.' ' ) .or. ( i.eq.len ) ) then
            stat = 2
            ir   = i
         else
            i    = i + 1
         endif
      enddo
 
c     Convert the substring into a number
      if ( stat.eq.2 ) then
         n = n + 1
         read(string(il:ir),*) subfloat(n)
         stat = 0
      endif
 
      if ( i.lt.len ) goto 100
 
 
c     -------- Get the region number
      
      iregion = nint(subfloat(1))
 
c     -------- Get the corners of the region
      
      if ( n.eq.5 ) then     ! lonw(2),lone(3),lats(4),latn(5)
 
         xcorner(1) = subfloat(2)
         ycorner(1) = subfloat(4)
 
         xcorner(2) = subfloat(3)
         ycorner(2) = subfloat(4)
 
         xcorner(3) = subfloat(3)
         ycorner(3) = subfloat(5)
         
         xcorner(4) = subfloat(2)
         ycorner(4) = subfloat(5)
        
      elseif ( n.eq.9 ) then     ! lon1,lat1,lon2,lat2,lon3,lon4,lat4
 
         xcorner(1) = subfloat(2)
         ycorner(1) = subfloat(3)
 
         xcorner(2) = subfloat(4)
         ycorner(2) = subfloat(5)
 
         xcorner(3) = subfloat(6)
         ycorner(3) = subfloat(7)
         
         xcorner(4) = subfloat(8)
         ycorner(4) = subfloat(9)
 
      else
         
         print*,' ERROR: invalid region specification '
         print*,'     ',trim(string)
         stop
         
      endif
         
 
      end
 
c     --------------------------------------------------------------------------
c     Decide whether lat/lon point is in or out of region
c     --------------------------------------------------------------------------
      
      integer function inregion (lon,lat,xcorner,ycorner)
      
c     Decide whether point (lon/lat) is in the region specified by <xcorner(1..4),
c     ycorner(1..4).
      
      implicit none
      
c     Declaration of subroutine parameters
      real    lon,lat
      real    xcorner(4),ycorner(4)
 
c     Local variables
      integer flag
      real    xmin,xmax,ymin,ymax
      integer i
 
c     Reset the flag
      flag = 0
 
c     Set some boundaries
      xmax = xcorner(1)
      xmin = xcorner(1)
      ymax = ycorner(1)
      ymin = ycorner(1)
      do i=2,4
        if (xcorner(i).lt.xmin) xmin = xcorner(i)
        if (xcorner(i).gt.xmax) xmax = xcorner(i)
        if (ycorner(i).lt.ymin) ymin = ycorner(i)
        if (ycorner(i).gt.ymax) ymax = ycorner(i)
      enddo
 
c     Do the tests - set flag=1 if all tests pased
      if (lon.lt.xmin) goto 970
      if (lon.gt.xmax) goto 970
      if (lat.lt.ymin) goto 970
      if (lat.gt.ymax) goto 970
      
      if ((lon-xcorner(1))*(ycorner(2)-ycorner(1))-
     >    (lat-ycorner(1))*(xcorner(2)-xcorner(1)).gt.0.) goto 970
      if ((lon-xcorner(2))*(ycorner(3)-ycorner(2))-
     >    (lat-ycorner(2))*(xcorner(3)-xcorner(2)).gt.0.) goto 970
      if ((lon-xcorner(3))*(ycorner(4)-ycorner(3))-
     >    (lat-ycorner(3))*(xcorner(4)-xcorner(3)).gt.0.) goto 970
      if ((lon-xcorner(4))*(ycorner(1)-ycorner(4))-
     >    (lat-ycorner(4))*(xcorner(1)-xcorner(4)).gt.0.) goto 970
 
      flag = 1
 
c     Return the value
 970  continue
      
      inregion = flag
      
      return
      
      end
 
   
c     --------------------------------------------------------------------------
c     Spherical distance between lat/lon points                                                       
c     --------------------------------------------------------------------------
 
      real function sdis(xp,yp,xq,yq)
c
c     calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
c
      real      re
      parameter (re=6370.)
      real      pi180
      parameter (pi180=3.14159/180.)
      real      xp,yp,xq,yq,arg
 
      arg=sin(pi180*yp)*sin(pi180*yq)+
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
      if (arg.lt.-1.) arg=-1.
      if (arg.gt.1.) arg=1.
 
      sdis=re*acos(arg)
 
      end
 
 
c     ****************************************************************
c     * Given some spherical polygon S and some point X known to be  *
c     * located inside S, these routines will determine if an arbit- *
c     * -rary point P lies inside S, outside S, or on its boundary.  *
c     * The calling program must first call DefSPolyBndry to define  *
c     * the boundary of S and the point X. Any subsequent call to    *
c     * subroutine LctPtRelBndry will determine if some point P lies *
c     * inside or outside S, or on its boundary. (Usually            *
c     * DefSPolyBndry is called once, then LctPrRelBndry is called   *
c     * many times).                                                 *
c     *                                                              * 
c     * REFERENCE:            Bevis, M. and Chatelain, J.-L. (1989)  * 
c     *                       Maflaematical Geology, vol 21.         *
c     * VERSION 1.0                                                  *
c     ****************************************************************
 
      Subroutine DefSPolyBndry(vlat,vlon,nv,xlat, xlon)
 
c     ****************************************************************
c     * This mmn entry point is used m define ~e spheric~ polygon S  *
c     * and the point X.                                             *
c     * ARGUMENTS:                                                   *
c     * vlat,vlon (sent) ... vectors containing the latitude and     * 
c     *                      longitude of each vertex of the         *
c     *                      spherical polygon S. The ith.vertex is  *
c     *                      located at [vlat(i),vlon(i)].           *
c     * nv        (sent) ... the number of vertices and sides in the *
c     *                      spherical polygon S                     *
c     * xlat,xlon (sent) ... latitude and longitude of some point X  *
c     *                      located inside S. X must not be located *
c     *                      on any great circle that includes two   *
c     *                      vertices of S.                          *
c     *                                                              *
c     * UNITS AND SIGN CONVENTION:                                   *
c     *  Latitudes and longitudes are specified in degrees.          *
c     *  Latitudes are positive to the north and negative to the     *
c     *  south.                                                      *
c     *  Longitudes are positive to the east and negative to the     *
c     *  west.                                                       *
c     *                                                              * 
c     * VERTEX ENUMERATION:                                          * 
c     * The vertices of S should be numbered sequentially around the *
c     * border of the spherical polygon. Vertex 1 lies between vertex*
c     * nv and vertex 2. Neighbouring vertices must be seperated by  *
c     * less than 180 degrees. (In order to generate a polygon side  *
c     * whose arc length equals or exceeds 180 degrees simply        *
c     * introduce an additional (pseudo)vertex). Having chosen       *
c     * vertex 1, the user may number the remaining vertices in      *
c     * either direction. However if the user wishes to use the      *
c     * subroutine SPA to determine the area of the polygon S (Bevis *
c     * & Cambareri, 1987, Math. Geol., v.19, p. 335-346) then he or *
c     * she must follow the convention whereby in moving around the  *
c     * polygon border in the direction of increasing vertex number  *
c     * clockwise bends occur at salient vertices. A vertex is       *
c     * salient if the interior angle is less than 180 degrees.      *
c     * (In the case of a convex polygon this convention implies     *
c     * that vertices are numbered in clockwise sequence).           *
c     ****************************************************************
 
      implicit none
      
      integer mxnv,nv
 
c     ----------------------------------------------------------------
c     Edit next statement to increase maximum number of vertices that 
c     may be used to define the spherical polygon S               
c     The value of parameter mxnv in subroutine LctPtRelBndry must match
c     that of parameter mxnv in this subroutine, as assigned above.
c     ----------------------------------------------------------------
      parameter (mxnv=2000)
 
      real  vlat(nv),vlon(nv),xlat,xlon,dellon
      real  tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
      integer i,ibndry,nv_c,ip
 
      data ibndry/0/
      
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
      if (nv.gt.mxnv) then
         print *,'nv exceeds maximum allowed value'
         print *,'adjust parameter mxnv in subroutine DefSPolyBndry'
         stop
      endif
 
      ibndry=1                  ! boundary defined at least once (flag)
      nv_c=nv                   ! copy for named common
      xlat_c=xlat               ! . . . .
      xlon_c=xlon               !
 
      do i=1,nv
         vlat_c(i)=vlat(i)      ! "
         vlon_c(i)=vlon(i)      !
 
         call TrnsfmLon(xlat,xlon,vlat(i),vlon(i),tlonv(i))
 
         if (i.gt.1) then
            ip=i-1
         else
            ip=nv
         endif
         
         if ((vlat(i).eq.vlat(ip)).and.(vlon(i).eq.vlon(ip))) then
            print *,'DefSPolyBndry detects user error:'
            print *,'vertices ',i,' and ',ip,' are not distinct'
            print*,'lat ',i,ip,vlat(i),vlat(ip)
            print*,'lon ',i,ip,vlon(i),vlon(ip)            
            stop
         endif
 
         if (tlonv(i).eq.tlonv(ip)) then
            print *,'DefSPolyBndry detects user error:'
            print *,'vertices ',i,' & ',ip,' on same gt. circle as X'
            stop
         endif
 
         if (vlat(i).eq.(-vlat(ip))) then
            dellon=vlon(i)-vlon(ip)
            if (dellon.gt.+180.) dellon=dellon-360.
            if (dellon.lt.-180.) dellon=dellon-360.
            if ((dellon.eq.+180.0).or.(dellon.eq.-180.0)) then
               print *,'DefSPolyBndry detects user error:'
               print *,'vertices ',i,' and ',ip,' are antipodal'
               stop
            endif
         endif
      enddo
 
      return
      
      end
 
 
c     ****************************************************************
 
      Subroutine LctPtRelBndry(plat,plon,location)
 
c     ****************************************************************
 
c     ****************************************************************
c     * This routine is used to see if some point P is located       *
c     * inside, outside or on the boundary of the spherical polygon  *
c     * S previously defined by a call to subroutine DefSPolyBndry.  *
c     * There is a single restriction on point P: it must not be     *
c     * antipodal to the point X defined in the call to DefSPolyBndry*
c     * (ie.P and X cannot be seperated by exactly 180 degrees).     *
c     * ARGUMENTS:                                                   *  
c     * plat,plon (sent)... the latitude and longitude of point P    *
c     * location (returned)... specifies the location of P:          *
c     *                        location=0 implies P is outside of S  *
c     *                        location=1 implies P is inside of S   *
c     *                        location=2 implies P on boundary of S *
c     *                        location=3 implies user error (P is   *
c     *                                     antipodal to X)          *
c     * UNFfS AND SIGN CONVENTION:                                   * 
c     *  Latitudes and longitudes are specified in degrees.          *
c     *  Latitudes are positive to the north and negative to the     *
c     *  south.                                                      *    
c     *  Longitudes are positive to the east and negative to the     *
c     *  west.                                                       *
c     ****************************************************************
      
      implicit none
      
      integer mxnv
 
c     ----------------------------------------------------------------
c     The statement below must match that in subroutine DefSPolyBndry
c     ----------------------------------------------------------------
 
      parameter (mxnv=2000)
 
      real tlonv(mxnv),vlat_c(mxnv),vlon_c(mxnv),xlat_c,xlon_c
      real plat,plon,vAlat,vAlon,vBlat,vBlon,tlonA,tlonB,tlonP
      real tlon_X,tlon_P,tlon_B,dellon
      integer i,ibndry,nv_c,location,icross,ibrngAB,ibrngAP,ibrngPB
      integer ibrng_BX,ibrng_BP,istrike
 
      common/spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry
 
      if (ibndry.eq.0) then     ! user has never defined the bndry
         print*,'Subroutine LctPtRelBndry detects user error:'
         print*,'Subroutine DefSPolyBndry must be called before'
         print*,'subroutine LctPtRelBndry can be called'
         stop
      endif
 
      if (plat.eq.(-xlat_c)) then
         dellon=plon-xlon_c
         if (dellon.lt.(-180.)) dellon=dellon+360.
         if (dellon.gt.+180.) dellon=dellon-360.
         if ((dellon.eq.+180.0).or.(dellon.eq.-180.)) then
            print*,'Warning: LctPtRelBndry detects case P antipodal
     >           to X'
            print*,'location of P relative to S is undetermined'
            location=3
            return
         endif
      endif 
 
      location=0                ! default ( P is outside S)
      icross=0                  ! initialize counter
 
      if ((plat.eq.xlat_c).and.(plon.eq.xlon_c)) then
         location=1
         return
      endif
 
      
      call TrnsfmLon (xlat_c,xlon_c,plat,plon,tlonP)
 
      do i=1,nv_c              ! start of loop over sides of S 
 
         vAlat=vlat_c(i)
         vAlon=vlon_c(i)
         tlonA=tlonv(i)
 
         if (i.lt.nv_c) then
            vBlat=vlat_c(i+1)
            vBlon=vlon_c(i+1)
            tlonB=tlonv(i+1)
         else
            vBlat=vlat_c(1)
            vBlon=vlon_c(1)
            tlonB=tlonv(1)
         endif
         
         istrike=0
         
         if (tlonP.eq.tlonA) then
            istrike=1
         else
            call EastOrWest(tlonA,tlonB,ibrngAB)
            call EastOrWest(tlonA,tlonP,ibrngAP)
            call EastOrWest(tlonP,tlonB,ibrngPB)
            
 
            if((ibrngAP.eq.ibrngAB).and.(ibrngPB.eq.ibrngAB)) istrike=1
         endif
 
         
         if (istrike.eq.1) then
 
            if ((plat.eq.vAlat).and.(plon.eq.vAlon)) then
               location=2       ! P lies on a vertex of S
               return
            endif
            call TrnsfmLon(vAlat,vAlon,xlat_c,xlon_c,tlon_X)
            call TrnsfmLon(vAlat,vAlon,vBlat,vBlon,tlon_B)
            call TrnsfmLon(vAlat,vAlon,plat,plon,tlon_P)
            
            if (tlon_P.eq.tlon_B) then
               location=2       ! P lies on side of S
               return 
            else
               call EastOrWest(tlon_B,tlon_X,ibrng_BX)
               call EastOrWest(tlon_B,tlon_P,ibrng_BP)
               if(ibrng_BX.eq.(-ibrng_BP)) icross=icross+1
            endif
            
         endif
      enddo                     ! end of loop over the sides of S
 
 
c     if the arc XP crosses the boundary S an even number of times then P
c     is in S
 
      if (mod(icross,2).eq.0) location=1
 
      return
 
      end
 
 
c     ****************************************************************
      
      subroutine TrnsfmLon(plat,plon,qlat,qlon,tranlon)
 
c     ****************************************************************
c     * This subroutine is required by subroutines DefSPolyBndry &   *
c     * LctPtRelBndry. It finds the 'longitude' of point Q in a      *
c     * geographic coordinate system for which point P acts as a     *
c     * 'north pole'. SENT: plat,plon,qlat,qlon, in degrees.         *
c     * RETURNED: tranlon, in degrees.                               *
c     ****************************************************************
 
      implicit none
 
      real pi,dtr,plat,plon,qlat,qlon,tranlon,t,b
      parameter (pi=3.141592654,dtr=pi/180.0)
 
      if (plat.eq.90.) then
         tranlon=qlon
      else
         t=sin((qlon-plon)*dtr)*cos(qlat*dtr)
         b=sin(dtr*qlat)*cos(plat*dtr)-cos(qlat*dtr)*sin(plat*dtr)
     >    *cos((qlon-plon)*dtr)
         tranlon=atan2(t,b)/dtr
      endif
 
      return
      end
 
c     ****************************************************************
 
      subroutine EastOrWest(clon,dlon,ibrng)
 
c     ****************************************************************
c     * This subroutine is required by subroutine LctPtRelBndry.     *
c     * This routine determines if in travelling the shortest path   *
c     * from point C (at longitude clon) to point D (at longitude    *
c     * dlon) one is heading east, west or neither.                  *
c     * SENT: clon,dlon; in degrees. RETURNED: ibrng                 *
c     * (1=east,-1=west, 0=neither).                                 *
c     ****************************************************************
 
      implicit none
      real clon,dlon,del
      integer ibrng
      del=dlon-clon
      if (del.gt.180.) del=del-360.
      if (del.lt.-180.) del=del+360.
      if ((del.gt.0.0).and.(del.ne.180.)) then
         ibrng=-1               ! (D is west of C)
      elseif ((del.lt.0.0).and.(del.ne.-180.)) then
         ibrng=+1               ! (D is east of C)
      else
         ibrng=0                ! (D north or south of C)
      endif
      return
      end

Rev 3	Rev 9
1	`PROGRAM select`	1	`PROGRAM select`
2		2
3	`c **************************************************************`	3	`c **************************************************************`
4	`c * Select trajectories from LSL file *`	4	`c * Select trajectories from LSL file *`
5	`c * Michael Sprenger / January, February 2008 *`	5	`c * Michael Sprenger / January, February 2008 *`
6	`c **************************************************************`	6	`c **************************************************************`
7		7
8	`implicit none`	8	`implicit none`
9		9
10	`c --------------------------------------------------------------`	10	`c --------------------------------------------------------------`
11	`c Declaration of parameters`	11	`c Declaration of parameters`
12	`c --------------------------------------------------------------`	12	`c --------------------------------------------------------------`
13		13
14	`c Maximum number of columns per trajectory`	14	`c Maximum number of columns per trajectory`
15	`integer maxcol`	15	`integer maxcol`
16	`parameter (maxcol=100)`	16	`parameter (maxcol=100)`
17		17
18	`c Maximum number of commands`	18	`c Maximum number of commands`
19	`integer maxcmd`	19	`integer maxcmd`
20	`parameter (maxcmd=10000)`	20	`parameter (maxcmd=10000)`
21		21
22	`c --------------------------------------------------------------`	22	`c --------------------------------------------------------------`
23	`c Declaration of variables`	23	`c Declaration of variables`
24	`c --------------------------------------------------------------`	24	`c --------------------------------------------------------------`
25		25
26	`c Input and output files`	26	`c Input and output files`
27	`character*120 inp_lslfile ! Input lsl file`	27	`character*120 inp_lslfile ! Input lsl file`
28	`character*120 out_lslfile ! Output lsl file`	28	`character*120 out_lslfile ! Output lsl file`
29	`character*120 inp_criteria ! Input file with criteria`	29	`character*120 inp_criteria ! Input file with criteria`
30	`integer inpmode ! Format of input file`	30	`integer inpmode ! Format of input file`
31	`integer outmode ! Format of output file`	31	`integer outmode ! Format of output file`
32	`character*80 outformat ! Trajectory/Boolean/Index of output`	32	`character*80 outformat ! Trajectory/Boolean/Index of output`
33	`character*80 regionf ! Name of the regionfile`	33	`character*80 regionf ! Name of the regionfile`
34		34
35	`c Trajectories`	35	`c Trajectories`
36	`integer ntim ! Number of times`	36	`integer ntim ! Number of times`
37	`integer ncol ! Number of columns (including time...)`	37	`integer ncol ! Number of columns (including time...)`
38	`integer ntrainp,ntraout ! Number of trajectories`	38	`integer ntrainp,ntraout ! Number of trajectories`
39	`character*80 vars(maxcol) ! Names of trajectory columns`	39	`character*80 vars(maxcol) ! Names of trajectory columns`
40	`integer basetime(6) ! Base time of trajectory (first line in lsl)`	40	`integer basetime(6) ! Base time of trajectory (first line in lsl)`
41	`real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories`	41	`real,allocatable, dimension (:,:,:) :: trainp ! Input trajectories`
42	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories`	42	`real,allocatable, dimension (:,:,:) :: traout ! Output trajectories`
43	`integer,allocatable,dimension (:) :: selflag ! Flag for selection`	43	`integer,allocatable,dimension (:) :: selflag ! Flag for selection`
44	`real,allocatable, dimension (:) :: time ! Times of the trajectory`	44	`real,allocatable, dimension (:) :: time ! Times of the trajectory`
45	`integer,allocatable,dimension (:,:) :: trigger ! Trigger column`	45	`integer,allocatable,dimension (:,:) :: trigger ! Trigger column`
46	`integer itrigger ! Column index for trigger`	46	`integer itrigger ! Column index for trigger`
47	`character*80 addtrigger ! Flag whether to add trigger column`	47	`character*80 addtrigger ! Flag whether to add trigger column`
48		48
49	`c Command stack`	49	`c Command stack`
50	`real cmd(maxcmd) ! Decoded slection criterion`	50	`real cmd(maxcmd) ! Decoded slection criterion`
51	`integer ncmd ! Number of commands`	51	`integer ncmd ! Number of commands`
52		52
53	`c Common block for initialisation of polygon check`	53	`c Common block for initialisation of polygon check`
54	`real tlonv(2000),vlat_c(2000),vlon_c(2000)`	54	`real tlonv(2000),vlat_c(2000),vlon_c(2000)`
55	`real xlat_c,xlon_c`	55	`real xlat_c,xlon_c`
56	`integer ibndry,nv_c`	56	`integer ibndry,nv_c`
57	`data ibndry/0/`	57	`data ibndry/0/`
58	`common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`	58	`common /spolybndry/vlat_c,vlon_c,nv_c,xlat_c,xlon_c,tlonv,ibndry`
59		59
60	`c Auxiliary variables`	60	`c Auxiliary variables`
61	`integer stat ! Logical (state) variable`	61	`integer stat ! Logical (state) variable`
62	`integer fid,fod,fcr ! File identifier for input and output`	62	`integer fid,fod,fcr ! File identifier for input and output`
63	`integer i,j,k ! Index counter`	63	`integer i,j,k ! Index counter`
64	`integer isok ! Flag for selected trajectory`	64	`integer isok ! Flag for selected trajectory`
65	`real param(1000) ! List of parameters`	65	`real param(1000) ! List of parameters`
66	`integer nparam ! Number of parameters`	66	`integer nparam ! Number of parameters`
67	`character*80 specialstr ! Name of special command`	67	`character*80 specialstr ! Name of special command`
68	`integer len`	68	`integer len`
69	`integer,allocatable,dimension (:) :: trigger1 ! Trigger column`	69	`integer,allocatable,dimension (:) :: trigger1 ! Trigger column`
70	`real,allocatable,dimension (:,:) :: trainp1 ! A single trajectory`	70	`real,allocatable,dimension (:,:) :: trainp1 ! A single trajectory`
71	`character ch`	71	`character ch`
72		72
73	`c --------------------------------------------------------------`	73	`c --------------------------------------------------------------`
74	`c Preparations`	74	`c Preparations`
75	`c --------------------------------------------------------------`	75	`c --------------------------------------------------------------`
76		76
77	`c Write start message`	77	`c Write start message`
78	`print*,'========================================================='`	78	`print*,'========================================================='`
79	`print,' START OF PROGRAM SELECT *'`	79	`print,' START OF PROGRAM SELECT *'`
80	`print*`	80	`print*`
81		81
82	`c Read parameter file`	82	`c Read parameter file`
83	`open(10,file='select.param')`	83	`open(10,file='select.param')`
84	`read(10,*) inp_lslfile`	84	`read(10,*) inp_lslfile`
85	`read(10,*) out_lslfile`	85	`read(10,*) out_lslfile`
86	`read(10,*) outformat`	86	`read(10,*) outformat`
87	`read(10,*) inp_criteria`	87	`read(10,*) inp_criteria`
88	`read(10,*) ntrainp`	88	`read(10,*) ntrainp`
89	`read(10,*) ntim`	89	`read(10,*) ntim`
90	`read(10,*) ncol`	90	`read(10,*) ncol`
91	`read(10,*) regionf`	91	`read(10,*) regionf`
92	`read(10,*) addtrigger`	92	`read(10,*) addtrigger`
93	`close(10)`	93	`close(10)`
94		94
95	`c Set the formats of the input and output files`	95	`c Set the formats of the input and output files`
96	`call mode_tra(inpmode,inp_lslfile)`	96	`call mode_tra(inpmode,inp_lslfile)`
97	`if (inpmode.eq.-1) inpmode=1`	97	`if (inpmode.eq.-1) inpmode=1`
98	`call mode_tra(outmode,out_lslfile)`	98	`call mode_tra(outmode,out_lslfile)`
99	`if ( (outmode.eq.-1).and.(outformat.ne.'startf') ) then`	99	`if ( (outmode.eq.-1).and.(outformat.ne.'startf') ) then`
100	`outmode=1`	100	`outmode=1`
101	`endif`	101	`endif`
102		102
103	`c Allocate memory for a single trajectory`	103	`c Allocate memory for a single trajectory`
104	`allocate(trainp(ntrainp,ntim,ncol),stat=stat)`	104	`allocate(trainp(ntrainp,ntim,ncol),stat=stat)`
105	`if (stat.ne.0) stop '* error allocating array trainp *'`	105	`if (stat.ne.0) stop '* error allocating array trainp *'`
106	`allocate(time(ntim),stat=stat)`	106	`allocate(time(ntim),stat=stat)`
107	`if (stat.ne.0) stop '* error allocating array time *'`	107	`if (stat.ne.0) stop '* error allocating array time *'`
108	`allocate(selflag(ntrainp),stat=stat)`	108	`allocate(selflag(ntrainp),stat=stat)`
109	`if (stat.ne.0) stop '* error allocating array selflag *'`	109	`if (stat.ne.0) stop '* error allocating array selflag *'`
110	`allocate(trigger(ntrainp,ntim),stat=stat)`	110	`allocate(trigger(ntrainp,ntim),stat=stat)`
111	`if (stat.ne.0) stop '* error allocating array trigger *'`	111	`if (stat.ne.0) stop '* error allocating array trigger *'`
112	`allocate(trigger1(ntim),stat=stat)`	112	`allocate(trigger1(ntim),stat=stat)`
113	`if (stat.ne.0) stop '* error allocating array trigger1 *'`	113	`if (stat.ne.0) stop '* error allocating array trigger1 *'`
114	`allocate(trainp1(ntim,ncol),stat=stat)`	114	`allocate(trainp1(ntim,ncol),stat=stat)`
115	`if (stat.ne.0) stop '* error allocating array trainp1 *'`	115	`if (stat.ne.0) stop '* error allocating array trainp1 *'`
116		116
117	`c Read the input trajectory file`	117	`c Read the input trajectory file`
118	`call ropen_tra(fid,inp_lslfile,ntrainp,ntim,ncol,`	118	`call ropen_tra(fid,inp_lslfile,ntrainp,ntim,ncol,`
119	`> basetime,vars,inpmode)`	119	`> basetime,vars,inpmode)`
120	`call read_tra (fid,trainp,ntrainp,ntim,ncol,inpmode)`	120	`call read_tra (fid,trainp,ntrainp,ntim,ncol,inpmode)`
121	`call close_tra(fid,inpmode)`	121	`call close_tra(fid,inpmode)`
122		122
123	`c Check that first four columns correspond to time,lon,lat,p`	123	`c Check that first four columns correspond to time,lon,lat,p`
124	`if ( (vars(1).ne.'time' ).or.`	124	`if ( (vars(1).ne.'time' ).or.`
125	`> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.`	125	`> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.`
126	`> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.`	126	`> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.`
127	`> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )`	127	`> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )`
128	`>then`	128	`>then`
129	`print*,' ERROR: problem with input trajectories ...'`	129	`print*,' ERROR: problem with input trajectories ...'`
130	`stop`	130	`stop`
131	`endif`	131	`endif`
132	`vars(1) = 'time'`	132	`vars(1) = 'time'`
133	`vars(2) = 'lon'`	133	`vars(2) = 'lon'`
134	`vars(3) = 'lat'`	134	`vars(3) = 'lat'`
135	`vars(4) = 'p'`	135	`vars(4) = 'p'`
136		136
137	`c Get the trajectory times from first trajectory`	137	`c Get the trajectory times from first trajectory`
138	`do i=1,ntim`	138	`do i=1,ntim`
139	`time(i)=trainp(1,i,1)`	139	`time(i)=trainp(1,i,1)`
140	`enddo`	140	`enddo`
141		141
142	`c Init the trigger field - first check whether it is already available`	142	`c Init the trigger field - first check whether it is already available`
143	`itrigger = 0`	143	`itrigger = 0`
144	`do i=1,ncol`	144	`do i=1,ncol`
145	`if ( vars(i).eq.'TRIGGER' ) itrigger = i`	145	`if ( vars(i).eq.'TRIGGER' ) itrigger = i`
146	`enddo`	146	`enddo`
147		147
148	`if ( itrigger.ne.0 ) then`	148	`if ( itrigger.ne.0 ) then`
149	`do i=1,ntrainp`	149	`do i=1,ntrainp`
150	`do j=1,ntim`	150	`do j=1,ntim`
151	`trigger(i,j) = nint( trainp(i,j,itrigger) )`	151	`trigger(i,j) = nint( trainp(i,j,itrigger) )`
152	`enddo`	152	`enddo`
153	`enddo`	153	`enddo`
154	`else`	154	`else`
155	`do i=1,ntrainp`	155	`do i=1,ntrainp`
156	`do j=1,ntim`	156	`do j=1,ntim`
157	`trigger(i,j) = 0`	157	`trigger(i,j) = 0`
158	`enddo`	158	`enddo`
159	`enddo`	159	`enddo`
160	`endif`	160	`endif`
161		161
162	`c Write some info about the trajectory`	162	`c Write some info about the trajectory`
163	`print*,'---- INPUT PARAMETERS -----------------------------------'`	163	`print*,'---- INPUT PARAMETERS -----------------------------------'`
164	`write(,)`	164	`write(,)`
165	`write(,) 'Input file : ',trim(inp_lslfile)`	165	`write(,) 'Input file : ',trim(inp_lslfile)`
166	`write(,) 'Output file : ',trim(out_lslfile)`	166	`write(,) 'Output file : ',trim(out_lslfile)`
167	`write(,) 'Output format : ',trim(outformat)`	167	`write(,) 'Output format : ',trim(outformat)`
168	`write(,) 'Criteria file : ',trim(inp_criteria)`	168	`write(,) 'Criteria file : ',trim(inp_criteria)`
169	`write(,) '# tra : ',ntrainp`	169	`write(,) '# tra : ',ntrainp`
170	`write(,) '# time : ',ntim`	170	`write(,) '# time : ',ntim`
171	`write(,) '# col : ',ncol`	171	`write(,) '# col : ',ncol`
172	`write(,) 'Region file : ',trim(regionf)`	172	`write(,) 'Region file : ',trim(regionf)`
173	`write(,) 'Add trigger : ',trim(addtrigger)`	173	`write(,) 'Add trigger : ',trim(addtrigger)`
174		174
175	`print*`	175	`print*`
176	`print*,'---- INPUT TRAJECTORY FILE ------------------------------'`	176	`print*,'---- INPUT TRAJECTORY FILE ------------------------------'`
177	`print*`	177	`print*`
178	`write(*,'(1x,a12,i4,a10)') 'Vars : ',1,trim(vars(1))`	178	`write(*,'(1x,a12,i4,a10)') 'Vars : ',1,trim(vars(1))`
179	`do i=2,ncol`	179	`do i=2,ncol`
180	`write(*,'(1x,a12,i4,a10)') ' ',i,trim(vars(i))`	180	`write(*,'(1x,a12,i4,a10)') ' ',i,trim(vars(i))`
181	`enddo`	181	`enddo`
182	`print*`	182	`print*`
183	`write(*,'(1x,a12,i4,f10.2)') 'Time : ',1,time(1)`	183	`write(*,'(1x,a12,i4,f10.2)') 'Time : ',1,time(1)`
184	`do i=2,ntim`	184	`do i=2,ntim`
185	`write(*,'(1x,a12,i4,f10.2)') ' ',i,time(i)`	185	`write(*,'(1x,a12,i4,f10.2)') ' ',i,time(i)`
186	`enddo`	186	`enddo`
187	`print*`	187	`print*`
188	`write(*,'(1x,a12,i4,i10)') 'Base date : ',1,basetime(1)`	188	`write(*,'(1x,a12,i4,i10)') 'Base date : ',1,basetime(1)`
189	`write(*,'(1x,a12,i4,i10)') ' ',2,basetime(2)`	189	`write(*,'(1x,a12,i4,i10)') ' ',2,basetime(2)`
190	`write(*,'(1x,a12,i4,i10)') ' ',3,basetime(3)`	190	`write(*,'(1x,a12,i4,i10)') ' ',3,basetime(3)`
191	`write(*,'(1x,a12,i4,i10)') ' ',4,basetime(4)`	191	`write(*,'(1x,a12,i4,i10)') ' ',4,basetime(4)`
192	`write(*,'(1x,a12,i4,i10)') ' ',5,basetime(5)`	192	`write(*,'(1x,a12,i4,i10)') ' ',5,basetime(5)`
193	`print*`	193	`print*`
194	`write(*,'(1x,a12,i4,i10)') 'Time range : ',6,basetime(6)`	194	`write(*,'(1x,a12,i4,i10)') 'Time range : ',6,basetime(6)`
195	`print*`	195	`print*`
196	`if ( itrigger.ne.0 ) then`	196	`if ( itrigger.ne.0 ) then`
197	`print*,'TRIGGER FIELD FOUND IN COLUMN ',itrigger`	197	`print*,'TRIGGER FIELD FOUND IN COLUMN ',itrigger`
198	`print*`	198	`print*`
199	`endif`	199	`endif`
200		200
201	`c Read and decode the selection criterion`	201	`c Read and decode the selection criterion`
202	`fcr = 10`	202	`fcr = 10`
203	`open(fcr,file=inp_criteria)`	203	`open(fcr,file=inp_criteria)`
204	`ncmd=maxcmd`	204	`ncmd=maxcmd`
205	`call decode(fcr,cmd,ncmd,vars,ncol,time,ntim,regionf)`	205	`call decode(fcr,cmd,ncmd,vars,ncol,time,ntim,regionf)`
206	`close(fcr)`	206	`close(fcr)`
207		207
208	`print*`	208	`print*`
209	`print*,'---- PSEUDO CODE FOR SELECTION --------------------------'`	209	`print*,'---- PSEUDO CODE FOR SELECTION --------------------------'`
210	`print*`	210	`print*`
211	`call dumpcode(cmd,ncmd,vars,ncol,time,ntim)`	211	`call dumpcode(cmd,ncmd,vars,ncol,time,ntim)`
212		212
213	`c --------------------------------------------------------------`	213	`c --------------------------------------------------------------`
214	`c Loop over all trajectories - selection`	214	`c Loop over all trajectories - selection`
215	`c --------------------------------------------------------------`	215	`c --------------------------------------------------------------`
216		216
217	`c Prepare string and parameters for SPECIAL commands`	217	`c Prepare string and parameters for SPECIAL commands`
218	`if ( cmd(1).eq.0 ) then`	218	`if ( cmd(1).eq.0 ) then`
219		219
220	`c Get command string`	220	`c Get command string`
221	`j = 2`	221	`j = 2`
222	`len = nint(cmd(j))`	222	`len = nint(cmd(j))`
223	`specialstr = ''`	223	`specialstr = ''`
224	`do k=1,len`	224	`do k=1,len`
225	`j = j + 1`	225	`j = j + 1`
226	`specialstr = trim(specialstr)//char(nint(cmd(j)))`	226	`specialstr = trim(specialstr)//char(nint(cmd(j)))`
227	`enddo`	227	`enddo`
228		228
229	`c Get paramters`	229	`c Get paramters`
230	`j = j + 1`	230	`j = j + 1`
231	`nparam = nint(cmd(j))`	231	`nparam = nint(cmd(j))`
232	`do k=1,nparam`	232	`do k=1,nparam`
233	`j = j + 1`	233	`j = j + 1`
234	`param(k) = cmd(j)`	234	`param(k) = cmd(j)`
235	`enddo`	235	`enddo`
236		236
237	`endif`	237	`endif`
238		238
239	`c Init the counter for selected trajectories`	239	`c Init the counter for selected trajectories`
240	`ntraout = 0`	240	`ntraout = 0`
241		241
242	`c Loop over all trajectories`	242	`c Loop over all trajectories`
243	`do i=1,ntrainp`	243	`do i=1,ntrainp`
244		244
245	`c Copy a single trajectory to <trainp1> and <trigger1>`	245	`c Copy a single trajectory to <trainp1> and <trigger1>`
246	`do j=1,ntim`	246	`do j=1,ntim`
247	`do k=1,ncol`	247	`do k=1,ncol`
248	`trainp1(j,k) = trainp(i,j,k)`	248	`trainp1(j,k) = trainp(i,j,k)`
249	`enddo`	249	`enddo`
250	`trigger1(j) = trigger(i,j)`	250	`trigger1(j) = trigger(i,j)`
251	`enddo`	251	`enddo`
252		252
253	`C Skip the trajectory if missing data are found for positions`	253	`C Skip the trajectory if missing data are found for positions`
254	`isok = 1`	254	`isok = 1`
255	`do j=1,ntim`	255	`do j=1,ntim`
256	`if ( trainp1(j,4).lt.0. ) isok = 0`	256	`if ( trainp1(j,4).lt.0. ) isok = 0`
257	`enddo`	257	`enddo`
258		258
259	`c Decide whether the trajectory is selected (handle SPECIAL commands)`	259	`c Decide whether the trajectory is selected (handle SPECIAL commands)`
260	`if ( isok.eq.1 ) then`	260	`if ( isok.eq.1 ) then`
261	`if (cmd(1).ne.0 ) then`	261	`if (cmd(1).ne.0 ) then`
262	`call select_tra (isok,cmd,ncmd,trainp1,trigger1,ntim,ncol)`	262	`call select_tra (isok,cmd,ncmd,trainp1,trigger1,ntim,ncol)`
263		263
264	`else`	264	`else`
265	`call special (isok,specialstr,trainp1,ntim,ncol,`	265	`call special (isok,specialstr,trainp1,ntim,ncol,`
266	`> vars,time,param,nparam)`	266	`> vars,time,param,nparam)`
267	`endif`	267	`endif`
268	`endif`	268	`endif`
269		269
270	`c The trigger might be changed in the selection - copy it`	270	`c The trigger might be changed in the selection - copy it`
271	`do j=1,ntim`	271	`do j=1,ntim`
272	`trigger(i,j) = trigger1(j)`	272	`trigger(i,j) = trigger1(j)`
273	`enddo`	273	`enddo`
274		274
275	`c Set flag for selected trajectories`	275	`c Set flag for selected trajectories`
276	`if (isok.eq.1) then`	276	`if (isok.eq.1) then`
277	`selflag(i) = 1`	277	`selflag(i) = 1`
278	`ntraout = ntraout + 1`	278	`ntraout = ntraout + 1`
279	`else`	279	`else`
280	`selflag(i) = 0`	280	`selflag(i) = 0`
281	`endif`	281	`endif`
282		282
283	`enddo`	283	`enddo`
284		284
285	`c --------------------------------------------------------------`	285	`c --------------------------------------------------------------`
286	`c Write output trajectories`	286	`c Write output trajectories`
287	`c --------------------------------------------------------------`	287	`c --------------------------------------------------------------`
288		288
289	`c ------ Write output trajectories -----------------------------`	289	`c ------ Write output trajectories -----------------------------`
290	`if ( outformat.eq.'trajectory' ) then`	290	`if ( outformat.eq.'trajectory' ) then`
291		291
292		292
293	`c Define the trigger field if it is not yet defined`	293	`c Define the trigger field if it is not yet defined`
294	`if ( ( addtrigger.eq.'-trigger' ).and.(itrigger.eq.0) ) then`	294	`if ( ( addtrigger.eq.'-trigger' ).and.(itrigger.eq.0) ) then`
295	`ncol = ncol + 1`	295	`ncol = ncol + 1`
296	`vars(ncol) = 'TRIGGER'`	296	`vars(ncol) = 'TRIGGER'`
297	`itrigger = ncol`	297	`itrigger = ncol`
298	`endif`	298	`endif`
299		299
300	`c Allocate memory for output trajectory`	300	`c Allocate memory for output trajectory`
301	`allocate(traout(ntraout,ntim,ncol),stat=stat)`	301	`allocate(traout(ntraout,ntim,ncol),stat=stat)`
302	`if (stat.ne.0) stop '* error allocating array apply *'`	302	`if (stat.ne.0) stop '* error allocating array apply *'`
303		303
304	`c Set output trajectories`	304	`c Set output trajectories`
305	`j = 0`	305	`j = 0`
306	`do i=1,ntrainp`	306	`do i=1,ntrainp`
307	`if (selflag(i).eq.1) then`	307	`if (selflag(i).eq.1) then`
308	`j = j + 1`	308	`j = j + 1`
309	`traout(j,1:ntim,1:ncol) = trainp(i,1:ntim,1:ncol)`	309	`traout(j,1:ntim,1:ncol) = trainp(i,1:ntim,1:ncol)`
310	`if ( itrigger.ne.0 ) then`	310	`if ( itrigger.ne.0 ) then`
311	`traout(j,1:ntim,itrigger) = real(trigger(i,1:ntim))`	311	`traout(j,1:ntim,itrigger) = real(trigger(i,1:ntim))`
312	`endif`	312	`endif`
313	`endif`	313	`endif`
314	`enddo`	314	`enddo`
315		315
316	`c Write trajectories`	316	`c Write trajectories`
317	`call wopen_tra(fod,out_lslfile,ntraout,ntim,ncol,`	317	`call wopen_tra(fod,out_lslfile,ntraout,ntim,ncol,`
318	`> basetime,vars,outmode)`	318	`> basetime,vars,outmode)`
319	`call write_tra(fod,traout,ntraout,ntim,ncol,outmode)`	319	`call write_tra(fod,traout,ntraout,ntim,ncol,outmode)`
320	`call close_tra(fod,outmode)`	320	`call close_tra(fod,outmode)`
321		321
322	`c ------ Write index list -------------------------------------`	322	`c ------ Write index list -------------------------------------`
323	`elseif ( outformat.eq.'index' ) then`	323	`elseif ( outformat.eq.'index' ) then`
324		324
325	`open(10,file=out_lslfile)`	325	`open(10,file=out_lslfile)`
326	`do i=1,ntrainp`	326	`do i=1,ntrainp`
327	`if ( selflag(i).eq.1) write(10,*) i`	327	`if ( selflag(i).eq.1) write(10,*) i`
328	`enddo`	328	`enddo`
329	`close(10)`	329	`close(10)`
330		330
331	`c ------ Write boolean list -----------------------------------`	331	`c ------ Write boolean list -----------------------------------`
332	`elseif ( outformat.eq.'boolean' ) then`	332	`elseif ( outformat.eq.'boolean' ) then`
333		333
334	`open(10,file=out_lslfile)`	334	`open(10,file=out_lslfile)`
335	`do i=1,ntrainp`	335	`do i=1,ntrainp`
336	`write(10,'(i1)') selflag(i)`	336	`write(10,'(i1)') selflag(i)`
337	`enddo`	337	`enddo`
338	`close(10)`	338	`close(10)`
339		339
340	`c ------ Write count -------------------------------------------`	340	`c ------ Write count -------------------------------------------`
341	`elseif ( outformat.eq.'count' ) then`	341	`elseif ( outformat.eq.'count' ) then`
342		342
343	`open(10,file=out_lslfile)`	343	`open(10,file=out_lslfile)`
344	`write(10,'(i7)') ntraout`	344	`write(10,'(i7)') ntraout`
345	`close(10)`	345	`close(10)`
346		346
347	`c ------ Write starting positions -----------------------------`	347	`c ------ Write starting positions -----------------------------`
348	`elseif ( outformat.eq.'startf' ) then`	348	`elseif ( outformat.eq.'startf' ) then`
349		349
350	`c Allocate memory for output trajectory`	350	`c Allocate memory for output trajectory`
351	`allocate(traout(ntraout,1,ncol),stat=stat)`	351	`allocate(traout(ntraout,1,ncol),stat=stat)`
352	`if (stat.ne.0) stop '* error allocating array apply *'`	352	`if (stat.ne.0) stop '* error allocating array apply *'`
353		353
354	`c Set output trajectories`	354	`c Set output trajectories`
355	`j = 0`	355	`j = 0`
356	`do i=1,ntrainp`	356	`do i=1,ntrainp`
357	`if (selflag(i).eq.1) then`	357	`if (selflag(i).eq.1) then`
358	`j = j + 1`	358	`j = j + 1`
359	`traout(j,1,:) = trainp(i,1,:)`	359	`traout(j,1,:) = trainp(i,1,:)`
360	`endif`	360	`endif`
361	`enddo`	361	`enddo`
362		362
363	`c Write trajectories`	363	`c Write trajectories`
364	`if (outmode.ne.-1) then`	364	`if (outmode.ne.-1) then`
365		365
366	`call wopen_tra(fod,out_lslfile,ntraout,1,ncol,`	366	`call wopen_tra(fod,out_lslfile,ntraout,1,ncol,`
367	`> basetime,vars,outmode)`	367	`> basetime,vars,outmode)`
368	`call write_tra(fod,traout,ntraout,1,ncol,outmode)`	368	`call write_tra(fod,traout,ntraout,1,ncol,outmode)`
369	`call close_tra(fod,outmode)`	369	`call close_tra(fod,outmode)`
370		370
371	`c Output as a triple list (corresponding to <startf> file)`	371	`c Output as a triple list (corresponding to <startf> file)`
372	`else`	372	`else`
373		373
374	`fid = 10`	374	`fid = 10`
375	`open(fid,file=out_lslfile)`	375	`open(fid,file=out_lslfile)`
376	`do i=1,ntraout`	376	`do i=1,ntraout`
377	`write(fid,'(3f10.3)') traout(i,1,2), ! longitude`	377	`write(fid,'(3f10.3)') traout(i,1,2), ! longitude`
378	`> traout(i,1,3), ! latitude`	378	`> traout(i,1,3), ! latitude`
379	`> traout(i,1,4) ! pressure`	379	`> traout(i,1,4) ! pressure`
380	`enddo`	380	`enddo`
381	`close(fid)`	381	`close(fid)`
382		382
383	`endif`	383	`endif`
384		384
385	`endif`	385	`endif`
386		386
387	`c Write some status information, and end of program message`	387	`c Write some status information, and end of program message`
388	`print*`	388	`print*`
389	`print*,'---- STATUS INFORMATION --------------------------------'`	389	`print*,'---- STATUS INFORMATION --------------------------------'`
390	`print*`	390	`print*`
391	`print*,' # input trajectories : ',ntrainp`	391	`print*,' # input trajectories : ',ntrainp`
392	`print*,' # output trajectories : ',ntraout`	392	`print*,' # output trajectories : ',ntraout`
393	`print*`	393	`print*`
394	`print,' END OF PROGRAM SELECT *'`	394	`print,' END OF PROGRAM SELECT *'`
395	`print*,'========================================================='`	395	`print*,'========================================================='`
396		396
397	`stop`	397	`stop`
398		398
399	`c Exception handling`	399	`c Exception handling`
400	`100 stop 'select: First column in input trajectory must be <time>'`	400	`100 stop 'select: First column in input trajectory must be <time>'`
401	`101 stop 'select: Input trajectory file is empty'`	401	`101 stop 'select: Input trajectory file is empty'`
402		402
403	`end`	403	`end`
404		404
405	`c --------------------------------------------------------------`	405	`c --------------------------------------------------------------`
406	`c Dump the command list`	406	`c Dump the command list`
407	`c --------------------------------------------------------------`	407	`c --------------------------------------------------------------`
408		408
409	`subroutine dumpcode(out,n,vars,nvars,times,ntimes)`	409	`subroutine dumpcode(out,n,vars,nvars,times,ntimes)`
410		410
411	`c Write the command list to screen. The command list is decoded`	411	`c Write the command list to screen. The command list is decoded`
412	`c by call to <decode>`	412	`c by call to <decode>`
413		413
414	`implicit none`	414	`implicit none`
415		415
416	`c Declaration of subroutine parameters`	416	`c Declaration of subroutine parameters`
417	`integer n`	417	`integer n`
418	`real out(n)`	418	`real out(n)`
419	`integer nvars`	419	`integer nvars`
420	`character*80 vars(nvars)`	420	`character*80 vars(nvars)`
421	`integer ntimes`	421	`integer ntimes`
422	`real times(ntimes)`	422	`real times(ntimes)`
423		423
424	`c A single command`	424	`c A single command`
425	`character*80 cmd`	425	`character*80 cmd`
426	`character*80 var,mode,strtim`	426	`character*80 var,mode,strtim`
427	`integer nval`	427	`integer nval`
428	`integer ntim`	428	`integer ntim`
429	`integer ivar,imode,icmd,itime`	429	`integer ivar,imode,icmd,itime`
430		430
431	`c Auxiliary variables`	431	`c Auxiliary variables`
432	`integer i,j`	432	`integer i,j`
433		433
434	`c Loop through the complete list`	434	`c Loop through the complete list`
435	`i=0`	435	`i=0`
436	`100 if (i.lt.n) then`	436	`100 if (i.lt.n) then`
437		437
438	`write(,) '---------------------------------------'`	438	`write(,) '---------------------------------------'`
439		439
440	`c Get command`	440	`c Get command`
441	`i=i+1`	441	`i=i+1`
442	`icmd=nint(out(i))`	442	`icmd=nint(out(i))`
443		443
444	`c Special handling of SPECIAL commands`	444	`c Special handling of SPECIAL commands`
445	`if ( icmd.eq.0 ) then`	445	`if ( icmd.eq.0 ) then`
446		446
447	`c Write 'header' for SPECIAL command`	447	`c Write 'header' for SPECIAL command`
448	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'SPECIAL'`	448	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'SPECIAL'`
449		449
450	`c Write command string`	450	`c Write command string`
451	`i = i + 1`	451	`i = i + 1`
452	`icmd = nint(out(i))`	452	`icmd = nint(out(i))`
453	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'LEN(CMD)'`	453	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'LEN(CMD)'`
454	`do j=1,icmd`	454	`do j=1,icmd`
455	`i = i + 1`	455	`i = i + 1`
456	`ivar = nint(out(i))`	456	`ivar = nint(out(i))`
457	`write(*,'(i5,f15.4,10x,a10)') i,out(i),char(ivar)`	457	`write(*,'(i5,f15.4,10x,a10)') i,out(i),char(ivar)`
458	`enddo`	458	`enddo`
459		459
460	`c Write parameters`	460	`c Write parameters`
461	`i = i + 1`	461	`i = i + 1`
462	`nval = nint(out(i))`	462	`nval = nint(out(i))`
463	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'`	463	`write(*,'(i5,f15.4,10x,a10)') i,out(i),'#PARAMETER'`
464	`do j=1,nval`	464	`do j=1,nval`
465	`i=i+1`	465	`i=i+1`
466	`if ( var.ne.'INPOLYGON' ) then`	466	`if ( var.ne.'INPOLYGON' ) then`
467	`write(*,'(i5,f15.4)') i,out(i)`	467	`write(*,'(i5,f15.4)') i,out(i)`
468	`else`	468	`else`
469	`write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))`	469	`write(*,'(i5,f15.4,a2)') i,out(i),char(nint(out(i)))`
470	`endif`	470	`endif`
471	`enddo`	471	`enddo`
472		472
473	`c Nothing else to do - exit`	473	`c Nothing else to do - exit`
474	`goto 200`	474	`goto 200`
475		475
476	`endif`	476	`endif`
477		477
478	`c Set the command`	478	`c Set the command`
479	`if (icmd.eq. 1) cmd='GT'`	479	`if (icmd.eq. 1) cmd='GT'`
480	`if (icmd.eq. 2) cmd='LT'`	480	`if (icmd.eq. 2) cmd='LT'`
481	`if (icmd.eq. 3) cmd='IN'`	481	`if (icmd.eq. 3) cmd='IN'`
482	`if (icmd.eq. 4) cmd='OUT'`	482	`if (icmd.eq. 4) cmd='OUT'`
483	`if (icmd.eq. 5) cmd='EQ'`	483	`if (icmd.eq. 5) cmd='EQ'`
484	`if (icmd.eq. 6) cmd='TRUE'`	484	`if (icmd.eq. 6) cmd='TRUE'`
485	`if (icmd.eq. 7) cmd='FALSE'`	485	`if (icmd.eq. 7) cmd='FALSE'`
486	`if (icmd.eq. 8) cmd='ALL'`	486	`if (icmd.eq. 8) cmd='ALL'`
487	`if (icmd.eq. 9) cmd='ANY'`	487	`if (icmd.eq. 9) cmd='ANY'`
488	`if (icmd.eq. 10) cmd='NONE'`	488	`if (icmd.eq. 10) cmd='NONE'`
489	`if (icmd.eq. -1) cmd='BEGIN'`	489	`if (icmd.eq. -1) cmd='BEGIN'`
490	`if (icmd.eq. -2) cmd='END'`	490	`if (icmd.eq. -2) cmd='END'`
491	`if (icmd.eq. -3) cmd='AND'`	491	`if (icmd.eq. -3) cmd='AND'`
492	`if (icmd.eq. -4) cmd='OR'`	492	`if (icmd.eq. -4) cmd='OR'`
493	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(cmd)`	493	`write(*,'(i5,f15.4,10x,a10)') i,out(i),trim(cmd)`
494	`if (icmd.lt.0) goto 100`	494	`if (icmd.lt.0) goto 100`
495		495
496	`c Get variable`	496	`c Get variable`
497	`i=i+1`	497	`i=i+1`
498	`ivar=nint(out(i))`	498	`ivar=nint(out(i))`
499	`if ( ivar.eq. -1 ) then`	499	`if ( ivar.eq. -1 ) then`
500	`var = 'DIST'`	500	`var = 'DIST'`

Subversion Repositories lagranto.ecmwf

(root)/trunk/select/select.f – Rev 3 → 9