Subversion Repositories lagranto.ecmwf

Rev

Rev 11 | Go to most recent revision | Details | Last modification | View Log | RSS feed

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