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