Subversion Repositories lagranto.wrf

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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