source: palm/trunk/SOURCE/spline_x.f90 @ 1

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 19.9 KB
Line 
1 SUBROUTINE spline_x( vad_in_out, ad_v, var_char )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: spline_x.f90,v $
11! Revision 1.8  2004/04/30 12:54:20  raasch
12! Names of transpose indices changed, enlarged transposition arrays introduced
13!
14! Revision 1.7  2003/03/16 09:47:17  raasch
15! Two underscores (_) are placed in front of all define-strings
16!
17! Revision 1.6  2001/03/30 07:53:19  raasch
18! Arrays r, vad and wrk_spline changed from 3D to 2D and removed from argument
19! list. Several loops over k/j (parallel/non-parallel) combined to one loop.
20! Application of long filter moved to this routine. All comments and identifiers
21! translated into English.
22!
23! Revision 1.5  2001/01/22 08:07:41  raasch
24! Module test_variables removed
25!
26! Revision 1.4  1999/03/25 07:32:31  raasch
27! Filterung der Ueberschwinger geschieht optional, ups_limit_e eingefuehrt,
28! Ueberschwinger werden in gewissen Grenzen erlaubt
29!
30! Revision 1.3  1999/02/26 17:53:26  schroeter
31! - Gradientenkontrolle fuer den nicht-parallelen Teil
32! - statistische Auswertung ueber den prozentualen Anteil des
33!   Upstream-Verfahrens an der Gesamtadvektion fuer nicht-
34!   parallelen Teil
35!
36! Revision 1.2  1999/02/17 09:30:19  raasch
37! Wertebegrenzung zur Verhinderung von Ueberschwingern
38!
39! Revision 1.1  1999/02/05 09:15:59  raasch
40! Initial revision
41!
42!
43! Description:
44! ------------
45! Upstream-spline advection along x
46!
47! Input/output parameters:
48! ad_v       = advecting wind speed component
49! vad_in_out = quantity to be advected, excluding ghost- or cyclic boundaries
50!              result is given to the calling routine in this array
51! var_char   = string which defines the quantity to be advected
52!
53! Internal arrays:
54! r          = 2D-working array (right hand side of linear equation, buffer for
55!              Long filter)
56! tf         = tendency field (2D), used for long filter
57! vad        = quantity to be advected (2D), including ghost- or cyclic
58!              boundarys along the direction of advection
59! wrk_long   = working array (long coefficients)
60! wrk_spline = working array (spline coefficients)
61!------------------------------------------------------------------------------!
62
63    USE advection
64    USE grid_variables
65    USE indices
66    USE statistics
67    USE control_parameters
68    USE transpose_indices
69
70    IMPLICIT NONE
71
72    CHARACTER (LEN=*) ::  var_char
73
74    INTEGER ::  component, i, j, k, sr
75    REAL    ::  overshoot_limit, sm_faktor, t1, t2, t3, ups_limit
76    REAL, DIMENSION(:,:),   ALLOCATABLE ::  r, tf, vad, wrk_spline
77    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  wrk_long
78
79#if defined( __parallel )
80    REAL ::  ad_v(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa),       &
81             vad_in_out(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa)
82#else
83    REAL ::  ad_v(nzb+1:nzt,nys:nyn,nxl:nxr),               &
84             vad_in_out(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
85#endif
86
87!
88!-- Set criteria for switching between upstream- and upstream-spline-method
89    IF ( var_char == 'u' )  THEN
90       overshoot_limit = overshoot_limit_u
91       ups_limit       = ups_limit_u
92       component       = 1
93    ELSEIF (  var_char == 'v' )  THEN
94       overshoot_limit = overshoot_limit_v
95       ups_limit       = ups_limit_v
96       component       = 2
97    ELSEIF (  var_char == 'w' )  THEN
98       overshoot_limit = overshoot_limit_w
99       ups_limit       = ups_limit_w
100       component       = 3
101    ELSEIF (  var_char == 'pt' )  THEN
102       overshoot_limit = overshoot_limit_pt
103       ups_limit       = ups_limit_pt
104       component       = 4
105    ELSEIF (  var_char == 'e' )  THEN
106       overshoot_limit = overshoot_limit_e
107       ups_limit       = ups_limit_e
108       component       = 5
109    ENDIF
110
111!
112!-- Initialize calculation of relative upstream fraction
113    sums_up_fraction_l(component,1,:) = 0.0
114
115#if defined( __parallel )
116
117!
118!-- Allocate working arrays
119    ALLOCATE( r(-1:nx+1,nys_x:nyn_x), vad(-1:nx+1,nys_x:nyn_x), &
120              wrk_spline(0:nx,nys_x:nyn_x) )
121    IF ( long_filter_factor /= 0.0 )  THEN
122       ALLOCATE( tf(0:nx,nys_x:nyn_x), wrk_long(0:nx,nys_x:nyn_x,1:3) )
123    ENDIF
124
125!
126!-- Loop over all gridpoints along z
127    DO  k = nzb_x, nzt_x
128
129!
130!--    Store array to be advected on work array and add cyclic boundary along x
131       vad(0:nx,nys_x:nyn_x) = vad_in_out(0:nx,nys_x:nyn_x,k)
132       vad(-1,:)   = vad(nx,:)
133       vad(nx+1,:) = vad(0,:)
134
135!
136!--    Calculate right hand side
137       DO  j = nys_x, nyn_x
138          DO  i = 0, nx
139             r(i,j) = 3.0 * (                                                  &
140                            spl_tri_x(2,i) * ( vad(i,j) - vad(i-1,j) ) * ddx + &
141                            spl_tri_x(3,i) * ( vad(i+1,j) - vad(i,j) ) * ddx   &
142                            )
143          ENDDO
144       ENDDO
145
146!
147!--    Forward substitution
148       DO  j = nys_x, nyn_x
149          wrk_spline(0,j) = r(0,j)
150          DO  i = 1, nx
151             wrk_spline(i,j) = r(i,j) - spl_tri_x(5,i) * wrk_spline(i-1,j)
152          ENDDO
153       ENDDO
154
155!
156!--    Backward substitution (Sherman-Morrison-formula)
157       DO  j = nys_x, nyn_x
158          r(nx,j) = wrk_spline(nx,j) / spl_tri_x(4,nx)
159          DO  i = nx-1, 0, -1
160             r(i,j) = ( wrk_spline(i,j) - spl_tri_x(3,i) * r(i+1,j) ) / &
161                      spl_tri_x(4,i)
162          ENDDO
163          sm_faktor = ( r(0,j) + 0.5 * r(nx,j) / spl_gamma_x ) / &
164                      ( 1.0 + spl_z_x(0) + 0.5 * spl_z_x(nx) / spl_gamma_x )
165          DO  i = 0, nx
166             r(i,j) = r(i,j) - sm_faktor * spl_z_x(i)
167          ENDDO
168       ENDDO
169
170!
171!--    Add cyclic boundary to right hand side
172       r(-1,:)   = r(nx,:)
173       r(nx+1,:) = r(0,:)
174
175!
176!--    Calculate advection along x
177       DO  j = nys_x, nyn_x
178          DO  i = 0, nx
179
180             IF ( ad_v(i,j,k) == 0.0 )  THEN
181
182                vad_in_out(i,j,k) = vad(i,j)
183
184             ELSEIF ( ad_v(i,j,k) > 0.0 )  THEN
185
186                IF ( ABS( vad(i,j) - vad(i-1,j) ) <= ups_limit )  THEN
187                   vad_in_out(i,j,k) = vad(i,j) - dt_3d * ad_v(i,j,k) * &
188                                       ( vad(i,j) - vad(i-1,j) ) * ddx
189!
190!--                Calculate upstream fraction in % (s. flow_statistics)
191                   DO  sr = 0, statistic_regions
192                      sums_up_fraction_l(component,1,sr) = &
193                               sums_up_fraction_l(component,1,sr) + 1.0
194                   ENDDO
195                ELSE
196                   t1 = ad_v(i,j,k) * dt_3d * ddx
197                   t2 = 3.0 * ( vad(i-1,j) - vad(i,j) ) + &
198                        ( 2.0 * r(i,j) + r(i-1,j) ) * dx
199                   t3 = 2.0 * ( vad(i-1,j) - vad(i,j) ) + &
200                        ( r(i,j) + r(i-1,j) ) * dx
201                   vad_in_out(i,j,k) = vad(i,j) - r(i,j) * t1 * dx + &
202                                       t2 * t1**2 - t3 * t1**3
203                   IF ( vad(i-1,j) == vad(i,j) )  THEN
204                      vad_in_out(i,j,k) = vad(i,j)
205                   ENDIF
206                ENDIF
207
208             ELSE
209
210                IF ( ABS( vad(i,j) - vad(i+1,j) ) <= ups_limit )  THEN
211                   vad_in_out(i,j,k) = vad(i,j) - dt_3d * ad_v(i,j,k) * &
212                                       ( vad(i+1,j) - vad(i,j) ) * ddx
213!
214!--                Calculate upstream fraction in % (s. flow_statistics)
215                   DO  sr = 0, statistic_regions
216                      sums_up_fraction_l(component,1,sr) = &
217                               sums_up_fraction_l(component,1,sr) + 1.0
218                   ENDDO
219                ELSE
220                   t1 = -ad_v(i,j,k) * dt_3d * ddx
221                   t2 = 3.0 * ( vad(i,j) - vad(i+1,j) ) + &
222                        ( 2.0 * r(i,j) + r(i+1,j) ) * dx
223                   t3 = 2.0 * ( vad(i,j) - vad(i+1,j) ) + &
224                        ( r(i,j) + r(i+1,j) ) * dx
225                   vad_in_out(i,j,k) = vad(i,j) + r(i,j) * t1 * dx - &
226                                       t2 * t1**2 + t3 * t1**3
227                   IF ( vad(i+1,j) == vad(i,j) )  THEN
228                      vad_in_out(i,j,k) = vad(i,j)
229                   ENDIF
230                ENDIF
231
232             ENDIF
233          ENDDO
234       ENDDO
235
236!
237!--    Limit values in order to prevent overshooting
238       IF ( cut_spline_overshoot )  THEN
239
240          DO  j = nys_x, nyn_x
241             DO  i = 0, nx
242                IF ( ad_v(i,j,k) > 0.0 ) THEN
243                   IF ( vad(i,j) > vad(i-1,j) )  THEN
244                      vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), &
245                                               vad(i,j)   + overshoot_limit )
246                      vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), &
247                                               vad(i-1,j) - overshoot_limit )
248                   ELSE
249                      vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), &
250                                               vad(i,j)   - overshoot_limit )
251                      vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), &
252                                               vad(i-1,j) + overshoot_limit )
253                   ENDIF
254                ELSE
255                   IF ( vad(i,j) > vad(i+1,j) )  THEN
256                      vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), &
257                                               vad(i,j)   + overshoot_limit )
258                      vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), &
259                                               vad(i+1,j) - overshoot_limit )
260                   ELSE
261                      vad_in_out(i,j,k) = MAX( vad_in_out(i,j,k), &
262                                               vad(i,j)   - overshoot_limit )
263                      vad_in_out(i,j,k) = MIN( vad_in_out(i,j,k), &
264                                               vad(i+1,j) + overshoot_limit )
265                   ENDIF
266                ENDIF
267             ENDDO
268          ENDDO
269
270       ENDIF
271
272!
273!--    Long-filter (acting on tendency only)
274       IF ( long_filter_factor /= 0.0 )  THEN
275
276!
277!--       Compute tendency
278          DO  j = nys_x, nyn_x
279             DO  i = 0, nx
280                tf(i,j) = vad_in_out(i,j,k) - vad(i,j)
281             ENDDO
282          ENDDO
283
284!
285!--       Apply the filter.
286          DO  j = nys_x, nyn_x
287             wrk_long(0,j,1) = 2.0 * ( 1.0 + long_filter_factor )
288             wrk_long(0,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(0,j,1)
289             wrk_long(0,j,3) = ( long_filter_factor * tf(nx,j) + &
290                                 2.0 * tf(0,j) + tf(1,j)         &
291                               ) / wrk_long(0,j,1)
292             DO  i = 1, nx-1
293                wrk_long(i,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
294                                  ( 1.0 - long_filter_factor ) * wrk_long(i-1,j,2)
295                wrk_long(i,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(i,j,1)
296                wrk_long(i,j,3) = ( tf(i-1,j) + 2.0 * tf(i,j)  +               &
297                                    tf(i+1,j) - ( 1.0 - long_filter_factor ) * &
298                                    wrk_long(i-1,j,3) ) / wrk_long(i,j,1)
299             ENDDO
300             wrk_long(nx,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
301                                ( 1.0 - long_filter_factor ) * wrk_long(nx-1,j,2)
302             wrk_long(nx,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(nx,j,1)
303             wrk_long(nx,j,3) = ( tf(nx-1,j) + 2.0 * tf(nx,j) +  &
304                                  long_filter_factor * tf(0,j) - &
305                                  ( 1.0 - long_filter_factor ) * &
306                                  wrk_long(nx-1,j,3)             &
307                                ) / wrk_long(nx,j,1)
308             r(nx,j) = wrk_long(nx,j,3)
309          ENDDO
310
311          DO  i = nx-1, 0, -1
312             DO  j = nys_x, nyn_x
313                r(i,j) = wrk_long(i,j,3) - wrk_long(i,j,2) * r(i+1,j)
314             ENDDO
315          ENDDO
316
317          DO  j = nys_x, nyn_x
318             DO  i = 0, nx
319                vad_in_out(i,j,k) = vad(i,j) + r(i,j)
320             ENDDO
321          ENDDO
322
323       ENDIF  ! Long filter
324
325    ENDDO
326
327#else
328
329!
330!-- Allocate working arrays
331    ALLOCATE( r(nzb+1:nzt,nxl-1:nxr+1), vad(nzb:nzt+1,nxl-1:nxr+1), &
332              wrk_spline(nzb+1:nzt,nxl-1:nxr+1) )
333    IF ( long_filter_factor /= 0.0 )  THEN
334       ALLOCATE( tf(nzb+1:nzt,nxl-1:nxr+1), wrk_long(nzb+1:nzt,0:nx,1:3) )
335    ENDIF
336
337!
338!-- Loop over all gridpoints along y
339    DO  j = nys, nyn
340   
341!
342!--    Store array to be advected on work array and add cyclic boundary along x
343       vad(:,:) = vad_in_out(:,j,:)
344       vad(:,-1)   = vad(:,nx)
345       vad(:,nx+1) = vad(:,0)
346
347!
348!--    Calculate right hand side
349       DO  i = 0, nx
350          DO  k = nzb+1, nzt
351             r(k,i) = 3.0 * (                                                  &
352                            spl_tri_x(2,i) * ( vad(k,i) - vad(k,i-1) ) * ddx + &
353                            spl_tri_x(3,i) * ( vad(k,i+1) - vad(k,i) ) * ddx   &
354                            )
355          ENDDO
356       ENDDO
357   
358!
359!--    Forward substitution
360       DO  k = nzb+1, nzt
361          wrk_spline(k,0) = r(k,0)
362       ENDDO
363
364       DO  i = 1, nx
365          DO  k = nzb+1, nzt
366             wrk_spline(k,i) = r(k,i) - spl_tri_x(5,i) * wrk_spline(k,i-1)
367          ENDDO
368       ENDDO
369
370!
371!--    Backward substitution (Sherman-Morrison-formula)
372       DO  k = nzb+1, nzt
373          r(k,nx) = wrk_spline(k,nx) / spl_tri_x(4,nx)
374       ENDDO
375
376       DO  k = nzb+1, nzt
377          DO  i = nx-1, 0, -1
378             r(k,i) = ( wrk_spline(k,i) - spl_tri_x(3,i) * r(k,i+1) ) / &
379                        spl_tri_x(4,i)
380          ENDDO
381          sm_faktor = ( r(k,0) + 0.5 * r(k,nx) / spl_gamma_x ) / &
382                      ( 1.0 + spl_z_x(0) + 0.5 * spl_z_x(nx) / spl_gamma_x )
383          DO  i = 0, nx
384             r(k,i) = r(k,i) - sm_faktor * spl_z_x(i)
385          ENDDO
386       ENDDO
387
388!
389!--    Add cyclic boundary to the right hand side
390       r(:,-1)   = r(:,nx)
391       r(:,nx+1) = r(:,0)
392
393!
394!--    Calculate advection along x
395       DO  i = 0, nx
396          DO  k = nzb+1, nzt
397
398             IF (ad_v(k,j,i) == 0.0)  THEN
399
400                vad_in_out(k,j,i) = vad(k,i)
401
402             ELSEIF ( ad_v(k,j,i) > 0.0)  THEN
403               
404                IF ( ABS( vad(k,i) - vad(k,i-1) ) <= ups_limit ) THEN
405                   vad_in_out(k,j,i) = vad(k,i) - dt_3d * ad_v(k,j,i) * &
406                                       ( vad(k,i) - vad(k,i-1) ) * ddx
407!
408!--                Calculate upstream fraction in % (s. flow_statistics)
409                   DO  sr = 0, statistic_regions
410                      sums_up_fraction_l(component,1,sr) = &
411                               sums_up_fraction_l(component,1,sr) + 1.0
412                   ENDDO
413                ELSE
414                   t1 = ad_v(k,j,i) * dt_3d * ddx
415                   t2 = 3.0 * ( vad(k,i-1) - vad(k,i) ) + &
416                        ( 2.0 * r(k,i) + r(k,i-1) ) * dx
417                   t3 = 2.0 * ( vad(k,i-1) - vad(k,i) ) + &
418                        ( r(k,i) + r(k,i-1) ) * dx
419                   vad_in_out(k,j,i) = vad(k,i) - r(k,i) * t1 * dx + &
420                                       t2 * t1**2 - t3 * t1**3
421                   IF ( vad(k,i-1) == vad(k,i) )  THEN
422                      vad_in_out(k,j,i) = vad(k,i)
423                   ENDIF
424                ENDIF
425
426             ELSE
427
428                IF ( ABS( vad(k,i) - vad(k,i+1) ) <= ups_limit ) THEN
429                   vad_in_out(k,j,i) = vad(k,i) - dt_3d * ad_v(k,j,i) * &
430                                       ( vad(k,i+1) - vad(k,i) ) * ddx
431!
432!--                Calculate upstream fraction in % (s. flow_statistics)
433                   DO  sr = 0, statistic_regions
434                      sums_up_fraction_l(component,1,sr) = &
435                               sums_up_fraction_l(component,1,sr) + 1.0
436                   ENDDO
437                ELSE
438                   t1 = -ad_v(k,j,i) * dt_3d * ddx
439                   t2 = 3.0 * ( vad(k,i) - vad(k,i+1) ) + &
440                        ( 2.0 * r(k,i) + r(k,i+1)) * dx
441                   t3 = 2.0 * ( vad(k,i) - vad(k,i+1) ) + &
442                        ( r(k,i) + r(k,i+1) ) * dx
443                   vad_in_out(k,j,i) = vad(k,i) + r(k,i) * t1 * dx - &
444                                       t2 * t1**2 + t3 * t1**3
445                   IF ( vad(k,i+1) == vad(k,i) )  THEN
446                      vad_in_out(k,j,i) = vad(k,i)
447                   ENDIF
448                ENDIF
449
450             ENDIF
451          ENDDO
452       ENDDO
453
454!
455!--    Limit values in order to prevent overshooting
456       IF ( cut_spline_overshoot )  THEN
457
458          DO  i = 0, nx
459             DO  k =  nzb+1, nzt
460                IF ( ad_v(k,j,i) > 0.0 ) THEN
461                   IF ( vad(k,i) > vad(k,i-1) )  THEN
462                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
463                                               vad(k,i)   + overshoot_limit )
464                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
465                                               vad(k,i-1) - overshoot_limit )
466                   ELSE
467                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
468                                               vad(k,i)   - overshoot_limit )
469                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
470                                               vad(k,i-1) + overshoot_limit )
471                   ENDIF
472                ELSE
473                   IF ( vad(k,i) > vad(k,i+1) )  THEN
474                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
475                                               vad(k,i)   + overshoot_limit )
476                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
477                                               vad(k,i+1) - overshoot_limit )
478                   ELSE
479                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
480                                               vad(k,i)   - overshoot_limit )
481                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
482                                               vad(k,i+1) + overshoot_limit )
483                   ENDIF
484                ENDIF
485             ENDDO
486          ENDDO
487
488       ENDIF
489
490!
491!--    Long filter (acting on tendency only)
492       IF ( long_filter_factor /= 0.0 )  THEN
493
494!
495!--       Compute tendency
496          DO  i = nxl, nxr
497             DO  k = nzb+1, nzt
498                tf(k,i) = vad_in_out(k,j,i) - vad(k,i)
499             ENDDO
500          ENDDO
501
502!
503!--       Apply the filter
504          wrk_long(:,0,1) = 2.0 * ( 1.0 + long_filter_factor )
505          wrk_long(:,0,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,0,1)
506          wrk_long(:,0,3) = ( long_filter_factor * tf(:,nx) + &
507                                   2.0 * tf(:,0) + tf(:,1) ) / wrk_long(:,0,1)
508
509          DO  i = 1, nx-1
510             DO  k = nzb+1, nzt
511                wrk_long(k,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
512                                  ( 1.0 - long_filter_factor ) *       &
513                                  wrk_long(k,i-1,2)
514                wrk_long(k,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(k,i,1)
515                wrk_long(k,i,3) = ( tf(k,i-1) + 2.0 * tf(k,i) +                &
516                                    tf(k,i+1) - ( 1.0 - long_filter_factor ) * &
517                                    wrk_long(k,i-1,3) ) / wrk_long(k,i,1)
518             ENDDO
519             wrk_long(:,nx,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
520                                      ( 1.0 - long_filter_factor ) * &
521                                        wrk_long(:,nx-1,2)
522             wrk_long(:,nx,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,nx,1)
523             wrk_long(:,nx,3) = ( tf(:,nx-1) + 2.0 * tf(:,nx) +    &
524                                    long_filter_factor * tf(:,0) - &
525                                    ( 1.0 - long_filter_factor ) * &
526                                    wrk_long(:,nx-1,3) ) / wrk_long(:,nx,1)
527             r(:,nx) = wrk_long(:,nx,3)
528          ENDDO
529          DO  i = nx-1, 0, -1
530             DO  k = nzb+1, nzt
531                r(k,i) = wrk_long(k,i,3) - wrk_long(k,i,2) * r(k,i+1)
532             ENDDO
533          ENDDO
534          DO  i = 0, nx
535             DO  k = nzb+1, nzt
536                vad_in_out(k,j,i) = vad(k,i) + r(k,i)
537             ENDDO
538          ENDDO
539
540       ENDIF  ! Long filter
541
542    ENDDO
543#endif
544
545    IF ( long_filter_factor /= 0.0 )  DEALLOCATE( tf, wrk_long )
546    DEALLOCATE( r, vad, wrk_spline )
547
548 END SUBROUTINE spline_x
Note: See TracBrowser for help on using the repository browser.