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

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

typo in file headers removed

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