Subversion Repositories lagranto.ecmwf

Rev

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