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

Subversion Repositories lagranto.ecmwf

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