Subversion Repositories lagranto.icon

Rev

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

Rev Author Line No. Line
3 michaesp 1
      PROGRAM mege
2
 
3
c     ***********************************************************************
4
c     * Merge two trajectory files                                          *
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     Merging mode
23
      integer                                mode                             
24
 
25
c     Input and output format for trajectories (see iotra.f)
26
      character*80                           inpfile1       ! Input filename 1
27
      character*80                           inpfile2       ! Input filename 2
28
      character*80                           outfile        ! Output filename
29
      integer                                inpmode1       ! Input format 1
30
      integer                                inpmode2       ! Input format 2
31
      integer                                outmode        ! Output format
32
 
33
c     Trajectories
34
      integer                                ntra_inp1      ! Number of trajectories
35
      integer                                ntim_inp1      ! Number of times
36
      integer                                ncol_inp1      ! Number of columns
37
      real,allocatable, dimension (:,:,:) :: tra_inp1       ! Trajectories (ntra,ntim,ncol)
38
      character*80                           vars_inp1(100) ! Variable names
39
      real                                   time_inp1(500) ! Times of input trajectory
40
      integer                                refdate1(6)    ! Reference date
41
 
42
      integer                                ntra_inp2      ! Number of trajectories
43
      integer                                ntim_inp2      ! Number of times
44
      integer                                ncol_inp2      ! Number of columns
45
      real,allocatable, dimension (:,:,:) :: tra_inp2       ! Trajectories (ntra,ntim,ncol)
46
      character*80                           vars_inp2(100) ! Variable names
47
      real                                   time_inp2(500) ! Times of input trajectory
48
      integer                                refdate2(6)    ! Reference date
49
 
50
      integer                                ntra_out      ! Number of trajectories
51
      integer                                ntim_out      ! Number of times
52
      integer                                ncol_out      ! Number of columns
53
      real,allocatable, dimension (:,:,:) :: tra_out       ! Trajectories (ntra,ntim,ncol)
54
      character*80                           vars_out(100) ! Variable names
55
 
56
      real                                   time_out(500) ! Times of output trajectory
57
      integer                                refdate(6)    ! Reference date
58
 
59
c     Auxiliary variables
60
      integer                                i,j,k
61
      integer                                ind(200)
62
      integer                                itr(200)
63
      integer                                isok
64
      integer                                stat
65
      integer                                fid
66
      real                                   rswap
67
      integer                                iswap
68
 
69
c     ----------------------------------------------------------------------
70
c     Read and handle parameters
71
c     ----------------------------------------------------------------------
72
 
73
c     Read parameters
74
      open(10,file='mergetra.param')
75
       read(10,*) inpfile1
76
       read(10,*) inpfile2
77
       read(10,*) outfile
78
       read(10,*) ntra_inp1,ntim_inp1,ncol_inp1
79
       read(10,*) ntra_inp2,ntim_inp2,ncol_inp2
80
      close(10)
81
 
82
c     Determine the formats
83
      call mode_tra(inpmode1,inpfile1)
84
      if (inpmode1.eq.-1) inpmode1=1
85
      call mode_tra(inpmode2,inpfile2)
86
      if (inpmode2.eq.-1) inpmode2=1
87
      call mode_tra(outmode,outfile)
88
      if (outmode.eq.-1) outmode=1
89
 
90
c     ----------------------------------------------------------------------
91
c     Read input trajectories
92
c     ----------------------------------------------------------------------
93
 
94
c     Allocate memory
95
      allocate(tra_inp1(ntra_inp1,ntim_inp1,ncol_inp1),stat=stat)
96
      if (stat.ne.0) print*,'*** error allocating array tra_inp1    ***' 
97
      allocate(tra_inp2(ntra_inp2,ntim_inp2,ncol_inp2),stat=stat)
98
      if (stat.ne.0) print*,'*** error allocating array tra_inp2    ***' 
99
 
100
c     Read inpufile
101
      call ropen_tra(fid,inpfile1,ntra_inp1,ntim_inp1,ncol_inp1,
102
     >               refdate1,vars_inp1,inpmode1)
103
      call read_tra (fid,tra_inp1,ntra_inp1,ntim_inp1,ncol_inp1,
104
     >               inpmode1)
105
      call close_tra(fid,inpmode1)
106
 
107
      call ropen_tra(fid,inpfile2,ntra_inp2,ntim_inp2,ncol_inp2,
108
     >               refdate2,vars_inp2,inpmode2)
109
      call read_tra (fid,tra_inp2,ntra_inp2,ntim_inp2,ncol_inp2,
110
     >               inpmode2)
111
      call close_tra(fid,inpmode2)
112
 
113
c     Get the times of the trajectories
114
      do i=1,ntim_inp1
115
         time_inp1(i) = tra_inp1(1,i,1)
116
      enddo
117
      do i=1,ntim_inp2
118
         time_inp2(i) = tra_inp2(1,i,1)
119
      enddo
120
 
121
c     Check format 
122
      if (  vars_inp1(1).ne.'time' ) goto 990
123
      if ( (vars_inp1(2).ne.'lon').and.
124
     >     (vars_inp1(2).ne.'xpos') ) goto 990
125
      if ( (vars_inp1(3).ne.'lat').and.
126
     >     (vars_inp1(3).ne.'ypos') ) goto 990
127
      if ( (vars_inp1(4).ne.'z'  ).and.
128
     >     (vars_inp1(4).ne.'zpos') ) goto 990
129
 
130
      if (  vars_inp2(1).ne.'time' ) goto 990
131
      if ( (vars_inp2(2).ne.'lon').and.
132
     >     (vars_inp2(2).ne.'xpos') ) goto 990
133
      if ( (vars_inp2(3).ne.'lat').and.
134
     >     (vars_inp2(3).ne.'ypos') ) goto 990
135
      if ( (vars_inp2(4).ne.'z'  ).and.
136
     >     (vars_inp2(4).ne.'zpos') ) goto 990
137
 
138
c     Check whether the reference dates are equal
139
      do i=1,5
140
         if ( refdate1(i).ne.refdate2(i) ) then
141
            print*,' ERROR; Reference dates must be the same... Stop'
142
            print*,(refdate1(j),j=1,5)
143
            print*,(refdate2(j),j=1,5)
144
            stop
145
         endif
146
      enddo
147
 
148
c     ----------------------------------------------------------------------
149
c     Decide what to do
150
c     ----------------------------------------------------------------------
151
 
152
c     Init the mode
153
      mode = 0
154
 
155
c     ---- Check whether only the columns should be combined (mode 1) -------
156
      print*,'Testing for mode 1 (combine columns)'
157
 
158
      if ( ntim_inp1.ne.ntim_inp2 ) goto 100
159
      print*,'  -> ntim       ok'
160
 
161
      if ( ntra_inp1.ne.ntra_inp2 ) goto 100
162
      print*,'  -> ntra       ok'
163
 
164
      do i=1,ntim_inp1
165
         if ( time_inp1(i).ne.time_inp2(i) ) goto 100
166
      enddo
167
      print*,'  -> times      ok'
168
 
169
      do i=1,ntra_inp1
170
         do j=1,ntim_inp1
171
            if ( tra_inp1(i,j,1).ne.tra_inp2(i,j,1) ) goto 100
172
            if ( tra_inp1(i,j,2).ne.tra_inp2(i,j,2) ) goto 100
173
            if ( tra_inp1(i,j,3).ne.tra_inp2(i,j,3) ) goto 100
174
            if ( tra_inp1(i,j,4).ne.tra_inp2(i,j,4) ) goto 100
175
         enddo
176
      enddo
177
      print*,'  -> lon,lat,p  ok     => Mode 1 accepted'
178
 
179
      mode = 1
180
      goto 130
181
 
182
 
183
 100  continue
184
 
185
c     ---- Check whether second file to appended to first one (mode 2) -----
186
      print*,'Testing for mode 2 (append file 2 to file 1)'
187
 
188
      if ( ntim_inp1.ne.ntim_inp2 ) goto 110
189
      print*,'  -> ntim       ok'
190
 
191
      if ( ncol_inp1.ne.ncol_inp2 ) goto 110
192
      print*,'  -> ncol       ok'
193
 
194
      do i=1,ntim_inp1
195
         if ( time_inp1(i).ne.time_inp2(i) ) goto 110
196
      enddo
197
      print*,'  -> times      ok'
198
 
199
      do i=1,ncol_inp1
200
         if ( vars_inp1(i).ne.vars_inp2(i) ) goto 110
201
      enddo
202
      print*,'  -> vars       ok    => Mode 2 accepted'
203
 
204
      mode = 2
205
      goto 130
206
 
207
 110  continue
208
 
209
c     ----- Check whether to combine different times (mode 3) --------------
210
      print*,'Testing for mode 3 (combining times)'
211
 
212
      if ( ntra_inp1.ne.ntra_inp2 ) goto 120
213
      print*,'  -> ntra       ok'
214
 
215
      if ( ncol_inp1.ne.ncol_inp2 ) goto 120
216
      print*,'  -> ncol       ok'
217
 
218
      do i=1,ncol_inp1
219
         if ( vars_inp1(i).ne.vars_inp2(i) ) goto 120
220
      enddo
221
      print*,'  -> vars       ok    => Mode 3 accepted'
222
 
223
      mode = 3
224
      goto 130
225
 
226
 120  continue
227
 
228
c     -----  Stop if no valid merging mode could be determined -------------
229
      print*,' ERROR: could not determine merging mode ...'
230
      stop
231
 
232
 
233
c     Exit point for mode decision
234
 130  continue       
235
 
236
c     ----------------------------------------------------------------------
237
c     Merging depending on mode
238
c     ----------------------------------------------------------------------
239
 
240
c     ------ Merging of columns (mode 1) -----------------------------------
241
      if ( mode.ne.1 ) goto 200
242
 
243
c     Get the complete list of additional field
244
      do i=1,ncol_inp1
245
         ncol_out           = ncol_out + 1
246
         vars_out(ncol_out) = vars_inp1(i)
247
         ind(ncol_out)      = i 
248
         itr(ncol_out)      = 1
249
      enddo
250
      do i=5,ncol_inp2
251
 
252
         isok=1
253
         do j=1,ncol_out
254
            if ( vars_inp2(i).eq.vars_out(j) ) isok=0
255
         enddo
256
 
257
         if ( isok.eq.1 ) then
258
            ncol_out           = ncol_out + 1
259
            vars_out(ncol_out) = vars_inp2(i)
260
            ind(ncol_out)      = i 
261
            itr(ncol_out)      = 2
262
         endif
263
 
264
      enddo      
265
 
266
c     Allocate memory for output trajectory
267
      ntra_out = ntra_inp1
268
      ntim_out = ntim_inp1
269
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
270
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
271
 
272
c     Save the trajectory
273
      do i=1,ntra_out
274
         do j=1,ntim_out
275
            do k=1,ncol_out
276
 
277
               if ( itr(k).eq.1) then
278
                  tra_out(i,j,k) = tra_inp1(i,j,ind(k))
279
               else
280
                  tra_out(i,j,k) = tra_inp2(i,j,ind(k))
281
               endif
282
 
283
            enddo
284
         enddo
285
      enddo   
286
 
287
 200  continue
288
 
289
c     ------ Appending files (mode 2) --------------------------------------
290
      if ( mode.ne.2 ) goto 210
291
 
292
c     Allocate memory for output trajectory
293
      ntra_out = ntra_inp1 + ntra_inp2
294
      ntim_out = ntim_inp1
295
      ncol_out = ncol_inp1
296
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
297
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
298
 
299
c     Set the field names for the output trajectory
300
      do i=1,ncol_out
301
         vars_out(i) = vars_inp1(i)
302
      enddo
303
 
304
c     Save the trajectory
305
      do i=1,ntra_out
306
         do j=1,ntim_out
307
            do k=1,ncol_out
308
 
309
               if ( i.le.ntra_inp1 ) then
310
                  tra_out(i,j,k) = tra_inp1(i          ,j,k)
311
               else
312
                  tra_out(i,j,k) = tra_inp2(i-ntra_inp1,j,k)
313
               endif
314
 
315
            enddo
316
         enddo
317
      enddo
318
 
319
 210  continue
320
 
321
c     ------ Combining times (mode 3) --------------------------------------
322
      if ( mode.ne.3 ) goto 220
323
 
324
c     Get a list of all output times
325
      ntim_out = 0
326
      do i=1,ntim_inp1
327
         isok = 1
328
         do j=1,ntim_out
329
            if ( time_inp1(i).eq.time_out(j) ) isok=0
330
         enddo
331
         if (isok.eq.1 ) then
332
            ntim_out           = ntim_out + 1
333
            time_out(ntim_out) = time_inp1(i)
334
            itr(ntim_out)      = 1
335
            ind(ntim_out)      = i
336
         endif
337
      enddo
338
      do i=1,ntim_inp2
339
         isok = 1
340
         do j=1,ntim_out
341
            if ( time_inp2(i).eq.time_out(j) ) isok=0
342
         enddo
343
         if (isok.eq.1 ) then
344
            ntim_out           = ntim_out + 1
345
            time_out(ntim_out) = time_inp2(i)
346
            itr(ntim_out)      = 2
347
            ind(ntim_out)      = i
348
         endif
349
      enddo
350
 
351
c     Sort the times
352
      do i=1,ntim_out
353
         do j=i+1,ntim_out
354
            if ( time_out(j).lt.time_out(i) ) then
355
 
356
               rswap       = time_out(i)
357
               time_out(i) = time_out(j)
358
               time_out(j) = rswap
359
 
360
               iswap       = itr(i)
361
               itr(i)      = itr(j)
362
               itr(j)      = iswap
363
 
364
               iswap       = ind(i)
365
               ind(i)      = ind(j)
366
               ind(j)      = iswap
367
 
368
            endif
369
         enddo
370
      enddo
371
 
372
c     Allocate memory for output trajectory
373
      ntra_out = ntra_inp1 
374
      ncol_out = ncol_inp1
375
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
376
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
377
 
378
c     Set the field names for the output trajectory
379
      do i=1,ncol_out
380
         vars_out(i) = vars_inp1(i)
381
      enddo
382
 
383
c     Save the trajectory
384
      do i=1,ntra_out
385
         do j=1,ntim_out
386
            do k=1,ncol_out
387
 
388
               if ( itr(j).eq.1 ) then
389
                  tra_out(i,j,k) = tra_inp1(i,ind(j),k)
390
               else
391
                  tra_out(i,j,k) = tra_inp2(i,ind(j),k)
392
               endif
393
 
394
            enddo
395
         enddo
396
      enddo
397
 
398
 
399
 
400
 220  continue
401
 
402
 
403
c     ----------------------------------------------------------------------
404
c     Write the output trajectory
405
c     ----------------------------------------------------------------------
406
 
407
      call wopen_tra(fid,outfile,ntra_out,ntim_out,ncol_out,
408
     >               refdate1,vars_out,outmode)
409
      call write_tra(fid,tra_out,ntra_out,ntim_out,ncol_out,outmode)
410
      call close_tra(fid,outmode)
411
 
412
c     ----------------------------------------------------------------------
413
c     Error handling
414
c     ----------------------------------------------------------------------
415
 
416
      stop
417
 
418
 990  print*,'First columns must be <time,lon,lat,z>... Stop'
419
      stop
420
 
421
 
422
      end
423
 
424
 
425
 
426
 
427
 
428
 
429