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

Subversion Repositories lagranto.ecmwf

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