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

Subversion Repositories lagranto.ecmwf

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