Subversion Repositories lagranto.ecmwf

Rev

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