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

Subversion Repositories lagranto.um

(root)/tags/1.2/goodies/mergetra.f – Rev 3 → 6