WebSVN – lagranto.um – Diff – /tags/2.0/goodies/difference.f


 
      PROGRAM difference
      
c     ***********************************************************************
c     * Get the difference between tow trajectory files                     *
c     * Michael Sprenger / Spring, summer, autumn 2010                      *
c     ***********************************************************************
 
      implicit none
      
c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------
 
c     Field name and mode for difference calculation
      character*80                           mode                     ! Difference mode
      character*80                           fieldname                ! Name of differencing
 
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
 
c     Input trajectories
      integer                                ntra1      ,ntra2        ! Number of trajectories
      integer                                ntim1      ,ntim2        ! Number of times
      integer                                ncol1      ,ncol2        ! Number of columns
      real,allocatable, dimension (:,:,:) :: trainp1    ,trainp2      ! Trajectories (ntra,ntim,ncol)
      character*80                           vars1(100) ,vars2(100)   ! Variable names
      integer                                refdate1(6),refdate2(6)  ! Reference date
      
c     Output/comparison trajectory
      integer                                ntra                    ! Number of trajectories
      integer                                ntim                    ! Number of times
      integer                                ncol                    ! Number of columns
      real,allocatable, dimension (:,:,:) :: traout                  ! Trajectories (ntra,ntim,ncol)
      character*80                           vars(100)               ! Variable names
      integer                                refdate(6)              ! Reference date
 
c     Auxiliary variables
      integer                                inpmode1,inpmode2
      integer                                outmode
      integer                                stat
      integer                                fid
      integer                                i,j,k
      integer                                ind1,ind2
      integer                                isok
      character                              ch
      real,allocatable, dimension (:) ::     diff
      integer                                outind
 
c     Externals 
      real,external                   ::     sdis
      
c     ----------------------------------------------------------------------
c     Preparations
c     ----------------------------------------------------------------------
 
c     Read parameters
      open(10,file='difference.param')
       read(10,*) inpfile1
       read(10,*) inpfile2
       read(10,*) outfile
       read(10,*) ntra1,ntim1,ncol1
       read(10,*) ntra2,ntim2,ncol2
       read(10,*) mode
       read(10,*) fieldname
      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     Set number of trajectories for output
      if ( ntra1.lt.ntra2) then
         ntra = ntra1
      else
         ntra = ntra2
      endif
 
c     Set number of times for output
      if ( mode.eq.'single' ) then
         ntim = ntim1
      else
         ntim = 1
      endif
 
c     Set the column names for output 
      if ( fieldname.eq.'LATLON') then         
         ncol    = 1 + 3 + 3     + 1
         vars(1) = 'time'
         vars(2) = 'lon[1]'
         vars(3) = 'lat[1]'
         vars(4) = 'p[1]'
         vars(5) = 'lon[2]'
         vars(6) = 'lat[2]'
         vars(7) = 'p[2]'
         vars(8) = 'SDIS'
      else
         ncol = 1 + 3 + 3 + 2 + 1
         vars( 1) = 'time'
         vars( 2) = 'lon[1]'
         vars( 3) = 'lat[1]'
         vars( 4) = 'p[1]'
         vars( 5) = 'lon[2]'
         vars( 6) = 'lat[2]'
         vars( 7) = 'p[2]'
         vars( 8) = trim(fieldname)//'[1]'
         vars( 9) = trim(fieldname)//'[2]'
         vars(10) = trim(fieldname)//'[1-2]'
      endif
 
c     Allocate memory
      allocate(trainp1(ntra1,ntim1,ncol1),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp1   ***' 
      allocate(trainp2(ntra2,ntim2,ncol2),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp2   ***' 
      allocate(traout(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array traout    ***' 
      allocate(diff(ntim1),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array diff      ***' 
 
c     Read inpufiles
      call ropen_tra(fid,inpfile1,ntra1,ntim1,ncol1,
     >                   refdate1,vars1,inpmode1)
      call read_tra (fid,trainp1,ntra1,ntim1,ncol1,inpmode1)
      call close_tra(fid,inpmode1)
 
      call ropen_tra(fid,inpfile2,ntra2,ntim2,ncol2,
     >                   refdate2,vars2,inpmode2)
      call read_tra (fid,trainp2,ntra2,ntim2,ncol2,inpmode2)
      call close_tra(fid,inpmode2)
 
c     Check dimensions of the two trajectory files (#tim hard check)
      if (ntim1.ne.ntim2) then
         print*,'Trajectoy files have different time dimensions... Stop'
         stop
      endif
 
c     Check dimensions of the two trajectory files (#tra soft check)
      if (ntra1.ne.ntra2) then
         print*,'Differing number of trajectories... proceed [y/n]'
         read*,ch
         if (ch.eq.'n') stop
      endif
 
c     Check whether difference field is available on both files 
      if ( fieldname.ne.'LATLON') then
         ind1 = 0
         ind2 = 0
         do i=1,ncol
            if ( fieldname.eq.vars1(i) ) ind1 = i
            if ( fieldname.eq.vars2(i) ) ind2 = i
         enddo
         if ( (ind1.eq.0).or.(ind2.eq.0) ) then
            print*,'Field ',trim(fieldname),' not available... Stop'
            stop
         endif
      endif
 
c     Check reference dates (soft check)
      isok = 1
      do i=1,6
         if ( refdate1(i).ne.refdate2(i) ) isok = 0
      enddo
      if ( isok.eq.0 ) then
         print*,'Warning: reference dates differ... proceed [y/n]'
         read*,ch
         if (ch.eq.'n') stop
      endif
 
c     Check trajectory times (soft check)
      isok = 1
      do i=1,ntim
         if ( trainp1(1,i,1).ne.trainp2(1,i,1) ) isok = 0
      enddo
      if ( isok.eq.0 ) then
         print*,'Warning: trajectory times differ... proceed [y/n]'
         read*,ch
         if (ch.eq.'n') stop
      endif
 
c     Copy reference date to output
      do i=1,6
         refdate(i) = refdate1(i)
      enddo
      
c     ----------------------------------------------------------------------
c     Calculate the difference (depending on mode)
c     ----------------------------------------------------------------------
 
c     Loop over all trajectories
      do i=1,ntra
 
c        Calculate difference for all times
         do j=1,ntim1
            
c           Calculate the difference (distance or absolute value)
            if (fieldname.eq.'LATLON') then
               diff(j) = sdis( trainp1(i,j,2),trainp1(i,j,3),
     >                         trainp2(i,j,2),trainp2(i,j,3) )         
            else
               diff(j) = abs(trainp1(i,j,ind1) - trainp2(i,j,ind2))
            endif
 
         enddo
 
c        Save output for each time
         if ( mode.eq.'single' ) then
 
            do j=1,ntim
            
               if ( fieldname.eq.'LATLON' ) then
                  traout(i,j, 1) = trainp1(i,j,1)    ! time
                  traout(i,j, 2) = trainp1(i,j,2)    ! lon[1]
                  traout(i,j, 3) = trainp1(i,j,3)    ! lat[1]
                  traout(i,j, 4) = trainp1(i,j,4)    ! p[1]
                  traout(i,j, 5) = trainp2(i,j,2)    ! lon[2]
                  traout(i,j, 6) = trainp2(i,j,3)    ! lat[2]
                  traout(i,j, 7) = trainp2(i,j,4)    ! p[2]
                  traout(i,j, 8) = diff(j)           ! SDIS(j)
               else
                  traout(i,j, 1) = trainp1(i,j,1)    ! time
                  traout(i,j, 2) = trainp1(i,j,2)    ! lon[1]
                  traout(i,j, 3) = trainp1(i,j,3)    ! lat[1]
                  traout(i,j, 4) = trainp1(i,j,4)    ! p[1]
                  traout(i,j, 5) = trainp2(i,j,2)    ! lon[2]
                  traout(i,j, 6) = trainp2(i,j,3)    ! lat[2]
                  traout(i,j, 7) = trainp2(i,j,4)    ! p[2]
                  traout(i,j, 8) = trainp1(i,j,ind1) ! field[1]
                  traout(i,j, 9) = trainp2(i,j,ind2) ! field[2]
                  traout(i,j,10) = diff(j)           ! SDIS(j)
               endif
 
            enddo
 
c        Save only maximum
         elseif ( mode.eq.'max') then
            
            outind  = 1
            do j=2,ntim1
               if ( diff(j).gt.diff(outind) ) outind = j
            enddo
 
            if ( fieldname.eq.'LATLON' ) then
               traout(i,1, 1) = trainp1(i,outind,1)     ! time
               traout(i,1, 2) = trainp1(i,outind,2)     ! lon[1]
               traout(i,1, 3) = trainp1(i,outind,3)     ! lat[1]
               traout(i,1, 4) = trainp1(i,outind,4)     ! p[1]
               traout(i,1, 5) = trainp2(i,outind,2)     ! lon[2]
               traout(i,1, 6) = trainp2(i,outind,3)     ! lat[2]
               traout(i,1, 7) = trainp2(i,outind,4)     ! p[2]
               traout(i,1, 8) = diff(outind)            ! SDIS
            else
               traout(i,1, 1) = trainp1(i,outind,1)     ! time
               traout(i,1, 2) = trainp1(i,outind,2)     ! lon[1]
               traout(i,1, 3) = trainp1(i,outind,3)     ! lat[1]
               traout(i,1, 4) = trainp1(i,outind,4)     ! p[1]
               traout(i,1, 5) = trainp2(i,outind,2)     ! lon[2]
               traout(i,1, 6) = trainp2(i,outind,3)     ! lat[2]
               traout(i,1, 7) = trainp2(i,outind,4)     ! p[2]
               traout(i,1, 8) = trainp1(i,outind,ind1)  ! field[1]
               traout(i,1, 9) = trainp2(i,outind,ind2)  ! field[2]
               traout(i,1,10) = diff(outind)            ! SDIS(j)
            endif
            
         endif
 
      enddo
    
c     ----------------------------------------------------------------------
c     Write output
c     ----------------------------------------------------------------------
 
c     Write output file
      call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)
      call write_tra(fid,traout,ntra,ntim,ncol,outmode)
      call close_tra(fid,outmode)
      
      end
 
 
 
 
c     ***********************************************************************
c     * Subroutines                                                         *
c     ***********************************************************************
 
c     ----------------------------------------------------------------------
c     Spherical distance
c     ----------------------------------------------------------------------
 
      real function sdis(xp,yp,xq,yq)
c
c     calculates spherical distance (in km) between two points given
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
c
      real      pi180
      parameter (pi180=3.14159/180.)
      real      re
      parameter (re=6370.)
      real      degkm
      parameter (degkm=111.1775)
      real      xp,yp,xq,yq,arg
 
      if ( (abs(xp-xq).gt.0.05).and.(abs(yp-yq).gt.0.05) ) then
         arg=sin(pi180*yp)*sin(pi180*yq)+
     >       cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
         if (arg.lt.-1.) arg=-1.
         if (arg.gt.1.) arg=1.
         sdis=re*acos(arg)
      else
         sdis= (yp-yq)**2 + ( (xp-xq) * cos( pi180*0.5*(yp+yq) ) )**2
         sdis = deg2km * sqrt(sdis)
      endif
 
      end
 
 
      

Rev 3	Rev 10
1		1
2	`PROGRAM difference`	2	`PROGRAM difference`
3		3
4	`c ***********************************************************************`	4	`c ***********************************************************************`
5	`c * Get the difference between tow trajectory files *`	5	`c * Get the difference between tow trajectory files *`
6	`c * Michael Sprenger / Spring, summer, autumn 2010 *`	6	`c * Michael Sprenger / Spring, summer, autumn 2010 *`
7	`c ***********************************************************************`	7	`c ***********************************************************************`
8		8
9	`implicit none`	9	`implicit none`
10		10
11	`c ----------------------------------------------------------------------`	11	`c ----------------------------------------------------------------------`
12	`c Declaration of variables`	12	`c Declaration of variables`
13	`c ----------------------------------------------------------------------`	13	`c ----------------------------------------------------------------------`
14		14
15	`c Field name and mode for difference calculation`	15	`c Field name and mode for difference calculation`
16	`character*80 mode ! Difference mode`	16	`character*80 mode ! Difference mode`
17	`character*80 fieldname ! Name of differencing`	17	`character*80 fieldname ! Name of differencing`
18		18
19	`c Input and output format for trajectories (see iotra.f)`	19	`c Input and output format for trajectories (see iotra.f)`
20	`character*80 inpfile1 ! Input filename 1`	20	`character*80 inpfile1 ! Input filename 1`
21	`character*80 inpfile2 ! Input filename 2`	21	`character*80 inpfile2 ! Input filename 2`
22	`character*80 outfile ! Output filename`	22	`character*80 outfile ! Output filename`
23		23
24	`c Input trajectories`	24	`c Input trajectories`
25	`integer ntra1 ,ntra2 ! Number of trajectories`	25	`integer ntra1 ,ntra2 ! Number of trajectories`
26	`integer ntim1 ,ntim2 ! Number of times`	26	`integer ntim1 ,ntim2 ! Number of times`
27	`integer ncol1 ,ncol2 ! Number of columns`	27	`integer ncol1 ,ncol2 ! Number of columns`
28	`real,allocatable, dimension (:,:,:) :: trainp1 ,trainp2 ! Trajectories (ntra,ntim,ncol)`	28	`real,allocatable, dimension (:,:,:) :: trainp1 ,trainp2 ! Trajectories (ntra,ntim,ncol)`
29	`character*80 vars1(100) ,vars2(100) ! Variable names`	29	`character*80 vars1(100) ,vars2(100) ! Variable names`
30	`integer refdate1(6),refdate2(6) ! Reference date`	30	`integer refdate1(6),refdate2(6) ! Reference date`
31		31
32	`c Output/comparison trajectory`	32	`c Output/comparison trajectory`
33	`integer ntra ! Number of trajectories`	33	`integer ntra ! Number of trajectories`
34	`integer ntim ! Number of times`	34	`integer ntim ! Number of times`
35	`integer ncol ! Number of columns`	35	`integer ncol ! Number of columns`
36	`real,allocatable, dimension (:,:,:) :: traout ! Trajectories (ntra,ntim,ncol)`	36	`real,allocatable, dimension (:,:,:) :: traout ! Trajectories (ntra,ntim,ncol)`
37	`character*80 vars(100) ! Variable names`	37	`character*80 vars(100) ! Variable names`
38	`integer refdate(6) ! Reference date`	38	`integer refdate(6) ! Reference date`
39		39
40	`c Auxiliary variables`	40	`c Auxiliary variables`
41	`integer inpmode1,inpmode2`	41	`integer inpmode1,inpmode2`
42	`integer outmode`	42	`integer outmode`
43	`integer stat`	43	`integer stat`
44	`integer fid`	44	`integer fid`
45	`integer i,j,k`	45	`integer i,j,k`
46	`integer ind1,ind2`	46	`integer ind1,ind2`
47	`integer isok`	47	`integer isok`
48	`character ch`	48	`character ch`
49	`real,allocatable, dimension (:) :: diff`	49	`real,allocatable, dimension (:) :: diff`
50	`integer outind`	50	`integer outind`
51		51
52	`c Externals`	52	`c Externals`
53	`real,external :: sdis`	53	`real,external :: sdis`
54		54
55	`c ----------------------------------------------------------------------`	55	`c ----------------------------------------------------------------------`
56	`c Preparations`	56	`c Preparations`
57	`c ----------------------------------------------------------------------`	57	`c ----------------------------------------------------------------------`
58		58
59	`c Read parameters`	59	`c Read parameters`
60	`open(10,file='difference.param')`	60	`open(10,file='difference.param')`
61	`read(10,*) inpfile1`	61	`read(10,*) inpfile1`
62	`read(10,*) inpfile2`	62	`read(10,*) inpfile2`
63	`read(10,*) outfile`	63	`read(10,*) outfile`
64	`read(10,*) ntra1,ntim1,ncol1`	64	`read(10,*) ntra1,ntim1,ncol1`
65	`read(10,*) ntra2,ntim2,ncol2`	65	`read(10,*) ntra2,ntim2,ncol2`
66	`read(10,*) mode`	66	`read(10,*) mode`
67	`read(10,*) fieldname`	67	`read(10,*) fieldname`
68	`close(10)`	68	`close(10)`
69		69
70	`c Determine the formats`	70	`c Determine the formats`
71	`call mode_tra(inpmode1,inpfile1)`	71	`call mode_tra(inpmode1,inpfile1)`
72	`if (inpmode1.eq.-1) inpmode1=1`	72	`if (inpmode1.eq.-1) inpmode1=1`
73		73
74	`call mode_tra(inpmode2,inpfile2)`	74	`call mode_tra(inpmode2,inpfile2)`
75	`if (inpmode2.eq.-1) inpmode2=1`	75	`if (inpmode2.eq.-1) inpmode2=1`
76		76
77	`call mode_tra(outmode,outfile)`	77	`call mode_tra(outmode,outfile)`
78	`if (outmode.eq.-1) outmode=1`	78	`if (outmode.eq.-1) outmode=1`
79		79
80	`c Set number of trajectories for output`	80	`c Set number of trajectories for output`
81	`if ( ntra1.lt.ntra2) then`	81	`if ( ntra1.lt.ntra2) then`
82	`ntra = ntra1`	82	`ntra = ntra1`
83	`else`	83	`else`
84	`ntra = ntra2`	84	`ntra = ntra2`
85	`endif`	85	`endif`
86		86
87	`c Set number of times for output`	87	`c Set number of times for output`
88	`if ( mode.eq.'single' ) then`	88	`if ( mode.eq.'single' ) then`
89	`ntim = ntim1`	89	`ntim = ntim1`
90	`else`	90	`else`
91	`ntim = 1`	91	`ntim = 1`
92	`endif`	92	`endif`
93		93
94	`c Set the column names for output`	94	`c Set the column names for output`
95	`if ( fieldname.eq.'LATLON') then`	95	`if ( fieldname.eq.'LATLON') then`
96	`ncol = 1 + 3 + 3 + 1`	96	`ncol = 1 + 3 + 3 + 1`
97	`vars(1) = 'time'`	97	`vars(1) = 'time'`
98	`vars(2) = 'lon[1]'`	98	`vars(2) = 'lon[1]'`
99	`vars(3) = 'lat[1]'`	99	`vars(3) = 'lat[1]'`
100	`vars(4) = 'p[1]'`	100	`vars(4) = 'p[1]'`
101	`vars(5) = 'lon[2]'`	101	`vars(5) = 'lon[2]'`
102	`vars(6) = 'lat[2]'`	102	`vars(6) = 'lat[2]'`
103	`vars(7) = 'p[2]'`	103	`vars(7) = 'p[2]'`
104	`vars(8) = 'SDIS'`	104	`vars(8) = 'SDIS'`
105	`else`	105	`else`
106	`ncol = 1 + 3 + 3 + 2 + 1`	106	`ncol = 1 + 3 + 3 + 2 + 1`
107	`vars( 1) = 'time'`	107	`vars( 1) = 'time'`
108	`vars( 2) = 'lon[1]'`	108	`vars( 2) = 'lon[1]'`
109	`vars( 3) = 'lat[1]'`	109	`vars( 3) = 'lat[1]'`
110	`vars( 4) = 'p[1]'`	110	`vars( 4) = 'p[1]'`
111	`vars( 5) = 'lon[2]'`	111	`vars( 5) = 'lon[2]'`
112	`vars( 6) = 'lat[2]'`	112	`vars( 6) = 'lat[2]'`
113	`vars( 7) = 'p[2]'`	113	`vars( 7) = 'p[2]'`
114	`vars( 8) = trim(fieldname)//'[1]'`	114	`vars( 8) = trim(fieldname)//'[1]'`
115	`vars( 9) = trim(fieldname)//'[2]'`	115	`vars( 9) = trim(fieldname)//'[2]'`
116	`vars(10) = trim(fieldname)//'[1-2]'`	116	`vars(10) = trim(fieldname)//'[1-2]'`
117	`endif`	117	`endif`
118		118
119	`c Allocate memory`	119	`c Allocate memory`
120	`allocate(trainp1(ntra1,ntim1,ncol1),stat=stat)`	120	`allocate(trainp1(ntra1,ntim1,ncol1),stat=stat)`
121	`if (stat.ne.0) print,' error allocating array trainp1 *'`	121	`if (stat.ne.0) print,' error allocating array trainp1 *'`
122	`allocate(trainp2(ntra2,ntim2,ncol2),stat=stat)`	122	`allocate(trainp2(ntra2,ntim2,ncol2),stat=stat)`
123	`if (stat.ne.0) print,' error allocating array trainp2 *'`	123	`if (stat.ne.0) print,' error allocating array trainp2 *'`
124	`allocate(traout(ntra,ntim,ncol),stat=stat)`	124	`allocate(traout(ntra,ntim,ncol),stat=stat)`
125	`if (stat.ne.0) print,' error allocating array traout *'`	125	`if (stat.ne.0) print,' error allocating array traout *'`
126	`allocate(diff(ntim1),stat=stat)`	126	`allocate(diff(ntim1),stat=stat)`
127	`if (stat.ne.0) print,' error allocating array diff *'`	127	`if (stat.ne.0) print,' error allocating array diff *'`
128		128
129	`c Read inpufiles`	129	`c Read inpufiles`
130	`call ropen_tra(fid,inpfile1,ntra1,ntim1,ncol1,`	130	`call ropen_tra(fid,inpfile1,ntra1,ntim1,ncol1,`
131	`> refdate1,vars1,inpmode1)`	131	`> refdate1,vars1,inpmode1)`
132	`call read_tra (fid,trainp1,ntra1,ntim1,ncol1,inpmode1)`	132	`call read_tra (fid,trainp1,ntra1,ntim1,ncol1,inpmode1)`
133	`call close_tra(fid,inpmode1)`	133	`call close_tra(fid,inpmode1)`
134		134
135	`call ropen_tra(fid,inpfile2,ntra2,ntim2,ncol2,`	135	`call ropen_tra(fid,inpfile2,ntra2,ntim2,ncol2,`
136	`> refdate2,vars2,inpmode2)`	136	`> refdate2,vars2,inpmode2)`
137	`call read_tra (fid,trainp2,ntra2,ntim2,ncol2,inpmode2)`	137	`call read_tra (fid,trainp2,ntra2,ntim2,ncol2,inpmode2)`
138	`call close_tra(fid,inpmode2)`	138	`call close_tra(fid,inpmode2)`
139		139
140	`c Check dimensions of the two trajectory files (#tim hard check)`	140	`c Check dimensions of the two trajectory files (#tim hard check)`
141	`if (ntim1.ne.ntim2) then`	141	`if (ntim1.ne.ntim2) then`
142	`print*,'Trajectoy files have different time dimensions... Stop'`	142	`print*,'Trajectoy files have different time dimensions... Stop'`
143	`stop`	143	`stop`
144	`endif`	144	`endif`
145		145
146	`c Check dimensions of the two trajectory files (#tra soft check)`	146	`c Check dimensions of the two trajectory files (#tra soft check)`
147	`if (ntra1.ne.ntra2) then`	147	`if (ntra1.ne.ntra2) then`
148	`print*,'Differing number of trajectories... proceed [y/n]'`	148	`print*,'Differing number of trajectories... proceed [y/n]'`
149	`read*,ch`	149	`read*,ch`
150	`if (ch.eq.'n') stop`	150	`if (ch.eq.'n') stop`
151	`endif`	151	`endif`
152		152
153	`c Check whether difference field is available on both files`	153	`c Check whether difference field is available on both files`
154	`if ( fieldname.ne.'LATLON') then`	154	`if ( fieldname.ne.'LATLON') then`
155	`ind1 = 0`	155	`ind1 = 0`
156	`ind2 = 0`	156	`ind2 = 0`
157	`do i=1,ncol`	157	`do i=1,ncol`
158	`if ( fieldname.eq.vars1(i) ) ind1 = i`	158	`if ( fieldname.eq.vars1(i) ) ind1 = i`
159	`if ( fieldname.eq.vars2(i) ) ind2 = i`	159	`if ( fieldname.eq.vars2(i) ) ind2 = i`
160	`enddo`	160	`enddo`
161	`if ( (ind1.eq.0).or.(ind2.eq.0) ) then`	161	`if ( (ind1.eq.0).or.(ind2.eq.0) ) then`
162	`print*,'Field ',trim(fieldname),' not available... Stop'`	162	`print*,'Field ',trim(fieldname),' not available... Stop'`
163	`stop`	163	`stop`
164	`endif`	164	`endif`
165	`endif`	165	`endif`
166		166
167	`c Check reference dates (soft check)`	167	`c Check reference dates (soft check)`
168	`isok = 1`	168	`isok = 1`
169	`do i=1,6`	169	`do i=1,6`
170	`if ( refdate1(i).ne.refdate2(i) ) isok = 0`	170	`if ( refdate1(i).ne.refdate2(i) ) isok = 0`
171	`enddo`	171	`enddo`
172	`if ( isok.eq.0 ) then`	172	`if ( isok.eq.0 ) then`
173	`print*,'Warning: reference dates differ... proceed [y/n]'`	173	`print*,'Warning: reference dates differ... proceed [y/n]'`
174	`read*,ch`	174	`read*,ch`
175	`if (ch.eq.'n') stop`	175	`if (ch.eq.'n') stop`
176	`endif`	176	`endif`
177		177
178	`c Check trajectory times (soft check)`	178	`c Check trajectory times (soft check)`
179	`isok = 1`	179	`isok = 1`
180	`do i=1,ntim`	180	`do i=1,ntim`
181	`if ( trainp1(1,i,1).ne.trainp2(1,i,1) ) isok = 0`	181	`if ( trainp1(1,i,1).ne.trainp2(1,i,1) ) isok = 0`
182	`enddo`	182	`enddo`
183	`if ( isok.eq.0 ) then`	183	`if ( isok.eq.0 ) then`
184	`print*,'Warning: trajectory times differ... proceed [y/n]'`	184	`print*,'Warning: trajectory times differ... proceed [y/n]'`
185	`read*,ch`	185	`read*,ch`
186	`if (ch.eq.'n') stop`	186	`if (ch.eq.'n') stop`
187	`endif`	187	`endif`
188		188
189	`c Copy reference date to output`	189	`c Copy reference date to output`
190	`do i=1,6`	190	`do i=1,6`
191	`refdate(i) = refdate1(i)`	191	`refdate(i) = refdate1(i)`
192	`enddo`	192	`enddo`
193		193
194	`c ----------------------------------------------------------------------`	194	`c ----------------------------------------------------------------------`
195	`c Calculate the difference (depending on mode)`	195	`c Calculate the difference (depending on mode)`
196	`c ----------------------------------------------------------------------`	196	`c ----------------------------------------------------------------------`
197		197
198	`c Loop over all trajectories`	198	`c Loop over all trajectories`
199	`do i=1,ntra`	199	`do i=1,ntra`
200		200
201	`c Calculate difference for all times`	201	`c Calculate difference for all times`
202	`do j=1,ntim1`	202	`do j=1,ntim1`
203		203
204	`c Calculate the difference (distance or absolute value)`	204	`c Calculate the difference (distance or absolute value)`
205	`if (fieldname.eq.'LATLON') then`	205	`if (fieldname.eq.'LATLON') then`
206	`diff(j) = sdis( trainp1(i,j,2),trainp1(i,j,3),`	206	`diff(j) = sdis( trainp1(i,j,2),trainp1(i,j,3),`
207	`> trainp2(i,j,2),trainp2(i,j,3) )`	207	`> trainp2(i,j,2),trainp2(i,j,3) )`
208	`else`	208	`else`
209	`diff(j) = abs(trainp1(i,j,ind1) - trainp2(i,j,ind2))`	209	`diff(j) = abs(trainp1(i,j,ind1) - trainp2(i,j,ind2))`
210	`endif`	210	`endif`
211		211
212	`enddo`	212	`enddo`
213		213
214	`c Save output for each time`	214	`c Save output for each time`
215	`if ( mode.eq.'single' ) then`	215	`if ( mode.eq.'single' ) then`
216		216
217	`do j=1,ntim`	217	`do j=1,ntim`
218		218
219	`if ( fieldname.eq.'LATLON' ) then`	219	`if ( fieldname.eq.'LATLON' ) then`
220	`traout(i,j, 1) = trainp1(i,j,1) ! time`	220	`traout(i,j, 1) = trainp1(i,j,1) ! time`
221	`traout(i,j, 2) = trainp1(i,j,2) ! lon[1]`	221	`traout(i,j, 2) = trainp1(i,j,2) ! lon[1]`
222	`traout(i,j, 3) = trainp1(i,j,3) ! lat[1]`	222	`traout(i,j, 3) = trainp1(i,j,3) ! lat[1]`
223	`traout(i,j, 4) = trainp1(i,j,4) ! p[1]`	223	`traout(i,j, 4) = trainp1(i,j,4) ! p[1]`
224	`traout(i,j, 5) = trainp2(i,j,2) ! lon[2]`	224	`traout(i,j, 5) = trainp2(i,j,2) ! lon[2]`
225	`traout(i,j, 6) = trainp2(i,j,3) ! lat[2]`	225	`traout(i,j, 6) = trainp2(i,j,3) ! lat[2]`
226	`traout(i,j, 7) = trainp2(i,j,4) ! p[2]`	226	`traout(i,j, 7) = trainp2(i,j,4) ! p[2]`
227	`traout(i,j, 8) = diff(j) ! SDIS(j)`	227	`traout(i,j, 8) = diff(j) ! SDIS(j)`
228	`else`	228	`else`
229	`traout(i,j, 1) = trainp1(i,j,1) ! time`	229	`traout(i,j, 1) = trainp1(i,j,1) ! time`
230	`traout(i,j, 2) = trainp1(i,j,2) ! lon[1]`	230	`traout(i,j, 2) = trainp1(i,j,2) ! lon[1]`
231	`traout(i,j, 3) = trainp1(i,j,3) ! lat[1]`	231	`traout(i,j, 3) = trainp1(i,j,3) ! lat[1]`
232	`traout(i,j, 4) = trainp1(i,j,4) ! p[1]`	232	`traout(i,j, 4) = trainp1(i,j,4) ! p[1]`
233	`traout(i,j, 5) = trainp2(i,j,2) ! lon[2]`	233	`traout(i,j, 5) = trainp2(i,j,2) ! lon[2]`
234	`traout(i,j, 6) = trainp2(i,j,3) ! lat[2]`	234	`traout(i,j, 6) = trainp2(i,j,3) ! lat[2]`
235	`traout(i,j, 7) = trainp2(i,j,4) ! p[2]`	235	`traout(i,j, 7) = trainp2(i,j,4) ! p[2]`
236	`traout(i,j, 8) = trainp1(i,j,ind1) ! field[1]`	236	`traout(i,j, 8) = trainp1(i,j,ind1) ! field[1]`
237	`traout(i,j, 9) = trainp2(i,j,ind2) ! field[2]`	237	`traout(i,j, 9) = trainp2(i,j,ind2) ! field[2]`
238	`traout(i,j,10) = diff(j) ! SDIS(j)`	238	`traout(i,j,10) = diff(j) ! SDIS(j)`
239	`endif`	239	`endif`
240		240
241	`enddo`	241	`enddo`
242		242
243	`c Save only maximum`	243	`c Save only maximum`
244	`elseif ( mode.eq.'max') then`	244	`elseif ( mode.eq.'max') then`
245		245
246	`outind = 1`	246	`outind = 1`
247	`do j=2,ntim1`	247	`do j=2,ntim1`
248	`if ( diff(j).gt.diff(outind) ) outind = j`	248	`if ( diff(j).gt.diff(outind) ) outind = j`
249	`enddo`	249	`enddo`
250		250
251	`if ( fieldname.eq.'LATLON' ) then`	251	`if ( fieldname.eq.'LATLON' ) then`
252	`traout(i,1, 1) = trainp1(i,outind,1) ! time`	252	`traout(i,1, 1) = trainp1(i,outind,1) ! time`
253	`traout(i,1, 2) = trainp1(i,outind,2) ! lon[1]`	253	`traout(i,1, 2) = trainp1(i,outind,2) ! lon[1]`
254	`traout(i,1, 3) = trainp1(i,outind,3) ! lat[1]`	254	`traout(i,1, 3) = trainp1(i,outind,3) ! lat[1]`
255	`traout(i,1, 4) = trainp1(i,outind,4) ! p[1]`	255	`traout(i,1, 4) = trainp1(i,outind,4) ! p[1]`
256	`traout(i,1, 5) = trainp2(i,outind,2) ! lon[2]`	256	`traout(i,1, 5) = trainp2(i,outind,2) ! lon[2]`
257	`traout(i,1, 6) = trainp2(i,outind,3) ! lat[2]`	257	`traout(i,1, 6) = trainp2(i,outind,3) ! lat[2]`
258	`traout(i,1, 7) = trainp2(i,outind,4) ! p[2]`	258	`traout(i,1, 7) = trainp2(i,outind,4) ! p[2]`
259	`traout(i,1, 8) = diff(outind) ! SDIS`	259	`traout(i,1, 8) = diff(outind) ! SDIS`
260	`else`	260	`else`
261	`traout(i,1, 1) = trainp1(i,outind,1) ! time`	261	`traout(i,1, 1) = trainp1(i,outind,1) ! time`
262	`traout(i,1, 2) = trainp1(i,outind,2) ! lon[1]`	262	`traout(i,1, 2) = trainp1(i,outind,2) ! lon[1]`
263	`traout(i,1, 3) = trainp1(i,outind,3) ! lat[1]`	263	`traout(i,1, 3) = trainp1(i,outind,3) ! lat[1]`
264	`traout(i,1, 4) = trainp1(i,outind,4) ! p[1]`	264	`traout(i,1, 4) = trainp1(i,outind,4) ! p[1]`
265	`traout(i,1, 5) = trainp2(i,outind,2) ! lon[2]`	265	`traout(i,1, 5) = trainp2(i,outind,2) ! lon[2]`
266	`traout(i,1, 6) = trainp2(i,outind,3) ! lat[2]`	266	`traout(i,1, 6) = trainp2(i,outind,3) ! lat[2]`
267	`traout(i,1, 7) = trainp2(i,outind,4) ! p[2]`	267	`traout(i,1, 7) = trainp2(i,outind,4) ! p[2]`
268	`traout(i,1, 8) = trainp1(i,outind,ind1) ! field[1]`	268	`traout(i,1, 8) = trainp1(i,outind,ind1) ! field[1]`
269	`traout(i,1, 9) = trainp2(i,outind,ind2) ! field[2]`	269	`traout(i,1, 9) = trainp2(i,outind,ind2) ! field[2]`
270	`traout(i,1,10) = diff(outind) ! SDIS(j)`	270	`traout(i,1,10) = diff(outind) ! SDIS(j)`
271	`endif`	271	`endif`
272		272
273	`endif`	273	`endif`
274		274
275	`enddo`	275	`enddo`
276		276
277	`c ----------------------------------------------------------------------`	277	`c ----------------------------------------------------------------------`
278	`c Write output`	278	`c Write output`
279	`c ----------------------------------------------------------------------`	279	`c ----------------------------------------------------------------------`
280		280
281	`c Write output file`	281	`c Write output file`
282	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`	282	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`
283	`call write_tra(fid,traout,ntra,ntim,ncol,outmode)`	283	`call write_tra(fid,traout,ntra,ntim,ncol,outmode)`
284	`call close_tra(fid,outmode)`	284	`call close_tra(fid,outmode)`
285		285
286	`end`	286	`end`
287		287
288		288
289		289
290		290
291	`c ***********************************************************************`	291	`c ***********************************************************************`
292	`c * Subroutines *`	292	`c * Subroutines *`
293	`c ***********************************************************************`	293	`c ***********************************************************************`
294		294
295	`c ----------------------------------------------------------------------`	295	`c ----------------------------------------------------------------------`
296	`c Spherical distance`	296	`c Spherical distance`
297	`c ----------------------------------------------------------------------`	297	`c ----------------------------------------------------------------------`
298		298
299	`real function sdis(xp,yp,xq,yq)`	299	`real function sdis(xp,yp,xq,yq)`
300	`c`	300	`c`
301	`c calculates spherical distance (in km) between two points given`	301	`c calculates spherical distance (in km) between two points given`
302	`c by their spherical coordinates (xp,yp) and (xq,yq), respectively.`	302	`c by their spherical coordinates (xp,yp) and (xq,yq), respectively.`
303	`c`	303	`c`
304	`real pi180`	304	`real pi180`
305	`parameter (pi180=3.14159/180.)`	305	`parameter (pi180=3.14159/180.)`
306	`real re`	306	`real re`
307	`parameter (re=6370.)`	307	`parameter (re=6370.)`
308	`real degkm`	308	`real degkm`
309	`parameter (degkm=111.1775)`	309	`parameter (degkm=111.1775)`
310	`real xp,yp,xq,yq,arg`	310	`real xp,yp,xq,yq,arg`
311		311
312	`if ( (abs(xp-xq).gt.0.05).and.(abs(yp-yq).gt.0.05) ) then`	312	`if ( (abs(xp-xq).gt.0.05).and.(abs(yp-yq).gt.0.05) ) then`
313	`arg=sin(pi180yp)sin(pi180*yq)+`	313	`arg=sin(pi180yp)sin(pi180*yq)+`
314	`> cos(pi180yp)cos(pi180yq)cos(pi180*(xp-xq))`	314	`> cos(pi180yp)cos(pi180yq)cos(pi180*(xp-xq))`
315	`if (arg.lt.-1.) arg=-1.`	315	`if (arg.lt.-1.) arg=-1.`
316	`if (arg.gt.1.) arg=1.`	316	`if (arg.gt.1.) arg=1.`
317	`sdis=re*acos(arg)`	317	`sdis=re*acos(arg)`
318	`else`	318	`else`
319	`sdis= (yp-yq)*2 + ( (xp-xq) cos( pi1800.5(yp+yq) ) )**2`	319	`sdis= (yp-yq)*2 + ( (xp-xq) cos( pi1800.5(yp+yq) ) )**2`
320	`sdis = deg2km * sqrt(sdis)`	320	`sdis = deg2km * sqrt(sdis)`
321	`endif`	321	`endif`
322		322
323	`end`	323	`end`
324		324
325		325
326		326

Subversion Repositories lagranto.um

(root)/tags/2.0/goodies/difference.f – Rev 3 → 10