Subversion Repositories lagranto.ecmwf

Rev

Rev 5 | Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
      PROGRAM extract
2
 
3
c     ***********************************************************************
4
c     * Extract columns of a trajectory file                                *
5
c     * Michael Sprenger / Spring, summer 2010                              *
6
c     ***********************************************************************
7
 
8
      implicit none
9
 
10
c     ----------------------------------------------------------------------
11
c     Declaration of parameters
12
c     ----------------------------------------------------------------------
13
 
14
c     Numerical epsilon
15
      real      eps
16
      parameter (eps=0.0001)
17
 
18
c     ----------------------------------------------------------------------
19
c     Declaration of variables
20
c     ----------------------------------------------------------------------
21
 
22
c     Extraction mode
23
      character*80                           mode                             
24
 
25
c     Input and output format for trajectories (see iotra.f)
26
      character*80                           inpfile       ! Input filename
27
      character*80                           outfile       ! Output filename
28
 
29
c     Trajectories
30
      integer                                ntra_inp      ! Number of trajectories
31
      integer                                ntim_inp      ! Number of times
32
      integer                                ncol_inp      ! Number of columns
33
      real,allocatable, dimension (:,:,:) :: tra_inp       ! Trajectories (ntra,ntim,ncol)
34
      character*80                           vars_inp(100) ! Variable names
35
 
36
      integer                                ntra_out      ! Number of trajectories
37
      integer                                ntim_out      ! Number of times
38
      integer                                ncol_out      ! Number of columns
39
      real,allocatable, dimension (:,:,:) :: tra_out       ! Trajectories (ntra,ntim,ncol)
40
      integer,allocatable, dimension (:)  :: ind           ! Index for selection
41
      character*80                           vars_out(100) ! Variable names
42
 
43
      real                                   time_inp(500) ! Times of input trajectory
44
      real                                   time_out(500) ! Times of output trajectory
45
      integer                                refdate(6)    ! Reference date
46
 
47
c     Auxiliary variables
48
      integer                                inpmode
49
      integer                                outmode
50
      integer                                stat
51
      integer                                fid
52
      integer                                i,j,k,n,j0,j1
53
      character*80                           str
54
      character*80                           split_str(100)
55
      integer                                split_n
56
      integer                                isstr,nvars,ileft,iright
57
      character*80                           vars(100)
58
      character                              ch
59
      real                                   tmp0,tmp1
60
      integer                                ind1
61
	  character*2000                         linestr
62
	  integer                                istr(100)
63
      integer                                nstr
64
      character*80                           strsplit(100)
65
      integer                                flag
66
 
67
c     ----------------------------------------------------------------------
68
c     Read and handle parameters
69
c     ----------------------------------------------------------------------
70
 
71
c     Read parameters
72
      open(10,file='extract.param')
73
       read(10,*) inpfile
74
       read(10,*) outfile
75
       read(10,*) mode
76
       read(10,*) ntra_inp,ntim_inp,ncol_inp
77
       read(10,*) str
78
      close(10)
79
 
80
c     Split the input string
81
      isstr   = 0
82
      split_n = 0
83
      do i=1,80
84
         ch = str(i:i)
85
         if ( (isstr.eq.0).and.(ch.ne.' ') ) then
86
            isstr=1
87
            ileft=i
88
         elseif ( (isstr.eq.1).and.(ch.eq.' ') ) then
89
            isstr              = 0
90
            iright             = i-1
91
            split_n            = split_n+1
92
            split_str(split_n) = str(ileft:iright)
93
         endif
94
      enddo
95
 
96
c     Determine the formats
97
      call mode_tra(inpmode,inpfile)
98
      if (inpmode.eq.-1) inpmode=1
99
      call mode_tra(outmode,outfile)
100
      if ( (mode.ne.'-startf').and.(outmode.eq.-1) ) then
101
         outmode=1
102
      endif
103
 
104
c     ----------------------------------------------------------------------
105
c     Read input trajectories
106
c     ----------------------------------------------------------------------
107
 
108
c     Allocate memory
109
      allocate(tra_inp(ntra_inp,ntim_inp,ncol_inp),stat=stat)
110
      if (stat.ne.0) print*,'*** error allocating array tra_inp    ***' 
111
      allocate(ind(ntra_inp),stat=stat)
112
      if (stat.ne.0) print*,'*** error allocating array ind        ***' 
113
 
114
c     Read inpufile
115
      fid = 10
116
      call ropen_tra(fid,inpfile,ntra_inp,ntim_inp,ncol_inp,
117
     >               refdate,vars_inp,inpmode)
118
      call read_tra (fid,tra_inp,ntra_inp,ntim_inp,ncol_inp,inpmode)
119
      call close_tra(fid,inpmode)
120
 
121
c     Check format 
122
      if (  vars_inp(1).ne.'time') goto 990
123
      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
125
      if ( (vars_inp(4).ne.'p'  ).and.(vars_inp(4).ne.'ppos') ) goto 990
126
 
127
c     ----------------------------------------------------------------------
128
c     Option -vars : Extract columns of variables
129
c     ----------------------------------------------------------------------
130
 
131
      if ( mode.ne.'-var' ) goto 100
132
 
133
c     Set the first for columns of the output
134
      ncol_out   = 4
135
      vars_out(1)='time'
136
      vars_out(2)='lon'
137
      vars_out(3)='lat'
138
      vars_out(4)='p'
139
      ind(1)     =1
140
      ind(2)     =2
141
      ind(3)     =3
142
      ind(4)     =4
143
 
144
c     Get final list of extraction columns (set number of columns)
145
      do i=1,split_n
146
 
147
         if (split_str(i).eq.'to') then
148
            do j=1,ncol_inp
149
               if ( vars_inp(j).eq.split_str(i-1) ) j0 = j + 1
150
            enddo
151
            do j=1,ncol_inp
152
               if ( vars_inp(j).eq.split_str(i+1) ) j1 = j - 1
153
            enddo
154
            do j=j0,j1
155
               ncol_out           = ncol_out + 1
156
               vars_out(ncol_out) = vars_inp(j)
157
            enddo
158
         else
159
            ncol_out           = ncol_out + 1
160
            vars_out(ncol_out) = split_str(i)
161
         endif
162
 
163
      enddo
164
 
165
c     Set the dimensions of the output trajectory
166
      ntra_out = ntra_inp
167
      ntim_out = ntim_inp
168
 
169
c     Allocate memory for output trajectory
170
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
171
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
172
 
173
c     Extract <time,lon,lat,p> columns
174
      do i=1,ntra_out
175
         do j=1,ntim_out
176
            tra_out(i,j,1) = tra_inp(i,j,1)
177
            tra_out(i,j,2) = tra_inp(i,j,2)
178
            tra_out(i,j,3) = tra_inp(i,j,3)
179
            tra_out(i,j,4) = tra_inp(i,j,4)
180
         enddo
181
      enddo
182
 
183
c     Get indices for new columns (1..4 are already ok)
184
      do i=5,ncol_out
185
         ind(i)=0
186
      enddo
187
      do i=5,ncol_out
188
         do j=1,ncol_inp
189
            if ( vars_inp(j).eq.vars_out(i) ) ind(i) = j
190
         enddo
191
      enddo
192
 
193
c     Check if all selected columns are available
194
      do i=1,ncol_out
195
         if ( ind(i).eq.0 ) then
196
            print*,'Invalid column in ',trim(str)
197
            stop
198
         endif
199
      enddo
200
 
201
c     Extract the column
202
      do i=1,ntra_out
203
         do j=1,ntim_out
204
            do k=1,ncol_out
205
               tra_out(i,j,k) = tra_inp(i,j,ind(k))
206
            enddo
207
         enddo
208
      enddo
209
 
210
 100  continue
211
 
212
c     ----------------------------------------------------------------------
213
c     Option -times : Extract times of trajectories
214
c     ----------------------------------------------------------------------
215
 
216
      if ( mode.ne.'-time' ) goto 110
217
 
218
c     Set the dimension of the output trajectory
219
      ntim_out = 0
220
 
221
c     Get the list of times for the input trajectory
222
      do i=1,ntim_inp
223
         time_inp(i) = tra_inp(1,i,1)
224
      enddo
225
 
226
c     Get final list of extraction times (set number of times)
227
      do i=1,split_n
228
 
229
         if (split_str(i).eq.'to') then
230
            read(split_str(i-1),*) tmp0
231
            do j=1,ntim_inp
232
               if ( time_inp(j).eq.tmp0 ) j0 = j + 1
233
            enddo
234
            read(split_str(i+1),*) tmp0
235
            do j=1,ntim_inp
236
               if ( time_inp(j).eq.tmp0 ) j1 = j - 1
237
            enddo
238
            do j=j0,j1
239
               ntim_out           = ntim_out + 1
240
               time_out(ntim_out) = time_inp(j)
241
            enddo
242
         elseif (split_str(i).eq.'first') then
243
            ntim_out           = ntim_out + 1
244
            time_out(ntim_out) = time_inp(1)
245
         elseif (split_str(i).eq.'last') then
246
            ntim_out           = ntim_out + 1
247
            time_out(ntim_out) = time_inp(ntim_inp)
248
         else
249
            ntim_out           = ntim_out + 1
250
            read(split_str(i),*) tmp0
251
            time_out(ntim_out) = tmp0
252
         endif
253
 
254
      enddo
255
 
256
c     Get the indices of the selected times
257
      do i=1,ntim_out
258
         ind(i) = 0
259
      enddo
260
      do i=1,ntim_out
261
         do j=1,ntim_inp
262
            if ( abs(time_out(i)-time_inp(j)).lt.eps) ind(i) = j
263
         enddo
264
      enddo
265
      do i=1,ntim_out
266
         if ( ind(i).eq.0) then
267
            print*,' Invalid time ',time_out(i)
268
            stop
269
         endif
270
      enddo
271
 
272
c     Set dimensions of output trajectory
273
      ntra_out = ntra_inp
274
      ncol_out = ncol_inp
275
 
276
c     Allocate memory for output trajectory
277
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
278
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
279
 
280
c     Copy the selected times to the output trajectory
281
      do i=1,ntra_out
282
         do j=1,ntim_out
283
            do k=1,ncol_out
284
               tra_out(i,j,k) = tra_inp(i,ind(j),k)
285
            enddo
286
         enddo
287
      enddo
288
 
289
c     Copy meta information
290
      do i=1,ncol_out
291
         vars_out(i) = vars_inp(i)
292
      enddo
293
 
294
 110  continue
295
 
296
c     ----------------------------------------------------------------------
297
c     Option -tra : Extract trajectories by number
298
c     ----------------------------------------------------------------------
299
 
300
      if ( mode.ne.'-tra' ) goto 120
301
 
302
c     Set the dimension of the output trajectory
303
      ntra_out = 0
304
 
305
c     Get final list of extraction times (set number of times)
306
      do i=1,split_n
307
 
308
         if (split_str(i).eq.'to') then
309
            read(split_str(i-1),*) tmp0
310
            read(split_str(i+1),*) tmp1
311
            do j=nint(tmp0)+1,nint(tmp1)-1
312
               ntra_out        = ntra_out + 1
313
               ind(ntra_out)   = j   
314
            enddo
315
         elseif (split_str(i).eq.'first') then
316
            ntra_out           = ntra_out + 1
317
            ind(ntra_out)      = 1
318
         elseif (split_str(i).eq.'last') then
319
            ntra_out           = ntra_out + 1
320
            ind(ntra_out)      = ntra_inp
321
         else
322
            ntra_out           = ntra_out + 1
323
            read(split_str(i),*) tmp0
324
            ind(ntra_out)      = nint(tmp0) 
325
         endif
326
 
327
      enddo
328
 
329
c     Check whether selected trajectories are ok
330
      do i=1,ntra_out
331
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
332
            print*,'Invalid trajectory selected ',ind(i)
333
            stop
334
         endif
335
      enddo
336
 
337
c     Set dimensions of output trajectory
338
      ntim_out = ntim_inp
339
      ncol_out = ncol_inp
340
 
341
c     Allocate memory for output trajectory
342
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
343
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
344
 
345
c     Copy the selected times to the output trajectory
346
      do i=1,ntra_out
347
         do j=1,ntim_out
348
            do k=1,ncol_out
349
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
350
            enddo
351
         enddo
352
      enddo
353
 
354
c     Copy meta information
355
      do i=1,ncol_out
356
         vars_out(i) = vars_inp(i)
357
      enddo
358
 
359
 120  continue
360
 
361
c     ----------------------------------------------------------------------
362
c     Option -startf : Extract starting positions for the trajectory file
363
c     ----------------------------------------------------------------------
364
 
365
      if ( mode.ne.'-startf' ) goto 130
366
 
367
c     Set the first for columns of the output
368
      ncol_out   = 4
369
      vars_out(1)='time'
370
      vars_out(2)='lon'
371
      vars_out(3)='lat'
372
      vars_out(4)='p'
373
      ind(1)     =1
374
      ind(2)     =2
375
      ind(3)     =3
376
      ind(4)     =4
377
 
378
c     Set dimensions of output trajectory
379
      ntim_out = 1
380
      ntra_out = ntra_inp
381
 
382
c     Allocate memory for output trajectory
383
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
384
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
385
 
386
c     Copy the selected times to the output trajectory
387
      do i=1,ntra_out
388
         do k=1,ncol_out
389
            tra_out(i,1,k) = tra_inp(i,1,k)
390
         enddo
391
      enddo
392
 
393
c     Copy meta information
394
      do i=1,ncol_out
395
         vars_out(i) = vars_inp(i)
396
      enddo
397
 
398
 130  continue
399
 
400
c     ----------------------------------------------------------------------
401
c     Option -index : Extract trajectories by index file
402
c     ----------------------------------------------------------------------
403
 
404
      if ( mode.ne.'-index' ) goto 140
405
 
406
c     Read the index file
407
      open(10,file=str)
408
 
409
        ntra_out = 1
410
 142    read(10,*,end=141) ind(ntra_out)
411
        ntra_out = ntra_out + 1
412
        goto 142 
413
 141    continue
414
        ntra_out = ntra_out - 1
415
 
416
      close(10)
417
 
418
c     Check whether selected trajectories are ok
419
      do i=1,ntra_out
420
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
421
            print*,'Invalid trajectory selected ',ind(i)
422
            stop
423
         endif
424
      enddo
425
 
426
c     Set dimensions of output trajectory
427
      ntim_out = ntim_inp
428
      ncol_out = ncol_inp
429
 
430
c     Allocate memory for output trajectory
431
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
432
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
433
 
434
c     Copy the selected times to the output trajectory
435
      do i=1,ntra_out
436
         do j=1,ntim_out
437
            do k=1,ncol_out
438
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
439
            enddo
440
         enddo
441
      enddo
442
 
443
c     Copy meta information
444
      do i=1,ncol_out
445
         vars_out(i) = vars_inp(i)
446
      enddo
447
 
448
 140  continue
449
 
450
c     ----------------------------------------------------------------------
451
c     Option -boolean : Extract trajectories by boolean file
452
c     ----------------------------------------------------------------------
453
 
454
      if ( mode.ne.'-boolean' ) goto 150
455
 
456
c     Read the index file
457
      open(10,file=str)
458
        ntra_out = 0
459
        do i=1,ntra_inp
460
           read(10,*) ind1
461
           if ( ind1.eq.1 ) then
462
              ntra_out      = ntra_out + 1
463
              ind(ntra_out) = i
464
           endif
465
        enddo
466
      close(10)
467
 
468
c     Check whether selected trajectories are ok
469
      do i=1,ntra_out
470
         if ( (ind(i).lt.1).or.(ind(i).gt.ntra_inp) ) then
471
            print*,'Invalid trajectory selected ',ind(i)
472
            stop
473
         endif
474
      enddo
475
 
476
c     Set dimensions of output trajectory
477
      ntim_out = ntim_inp
478
      ncol_out = ncol_inp
479
 
480
c     Allocate memory for output trajectory
481
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
482
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
483
 
484
c     Copy the selected times to the output trajectory
485
      do i=1,ntra_out
486
         do j=1,ntim_out
487
            do k=1,ncol_out
488
               tra_out(i,j,k) = tra_inp(ind(i),j,k)
489
            enddo
490
         enddo
491
      enddo
492
 
493
c     Copy meta information
494
      do i=1,ncol_out
495
         vars_out(i) = vars_inp(i)
496
      enddo
497
 
498
 150  continue
499
 
500
c     ----------------------------------------------------------------------
501
c     Option -pattern : Extract trajectories which match a regular expression
502
c     ----------------------------------------------------------------------
503
 
504
	  if ( mode.ne.'-pattern' ) goto 160
505
 
506
c	  All times and columns are extracted
507
      ncol_out = ncol_inp
508
      ntim_out = ntim_inp
509
      ntra_out = 0
510
 
511
c     Split the search string
512
	  nstr   = 0
513
	  ileft  = 0
514
	  iright = 0
515
      do i=1,len_trim(str)
516
      	if ( (str(i:i).eq.' ').and.(ileft.eq.0) ) then
517
      	   ileft = ileft + 1
518
        elseif ( (str(i:i).ne.' ').and.(ileft.eq.0) ) then
519
           ileft  = i
520
           iright = 0
521
        elseif ( (str(i:i).ne.' ').and.(ileft.ne.0) ) then
522
           iright = i
523
        elseif ( (str(i:i).eq.' ').and.(ileft.ne.0) ) then
524
           nstr           = nstr + 1
525
           strsplit(nstr) = str(ileft:iright)
526
           ileft          = 0
527
           iright         = 0
528
        endif
529
      enddo
530
      if ( (ileft.ne.0).and.(iright.ne.0) ) then
531
           nstr           = nstr + 1
532
           strsplit(nstr) = str(ileft:iright)
533
           ileft          = 0
534
      endif
535
 
536
c     Loop over the trajectories - check for matching pattern
537
	  do n=1,ntra_inp
538
 
539
	    ind(n) = 0
540
	  	do i=1,ntim_inp
541
 
542
           write(linestr,'(1f7.2,f9.2,f8.2,i6,100f10.3)')
543
     >                   (tra_inp(n,i,j),j=1,3),             ! time, lon, lat
544
     >                   nint(tra_inp(n,i,4)),               ! p
545
     >                   (tra_inp(n,i,j),j=5,ncol_inp)       ! fields
546
 
547
           flag = 1
548
           do k=1,nstr
549
             istr(k) = index(trim(linestr),trim(strsplit(k)))
550
             if ( istr(k).eq.0 ) flag = 0
551
           enddo
552
	       if ( flag.eq.1 ) ind(n) = 1
553
 
554
	  	enddo
555
	  	if ( ind(n).eq.1 ) ntra_out = ntra_out + 1
556
 
557
	  enddo
558
 
559
c     Allocate memory for output trajectory
560
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
561
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***'
562
 
563
c     Copy the selected times to the output trajectory
564
      ntra_out = 0
565
      do i=1,ntra_inp
566
      	if (ind(i).eq.1) then
567
      	    ntra_out = ntra_out + 1
568
         	do j=1,ntim_out
569
            	do k=1,ncol_out
570
               		tra_out(ntra_out,j,k) = tra_inp(i,j,k)
571
            	enddo
572
         	enddo
573
         endif
574
      enddo
575
 
576
c     Copy meta information
577
      do i=1,ncol_out
578
         vars_out(i) = vars_inp(i)
579
      enddo
580
 
581
 160  continue
582
 
583
 
584
c     ----------------------------------------------------------------------
585
c     Write output trajectories
586
c     ----------------------------------------------------------------------
587
 
588
c     Write output as trajectory file
589
      if (outmode.ge.1) then
590
 
591
         call wopen_tra(fid,outfile,ntra_out,ntim_out,ncol_out,
592
     >                  refdate,vars_out,outmode)
593
         call write_tra(fid,tra_out,ntra_out,ntim_out,ncol_out,outmode)
594
         call close_tra(fid,outmode)
595
 
596
c     Write output as (lon, lat, p)-list
597
      else
598
 
599
         open(10,file=outfile)
600
         do i=1,ntra_out
601
            write(10,'(3f10.2)') tra_out(i,1,2),    ! lon
602
     >                           tra_out(i,1,3),    ! lat
603
     >                           tra_out(i,1,4)     ! p
604
         enddo
605
         close(10)
606
 
607
      endif
608
 
609
 
610
 
611
!c     ----------------------------------------------------------------------
612
c     Error handling
613
c     ----------------------------------------------------------------------
614
 
615
      stop
616
 
617
 990  print*,'First columns must be <time,lon,lat,p>... Stop'
618
      stop
619
 
620
 
621
      end
622
 
623
 
624
 
625