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