WebSVN – lagranto.ecmwf – Diff – /trunk/goodies/mergetra.f


      PROGRAM mege
      
c     ***********************************************************************
c     * Merge two trajectory files                                          *
c     * Michael Sprenger / Spring, summer 2010                              *
c     ***********************************************************************
 
      implicit none
 
c     ----------------------------------------------------------------------
c     Declaration of parameters
c     ----------------------------------------------------------------------
 
c     Numerical epsilon
      real      eps
      parameter (eps=0.0001)
 
c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------
 
c     Merging mode
      integer                                mode
      character*20                           datecheck
 
c     Input and output format for trajectories (see iotra.f)
      character*80                           inpfile1       ! Input filename 1
      character*80                           inpfile2       ! Input filename 2
      character*80                           outfile        ! Output filename
      integer                                inpmode1       ! Input format 1
      integer                                inpmode2       ! Input format 2
      integer                                outmode        ! Output format
      
c     Trajectories
      integer                                ntra_inp1       ! Number of trajectories
      integer                                ntim_inp1       ! Number of times
      integer                                ncol_inp1       ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra_inp1        ! Trajectories (ntra,ntim,ncol)
      character*80                           vars_inp1(100)  ! Variable names
      real                                   time_inp1(10000) ! Times of input trajectory
      integer                                refdate1(6)     ! Reference date
 
      integer                                ntra_inp2       ! Number of trajectories
      integer                                ntim_inp2       ! Number of times
      integer                                ncol_inp2       ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra_inp2        ! Trajectories (ntra,ntim,ncol)
      character*80                           vars_inp2(100)  ! Variable names
      real                                   time_inp2(10000) ! Times of input trajectory
      integer                                refdate2(6)     ! Reference date
 
      integer                                ntra_out       ! Number of trajectories
      integer                                ntim_out       ! Number of times
      integer                                ncol_out       ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra_out        ! Trajectories (ntra,ntim,ncol)
      character*80                           vars_out(100)  ! Variable names
  
      real                                   time_out(20000)! Times of output trajectory
      integer                                refdate(6)     ! Reference date
 
c     Auxiliary variables
      integer                                i,j,k
      integer                                ind(20000)
      integer                                itr(20000)
      integer                                isok
      integer                                stat
      integer                                fid
      real                                   rswap
      integer                                iswap
 
c     ----------------------------------------------------------------------
c     Read and handle parameters
c     ----------------------------------------------------------------------
 
c     Read parameters
      open(10,file='mergetra.param')
       read(10,*) inpfile1
       read(10,*) inpfile2
       read(10,*) outfile
       read(10,*) ntra_inp1,ntim_inp1,ncol_inp1
       read(10,*) ntra_inp2,ntim_inp2,ncol_inp2
       read(10,*) datecheck
      close(10)
 
c     Determine the formats
      call mode_tra(inpmode1,inpfile1)
      if (inpmode1.eq.-1) inpmode1=1
      call mode_tra(inpmode2,inpfile2)
      if (inpmode2.eq.-1) inpmode2=1
      call mode_tra(outmode,outfile)
      if (outmode.eq.-1) outmode=1
 
c     ----------------------------------------------------------------------
c     Read input trajectories
c     ----------------------------------------------------------------------
 
c     Allocate memory
      allocate(tra_inp1(ntra_inp1,ntim_inp1,ncol_inp1),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_inp1    ***' 
      allocate(tra_inp2(ntra_inp2,ntim_inp2,ncol_inp2),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_inp2    ***' 
 
c     Read inpufile
      call ropen_tra(fid,inpfile1,ntra_inp1,ntim_inp1,ncol_inp1,
     >               refdate1,vars_inp1,inpmode1)
      call read_tra (fid,tra_inp1,ntra_inp1,ntim_inp1,ncol_inp1,
     >               inpmode1)
      call close_tra(fid,inpmode1)
 
      call ropen_tra(fid,inpfile2,ntra_inp2,ntim_inp2,ncol_inp2,
     >               refdate2,vars_inp2,inpmode2)
      call read_tra (fid,tra_inp2,ntra_inp2,ntim_inp2,ncol_inp2,
     >               inpmode2)
      call close_tra(fid,inpmode2)
 
c     Get the times of the trajectories
      do i=1,ntim_inp1
         time_inp1(i) = tra_inp1(1,i,1)
      enddo
      do i=1,ntim_inp2
         time_inp2(i) = tra_inp2(1,i,1)
      enddo
 
c     Check format 
      if (  vars_inp1(1).ne.'time' ) goto 990
      if ( (vars_inp1(2).ne.'lon').and.
     >     (vars_inp1(2).ne.'xpos') ) goto 990
      if ( (vars_inp1(3).ne.'lat').and.
     >     (vars_inp1(3).ne.'ypos') ) goto 990
      if ( (vars_inp1(4).ne.'p'  ).and.
     >     (vars_inp1(4).ne.'ppos') ) goto 990
 
      if (  vars_inp2(1).ne.'time' ) goto 990
      if ( (vars_inp2(2).ne.'lon').and.
     >     (vars_inp2(2).ne.'xpos') ) goto 990
      if ( (vars_inp2(3).ne.'lat').and.
     >     (vars_inp2(3).ne.'ypos') ) goto 990
      if ( (vars_inp2(4).ne.'p'  ).and.
     >     (vars_inp2(4).ne.'ppos') ) goto 990
 
c     Check whether the reference dates are equal
      if ( datecheck.eq.'datecheck' ) then
        do i=1,5
          if ( refdate1(i).ne.refdate2(i) ) then
            print*,' ERROR; Reference dates must be the same... Stop'
            print*,(refdate1(j),j=1,5)
            print*,(refdate2(j),j=1,5)
            stop
          endif
        enddo
      endif
            
c     ----------------------------------------------------------------------
c     Decide what to do
c     ----------------------------------------------------------------------
 
c     Init the mode
      mode = 0
 
c     ---- Check whether only the columns should be combined (mode 1) -------
      print*,'Testing for mode 1 (combine columns)'
 
      if ( ntim_inp1.ne.ntim_inp2 ) goto 100
      print*,'  -> ntim       ok'
 
      if ( ntra_inp1.ne.ntra_inp2 ) goto 100
      print*,'  -> ntra       ok'
 
      do i=1,ntim_inp1
         if ( time_inp1(i).ne.time_inp2(i) ) goto 100
      enddo
      print*,'  -> times      ok'
 
      do i=1,ntra_inp1
         do j=1,ntim_inp1
            if ( tra_inp1(i,j,1).ne.tra_inp2(i,j,1) ) goto 100
            if ( tra_inp1(i,j,2).ne.tra_inp2(i,j,2) ) goto 100
            if ( tra_inp1(i,j,3).ne.tra_inp2(i,j,3) ) goto 100
            if ( tra_inp1(i,j,4).ne.tra_inp2(i,j,4) ) goto 100
         enddo
      enddo
      print*,'  -> lon,lat,p  ok     => Mode 1 accepted'
      
      mode = 1
      goto 130
 
 
 100  continue
 
c     ---- Check whether second file to appended to first one (mode 2) -----
      print*,'Testing for mode 2 (append file 2 to file 1)'
      
      if ( ntim_inp1.ne.ntim_inp2 ) goto 110
      print*,'  -> ntim       ok'
 
      if ( ncol_inp1.ne.ncol_inp2 ) goto 110
      print*,'  -> ncol       ok'
 
      do i=1,ntim_inp1
         if ( time_inp1(i).ne.time_inp2(i) ) goto 110
      enddo
      print*,'  -> times      ok'
 
      do i=1,ncol_inp1
         if ( vars_inp1(i).ne.vars_inp2(i) ) goto 110
      enddo
      print*,'  -> vars       ok    => Mode 2 accepted'
      
      mode = 2
      goto 130
 
 110  continue
 
c     ----- Check whether to combine different times (mode 3) --------------
      print*,'Testing for mode 3 (combining times)'
      
      if ( ntra_inp1.ne.ntra_inp2 ) goto 120
      print*,'  -> ntra       ok'
 
      if ( ncol_inp1.ne.ncol_inp2 ) goto 120
      print*,'  -> ncol       ok'
 
      do i=1,ncol_inp1
         if ( vars_inp1(i).ne.vars_inp2(i) ) goto 120
      enddo
      print*,'  -> vars       ok    => Mode 3 accepted'
 
      mode = 3
      goto 130
      
 120  continue
 
c     -----  Stop if no valid merging mode could be determined -------------
      print*,' ERROR: could not determine merging mode ...'
      stop
 
 
c     Exit point for mode decision
 130  continue       
 
c     ----------------------------------------------------------------------
c     Merging depending on mode
c     ----------------------------------------------------------------------
 
c     ------ Merging of columns (mode 1) -----------------------------------
      if ( mode.ne.1 ) goto 200
 
c     Get the complete list of additional field
      do i=1,ncol_inp1
         ncol_out           = ncol_out + 1
         vars_out(ncol_out) = vars_inp1(i)
         ind(ncol_out)      = i 
         itr(ncol_out)      = 1
      enddo
      do i=5,ncol_inp2
         
         isok=1
         do j=1,ncol_out
            if ( vars_inp2(i).eq.vars_out(j) ) isok=0
         enddo
         
         if ( isok.eq.1 ) then
            ncol_out           = ncol_out + 1
            vars_out(ncol_out) = vars_inp2(i)
            ind(ncol_out)      = i 
            itr(ncol_out)      = 2
         endif
 
      enddo      
 
c     Allocate memory for output trajectory
      ntra_out = ntra_inp1
      ntim_out = ntim_inp1
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Save the trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               
               if ( itr(k).eq.1) then
                  tra_out(i,j,k) = tra_inp1(i,j,ind(k))
               else
                  tra_out(i,j,k) = tra_inp2(i,j,ind(k))
               endif
 
            enddo
         enddo
      enddo   
 
 200  continue
 
c     ------ Appending files (mode 2) --------------------------------------
      if ( mode.ne.2 ) goto 210
 
c     Allocate memory for output trajectory
      ntra_out = ntra_inp1 + ntra_inp2
      ntim_out = ntim_inp1
      ncol_out = ncol_inp1
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Set the field names for the output trajectory
      do i=1,ncol_out
         vars_out(i) = vars_inp1(i)
      enddo
 
c     Save the trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               
               if ( i.le.ntra_inp1 ) then
                  tra_out(i,j,k) = tra_inp1(i          ,j,k)
               else
                  tra_out(i,j,k) = tra_inp2(i-ntra_inp1,j,k)
               endif
 
            enddo
         enddo
      enddo
 
 210  continue
 
c     ------ Combining times (mode 3) --------------------------------------
      if ( mode.ne.3 ) goto 220
 
c     Get a list of all output times
      ntim_out = 0
      do i=1,ntim_inp1
         isok = 1
         do j=1,ntim_out
            if ( time_inp1(i).eq.time_out(j) ) isok=0
         enddo
         if (isok.eq.1 ) then
            ntim_out           = ntim_out + 1
            time_out(ntim_out) = time_inp1(i)
            itr(ntim_out)      = 1
            ind(ntim_out)      = i
         endif
      enddo
      do i=1,ntim_inp2
         isok = 1
         do j=1,ntim_out
            if ( time_inp2(i).eq.time_out(j) ) isok=0
         enddo
         if (isok.eq.1 ) then
            ntim_out           = ntim_out + 1
            time_out(ntim_out) = time_inp2(i)
            itr(ntim_out)      = 2
            ind(ntim_out)      = i
         endif
      enddo
 
c     Sort the times
      do i=1,ntim_out
         do j=i+1,ntim_out
            if ( time_out(j).lt.time_out(i) ) then
 
               rswap       = time_out(i)
               time_out(i) = time_out(j)
               time_out(j) = rswap
 
               iswap       = itr(i)
               itr(i)      = itr(j)
               itr(j)      = iswap
 
               iswap       = ind(i)
               ind(i)      = ind(j)
               ind(j)      = iswap
 
            endif
         enddo
      enddo
 
c     Allocate memory for output trajectory
      ntra_out = ntra_inp1 
      ncol_out = ncol_inp1
      allocate(tra_out(ntra_out,ntim_out,ncol_out),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra_out    ***' 
 
c     Set the field names for the output trajectory
      do i=1,ncol_out
         vars_out(i) = vars_inp1(i)
      enddo
 
c     Save the trajectory
      do i=1,ntra_out
         do j=1,ntim_out
            do k=1,ncol_out
               
               if ( itr(j).eq.1 ) then
                  tra_out(i,j,k) = tra_inp1(i,ind(j),k)
               else
                  tra_out(i,j,k) = tra_inp2(i,ind(j),k)
               endif
 
            enddo
         enddo
      enddo
 
 
 
 220  continue
 
 
c     ----------------------------------------------------------------------
c     Write the output trajectory
c     ----------------------------------------------------------------------
 
      call wopen_tra(fid,outfile,ntra_out,ntim_out,ncol_out,
     >               refdate1,vars_out,outmode)
      call write_tra(fid,tra_out,ntra_out,ntim_out,ncol_out,outmode)
      call close_tra(fid,outmode)
 
c     ----------------------------------------------------------------------
c     Error handling
c     ----------------------------------------------------------------------
      
      stop
   
 990  print*,'First columns must be <time,lon,lat,p>... Stop'
      stop
      
 
      end
 
 
 
 
 
 
 

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

Subversion Repositories lagranto.ecmwf

(root)/trunk/goodies/mergetra.f – Rev 27 → 46