WebSVN – lagranto.ecmwf – Diff – /trunk/goodies/extract.f


      PROGRAM extract
      
c     ***********************************************************************
c     * Extract columns of a trajectory file                                *
c     * Michael Sprenger / Spring, summer 2010                              *
c     ***********************************************************************
 
      implicit none
 
c     ----------------------------------------------------------------------
c     Declaration of parameters
c     ----------------------------------------------------------------------
 
c     Numerical epsiloN
      real      eps
      parameter (eps=0.0001)
 
c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------
 
c     Extraction mode
      character*80                           mode                             
 
c     Input and output format for trajectories (see iotra.f)
      character*80                           inpfile       ! Input filename
      character*80                           outfile       ! Output filename
 
c     Trajectories
      integer                                ntra_inp      ! Number of trajectories
      integer                                ntim_inp      ! Number of times
      integer                                ncol_inp      ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra_inp       ! Trajectories (ntra,ntim,ncol)
      character*80                           vars_inp(100) ! Variable names
 
      integer                                ntra_out      ! Number of trajectories
      integer                                ntim_out      ! Number of times
      integer                                ncol_out      ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra_out       ! Trajectories (ntra,ntim,ncol)
      integer,allocatable, dimension (:)  :: ind           ! Index for selection
      integer,allocatable, dimension (:)  :: isok          ! Index for selection
      character*80                           vars_out(100) ! Variable names
  
      real                                   time_inp(10000) ! Times of input trajectory
      real                                   time_out(10000) ! Times of output trajectory
      integer                                refdate(6)    ! Reference date
      integer                                ind_time(10000) ! Index for time selection
 
c     Auxiliary variables
      integer                                inpmode
      integer                                outmode
      integer                                stat
      integer                                fid
      integer                                i,j,k,n,j0,j1
      character*80                           str
      character*80                           split_str(100)
      integer                                split_n
      integer                                isstr,nvars,ileft,iright
      character*80                           vars(100)
      character                              ch
      real                                   tmp0,tmp1
      integer                                ind1
	  character*2000                         linestr
	  integer                                istr(100)
      integer                                nstr
      character*80                           strsplit(100)
      integer                                flag
 
c     ----------------------------------------------------------------------
c     Read and handle parameters
c     ----------------------------------------------------------------------
 
c     Read parameters
      open(10,file='extract.param')
       read(10,*) inpfile
       read(10,*) outfile
       read(10,*) mode
       read(10,*) ntra_inp,ntim_inp,ncol_inp
       read(10,*) str
      close(10)
 
c     Split the input string
      isstr   = 0
      split_n = 0
      do i=1,80
         ch = str(i:i)
         if ( (isstr.eq.0).and.(ch.ne.' ') ) then
            isstr=1
            ileft=i
         elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
            isstr              = 0
            iright             = i-1
            split_n            = split_n+1
            split_str(split_n) = str(ileft:iright)
         endif
      enddo
 
c     Determine the formats
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1
      call mode_tra(outmode,outfile)
      if ( (mode.ne.'-startf').and.(outmode.eq.-1) ) then
         outmode=1
      endif
 
c     ----------------------------------------------------------------------
c     Read input trajectories
c     ----------------------------------------------------------------------
 
c     Allocate memory
      allocate(tra_inp(ntra_inp,ntim_inp,ncol_inp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_inp    ***' 
      allocate(ind(ntra_inp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ind        ***' 
      allocate(isok(ntra_inp),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array isok       ***' 
 
c     Read inpufile
      fid = 10
      call ropen_tra(fid,inpfile,ntra_inp,ntim_inp,ncol_inp,
     >               refdate,vars_inp,inpmode)
      call read_tra (fid,tra_inp,ntra_inp,ntim_inp,ncol_inp,inpmode)
      call close_tra(fid,inpmode)
 
c     Check format 
      if (  vars_inp(1).ne.'time') goto 990
      if ( (vars_inp(2).ne.'lon').and.(vars_inp(2).ne.'xpos') ) goto 990
      if ( (vars_inp(3).ne.'lat').and.(vars_inp(3).ne.'ypos') ) goto 990
      if ( (vars_inp(4).ne.'p'  ).and.(vars_inp(4).ne.'ppos') ) goto 990
 
c     ----------------------------------------------------------------------
c     Option -vars : Extract columns of variables
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-var' ) goto 100
 
c     Set the first for columns of the output
      ncol_out   = 4
      vars_out(1)='time'
      vars_out(2)='lon'
      vars_out(3)='lat'
      vars_out(4)='p'
      ind(1)     =1
      ind(2)     =2
      ind(3)     =3
      ind(4)     =4
 
c     Get final list of extraction columns (set number of columns)
      do i=1,split_n
 
         if (split_str(i).eq.'to') then
            do j=1,ncol_inp
               if ( vars_inp(j).eq.split_str(i-1) ) j0 = j + 1
            enddo
            do j=1,ncol_inp
               if ( vars_inp(j).eq.split_str(i+1) ) j1 = j - 1
            enddo
            do j=j0,j1
               ncol_out           = ncol_out + 1
               vars_out(ncol_out) = vars_inp(j)
            enddo
         else
            ncol_out           = ncol_out + 1
            vars_out(ncol_out) = split_str(i)
         endif
 
      enddo
 
c     Set the dimensions of the output trajectory
      ntra_out = ntra_inp
      ntim_out = ntim_inp
 
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Extract <time,lon,lat,p> columns
      do i=1,ntra_out
         do j=1,ntim_out
            tra_out(i,j,1) = tra_inp(i,j,1)
            tra_out(i,j,2) = tra_inp(i,j,2)
            tra_out(i,j,3) = tra_inp(i,j,3)
            tra_out(i,j,4) = tra_inp(i,j,4)
         enddo
      enddo
 
c     Get indices for new columns (1..4 are already ok)
      do i=5,ncol_out
         ind(i)=0
      enddo
      do i=5,ncol_out
         do j=1,ncol_inp
            if ( vars_inp(j).eq.vars_out(i) ) ind(i) = j
         enddo
      enddo
 
c     Check if all selected columns are available
      do i=1,ncol_out
         if ( ind(i).eq.0 ) then
            print*,'Invalid column in ',trim(str)
            stop
         endif
      enddo
 
c     Extract the column
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               tra_out(i,j,k) = tra_inp(i,j,ind(k))
            enddo
         enddo
      enddo
 
 100  continue
 
c     ----------------------------------------------------------------------
c     Option -times : Extract times of trajectories
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-time' ) goto 110
 
c     Set the dimension of the output trajectory
      ntim_out = 0
 
c     Get the list of times for the input trajectory
      do i=1,ntim_inp
         time_inp(i) = tra_inp(1,i,1)
      enddo
 
c     Get final list of extraction times (set number of times)
      do i=1,split_n
 
         if (split_str(i).eq.'to') then
            read(split_str(i-1),*) tmp0
            do j=1,ntim_inp
               if ( time_inp(j).eq.tmp0 ) j0 = j + 1
            enddo
            read(split_str(i+1),*) tmp0
            do j=1,ntim_inp
               if ( time_inp(j).eq.tmp0 ) j1 = j - 1
            enddo
            do j=j0,j1
               ntim_out           = ntim_out + 1
               time_out(ntim_out) = time_inp(j)
            enddo
         elseif (split_str(i).eq.'first') then
            ntim_out           = ntim_out + 1
            time_out(ntim_out) = time_inp(1)
         elseif (split_str(i).eq.'last') then
            ntim_out           = ntim_out + 1
            time_out(ntim_out) = time_inp(ntim_inp)
         else
            ntim_out           = ntim_out + 1
            read(split_str(i),*) tmp0
            time_out(ntim_out) = tmp0
         endif
 
      enddo
 
c     Get the indices of the selected times
      do i=1,ntim_out
         ind_time(i) = 0
      enddo
      do i=1,ntim_out
         do j=1,ntim_inp
            if ( abs(time_out(i)-time_inp(j)).lt.eps) ind_time(i) = j
         enddo
      enddo
      do i=1,ntim_out
         if ( ind_time(i).eq.0) then
            print*,' Invalid time ',time_out(i)
            stop
         endif
      enddo
      
c     Set dimensions of output trajectory
      ntra_out = ntra_inp
      ncol_out = ncol_inp
 
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected times to the output trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               tra_out(i,j,k) = tra_inp(i,ind_time(j),k)
            enddo
         enddo
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 110  continue
 
c     ----------------------------------------------------------------------
c     Option -tra : Extract trajectories by number
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-tra' ) goto 120
 
c     Set the dimension of the output trajectory
      ntra_out = 0
 
c     Get final list of extraction times (set number of times)
      do i=1,split_n
 
         if (split_str(i).eq.'to') then
            read(split_str(i-1),*) tmp0
            read(split_str(i+1),*) tmp1
            do j=nint(tmp0)+1,nint(tmp1)-1
               ntra_out        = ntra_out + 1
               ind(ntra_out)   = j   
            enddo
         elseif (split_str(i).eq.'first') then
            ntra_out           = ntra_out + 1
            ind(ntra_out)      = 1
         elseif (split_str(i).eq.'last') then
            ntra_out           = ntra_out + 1
            ind(ntra_out)      = ntra_inp
         else
            ntra_out           = ntra_out + 1
            read(split_str(i),*) tmp0
            ind(ntra_out)      = nint(tmp0) 
         endif
 
      enddo
 
c     Check whether selected trajectories are ok
      do i=1,ntra_out
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
            print*,'Invalid trajectory selected ',ind(i)
            stop
         endif
      enddo
 
c     Set dimensions of output trajectory
      ntim_out = ntim_inp
      ncol_out = ncol_inp
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected times to the output trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
            enddo
         enddo
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 120  continue
 
c     ----------------------------------------------------------------------
c     Option -startf : Extract starting positions for the trajectory file
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-startf' ) goto 130
 
c     Set the first for columns of the output
      ncol_out   = 4
      vars_out(1)='time'
      vars_out(2)='lon'
      vars_out(3)='lat'
      vars_out(4)='p'
      ind(1)     =1
      ind(2)     =2
      ind(3)     =3
      ind(4)     =4
 
c     Set dimensions of output trajectory
      ntim_out = 1
      ntra_out = ntra_inp
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected times to the output trajectory
      do i=1,ntra_out
         do k=1,ncol_out
            tra_out(i,1,k) = tra_inp(i,1,k)
         enddo
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 130  continue
 
c     ----------------------------------------------------------------------
c     Option -index : Extract trajectories by index file
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-index' ) goto 140
 
c     Read the index file
      open(10,file=str)
 
        ntra_out = 1
 142    read(10,*,end=141) ind(ntra_out)
        ntra_out = ntra_out + 1
        goto 142 
 141    continue
        ntra_out = ntra_out - 1
 
      close(10)
 
c     Check whether selected trajectories are ok
      do i=1,ntra_out
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
            print*,'Invalid trajectory selected ',ind(i)
            stop
         endif
      enddo
 
c     Set dimensions of output trajectory
      ntim_out = ntim_inp
      ncol_out = ncol_inp
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected times to the output trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
            enddo
         enddo
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 140  continue
 
c     ----------------------------------------------------------------------
c     Option -boolean : Extract trajectories by boolean file
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-boolean' ) goto 150
 
c     Read the index file
      open(10,file=str)
        ntra_out = 0
        do i=1,ntra_inp
           read(10,*) ind1
           if ( ind1.eq.1 ) then
              ntra_out      = ntra_out + 1
              ind(ntra_out) = i
           endif
        enddo
      close(10)
 
c     Check whether selected trajectories are ok
      do i=1,ntra_out
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
            print*,'Invalid trajectory selected ',ind(i)
            stop
         endif
      enddo
 
c     Set dimensions of output trajectory
      ntim_out = ntim_inp
      ncol_out = ncol_inp
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected times to the output trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
            enddo
         enddo
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 150  continue
 
c     ----------------------------------------------------------------------
c     Option -pattern : Extract trajectories which match a regular expression
c     ----------------------------------------------------------------------
 
	  if ( mode.ne.'-pattern' ) goto 160
 
c	  All times and columns are extracted
      ncol_out = ncol_inp
      ntim_out = ntim_inp
      ntra_out = 0
 
c     Split the search string
	  nstr   = 0
	  ileft  = 0
	  iright = 0
      do i=1,len_trim(str)
      	if ( (str(i:i).eq.' ').and.(ileft.eq.0) ) then
      	   ileft = ileft + 1
        elseif ( (str(i:i).ne.' ').and.(ileft.eq.0) ) then
           ileft  = i
           iright = 0
        elseif ( (str(i:i).ne.' ').and.(ileft.ne.0) ) then
           iright = i
        elseif ( (str(i:i).eq.' ').and.(ileft.ne.0) ) then
           nstr           = nstr + 1
           strsplit(nstr) = str(ileft:iright)
           ileft          = 0
           iright         = 0
        endif
      enddo
      if ( (ileft.ne.0).and.(iright.ne.0) ) then
           nstr           = nstr + 1
           strsplit(nstr) = str(ileft:iright)
           ileft          = 0
      endif
 
c     Loop over the trajectories - check for matching pattern
	  do n=1,ntra_inp
 
	    ind(n) = 0
	  	do i=1,ntim_inp
 
           write(linestr,'(1f7.2,f9.2,f8.2,i6,100f10.3)')
     >                   (tra_inp(n,i,j),j=1,3),             ! time, lon, lat
     >                   nint(tra_inp(n,i,4)),               ! p
     >                   (tra_inp(n,i,j),j=5,ncol_inp)       ! fields
 
           flag = 1
           do k=1,nstr
             istr(k) = index(trim(linestr),trim(strsplit(k)))
             if ( istr(k).eq.0 ) flag = 0
           enddo
	       if ( flag.eq.1 ) ind(n) = 1
 
	  	enddo
	  	if ( ind(n).eq.1 ) ntra_out = ntra_out + 1
 
	  enddo
 
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***'
 
c     Copy the selected times to the output trajectory
      ntra_out = 0
      do i=1,ntra_inp
      	if (ind(i).eq.1) then
      	    ntra_out = ntra_out + 1
         	do j=1,ntim_out
            	do k=1,ncol_out
               		tra_out(ntra_out,j,k) = tra_inp(i,j,k)
            	enddo
         	enddo
         endif
      enddo
 
c     Copy meta information
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
 160  continue
 
c     ----------------------------------------------------------------------
c     Option -leaving : Extract all trajectories which leave domain
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-leaving' ) goto 170
 
c     Set dimensions of output trajectory
      ntim_out = ntim_inp
      ncol_out = ncol_inp
      ntra_out = 0
 
c     Copy the meta data
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
c     Determine the number of trajectories leaving domain
      do i=1,ntra_inp
         isok(i) = 1
         do j=1,ntim_inp
            if ( tra_inp(i,j,4).lt.0. ) isok(i) = 0
         enddo         
         if ( isok(i).eq.0 ) then
            ntra_out = ntra_out + 1
         endif
      enddo
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected trajectories to the output trajectory
      ntra_out = 0
      do i=1,ntra_inp
         if ( isok(i).eq.0 ) then
            ntra_out = ntra_out + 1
            do j=1,ntim_inp
               do k=1,ncol_out
                  tra_out(ntra_out,j,k) = tra_inp(i,j,k)
               enddo
            enddo
         endif
      enddo
         
c     Copy meta information
 
 170  continue
 
c     ----------------------------------------------------------------------
c     Option -staying : Extract all trajectories which stay in domain
c     ----------------------------------------------------------------------
 
      if ( mode.ne.'-staying' ) goto 180
 
c     Set dimensions of output trajectory
      ntim_out = ntim_inp
      ncol_out = ncol_inp
      ntra_out = 0
 
c     Copy the meta data
      do i=1,ncol_out
         vars_out(i) = vars_inp(i)
      enddo
 
c     Determine the number of trajectories staying in domain
      do i=1,ntra_inp
         isok(i) = 1
         do j=1,ntim_inp
            if ( tra_inp(i,j,4).lt.0. ) isok(i) = 0
         enddo         
         if ( isok(i).eq.1 ) then
            ntra_out = ntra_out + 1
         endif
      enddo
     
c     Allocate memory for output trajectory
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Copy the selected trajectories to the output trajectory
      ntra_out = 0
      do i=1,ntra_inp
         if ( isok(i).eq.1 ) then
            ntra_out = ntra_out + 1
            do j=1,ntim_inp
               do k=1,ncol_out
                  tra_out(ntra_out,j,k) = tra_inp(i,j,k)
               enddo
            enddo
         endif
      enddo
         
c     Copy meta information
 
 180  continue
 
c     ----------------------------------------------------------------------
c     Write output trajectories
c     ----------------------------------------------------------------------
 
c     Write output as trajectory file
      if (outmode.ge.1) then
 
         call wopen_tra(fid,outfile,ntra_out,ntim_out,ncol_out,
     >                  refdate,vars_out,outmode)
         call write_tra(fid,tra_out,ntra_out,ntim_out,ncol_out,outmode)
         call close_tra(fid,outmode)
 
c     Write output as (lon, lat, p)-list
      else
 
         open(10,file=outfile)
         do i=1,ntra_out
            write(10,'(3f10.2)') tra_out(i,1,2),    ! lon
     >                           tra_out(i,1,3),    ! lat
     >                           tra_out(i,1,4)     ! p
         enddo
         close(10)
 
      endif
 
 
 
!c     ----------------------------------------------------------------------
c     Error handling
c     ----------------------------------------------------------------------
      
      stop
   
 990  print*,'First columns must be <time,lon,lat,p>... Stop'
      stop
      
 
      end
 
      
 
      

Rev 5	Rev 27
1	`PROGRAM extract`	1	`PROGRAM extract`
2		2
3	`c ***********************************************************************`	3	`c ***********************************************************************`
4	`c * Extract columns of a trajectory file *`	4	`c * Extract columns of a trajectory file *`
5	`c * Michael Sprenger / Spring, summer 2010 *`	5	`c * Michael Sprenger / Spring, summer 2010 *`
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 Numerical epsilon`	14	`c Numerical epsiloN`
15	`real eps`	15	`real eps`
16	`parameter (eps=0.0001)`	16	`parameter (eps=0.0001)`
17		17
18	`c ----------------------------------------------------------------------`	18	`c ----------------------------------------------------------------------`
19	`c Declaration of variables`	19	`c Declaration of variables`
20	`c ----------------------------------------------------------------------`	20	`c ----------------------------------------------------------------------`
21		21
22	`c Extraction mode`	22	`c Extraction mode`
23	`character*80 mode`	23	`character*80 mode`
24		24
25	`c Input and output format for trajectories (see iotra.f)`	25	`c Input and output format for trajectories (see iotra.f)`
26	`character*80 inpfile ! Input filename`	26	`character*80 inpfile ! Input filename`
27	`character*80 outfile ! Output filename`	27	`character*80 outfile ! Output filename`
28		28
29	`c Trajectories`	29	`c Trajectories`
30	`integer ntra_inp ! Number of trajectories`	30	`integer ntra_inp ! Number of trajectories`
31	`integer ntim_inp ! Number of times`	31	`integer ntim_inp ! Number of times`
32	`integer ncol_inp ! Number of columns`	32	`integer ncol_inp ! Number of columns`
33	`real,allocatable, dimension (:,:,:) :: tra_inp ! Trajectories (ntra,ntim,ncol)`	33	`real,allocatable, dimension (:,:,:) :: tra_inp ! Trajectories (ntra,ntim,ncol)`
34	`character*80 vars_inp(100) ! Variable names`	34	`character*80 vars_inp(100) ! Variable names`
35		35
36	`integer ntra_out ! Number of trajectories`	36	`integer ntra_out ! Number of trajectories`
37	`integer ntim_out ! Number of times`	37	`integer ntim_out ! Number of times`
38	`integer ncol_out ! Number of columns`	38	`integer ncol_out ! Number of columns`
39	`real,allocatable, dimension (:,:,:) :: tra_out ! Trajectories (ntra,ntim,ncol)`	39	`real,allocatable, dimension (:,:,:) :: tra_out ! Trajectories (ntra,ntim,ncol)`
40	`integer,allocatable, dimension (:) :: ind ! Index for selection`	40	`integer,allocatable, dimension (:) :: ind ! Index for selection`
41	`integer,allocatable, dimension (:) :: isok ! Index for selection`	41	`integer,allocatable, dimension (:) :: isok ! Index for selection`
42	`character*80 vars_out(100) ! Variable names`	42	`character*80 vars_out(100) ! Variable names`
43		43
44	`real time_inp(500) ! Times of input trajectory`	44	`real time_inp(10000) ! Times of input trajectory`
45	`real time_out(500) ! Times of output trajectory`	45	`real time_out(10000) ! Times of output trajectory`
46	`integer refdate(6) ! Reference date`	46	`integer refdate(6) ! Reference date`
47	`integer ind_time(500) ! Index for time selection`	47	`integer ind_time(10000) ! Index for time selection`
48		48
49	`c Auxiliary variables`	49	`c Auxiliary variables`
50	`integer inpmode`	50	`integer inpmode`
51	`integer outmode`	51	`integer outmode`
52	`integer stat`	52	`integer stat`
53	`integer fid`	53	`integer fid`
54	`integer i,j,k,n,j0,j1`	54	`integer i,j,k,n,j0,j1`
55	`character*80 str`	55	`character*80 str`
56	`character*80 split_str(100)`	56	`character*80 split_str(100)`
57	`integer split_n`	57	`integer split_n`
58	`integer isstr,nvars,ileft,iright`	58	`integer isstr,nvars,ileft,iright`
59	`character*80 vars(100)`	59	`character*80 vars(100)`
60	`character ch`	60	`character ch`
61	`real tmp0,tmp1`	61	`real tmp0,tmp1`
62	`integer ind1`	62	`integer ind1`
63	`character*2000 linestr`	63	`character*2000 linestr`
64	`integer istr(100)`	64	`integer istr(100)`
65	`integer nstr`	65	`integer nstr`
66	`character*80 strsplit(100)`	66	`character*80 strsplit(100)`
67	`integer flag`	67	`integer flag`
68		68
69	`c ----------------------------------------------------------------------`	69	`c ----------------------------------------------------------------------`
70	`c Read and handle parameters`	70	`c Read and handle parameters`
71	`c ----------------------------------------------------------------------`	71	`c ----------------------------------------------------------------------`
72		72
73	`c Read parameters`	73	`c Read parameters`
74	`open(10,file='extract.param')`	74	`open(10,file='extract.param')`
75	`read(10,*) inpfile`	75	`read(10,*) inpfile`
76	`read(10,*) outfile`	76	`read(10,*) outfile`
77	`read(10,*) mode`	77	`read(10,*) mode`
78	`read(10,*) ntra_inp,ntim_inp,ncol_inp`	78	`read(10,*) ntra_inp,ntim_inp,ncol_inp`
79	`read(10,*) str`	79	`read(10,*) str`
80	`close(10)`	80	`close(10)`
81		81
82	`c Split the input string`	82	`c Split the input string`
83	`isstr = 0`	83	`isstr = 0`
84	`split_n = 0`	84	`split_n = 0`
85	`do i=1,80`	85	`do i=1,80`
86	`ch = str(i:i)`	86	`ch = str(i:i)`
87	`if ( (isstr.eq.0).and.(ch.ne.' ') ) then`	87	`if ( (isstr.eq.0).and.(ch.ne.' ') ) then`
88	`isstr=1`	88	`isstr=1`
89	`ileft=i`	89	`ileft=i`
90	`elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then`	90	`elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then`
91	`isstr = 0`	91	`isstr = 0`
92	`iright = i-1`	92	`iright = i-1`
93	`split_n = split_n+1`	93	`split_n = split_n+1`
94	`split_str(split_n) = str(ileft:iright)`	94	`split_str(split_n) = str(ileft:iright)`
95	`endif`	95	`endif`
96	`enddo`	96	`enddo`
97		97
98	`c Determine the formats`	98	`c Determine the formats`
99	`call mode_tra(inpmode,inpfile)`	99	`call mode_tra(inpmode,inpfile)`
100	`if (inpmode.eq.-1) inpmode=1`	100	`if (inpmode.eq.-1) inpmode=1`
101	`call mode_tra(outmode,outfile)`	101	`call mode_tra(outmode,outfile)`
102	`if ( (mode.ne.'-startf').and.(outmode.eq.-1) ) then`	102	`if ( (mode.ne.'-startf').and.(outmode.eq.-1) ) then`
103	`outmode=1`	103	`outmode=1`
104	`endif`	104	`endif`
105		105
106	`c ----------------------------------------------------------------------`	106	`c ----------------------------------------------------------------------`
107	`c Read input trajectories`	107	`c Read input trajectories`
108	`c ----------------------------------------------------------------------`	108	`c ----------------------------------------------------------------------`
109		109
110	`c Allocate memory`	110	`c Allocate memory`
111	`allocate(tra_inp(ntra_inp,ntim_inp,ncol_inp),stat=stat)`	111	`allocate(tra_inp(ntra_inp,ntim_inp,ncol_inp),stat=stat)`
112	`if (stat.ne.0) print,' error allocating array tra_inp *'`	112	`if (stat.ne.0) print,' error allocating array tra_inp *'`
113	`allocate(ind(ntra_inp),stat=stat)`	113	`allocate(ind(ntra_inp),stat=stat)`
114	`if (stat.ne.0) print,' error allocating array ind *'`	114	`if (stat.ne.0) print,' error allocating array ind *'`
115	`allocate(isok(ntra_inp),stat=stat)`	115	`allocate(isok(ntra_inp),stat=stat)`
116	`if (stat.ne.0) print,' error allocating array isok *'`	116	`if (stat.ne.0) print,' error allocating array isok *'`
117		117
118	`c Read inpufile`	118	`c Read inpufile`
119	`fid = 10`	119	`fid = 10`
120	`call ropen_tra(fid,inpfile,ntra_inp,ntim_inp,ncol_inp,`	120	`call ropen_tra(fid,inpfile,ntra_inp,ntim_inp,ncol_inp,`
121	`> refdate,vars_inp,inpmode)`	121	`> refdate,vars_inp,inpmode)`
122	`call read_tra (fid,tra_inp,ntra_inp,ntim_inp,ncol_inp,inpmode)`	122	`call read_tra (fid,tra_inp,ntra_inp,ntim_inp,ncol_inp,inpmode)`
123	`call close_tra(fid,inpmode)`	123	`call close_tra(fid,inpmode)`
124		124
125	`c Check format`	125	`c Check format`
126	`if ( vars_inp(1).ne.'time') goto 990`	126	`if ( vars_inp(1).ne.'time') goto 990`
127	`if ( (vars_inp(2).ne.'lon').and.(vars_inp(2).ne.'xpos') ) goto 990`	127	`if ( (vars_inp(2).ne.'lon').and.(vars_inp(2).ne.'xpos') ) goto 990`
128	`if ( (vars_inp(3).ne.'lat').and.(vars_inp(3).ne.'ypos') ) goto 990`	128	`if ( (vars_inp(3).ne.'lat').and.(vars_inp(3).ne.'ypos') ) goto 990`
129	`if ( (vars_inp(4).ne.'p' ).and.(vars_inp(4).ne.'ppos') ) goto 990`	129	`if ( (vars_inp(4).ne.'p' ).and.(vars_inp(4).ne.'ppos') ) goto 990`
130		130
131	`c ----------------------------------------------------------------------`	131	`c ----------------------------------------------------------------------`
132	`c Option -vars : Extract columns of variables`	132	`c Option -vars : Extract columns of variables`
133	`c ----------------------------------------------------------------------`	133	`c ----------------------------------------------------------------------`
134		134
135	`if ( mode.ne.'-var' ) goto 100`	135	`if ( mode.ne.'-var' ) goto 100`
136		136
137	`c Set the first for columns of the output`	137	`c Set the first for columns of the output`
138	`ncol_out = 4`	138	`ncol_out = 4`
139	`vars_out(1)='time'`	139	`vars_out(1)='time'`
140	`vars_out(2)='lon'`	140	`vars_out(2)='lon'`
141	`vars_out(3)='lat'`	141	`vars_out(3)='lat'`
142	`vars_out(4)='p'`	142	`vars_out(4)='p'`
143	`ind(1) =1`	143	`ind(1) =1`
144	`ind(2) =2`	144	`ind(2) =2`
145	`ind(3) =3`	145	`ind(3) =3`
146	`ind(4) =4`	146	`ind(4) =4`
147		147
148	`c Get final list of extraction columns (set number of columns)`	148	`c Get final list of extraction columns (set number of columns)`
149	`do i=1,split_n`	149	`do i=1,split_n`
150		150
151	`if (split_str(i).eq.'to') then`	151	`if (split_str(i).eq.'to') then`
152	`do j=1,ncol_inp`	152	`do j=1,ncol_inp`
153	`if ( vars_inp(j).eq.split_str(i-1) ) j0 = j + 1`	153	`if ( vars_inp(j).eq.split_str(i-1) ) j0 = j + 1`
154	`enddo`	154	`enddo`
155	`do j=1,ncol_inp`	155	`do j=1,ncol_inp`
156	`if ( vars_inp(j).eq.split_str(i+1) ) j1 = j - 1`	156	`if ( vars_inp(j).eq.split_str(i+1) ) j1 = j - 1`
157	`enddo`	157	`enddo`
158	`do j=j0,j1`	158	`do j=j0,j1`
159	`ncol_out = ncol_out + 1`	159	`ncol_out = ncol_out + 1`
160	`vars_out(ncol_out) = vars_inp(j)`	160	`vars_out(ncol_out) = vars_inp(j)`
161	`enddo`	161	`enddo`
162	`else`	162	`else`
163	`ncol_out = ncol_out + 1`	163	`ncol_out = ncol_out + 1`
164	`vars_out(ncol_out) = split_str(i)`	164	`vars_out(ncol_out) = split_str(i)`
165	`endif`	165	`endif`
166		166
167	`enddo`	167	`enddo`
168		168
169	`c Set the dimensions of the output trajectory`	169	`c Set the dimensions of the output trajectory`
170	`ntra_out = ntra_inp`	170	`ntra_out = ntra_inp`
171	`ntim_out = ntim_inp`	171	`ntim_out = ntim_inp`
172		172
173	`c Allocate memory for output trajectory`	173	`c Allocate memory for output trajectory`
174	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	174	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
175	`if (stat.ne.0) print,' error allocating array tra_out *'`	175	`if (stat.ne.0) print,' error allocating array tra_out *'`
176		176
177	`c Extract <time,lon,lat,p> columns`	177	`c Extract <time,lon,lat,p> columns`
178	`do i=1,ntra_out`	178	`do i=1,ntra_out`
179	`do j=1,ntim_out`	179	`do j=1,ntim_out`
180	`tra_out(i,j,1) = tra_inp(i,j,1)`	180	`tra_out(i,j,1) = tra_inp(i,j,1)`
181	`tra_out(i,j,2) = tra_inp(i,j,2)`	181	`tra_out(i,j,2) = tra_inp(i,j,2)`
182	`tra_out(i,j,3) = tra_inp(i,j,3)`	182	`tra_out(i,j,3) = tra_inp(i,j,3)`
183	`tra_out(i,j,4) = tra_inp(i,j,4)`	183	`tra_out(i,j,4) = tra_inp(i,j,4)`
184	`enddo`	184	`enddo`
185	`enddo`	185	`enddo`
186		186
187	`c Get indices for new columns (1..4 are already ok)`	187	`c Get indices for new columns (1..4 are already ok)`
188	`do i=5,ncol_out`	188	`do i=5,ncol_out`
189	`ind(i)=0`	189	`ind(i)=0`
190	`enddo`	190	`enddo`
191	`do i=5,ncol_out`	191	`do i=5,ncol_out`
192	`do j=1,ncol_inp`	192	`do j=1,ncol_inp`
193	`if ( vars_inp(j).eq.vars_out(i) ) ind(i) = j`	193	`if ( vars_inp(j).eq.vars_out(i) ) ind(i) = j`
194	`enddo`	194	`enddo`
195	`enddo`	195	`enddo`
196		196
197	`c Check if all selected columns are available`	197	`c Check if all selected columns are available`
198	`do i=1,ncol_out`	198	`do i=1,ncol_out`
199	`if ( ind(i).eq.0 ) then`	199	`if ( ind(i).eq.0 ) then`
200	`print*,'Invalid column in ',trim(str)`	200	`print*,'Invalid column in ',trim(str)`
201	`stop`	201	`stop`
202	`endif`	202	`endif`
203	`enddo`	203	`enddo`
204		204
205	`c Extract the column`	205	`c Extract the column`
206	`do i=1,ntra_out`	206	`do i=1,ntra_out`
207	`do j=1,ntim_out`	207	`do j=1,ntim_out`
208	`do k=1,ncol_out`	208	`do k=1,ncol_out`
209	`tra_out(i,j,k) = tra_inp(i,j,ind(k))`	209	`tra_out(i,j,k) = tra_inp(i,j,ind(k))`
210	`enddo`	210	`enddo`
211	`enddo`	211	`enddo`
212	`enddo`	212	`enddo`
213		213
214	`100 continue`	214	`100 continue`
215		215
216	`c ----------------------------------------------------------------------`	216	`c ----------------------------------------------------------------------`
217	`c Option -times : Extract times of trajectories`	217	`c Option -times : Extract times of trajectories`
218	`c ----------------------------------------------------------------------`	218	`c ----------------------------------------------------------------------`
219		219
220	`if ( mode.ne.'-time' ) goto 110`	220	`if ( mode.ne.'-time' ) goto 110`
221		221
222	`c Set the dimension of the output trajectory`	222	`c Set the dimension of the output trajectory`
223	`ntim_out = 0`	223	`ntim_out = 0`
224		224
225	`c Get the list of times for the input trajectory`	225	`c Get the list of times for the input trajectory`
226	`do i=1,ntim_inp`	226	`do i=1,ntim_inp`
227	`time_inp(i) = tra_inp(1,i,1)`	227	`time_inp(i) = tra_inp(1,i,1)`
228	`enddo`	228	`enddo`
229		229
230	`c Get final list of extraction times (set number of times)`	230	`c Get final list of extraction times (set number of times)`
231	`do i=1,split_n`	231	`do i=1,split_n`
232		232
233	`if (split_str(i).eq.'to') then`	233	`if (split_str(i).eq.'to') then`
234	`read(split_str(i-1),*) tmp0`	234	`read(split_str(i-1),*) tmp0`
235	`do j=1,ntim_inp`	235	`do j=1,ntim_inp`
236	`if ( time_inp(j).eq.tmp0 ) j0 = j + 1`	236	`if ( time_inp(j).eq.tmp0 ) j0 = j + 1`
237	`enddo`	237	`enddo`
238	`read(split_str(i+1),*) tmp0`	238	`read(split_str(i+1),*) tmp0`
239	`do j=1,ntim_inp`	239	`do j=1,ntim_inp`
240	`if ( time_inp(j).eq.tmp0 ) j1 = j - 1`	240	`if ( time_inp(j).eq.tmp0 ) j1 = j - 1`
241	`enddo`	241	`enddo`
242	`do j=j0,j1`	242	`do j=j0,j1`
243	`ntim_out = ntim_out + 1`	243	`ntim_out = ntim_out + 1`
244	`time_out(ntim_out) = time_inp(j)`	244	`time_out(ntim_out) = time_inp(j)`
245	`enddo`	245	`enddo`
246	`elseif (split_str(i).eq.'first') then`	246	`elseif (split_str(i).eq.'first') then`
247	`ntim_out = ntim_out + 1`	247	`ntim_out = ntim_out + 1`
248	`time_out(ntim_out) = time_inp(1)`	248	`time_out(ntim_out) = time_inp(1)`
249	`elseif (split_str(i).eq.'last') then`	249	`elseif (split_str(i).eq.'last') then`
250	`ntim_out = ntim_out + 1`	250	`ntim_out = ntim_out + 1`
251	`time_out(ntim_out) = time_inp(ntim_inp)`	251	`time_out(ntim_out) = time_inp(ntim_inp)`
252	`else`	252	`else`
253	`ntim_out = ntim_out + 1`	253	`ntim_out = ntim_out + 1`
254	`read(split_str(i),*) tmp0`	254	`read(split_str(i),*) tmp0`
255	`time_out(ntim_out) = tmp0`	255	`time_out(ntim_out) = tmp0`
256	`endif`	256	`endif`
257		257
258	`enddo`	258	`enddo`
259		259
260	`c Get the indices of the selected times`	260	`c Get the indices of the selected times`
261	`do i=1,ntim_out`	261	`do i=1,ntim_out`
262	`ind_time(i) = 0`	262	`ind_time(i) = 0`
263	`enddo`	263	`enddo`
264	`do i=1,ntim_out`	264	`do i=1,ntim_out`
265	`do j=1,ntim_inp`	265	`do j=1,ntim_inp`
266	`if ( abs(time_out(i)-time_inp(j)).lt.eps) ind_time(i) = j`	266	`if ( abs(time_out(i)-time_inp(j)).lt.eps) ind_time(i) = j`
267	`enddo`	267	`enddo`
268	`enddo`	268	`enddo`
269	`do i=1,ntim_out`	269	`do i=1,ntim_out`
270	`if ( ind_time(i).eq.0) then`	270	`if ( ind_time(i).eq.0) then`
271	`print*,' Invalid time ',time_out(i)`	271	`print*,' Invalid time ',time_out(i)`
272	`stop`	272	`stop`
273	`endif`	273	`endif`
274	`enddo`	274	`enddo`
275		275
276	`c Set dimensions of output trajectory`	276	`c Set dimensions of output trajectory`
277	`ntra_out = ntra_inp`	277	`ntra_out = ntra_inp`
278	`ncol_out = ncol_inp`	278	`ncol_out = ncol_inp`
279		279
280	`c Allocate memory for output trajectory`	280	`c Allocate memory for output trajectory`
281	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	281	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
282	`if (stat.ne.0) print,' error allocating array tra_out *'`	282	`if (stat.ne.0) print,' error allocating array tra_out *'`
283		283
284	`c Copy the selected times to the output trajectory`	284	`c Copy the selected times to the output trajectory`
285	`do i=1,ntra_out`	285	`do i=1,ntra_out`
286	`do j=1,ntim_out`	286	`do j=1,ntim_out`
287	`do k=1,ncol_out`	287	`do k=1,ncol_out`
288	`tra_out(i,j,k) = tra_inp(i,ind_time(j),k)`	288	`tra_out(i,j,k) = tra_inp(i,ind_time(j),k)`
289	`enddo`	289	`enddo`
290	`enddo`	290	`enddo`
291	`enddo`	291	`enddo`
292		292
293	`c Copy meta information`	293	`c Copy meta information`
294	`do i=1,ncol_out`	294	`do i=1,ncol_out`
295	`vars_out(i) = vars_inp(i)`	295	`vars_out(i) = vars_inp(i)`
296	`enddo`	296	`enddo`
297		297
298	`110 continue`	298	`110 continue`
299		299
300	`c ----------------------------------------------------------------------`	300	`c ----------------------------------------------------------------------`
301	`c Option -tra : Extract trajectories by number`	301	`c Option -tra : Extract trajectories by number`
302	`c ----------------------------------------------------------------------`	302	`c ----------------------------------------------------------------------`
303		303
304	`if ( mode.ne.'-tra' ) goto 120`	304	`if ( mode.ne.'-tra' ) goto 120`
305		305
306	`c Set the dimension of the output trajectory`	306	`c Set the dimension of the output trajectory`
307	`ntra_out = 0`	307	`ntra_out = 0`
308		308
309	`c Get final list of extraction times (set number of times)`	309	`c Get final list of extraction times (set number of times)`
310	`do i=1,split_n`	310	`do i=1,split_n`
311		311
312	`if (split_str(i).eq.'to') then`	312	`if (split_str(i).eq.'to') then`
313	`read(split_str(i-1),*) tmp0`	313	`read(split_str(i-1),*) tmp0`
314	`read(split_str(i+1),*) tmp1`	314	`read(split_str(i+1),*) tmp1`
315	`do j=nint(tmp0)+1,nint(tmp1)-1`	315	`do j=nint(tmp0)+1,nint(tmp1)-1`
316	`ntra_out = ntra_out + 1`	316	`ntra_out = ntra_out + 1`
317	`ind(ntra_out) = j`	317	`ind(ntra_out) = j`
318	`enddo`	318	`enddo`
319	`elseif (split_str(i).eq.'first') then`	319	`elseif (split_str(i).eq.'first') then`
320	`ntra_out = ntra_out + 1`	320	`ntra_out = ntra_out + 1`
321	`ind(ntra_out) = 1`	321	`ind(ntra_out) = 1`
322	`elseif (split_str(i).eq.'last') then`	322	`elseif (split_str(i).eq.'last') then`
323	`ntra_out = ntra_out + 1`	323	`ntra_out = ntra_out + 1`
324	`ind(ntra_out) = ntra_inp`	324	`ind(ntra_out) = ntra_inp`
325	`else`	325	`else`
326	`ntra_out = ntra_out + 1`	326	`ntra_out = ntra_out + 1`
327	`read(split_str(i),*) tmp0`	327	`read(split_str(i),*) tmp0`
328	`ind(ntra_out) = nint(tmp0)`	328	`ind(ntra_out) = nint(tmp0)`
329	`endif`	329	`endif`
330		330
331	`enddo`	331	`enddo`
332		332
333	`c Check whether selected trajectories are ok`	333	`c Check whether selected trajectories are ok`
334	`do i=1,ntra_out`	334	`do i=1,ntra_out`
335	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`	335	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`
336	`print*,'Invalid trajectory selected ',ind(i)`	336	`print*,'Invalid trajectory selected ',ind(i)`
337	`stop`	337	`stop`
338	`endif`	338	`endif`
339	`enddo`	339	`enddo`
340		340
341	`c Set dimensions of output trajectory`	341	`c Set dimensions of output trajectory`
342	`ntim_out = ntim_inp`	342	`ntim_out = ntim_inp`
343	`ncol_out = ncol_inp`	343	`ncol_out = ncol_inp`
344		344
345	`c Allocate memory for output trajectory`	345	`c Allocate memory for output trajectory`
346	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	346	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
347	`if (stat.ne.0) print,' error allocating array tra_out *'`	347	`if (stat.ne.0) print,' error allocating array tra_out *'`
348		348
349	`c Copy the selected times to the output trajectory`	349	`c Copy the selected times to the output trajectory`
350	`do i=1,ntra_out`	350	`do i=1,ntra_out`
351	`do j=1,ntim_out`	351	`do j=1,ntim_out`
352	`do k=1,ncol_out`	352	`do k=1,ncol_out`
353	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`	353	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`
354	`enddo`	354	`enddo`
355	`enddo`	355	`enddo`
356	`enddo`	356	`enddo`
357		357
358	`c Copy meta information`	358	`c Copy meta information`
359	`do i=1,ncol_out`	359	`do i=1,ncol_out`
360	`vars_out(i) = vars_inp(i)`	360	`vars_out(i) = vars_inp(i)`
361	`enddo`	361	`enddo`
362		362
363	`120 continue`	363	`120 continue`
364		364
365	`c ----------------------------------------------------------------------`	365	`c ----------------------------------------------------------------------`
366	`c Option -startf : Extract starting positions for the trajectory file`	366	`c Option -startf : Extract starting positions for the trajectory file`
367	`c ----------------------------------------------------------------------`	367	`c ----------------------------------------------------------------------`
368		368
369	`if ( mode.ne.'-startf' ) goto 130`	369	`if ( mode.ne.'-startf' ) goto 130`
370		370
371	`c Set the first for columns of the output`	371	`c Set the first for columns of the output`
372	`ncol_out = 4`	372	`ncol_out = 4`
373	`vars_out(1)='time'`	373	`vars_out(1)='time'`
374	`vars_out(2)='lon'`	374	`vars_out(2)='lon'`
375	`vars_out(3)='lat'`	375	`vars_out(3)='lat'`
376	`vars_out(4)='p'`	376	`vars_out(4)='p'`
377	`ind(1) =1`	377	`ind(1) =1`
378	`ind(2) =2`	378	`ind(2) =2`
379	`ind(3) =3`	379	`ind(3) =3`
380	`ind(4) =4`	380	`ind(4) =4`
381		381
382	`c Set dimensions of output trajectory`	382	`c Set dimensions of output trajectory`
383	`ntim_out = 1`	383	`ntim_out = 1`
384	`ntra_out = ntra_inp`	384	`ntra_out = ntra_inp`
385		385
386	`c Allocate memory for output trajectory`	386	`c Allocate memory for output trajectory`
387	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	387	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
388	`if (stat.ne.0) print,' error allocating array tra_out *'`	388	`if (stat.ne.0) print,' error allocating array tra_out *'`
389		389
390	`c Copy the selected times to the output trajectory`	390	`c Copy the selected times to the output trajectory`
391	`do i=1,ntra_out`	391	`do i=1,ntra_out`
392	`do k=1,ncol_out`	392	`do k=1,ncol_out`
393	`tra_out(i,1,k) = tra_inp(i,1,k)`	393	`tra_out(i,1,k) = tra_inp(i,1,k)`
394	`enddo`	394	`enddo`
395	`enddo`	395	`enddo`
396		396
397	`c Copy meta information`	397	`c Copy meta information`
398	`do i=1,ncol_out`	398	`do i=1,ncol_out`
399	`vars_out(i) = vars_inp(i)`	399	`vars_out(i) = vars_inp(i)`
400	`enddo`	400	`enddo`
401		401
402	`130 continue`	402	`130 continue`
403		403
404	`c ----------------------------------------------------------------------`	404	`c ----------------------------------------------------------------------`
405	`c Option -index : Extract trajectories by index file`	405	`c Option -index : Extract trajectories by index file`
406	`c ----------------------------------------------------------------------`	406	`c ----------------------------------------------------------------------`
407		407
408	`if ( mode.ne.'-index' ) goto 140`	408	`if ( mode.ne.'-index' ) goto 140`
409		409
410	`c Read the index file`	410	`c Read the index file`
411	`open(10,file=str)`	411	`open(10,file=str)`
412		412
413	`ntra_out = 1`	413	`ntra_out = 1`
414	`142 read(10,*,end=141) ind(ntra_out)`	414	`142 read(10,*,end=141) ind(ntra_out)`
415	`ntra_out = ntra_out + 1`	415	`ntra_out = ntra_out + 1`
416	`goto 142`	416	`goto 142`
417	`141 continue`	417	`141 continue`
418	`ntra_out = ntra_out - 1`	418	`ntra_out = ntra_out - 1`
419		419
420	`close(10)`	420	`close(10)`
421		421
422	`c Check whether selected trajectories are ok`	422	`c Check whether selected trajectories are ok`
423	`do i=1,ntra_out`	423	`do i=1,ntra_out`
424	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`	424	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`
425	`print*,'Invalid trajectory selected ',ind(i)`	425	`print*,'Invalid trajectory selected ',ind(i)`
426	`stop`	426	`stop`
427	`endif`	427	`endif`
428	`enddo`	428	`enddo`
429		429
430	`c Set dimensions of output trajectory`	430	`c Set dimensions of output trajectory`
431	`ntim_out = ntim_inp`	431	`ntim_out = ntim_inp`
432	`ncol_out = ncol_inp`	432	`ncol_out = ncol_inp`
433		433
434	`c Allocate memory for output trajectory`	434	`c Allocate memory for output trajectory`
435	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	435	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
436	`if (stat.ne.0) print,' error allocating array tra_out *'`	436	`if (stat.ne.0) print,' error allocating array tra_out *'`
437		437
438	`c Copy the selected times to the output trajectory`	438	`c Copy the selected times to the output trajectory`
439	`do i=1,ntra_out`	439	`do i=1,ntra_out`
440	`do j=1,ntim_out`	440	`do j=1,ntim_out`
441	`do k=1,ncol_out`	441	`do k=1,ncol_out`
442	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`	442	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`
443	`enddo`	443	`enddo`
444	`enddo`	444	`enddo`
445	`enddo`	445	`enddo`
446		446
447	`c Copy meta information`	447	`c Copy meta information`
448	`do i=1,ncol_out`	448	`do i=1,ncol_out`
449	`vars_out(i) = vars_inp(i)`	449	`vars_out(i) = vars_inp(i)`
450	`enddo`	450	`enddo`
451		451
452	`140 continue`	452	`140 continue`
453		453
454	`c ----------------------------------------------------------------------`	454	`c ----------------------------------------------------------------------`
455	`c Option -boolean : Extract trajectories by boolean file`	455	`c Option -boolean : Extract trajectories by boolean file`
456	`c ----------------------------------------------------------------------`	456	`c ----------------------------------------------------------------------`
457		457
458	`if ( mode.ne.'-boolean' ) goto 150`	458	`if ( mode.ne.'-boolean' ) goto 150`
459		459
460	`c Read the index file`	460	`c Read the index file`
461	`open(10,file=str)`	461	`open(10,file=str)`
462	`ntra_out = 0`	462	`ntra_out = 0`
463	`do i=1,ntra_inp`	463	`do i=1,ntra_inp`
464	`read(10,*) ind1`	464	`read(10,*) ind1`
465	`if ( ind1.eq.1 ) then`	465	`if ( ind1.eq.1 ) then`
466	`ntra_out = ntra_out + 1`	466	`ntra_out = ntra_out + 1`
467	`ind(ntra_out) = i`	467	`ind(ntra_out) = i`
468	`endif`	468	`endif`
469	`enddo`	469	`enddo`
470	`close(10)`	470	`close(10)`
471		471
472	`c Check whether selected trajectories are ok`	472	`c Check whether selected trajectories are ok`
473	`do i=1,ntra_out`	473	`do i=1,ntra_out`
474	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`	474	`if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then`
475	`print*,'Invalid trajectory selected ',ind(i)`	475	`print*,'Invalid trajectory selected ',ind(i)`
476	`stop`	476	`stop`
477	`endif`	477	`endif`
478	`enddo`	478	`enddo`
479		479
480	`c Set dimensions of output trajectory`	480	`c Set dimensions of output trajectory`
481	`ntim_out = ntim_inp`	481	`ntim_out = ntim_inp`
482	`ncol_out = ncol_inp`	482	`ncol_out = ncol_inp`
483		483
484	`c Allocate memory for output trajectory`	484	`c Allocate memory for output trajectory`
485	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`	485	`allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)`
486	`if (stat.ne.0) print,' error allocating array tra_out *'`	486	`if (stat.ne.0) print,' error allocating array tra_out *'`
487		487
488	`c Copy the selected times to the output trajectory`	488	`c Copy the selected times to the output trajectory`
489	`do i=1,ntra_out`	489	`do i=1,ntra_out`
490	`do j=1,ntim_out`	490	`do j=1,ntim_out`
491	`do k=1,ncol_out`	491	`do k=1,ncol_out`
492	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`	492	`tra_out(i,j,k) = tra_inp(ind(i),j,k)`
493	`enddo`	493	`enddo`
494	`enddo`	494	`enddo`
495	`enddo`	495	`enddo`
496		496
497	`c Copy meta information`	497	`c Copy meta information`
498	`do i=1,ncol_out`	498	`do i=1,ncol_out`
499	`vars_out(i) = vars_inp(i)`	499	`vars_out(i) = vars_inp(i)`
500	`enddo`	500	`enddo`

Subversion Repositories lagranto.ecmwf

(root)/trunk/goodies/extract.f – Rev 5 → 27