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


      PROGRAM timeresolution
      
c     ***********************************************************************
c     * Change time resolution of of a trajectory file                      *
c     * Michael Sprenger / Winter 2010                                      *
c     ***********************************************************************
 
      implicit none
      
c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------
 
c     Parameters
      character*80                           inpfile     
      character*80                           outfile     
      integer                                ntra,otim,ncol
      real                                   timeres
      character*80                           unit
      character*80                           mode
 
c     Trajectories
      character*80                           vars(100)   
      integer                                refdate(6)  
      integer                                ntim       
      real,allocatable, dimension (:,:,:) :: trainp        
      real,allocatable, dimension (:,:,:) :: traout        
      real,allocatable, dimension (:)     :: timold,timnew
      real,allocatable, dimension (:)     :: fldold,fldnew
 
c     Numerical constants
      real                                   eps
      parameter                              (eps=0.001)
 
c     Auxiliary variables
      integer                                inpmode
      integer                                outmode
      integer                                stat
      integer                                fid
      integer                                i,j,k
      real                                   hhmm,tfrac
      real                                   range
 
c     ----------------------------------------------------------------------
c     Parameter handling
c     ----------------------------------------------------------------------
 
c     Read parameters
      open(10,file='timeres.param')
       read(10,*) inpfile
       read(10,*) outfile
       read(10,*) ntra,otim,ncol
       read(10,*) timeres
       read(10,*) unit
       read(10,*) mode
      close(10)
 
c     Change unit to minutes
      if ( unit.eq.'min') then
         timeres = 1./60. * timeres
         unit    = 'h'
      endif
 
c     Determine the formats
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1
      call mode_tra(outmode,outfile)
      if (outmode.eq.-1) outmode=1
 
c     ----------------------------------------------------------------------
c     Read input trajectory and allocate memory
c     ----------------------------------------------------------------------
 
c     Allocate memory for input trajectories
      allocate(trainp(ntra,otim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp   ***' 
      allocate(timold(otim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array timold   ***' 
      allocate(fldold(otim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array fldold   ***' 
 
c     Read inpufile
      call ropen_tra(fid,inpfile,ntra,otim,ncol,refdate,vars,inpmode)
      call read_tra (fid,trainp,ntra,otim,ncol,inpmode)
      call close_tra(fid,inpmode)
 
c     Check that first four columns correspond to time,lon,lat,p
      if ( (vars(1).ne.'time' ).or.
     >     (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.
     >     (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.
     >     (vars(4).ne.'ppos' ).and.(vars(4).ne.'p'   ) )
     >then
         print*,' ERROR: problem with input trajectories ...'
         stop
      endif
 
c     Convert all times from hhmm to fractional time
      do i=1,ntra
         do j=1,otim
            hhmm = trainp(i,j,1)
            call hhmm2frac(hhmm,tfrac)
            trainp(i,j,1) = tfrac
         enddo
      enddo
 
c     Get the time range in minutes
      range = ( trainp(1,otim,1) - trainp(1,1,1) ) 
 
c     Determine the new number of times
      ntim = nint( abs( range ) / timeres ) + 1 
 
c     Check that the time range and new time resolution are consistent
      if ( abs( real(ntim-1) * timeres - abs(range) ).gt.eps ) then
         print*,' ERROR: time range and resolution are not compatible'
         print*,'   range              = ',range
         print*,'   (ntim-1) * timeres = ',real(ntim-1) * timeres
         stop
      endif
 
c     Allocate memory for output trajectories
      allocate(traout(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array trainp   ***' 
      allocate(timnew(ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array timnew   ***' 
      allocate(fldnew(ntim),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array fldnew   ***' 
 
c     Define the old and new times
      do i=1,otim
         timold(i) = trainp(1,i,1)
      enddo
      do i=1,ntim
         timnew(i) = timold(1) + real(i-1) * timeres
      enddo
 
c     ----------------------------------------------------------------------
c     Change time resolution
c     ----------------------------------------------------------------------
 
      do i=1,ntra
      do k=2,ncol
 
c        Copy old field
         do j=1,otim
            fldold(j) = trainp(i,j,k)
         enddo
            
c        Exception: Handle date line problem for longitude
         if ( k.eq.2 ) then 
            do j=2,otim
               if ( (fldold(j-1)-fldold(j)).gt.180. ) then
                  fldold(j) = fldold(j) + 360.
               else if ( (fldold(j-1)-fldold(j)).lt.-180. ) then
                  fldold(j) = fldold(j) - 360.
               endif
            enddo
         endif
         
c        Cubic spline fitting
         if ( mode.eq.'cubic' ) then
            call cubicfit (timold,fldold,otim,timnew,fldnew,ntim)
         else if (mode.eq.'linear' ) then
            call linearfit(timold,fldold,otim,timnew,fldnew,ntim)
         endif
 
c        Exception: Reverse date line handling for longitude
         if ( k.eq.2 ) then
            do j=1,ntim
               if ( fldnew(j).gt.180. ) then
                  fldnew(j) = fldnew(j) -360.
               else if ( fldnew(j).lt.-180. ) then
                  fldnew(j) = fldnew(j) +360.
               endif
            enddo
         endif
 
c        Save the new field in the output trajectory
         do j=1,ntim
            traout(i,j,1) = timnew(j)
            traout(i,j,k) = fldnew(j)
         enddo
 
      enddo
      enddo
 
c     ----------------------------------------------------------------------
c     Write output trajectory
c     ----------------------------------------------------------------------
 
c     Convert all times from fractional to hhmm time
      do i=1,ntra
         do j=1,ntim
            tfrac = traout(i,j,1)
            call frac2hhmm(tfrac,hhmm)
            traout(i,j,1) = hhmm
         enddo
      enddo
    
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     * REPARAMETERIZATION SUBROUTINES                                   *
c     ********************************************************************
 
c     -------------------------------------------------------------
c     Interpolation of the trajectory with linear interpolation
c     -------------------------------------------------------------
 
      SUBROUTINE linearfit (time,lon,n,sptime,splon,spn)
 
c     Given the curve <time,lon> with <n> data points, fit a
c     linear fit to this curve. The new curve is returned in 
c     <sptime,splon,spn> with <spn> data points. The parameter
c     <spn> specifies on entry the number of interpolated points
c     along the curve.
      
      implicit none
 
c     Declaration of subroutine parameters
      integer n
      real    time(n),lon(n)
      integer spn
      real    sptime(spn),splon(spn)
 
c     Auxiliary variables
      real    dt
      real    s
      integer i,j,iold
      real    order
 
c     Determine whether the input array is ascending or descending
      if (time(1).gt.time(n)) then
         order=-1.
      else
         order= 1.
      endif
 
c     Bring the time array into ascending order
      do i=1,n
         time(i)=order*time(i)
      enddo
 
c     Prepare the linear interpolation: define the new times
      dt=(time(n)-time(1))/real(spn-1)
      do i=1,spn
         sptime(i)=time(1)+real(i-1)*dt
      enddo
      
c     Do the interpolation
      iold = 1
      do i=1,spn
 
c        Decide which interval of the old time series must be taken
         do j=iold,n-1
            if ( ( sptime(i).ge.time(j  ) ).and.
     >           ( sptime(i).lt.time(j+1) ) ) 
     >      then
               iold = j
               exit
            endif
         enddo
         
c        Do the linear interpolation
         splon(i) = lon(iold) + ( lon(iold+1) - lon(iold) ) * 
     >       ( sptime(i) - time(iold) ) / ( time(iold+1) - time(iold) ) 
 
      enddo
 
c     Change the time arrays back: original order
      do i=1,spn
         sptime(i)=order*sptime(i)
      enddo
      do i=1,n
         time(i)=order*time(i)
      enddo
 
      return
      end
 
 
c     -------------------------------------------------------------
c     Interpolation of the trajectory with a natural cubic spline
c     -------------------------------------------------------------
 
      SUBROUTINE cubicfit (time,lon,n,sptime,splon,spn)
 
c     Given the curve <time,lon> with <n> data points, fit a
c     cubic spline to this curve. The new curve is returned in 
c     <sptime,splon,spn> with <spn> data points. The parameter
c     <spn> specifies on entry the number of spline interpolated points
c     along the curve.
      
      implicit none
 
c     Declaration of subroutine parameters
      integer n
      real    time(n),lon(n)
      integer spn
      real    sptime(spn),splon(spn)
 
c     Auxiliary variables
      real    y2ax(n)
      real    dt
      real    s
      integer i
      real    order
 
c     Determine whether the input array is ascending or descending
      if (time(1).gt.time(n)) then
         order=-1.
      else
         order= 1.
      endif
 
c     Bring the time array into ascending order
      do i=1,n
         time(i)=order*time(i)
      enddo
 
c     Prepare the (natural) cubic spline interpolation
      call spline (time,lon,n,1.e30,1.e30,y2ax)
      dt=(time(n)-time(1))/real(spn-1)
      do i=1,spn
         sptime(i)=time(1)+real(i-1)*dt
      enddo
      
c     Do the spline interpolation
      do i=1,spn
         call splint(time,lon,y2ax,n,sptime(i),s)
         splon(i)=s
      enddo
 
c     Change the time arrays back
      do i=1,spn
         sptime(i)=order*sptime(i)
      enddo
      do i=1,n
         time(i)=order*time(i)
      enddo
 
      return
      end
 
c     -------------------------------------------------------------
c     Basic routines for spline interpolation (Numerical Recipes)
c     -------------------------------------------------------------
 
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
      INTEGER n,NMAX
      REAL yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      REAL p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
 
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
      INTEGER n
      REAL x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      REAL a,b,h
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) then
         print*, 'bad xa input in splint'
         stop
      endif
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *2)/6.
      return
      END
    
 
      

Rev 11	Rev 13
1	`PROGRAM timeresolution`	1	`PROGRAM timeresolution`
2		2
3	`c ***********************************************************************`	3	`c ***********************************************************************`
4	`c * Change time resolution of of a trajectory file *`	4	`c * Change time resolution of of a trajectory file *`
5	`c * Michael Sprenger / Winter 2010 *`	5	`c * Michael Sprenger / Winter 2010 *`
6	`c ***********************************************************************`	6	`c ***********************************************************************`
7		7
8	`implicit none`	8	`implicit none`
9		9
10	`c ----------------------------------------------------------------------`	10	`c ----------------------------------------------------------------------`
11	`c Declaration of variables`	11	`c Declaration of variables`
12	`c ----------------------------------------------------------------------`	12	`c ----------------------------------------------------------------------`
13		13
14	`c Parameters`	14	`c Parameters`
15	`character*80 inpfile`	15	`character*80 inpfile`
16	`character*80 outfile`	16	`character*80 outfile`
17	`integer ntra,otim,ncol`	17	`integer ntra,otim,ncol`
18	`real timeres`	18	`real timeres`
19	`character*80 unit`	19	`character*80 unit`
20	`character*80 mode`	20	`character*80 mode`
21		21
22	`c Trajectories`	22	`c Trajectories`
23	`character*80 vars(100)`	23	`character*80 vars(100)`
24	`integer refdate(6)`	24	`integer refdate(6)`
25	`integer ntim`	25	`integer ntim`
26	`real,allocatable, dimension (:,:,:) :: trainp`	26	`real,allocatable, dimension (:,:,:) :: trainp`
27	`real,allocatable, dimension (:,:,:) :: traout`	27	`real,allocatable, dimension (:,:,:) :: traout`
28	`real,allocatable, dimension (:) :: timold,timnew`	28	`real,allocatable, dimension (:) :: timold,timnew`
29	`real,allocatable, dimension (:) :: fldold,fldnew`	29	`real,allocatable, dimension (:) :: fldold,fldnew`
30		30
31	`c Numerical constants`	31	`c Numerical constants`
32	`real eps`	32	`real eps`
33	`parameter (eps=0.001)`	33	`parameter (eps=0.001)`
34		34
35	`c Auxiliary variables`	35	`c Auxiliary variables`
36	`integer inpmode`	36	`integer inpmode`
37	`integer outmode`	37	`integer outmode`
38	`integer stat`	38	`integer stat`
39	`integer fid`	39	`integer fid`
40	`integer i,j,k`	40	`integer i,j,k`
41	`real hhmm,tfrac`	41	`real hhmm,tfrac`
42	`real range`	42	`real range`
43		43
44	`c ----------------------------------------------------------------------`	44	`c ----------------------------------------------------------------------`
45	`c Parameter handling`	45	`c Parameter handling`
46	`c ----------------------------------------------------------------------`	46	`c ----------------------------------------------------------------------`
47		47
48	`c Read parameters`	48	`c Read parameters`
49	`open(10,file='timeres.param')`	49	`open(10,file='timeres.param')`
50	`read(10,*) inpfile`	50	`read(10,*) inpfile`
51	`read(10,*) outfile`	51	`read(10,*) outfile`
52	`read(10,*) ntra,otim,ncol`	52	`read(10,*) ntra,otim,ncol`
53	`read(10,*) timeres`	53	`read(10,*) timeres`
54	`read(10,*) unit`	54	`read(10,*) unit`
55	`read(10,*) mode`	55	`read(10,*) mode`
56	`close(10)`	56	`close(10)`
57		57
58	`c Change unit to minutes`	58	`c Change unit to minutes`
59	`if ( unit.eq.'min') then`	59	`if ( unit.eq.'min') then`
60	`timeres = 1./60. * timeres`	60	`timeres = 1./60. * timeres`
61	`unit = 'h'`	61	`unit = 'h'`
62	`endif`	62	`endif`
63		63
64	`c Determine the formats`	64	`c Determine the formats`
65	`call mode_tra(inpmode,inpfile)`	65	`call mode_tra(inpmode,inpfile)`
66	`if (inpmode.eq.-1) inpmode=1`	66	`if (inpmode.eq.-1) inpmode=1`
67	`call mode_tra(outmode,outfile)`	67	`call mode_tra(outmode,outfile)`
68	`if (outmode.eq.-1) outmode=1`	68	`if (outmode.eq.-1) outmode=1`
69		69
70	`c ----------------------------------------------------------------------`	70	`c ----------------------------------------------------------------------`
71	`c Read input trajectory and allocate memory`	71	`c Read input trajectory and allocate memory`
72	`c ----------------------------------------------------------------------`	72	`c ----------------------------------------------------------------------`
73		73
74	`c Allocate memory for input trajectories`	74	`c Allocate memory for input trajectories`
75	`allocate(trainp(ntra,otim,ncol),stat=stat)`	75	`allocate(trainp(ntra,otim,ncol),stat=stat)`
76	`if (stat.ne.0) print,' error allocating array trainp *'`	76	`if (stat.ne.0) print,' error allocating array trainp *'`
77	`allocate(timold(otim),stat=stat)`	77	`allocate(timold(otim),stat=stat)`
78	`if (stat.ne.0) print,' error allocating array timold *'`	78	`if (stat.ne.0) print,' error allocating array timold *'`
79	`allocate(fldold(otim),stat=stat)`	79	`allocate(fldold(otim),stat=stat)`
80	`if (stat.ne.0) print,' error allocating array fldold *'`	80	`if (stat.ne.0) print,' error allocating array fldold *'`
81		81
82	`c Read inpufile`	82	`c Read inpufile`
83	`call ropen_tra(fid,inpfile,ntra,otim,ncol,refdate,vars,inpmode)`	83	`call ropen_tra(fid,inpfile,ntra,otim,ncol,refdate,vars,inpmode)`
84	`call read_tra (fid,trainp,ntra,otim,ncol,inpmode)`	84	`call read_tra (fid,trainp,ntra,otim,ncol,inpmode)`
85	`call close_tra(fid,inpmode)`	85	`call close_tra(fid,inpmode)`
86		86
87	`c Check that first four columns correspond to time,lon,lat,p`	87	`c Check that first four columns correspond to time,lon,lat,p`
88	`if ( (vars(1).ne.'time' ).or.`	88	`if ( (vars(1).ne.'time' ).or.`
89	`> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.`	89	`> (vars(2).ne.'xpos' ).and.(vars(2).ne.'lon' ).or.`
90	`> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.`	90	`> (vars(3).ne.'ypos' ).and.(vars(3).ne.'lat' ).or.`
91	`> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )`	91	`> (vars(4).ne.'ppos' ).and.(vars(4).ne.'p' ) )`
92	`>then`	92	`>then`
93	`print*,' ERROR: problem with input trajectories ...'`	93	`print*,' ERROR: problem with input trajectories ...'`
94	`stop`	94	`stop`
95	`endif`	95	`endif`
96		96
97	`c Convert all times from hhmm to fractional time`	97	`c Convert all times from hhmm to fractional time`
98	`do i=1,ntra`	98	`do i=1,ntra`
99	`do j=1,otim`	99	`do j=1,otim`
100	`hhmm = trainp(i,j,1)`	100	`hhmm = trainp(i,j,1)`
101	`call hhmm2frac(hhmm,tfrac)`	101	`call hhmm2frac(hhmm,tfrac)`
102	`trainp(i,j,1) = tfrac`	102	`trainp(i,j,1) = tfrac`
103	`enddo`	103	`enddo`
104	`enddo`	104	`enddo`
105		105
106	`c Get the time range in minutes`	106	`c Get the time range in minutes`
107	`range = ( trainp(1,otim,1) - trainp(1,1,1) )`	107	`range = ( trainp(1,otim,1) - trainp(1,1,1) )`
108		108
109	`c Determine the new number of times`	109	`c Determine the new number of times`
110	`ntim = nint( abs( range ) / timeres ) + 1`	110	`ntim = nint( abs( range ) / timeres ) + 1`
111		111
112	`c Check that the time range and new time resolution are consistent`	112	`c Check that the time range and new time resolution are consistent`
113	`if ( abs( real(ntim-1) * timeres -range ).gt.eps ) then`	113	`if ( abs( real(ntim-1) * timeres - abs(range) ).gt.eps ) then`
114	`print*,' ERROR: time range and resolution are not compatible'`	114	`print*,' ERROR: time range and resolution are not compatible'`
115	`print*,' range = ',range`	115	`print*,' range = ',range`
116	`print,' (ntim-1) timeres = ',real(ntim-1) * timeres`	116	`print,' (ntim-1) timeres = ',real(ntim-1) * timeres`
117	`stop`	117	`stop`
118	`endif`	118	`endif`
119		119
120	`c Allocate memory for output trajectories`	120	`c Allocate memory for output trajectories`
121	`allocate(traout(ntra,ntim,ncol),stat=stat)`	121	`allocate(traout(ntra,ntim,ncol),stat=stat)`
122	`if (stat.ne.0) print,' error allocating array trainp *'`	122	`if (stat.ne.0) print,' error allocating array trainp *'`
123	`allocate(timnew(ntim),stat=stat)`	123	`allocate(timnew(ntim),stat=stat)`
124	`if (stat.ne.0) print,' error allocating array timnew *'`	124	`if (stat.ne.0) print,' error allocating array timnew *'`
125	`allocate(fldnew(ntim),stat=stat)`	125	`allocate(fldnew(ntim),stat=stat)`
126	`if (stat.ne.0) print,' error allocating array fldnew *'`	126	`if (stat.ne.0) print,' error allocating array fldnew *'`
127		127
128	`c Define the old and new times`	128	`c Define the old and new times`
129	`do i=1,otim`	129	`do i=1,otim`
130	`timold(i) = trainp(1,i,1)`	130	`timold(i) = trainp(1,i,1)`
131	`enddo`	131	`enddo`
132	`do i=1,ntim`	132	`do i=1,ntim`
133	`timnew(i) = timold(1) + real(i-1) * timeres`	133	`timnew(i) = timold(1) + real(i-1) * timeres`
134	`enddo`	134	`enddo`
135		135
136	`c ----------------------------------------------------------------------`	136	`c ----------------------------------------------------------------------`
137	`c Change time resolution`	137	`c Change time resolution`
138	`c ----------------------------------------------------------------------`	138	`c ----------------------------------------------------------------------`
139		139
140	`do i=1,ntra`	140	`do i=1,ntra`
141	`do k=2,ncol`	141	`do k=2,ncol`
142		142
143	`c Copy old field`	143	`c Copy old field`
144	`do j=1,otim`	144	`do j=1,otim`
145	`fldold(j) = trainp(i,j,k)`	145	`fldold(j) = trainp(i,j,k)`
146	`enddo`	146	`enddo`
147		147
148	`c Exception: Handle date line problem for longitude`	148	`c Exception: Handle date line problem for longitude`
149	`if ( k.eq.2 ) then`	149	`if ( k.eq.2 ) then`
150	`do j=2,otim`	150	`do j=2,otim`
151	`if ( (fldold(j-1)-fldold(j)).gt.180. ) then`	151	`if ( (fldold(j-1)-fldold(j)).gt.180. ) then`
152	`fldold(j) = fldold(j) + 360.`	152	`fldold(j) = fldold(j) + 360.`
153	`else if ( (fldold(j-1)-fldold(j)).lt.-180. ) then`	153	`else if ( (fldold(j-1)-fldold(j)).lt.-180. ) then`
154	`fldold(j) = fldold(j) - 360.`	154	`fldold(j) = fldold(j) - 360.`
155	`endif`	155	`endif`
156	`enddo`	156	`enddo`
157	`endif`	157	`endif`
158		158
159	`c Cubic spline fitting`	159	`c Cubic spline fitting`
160	`if ( mode.eq.'cubic' ) then`	160	`if ( mode.eq.'cubic' ) then`
161	`call cubicfit (timold,fldold,otim,timnew,fldnew,ntim)`	161	`call cubicfit (timold,fldold,otim,timnew,fldnew,ntim)`
162	`else if (mode.eq.'linear' ) then`	162	`else if (mode.eq.'linear' ) then`
163	`call linearfit(timold,fldold,otim,timnew,fldnew,ntim)`	163	`call linearfit(timold,fldold,otim,timnew,fldnew,ntim)`
164	`endif`	164	`endif`
165		165
166	`c Exception: Reverse date line handling for longitude`	166	`c Exception: Reverse date line handling for longitude`
167	`if ( k.eq.2 ) then`	167	`if ( k.eq.2 ) then`
168	`do j=1,ntim`	168	`do j=1,ntim`
169	`if ( fldnew(j).gt.180. ) then`	169	`if ( fldnew(j).gt.180. ) then`
170	`fldnew(j) = fldnew(j) -360.`	170	`fldnew(j) = fldnew(j) -360.`
171	`else if ( fldnew(j).lt.-180. ) then`	171	`else if ( fldnew(j).lt.-180. ) then`
172	`fldnew(j) = fldnew(j) +360.`	172	`fldnew(j) = fldnew(j) +360.`
173	`endif`	173	`endif`
174	`enddo`	174	`enddo`
175	`endif`	175	`endif`
176		176
177	`c Save the new field in the output trajectory`	177	`c Save the new field in the output trajectory`
178	`do j=1,ntim`	178	`do j=1,ntim`
179	`traout(i,j,1) = timnew(j)`	179	`traout(i,j,1) = timnew(j)`
180	`traout(i,j,k) = fldnew(j)`	180	`traout(i,j,k) = fldnew(j)`
181	`enddo`	181	`enddo`
182		182
183	`enddo`	183	`enddo`
184	`enddo`	184	`enddo`
185		185
186	`c ----------------------------------------------------------------------`	186	`c ----------------------------------------------------------------------`
187	`c Write output trajectory`	187	`c Write output trajectory`
188	`c ----------------------------------------------------------------------`	188	`c ----------------------------------------------------------------------`
189		189
190	`c Convert all times from fractional to hhmm time`	190	`c Convert all times from fractional to hhmm time`
191	`do i=1,ntra`	191	`do i=1,ntra`
192	`do j=1,ntim`	192	`do j=1,ntim`
193	`tfrac = traout(i,j,1)`	193	`tfrac = traout(i,j,1)`
194	`call frac2hhmm(tfrac,hhmm)`	194	`call frac2hhmm(tfrac,hhmm)`
195	`traout(i,j,1) = hhmm`	195	`traout(i,j,1) = hhmm`
196	`enddo`	196	`enddo`
197	`enddo`	197	`enddo`
198		198
199	`c Write output file`	199	`c Write output file`
200	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`	200	`call wopen_tra(fid,outfile,ntra,ntim,ncol,refdate,vars,outmode)`
201	`call write_tra(fid,traout,ntra,ntim,ncol,outmode)`	201	`call write_tra(fid,traout,ntra,ntim,ncol,outmode)`
202	`call close_tra(fid,outmode)`	202	`call close_tra(fid,outmode)`
203		203
204	`end`	204	`end`
205		205
206		206
207	`c ********************************************************************`	207	`c ********************************************************************`
208	`c * REPARAMETERIZATION SUBROUTINES *`	208	`c * REPARAMETERIZATION SUBROUTINES *`
209	`c ********************************************************************`	209	`c ********************************************************************`
210		210
211	`c -------------------------------------------------------------`	211	`c -------------------------------------------------------------`
212	`c Interpolation of the trajectory with linear interpolation`	212	`c Interpolation of the trajectory with linear interpolation`
213	`c -------------------------------------------------------------`	213	`c -------------------------------------------------------------`
214		214
215	`SUBROUTINE linearfit (time,lon,n,sptime,splon,spn)`	215	`SUBROUTINE linearfit (time,lon,n,sptime,splon,spn)`
216		216
217	`c Given the curve <time,lon> with <n> data points, fit a`	217	`c Given the curve <time,lon> with <n> data points, fit a`
218	`c linear fit to this curve. The new curve is returned in`	218	`c linear fit to this curve. The new curve is returned in`
219	`c <sptime,splon,spn> with <spn> data points. The parameter`	219	`c <sptime,splon,spn> with <spn> data points. The parameter`
220	`c <spn> specifies on entry the number of interpolated points`	220	`c <spn> specifies on entry the number of interpolated points`
221	`c along the curve.`	221	`c along the curve.`
222		222
223	`implicit none`	223	`implicit none`
224		224
225	`c Declaration of subroutine parameters`	225	`c Declaration of subroutine parameters`
226	`integer n`	226	`integer n`
227	`real time(n),lon(n)`	227	`real time(n),lon(n)`
228	`integer spn`	228	`integer spn`
229	`real sptime(spn),splon(spn)`	229	`real sptime(spn),splon(spn)`
230		230
231	`c Auxiliary variables`	231	`c Auxiliary variables`
232	`real dt`	232	`real dt`
233	`real s`	233	`real s`
234	`integer i,j,iold`	234	`integer i,j,iold`
235	`real order`	235	`real order`
236		236
237	`c Determine whether the input array is ascending or descending`	237	`c Determine whether the input array is ascending or descending`
238	`if (time(1).gt.time(n)) then`	238	`if (time(1).gt.time(n)) then`
239	`order=-1.`	239	`order=-1.`
240	`else`	240	`else`
241	`order= 1.`	241	`order= 1.`
242	`endif`	242	`endif`
243		243
244	`c Bring the time array into ascending order`	244	`c Bring the time array into ascending order`
245	`do i=1,n`	245	`do i=1,n`
246	`time(i)=order*time(i)`	246	`time(i)=order*time(i)`
247	`enddo`	247	`enddo`
248		248
249	`c Prepare the linear interpolation: define the new times`	249	`c Prepare the linear interpolation: define the new times`
250	`dt=(time(n)-time(1))/real(spn-1)`	250	`dt=(time(n)-time(1))/real(spn-1)`
251	`do i=1,spn`	251	`do i=1,spn`
252	`sptime(i)=time(1)+real(i-1)*dt`	252	`sptime(i)=time(1)+real(i-1)*dt`
253	`enddo`	253	`enddo`
254		254
255	`c Do the interpolation`	255	`c Do the interpolation`
256	`iold = 1`	256	`iold = 1`
257	`do i=1,spn`	257	`do i=1,spn`
258		258
259	`c Decide which interval of the old time series must be taken`	259	`c Decide which interval of the old time series must be taken`
260	`do j=iold,n-1`	260	`do j=iold,n-1`
261	`if ( ( sptime(i).ge.time(j ) ).and.`	261	`if ( ( sptime(i).ge.time(j ) ).and.`
262	`> ( sptime(i).lt.time(j+1) ) )`	262	`> ( sptime(i).lt.time(j+1) ) )`
263	`> then`	263	`> then`
264	`iold = j`	264	`iold = j`
265	`exit`	265	`exit`
266	`endif`	266	`endif`
267	`enddo`	267	`enddo`
268		268
269	`c Do the linear interpolation`	269	`c Do the linear interpolation`
270	`splon(i) = lon(iold) + ( lon(iold+1) - lon(iold) ) *`	270	`splon(i) = lon(iold) + ( lon(iold+1) - lon(iold) ) *`
271	`> ( sptime(i) - time(iold) ) / ( time(iold+1) - time(iold) )`	271	`> ( sptime(i) - time(iold) ) / ( time(iold+1) - time(iold) )`
272		272
273	`enddo`	273	`enddo`
274		274
275	`c Change the time arrays back: original order`	275	`c Change the time arrays back: original order`
276	`do i=1,spn`	276	`do i=1,spn`
277	`sptime(i)=order*sptime(i)`	277	`sptime(i)=order*sptime(i)`
278	`enddo`	278	`enddo`
279	`do i=1,n`	279	`do i=1,n`
280	`time(i)=order*time(i)`	280	`time(i)=order*time(i)`
281	`enddo`	281	`enddo`
282		282
283	`return`	283	`return`
284	`end`	284	`end`
285		285
286		286
287	`c -------------------------------------------------------------`	287	`c -------------------------------------------------------------`
288	`c Interpolation of the trajectory with a natural cubic spline`	288	`c Interpolation of the trajectory with a natural cubic spline`
289	`c -------------------------------------------------------------`	289	`c -------------------------------------------------------------`
290		290
291	`SUBROUTINE cubicfit (time,lon,n,sptime,splon,spn)`	291	`SUBROUTINE cubicfit (time,lon,n,sptime,splon,spn)`
292		292
293	`c Given the curve <time,lon> with <n> data points, fit a`	293	`c Given the curve <time,lon> with <n> data points, fit a`
294	`c cubic spline to this curve. The new curve is returned in`	294	`c cubic spline to this curve. The new curve is returned in`
295	`c <sptime,splon,spn> with <spn> data points. The parameter`	295	`c <sptime,splon,spn> with <spn> data points. The parameter`
296	`c <spn> specifies on entry the number of spline interpolated points`	296	`c <spn> specifies on entry the number of spline interpolated points`
297	`c along the curve.`	297	`c along the curve.`
298		298
299	`implicit none`	299	`implicit none`
300		300
301	`c Declaration of subroutine parameters`	301	`c Declaration of subroutine parameters`
302	`integer n`	302	`integer n`
303	`real time(n),lon(n)`	303	`real time(n),lon(n)`
304	`integer spn`	304	`integer spn`
305	`real sptime(spn),splon(spn)`	305	`real sptime(spn),splon(spn)`
306		306
307	`c Auxiliary variables`	307	`c Auxiliary variables`
308	`real y2ax(n)`	308	`real y2ax(n)`
309	`real dt`	309	`real dt`
310	`real s`	310	`real s`
311	`integer i`	311	`integer i`
312	`real order`	312	`real order`
313		313
314	`c Determine whether the input array is ascending or descending`	314	`c Determine whether the input array is ascending or descending`
315	`if (time(1).gt.time(n)) then`	315	`if (time(1).gt.time(n)) then`
316	`order=-1.`	316	`order=-1.`
317	`else`	317	`else`
318	`order= 1.`	318	`order= 1.`
319	`endif`	319	`endif`
320		320
321	`c Bring the time array into ascending order`	321	`c Bring the time array into ascending order`
322	`do i=1,n`	322	`do i=1,n`
323	`time(i)=order*time(i)`	323	`time(i)=order*time(i)`
324	`enddo`	324	`enddo`
325		325
326	`c Prepare the (natural) cubic spline interpolation`	326	`c Prepare the (natural) cubic spline interpolation`
327	`call spline (time,lon,n,1.e30,1.e30,y2ax)`	327	`call spline (time,lon,n,1.e30,1.e30,y2ax)`
328	`dt=(time(n)-time(1))/real(spn-1)`	328	`dt=(time(n)-time(1))/real(spn-1)`
329	`do i=1,spn`	329	`do i=1,spn`
330	`sptime(i)=time(1)+real(i-1)*dt`	330	`sptime(i)=time(1)+real(i-1)*dt`
331	`enddo`	331	`enddo`
332		332
333	`c Do the spline interpolation`	333	`c Do the spline interpolation`
334	`do i=1,spn`	334	`do i=1,spn`
335	`call splint(time,lon,y2ax,n,sptime(i),s)`	335	`call splint(time,lon,y2ax,n,sptime(i),s)`
336	`splon(i)=s`	336	`splon(i)=s`
337	`enddo`	337	`enddo`
338		338
339	`c Change the time arrays back`	339	`c Change the time arrays back`
340	`do i=1,spn`	340	`do i=1,spn`
341	`sptime(i)=order*sptime(i)`	341	`sptime(i)=order*sptime(i)`
342	`enddo`	342	`enddo`
343	`do i=1,n`	343	`do i=1,n`
344	`time(i)=order*time(i)`	344	`time(i)=order*time(i)`
345	`enddo`	345	`enddo`
346		346
347	`return`	347	`return`
348	`end`	348	`end`
349		349
350	`c -------------------------------------------------------------`	350	`c -------------------------------------------------------------`
351	`c Basic routines for spline interpolation (Numerical Recipes)`	351	`c Basic routines for spline interpolation (Numerical Recipes)`
352	`c -------------------------------------------------------------`	352	`c -------------------------------------------------------------`
353		353
354	`SUBROUTINE spline(x,y,n,yp1,ypn,y2)`	354	`SUBROUTINE spline(x,y,n,yp1,ypn,y2)`
355	`INTEGER n,NMAX`	355	`INTEGER n,NMAX`
356	`REAL yp1,ypn,x(n),y(n),y2(n)`	356	`REAL yp1,ypn,x(n),y(n),y2(n)`
357	`PARAMETER (NMAX=500)`	357	`PARAMETER (NMAX=500)`
358	`INTEGER i,k`	358	`INTEGER i,k`
359	`REAL p,qn,sig,un,u(NMAX)`	359	`REAL p,qn,sig,un,u(NMAX)`
360	`if (yp1.gt..99e30) then`	360	`if (yp1.gt..99e30) then`
361	`y2(1)=0.`	361	`y2(1)=0.`
362	`u(1)=0.`	362	`u(1)=0.`
363	`else`	363	`else`
364	`y2(1)=-0.5`	364	`y2(1)=-0.5`
365	`u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)`	365	`u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)`
366	`endif`	366	`endif`
367	`do 11 i=2,n-1`	367	`do 11 i=2,n-1`
368	`sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))`	368	`sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))`
369	`p=sig*y2(i-1)+2.`	369	`p=sig*y2(i-1)+2.`
370	`y2(i)=(sig-1.)/p`	370	`y2(i)=(sig-1.)/p`
371	`u(i)=(6.*((y(i+1)-y(i))/(x(i+`	371	`u(i)=(6.*((y(i+1)-y(i))/(x(i+`
372	`1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig`	372	`1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig`
373	`*u(i-1))/p`	373	`*u(i-1))/p`
374	`11 continue`	374	`11 continue`
375	`if (ypn.gt..99e30) then`	375	`if (ypn.gt..99e30) then`
376	`qn=0.`	376	`qn=0.`
377	`un=0.`	377	`un=0.`
378	`else`	378	`else`
379	`qn=0.5`	379	`qn=0.5`
380	`un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))`	380	`un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))`
381	`endif`	381	`endif`
382	`y2(n)=(un-qnu(n-1))/(qny2(n-1)+1.)`	382	`y2(n)=(un-qnu(n-1))/(qny2(n-1)+1.)`
383	`do 12 k=n-1,1,-1`	383	`do 12 k=n-1,1,-1`
384	`y2(k)=y2(k)*y2(k+1)+u(k)`	384	`y2(k)=y2(k)*y2(k+1)+u(k)`
385	`12 continue`	385	`12 continue`
386	`return`	386	`return`
387	`END`	387	`END`
388		388
389	`SUBROUTINE splint(xa,ya,y2a,n,x,y)`	389	`SUBROUTINE splint(xa,ya,y2a,n,x,y)`
390	`INTEGER n`	390	`INTEGER n`
391	`REAL x,y,xa(n),y2a(n),ya(n)`	391	`REAL x,y,xa(n),y2a(n),ya(n)`
392	`INTEGER k,khi,klo`	392	`INTEGER k,khi,klo`
393	`REAL a,b,h`	393	`REAL a,b,h`
394	`klo=1`	394	`klo=1`
395	`khi=n`	395	`khi=n`
396	`1 if (khi-klo.gt.1) then`	396	`1 if (khi-klo.gt.1) then`
397	`k=(khi+klo)/2`	397	`k=(khi+klo)/2`
398	`if(xa(k).gt.x)then`	398	`if(xa(k).gt.x)then`
399	`khi=k`	399	`khi=k`
400	`else`	400	`else`
401	`klo=k`	401	`klo=k`
402	`endif`	402	`endif`
403	`goto 1`	403	`goto 1`
404	`endif`	404	`endif`
405	`h=xa(khi)-xa(klo)`	405	`h=xa(khi)-xa(klo)`
406	`if (h.eq.0.) then`	406	`if (h.eq.0.) then`
407	`print*, 'bad xa input in splint'`	407	`print*, 'bad xa input in splint'`
408	`stop`	408	`stop`
409	`endif`	409	`endif`
410	`a=(xa(khi)-x)/h`	410	`a=(xa(khi)-x)/h`
411	`b=(x-xa(klo))/h`	411	`b=(x-xa(klo))/h`
412	`y=aya(klo)+bya(khi)+((a*3-a)y2a(klo)+(b*3-b)y2a(khi))(h*`	412	`y=aya(klo)+bya(khi)+((a*3-a)y2a(klo)+(b*3-b)y2a(khi))(h*`
413	`*2)/6.`	413	`*2)/6.`
414	`return`	414	`return`
415	`END`	415	`END`
416		416
417		417
418		418

Subversion Repositories lagranto.ecmwf

(root)/trunk/goodies/timeres.f – Rev 11 → 13