Subversion Repositories lagranto.um

Rev

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