source: palm/trunk/SOURCE/spline_y.f90 @ 791

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

typo in file headers removed

  • Property svn:keywords set to Id
File size: 19.0 KB
Line 
1 SUBROUTINE spline_y( vad_in_out, ad_v, var_char )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: spline_y.f90 484 2010-02-05 07:36:54Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.9  2004/04/30 12:54:37  raasch
14! Names of transpose indices changed, enlarged transposition arrays introduced
15!
16! Revision 1.1  1999/02/05 09:16:31  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:nya,nxl_y:nxr_ya,nzb_y:nzt_ya), &
58             vad_in_out(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya)
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,2,:) = 0.0
91
92#if defined( __parallel )
93
94!
95!-- Allocate working arrays
96    ALLOCATE( r(-1:ny+1,nxl_y:nxr_y),   &
97              vad(-1:ny+1,nxl_y:nxr_y), &
98              wrk_spline(0:ny,nxl_y:nxr_y) )
99    IF ( long_filter_factor /= 0.0 )  THEN
100       ALLOCATE( tf(0:ny,nxl_y:nxr_y), &
101                 wrk_long(0:ny,nxl_y:nxr_y,1:3) )
102    ENDIF
103
104!
105!-- Loop over all gridpoints along z
106    DO  k = nzb_y, nzt_y
107
108!
109!--    Store array to be advected on work array and add cyclic boundary along y
110       vad(0:ny,nxl_y:nxr_y) = vad_in_out(0:ny,nxl_y:nxr_y,k)
111       vad(-1,:)   = vad(ny,:)
112       vad(ny+1,:) = vad(0,:)
113
114!
115!--    Calculate right hand side
116       DO  i = nxl_y, nxr_y
117          DO  j = 0, ny
118             r(j,i) = 3.0 * (                                                  &
119                            spl_tri_y(2,j) * ( vad(j,i) - vad(j-1,i) ) * ddy + &
120                            spl_tri_y(3,j) * ( vad(j+1,i) - vad(j,i) ) * ddy   &
121                            )
122          ENDDO
123       ENDDO
124   
125!
126!--    Forward substitution
127       DO  i = nxl_y, nxr_y
128          wrk_spline(0,i) = r(0,i)
129          DO  j = 1, ny
130             wrk_spline(j,i) = r(j,i) - spl_tri_y(5,j) * wrk_spline(j-1,i)
131          ENDDO
132       ENDDO
133   
134!
135!--    Backward substitution (sherman-Morrison-formula)
136       DO  i = nxl_y, nxr_y
137          r(ny,i) = wrk_spline(ny,i) / spl_tri_y(4,ny)
138          DO  j = ny-1, 0, -1
139             r(j,i) = ( wrk_spline(j,i) - spl_tri_y(3,j) * r(j+1,i) ) / &
140                        spl_tri_y(4,j)
141          ENDDO
142          sm_faktor = ( r(0,i) + 0.5 * r(ny,i) / spl_gamma_y ) / &
143                      ( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y )
144          DO  j = 0, ny
145             r(j,i) = r(j,i) - sm_faktor * spl_z_y(j)
146          ENDDO
147       ENDDO
148
149!
150!--    Add cyclic boundary conditions to right hand side
151       r(-1,:)   = r(ny,:)
152       r(ny+1,:) = r(0,:)
153
154!
155!--    Calculate advection along y
156       DO  i = nxl_y, nxr_y
157          DO  j = 0, ny
158
159             IF ( ad_v(j,i,k) == 0.0 )  THEN
160
161                vad_in_out(j,i,k) = vad(j,i)
162
163             ELSEIF ( ad_v(j,i,k) > 0.0 )  THEN
164
165                IF ( ABS( vad(j,i) - vad(j-1,i) ) <= ups_limit )  THEN
166                   vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * &
167                                       ( vad(j,i) - vad(j-1,i) ) * ddy
168!
169!--                Calculate upstream fraction in % (s. flow_statistics)
170                   DO  sr = 0, statistic_regions
171                      sums_up_fraction_l(component,2,sr) = &
172                               sums_up_fraction_l(component,2,sr) + 1.0
173                   ENDDO
174                ELSE
175                   t1 = ad_v(j,i,k) * dt_3d * ddy
176                   t2 = 3.0 * ( vad(j-1,i) - vad(j,i) ) + &
177                        ( 2.0 * r(j,i) + r(j-1,i) ) * dy
178                   t3 = 2.0 * ( vad(j-1,i) - vad(j,i) ) + &
179                              ( r(j,i) + r(j-1,i) ) * dy
180                   vad_in_out(j,i,k) = vad(j,i) - r(j,i) * t1 * dy + &
181                                       t2 * t1**2 - t3 * t1**3
182                   IF ( vad(j-1,i) == vad(j,i) )  THEN
183                      vad_in_out(j,i,k) = vad(j,i)
184                   ENDIF
185                ENDIF
186
187             ELSE
188
189                IF ( ABS( vad(j,i) - vad(j+1,i) ) <= ups_limit )  THEN
190                   vad_in_out(j,i,k) = vad(j,i) - dt_3d * ad_v(j,i,k) * &
191                                       ( vad(j+1,i) - vad(j,i) ) * ddy
192!
193!--                Calculate upstream fraction in % (s. flow_statistics)
194                   DO  sr = 0, statistic_regions
195                      sums_up_fraction_l(component,2,sr) = &
196                               sums_up_fraction_l(component,2,sr) + 1.0
197                   ENDDO
198                ELSE
199                   t1 = -ad_v(j,i,k) * dt_3d * ddy
200                   t2 = 3.0 * ( vad(j,i) - vad(j+1,i) ) + &
201                        ( 2.0 * r(j,i) + r(j+1,i) ) * dy
202                   t3 = 2.0 * ( vad(j,i) - vad(j+1,i) ) + &
203                        ( r(j,i) + r(j+1,i) ) * dy
204                   vad_in_out(j,i,k) = vad(j,i) + r(j,i) * t1 * dy - &
205                                       t2 * t1**2 + t3 * t1**3
206                   IF ( vad(j+1,i) == vad(j,i) )  THEN
207                      vad_in_out(j,i,k) = vad(j,i)
208                   ENDIF
209                ENDIF
210
211             ENDIF
212          ENDDO
213       ENDDO
214
215!
216!--    Limit values in order to prevent overshooting
217       IF ( cut_spline_overshoot )  THEN
218
219          DO  i = nxl_y, nxr_y
220             DO  j = 0, ny
221                IF ( ad_v(j,i,k) > 0.0 ) THEN
222                   IF ( vad(j,i) > vad(j-1,i) )  THEN
223                      vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
224                                               vad(j,i)   + overshoot_limit )
225                      vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
226                                               vad(j-1,i) - overshoot_limit )
227                   ELSE
228                      vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
229                                               vad(j,i)   - overshoot_limit )
230                      vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
231                                               vad(j-1,i) + overshoot_limit )
232                   ENDIF
233                ELSE
234                   IF ( vad(j,i) > vad(j+1,i) )  THEN
235                      vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
236                                               vad(j,i)   + overshoot_limit )
237                      vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
238                                               vad(j+1,i) - overshoot_limit )
239                   ELSE
240                      vad_in_out(j,i,k) = MAX( vad_in_out(j,i,k), &
241                                               vad(j,i)   - overshoot_limit )
242                      vad_in_out(j,i,k) = MIN( vad_in_out(j,i,k), &
243                                               vad(j+1,i) + overshoot_limit )
244                   ENDIF
245                ENDIF
246             ENDDO
247          ENDDO
248
249       ENDIF
250
251!
252!--    Long-filter (acting on tendency only)
253       IF ( long_filter_factor /= 0.0 )  THEN
254
255!
256!--       Compute tendency. Filter only acts on this quantity.
257          DO  i = nxl_y, nxr_y
258             DO  j = 0, ny
259                tf(j,i) = vad_in_out(j,i,k) - vad(j,i)
260             ENDDO
261          ENDDO
262
263!
264!--       Apply the filter.
265          DO  i = nxl_y, nxr_y
266             wrk_long(0,i,1) = 2.0 * ( 1.0 + long_filter_factor )
267             wrk_long(0,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(0,i,1)
268             wrk_long(0,i,3) = ( long_filter_factor * tf(ny,i) + &
269                                2.0 * tf(0,i) + tf(1,i)          &
270                               ) / wrk_long(0,i,1)
271             DO  j = 1, ny-1
272                wrk_long(j,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
273                                  ( 1.0 - long_filter_factor ) * wrk_long(j-1,i,2)
274                wrk_long(j,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(j,i,1)
275                wrk_long(j,i,3) = ( tf(j-1,i) + 2.0 * tf(j,i) +                &
276                                    tf(j+1,i) - ( 1.0 - long_filter_factor ) * &
277                                    wrk_long(j-1,i,3) ) / wrk_long(j,i,1)
278             ENDDO
279             wrk_long(ny,i,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
280                                ( 1.0 - long_filter_factor ) * wrk_long(ny-1,i,2)
281             wrk_long(ny,i,2) = ( 1.0 - long_filter_factor ) / wrk_long(ny,i,1)
282             wrk_long(ny,i,3) = ( tf(ny-1,i) + 2.0 * tf(ny,i) +  &
283                                  long_filter_factor * tf(0,i) - &
284                                  ( 1.0 - long_filter_factor ) * &
285                                  wrk_long(ny-1,i,3)             &
286                                ) / wrk_long(ny,i,1)
287             r(ny,i) = wrk_long(ny,i,3)
288          ENDDO
289
290          DO  j = ny-1, 0, -1
291             DO  i = nxl_y, nxr_y
292                r(j,i) = wrk_long(j,i,3) - wrk_long(j,i,2) * r(j+1,i)
293             ENDDO
294          ENDDO
295
296          DO  i = nxl_y, nxr_y
297             DO  j = 0, ny
298                vad_in_out(j,i,k) = vad(j,i) + r(j,i)
299             ENDDO
300          ENDDO
301
302       ENDIF  ! Long filter
303
304    ENDDO
305
306#else
307
308!
309!-- Allocate working arrays
310    ALLOCATE( r(nzb+1:nzt,nys-1:nyn+1), vad(nzb:nzt+1,nys-1:nyn+1), &
311              wrk_spline(nzb+1:nzt,nys-1:nyn+1) )
312    IF ( long_filter_factor /= 0.0 )  THEN
313       ALLOCATE( tf(nzb+1:nzt,nys-1:nyn+1), wrk_long(nzb+1:nzt,0:ny,1:3) )
314    ENDIF
315
316!
317!-- Loop over all gridpoints along x
318    DO  i = nxl, nxr
319
320!
321!--    Store array to be advected on work array and add cyclic boundary along x
322       vad(:,:)    = vad_in_out(:,:,i)
323       vad(:,-1)   = vad(:,ny)
324       vad(:,ny+1) = vad(:,0)
325
326!
327!--    Calculate right hand side
328       DO  j = 0, ny
329          DO  k = nzb+1, nzt
330             r(k,j) = 3.0 *  (                                                  &
331                             spl_tri_y(2,j) * ( vad(k,j) - vad(k,j-1) ) * ddy + &
332                             spl_tri_y(3,j) * ( vad(k,j+1) - vad(k,j) ) * ddy   &
333                             )
334          ENDDO
335       ENDDO
336   
337!
338!--    Forward substitution
339       DO  k = nzb+1, nzt
340          wrk_spline(k,0) = r(k,0)
341       ENDDO
342
343       DO  j = 1, ny
344          DO  k = nzb+1, nzt
345             wrk_spline(k,j) = r(k,j) - spl_tri_y(5,j) * wrk_spline(k,j-1)
346          ENDDO
347       ENDDO
348
349!
350!--    Backward substitution (Sherman-Morrison-formula)
351       DO  k = nzb+1, nzt
352          r(k,ny) = wrk_spline(k,ny) / spl_tri_y(4,ny)
353       ENDDO
354
355       DO  k = nzb+1, nzt
356          DO  j = ny-1, 0, -1
357             r(k,j) = ( wrk_spline(k,j) - spl_tri_y(3,j) * r(k,j+1) ) / &
358                        spl_tri_y(4,j)
359          ENDDO
360          sm_faktor = ( r(k,0) + 0.5 * r(k,ny) / spl_gamma_y ) / &
361                      ( 1.0 + spl_z_y(0) + 0.5 * spl_z_y(ny) / spl_gamma_y )
362          DO  j = 0, ny
363             r(k,j) = r(k,j) - sm_faktor * spl_z_y(j)
364          ENDDO
365       ENDDO
366
367!
368!--    Add cyclic boundary to the right hand side
369       r(:,-1)   = r(:,ny)
370       r(:,ny+1) = r(:,0)
371   
372!
373!--    Calculate advection along y
374       DO  j = 0, ny
375          DO  k = nzb+1, nzt
376
377             IF ( ad_v(k,j,i) == 0.0 )  THEN
378
379                vad_in_out(k,j,i) = vad(k,j)
380
381             ELSEIF ( ad_v(k,j,i) > 0.0 )  THEN
382
383                IF ( ABS( vad(k,j) - vad(k,j-1) ) <= ups_limit ) THEN
384                   vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
385                                       ( vad(k,j) - vad(k,j-1) ) * ddy
386!
387!--                Calculate upstream fraction in % (s. flow_statistics)
388                   DO  sr = 0, statistic_regions
389                      sums_up_fraction_l(component,2,sr) = &
390                               sums_up_fraction_l(component,2,sr) + 1.0
391                   ENDDO
392                ELSE
393                   t1 = ad_v(k,j,i) * dt_3d * ddy
394                   t2 = 3.0 * ( vad(k,j-1) - vad(k,j) ) + &
395                        ( 2.0 * r(k,j) + r(k,j-1) ) * dy
396                   t3 = 2.0 * ( vad(k,j-1) - vad(k,j) ) + &
397                        ( r(k,j) + r(k,j-1) ) * dy
398                   vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1 * dy + &
399                                       t2 * t1**2 - t3 * t1**3
400                   IF ( vad(k,j-1) == vad(k,j) )  THEN
401                      vad_in_out(k,j,i) = vad(k,j)
402                   ENDIF
403                ENDIF
404
405             ELSE
406               
407                IF ( ABS( vad(k,j) - vad(k,j+1) ) <= ups_limit ) THEN
408                   vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
409                                       ( vad(k,j+1) - vad(k,j) ) * ddy
410!
411!--                Calculate upstream fraction in % (s. flow_statistics)
412                   DO  sr = 0, statistic_regions
413                      sums_up_fraction_l(component,2,sr) = &
414                               sums_up_fraction_l(component,2,sr) + 1.0
415                   ENDDO
416                ELSE
417                   t1 = -ad_v(k,j,i) * dt_3d * ddy
418                   t2 = 3.0 * ( vad(k,j) - vad(k,j+1) ) + &
419                        ( 2.0 * r(k,j) + r(k,j+1) ) * dy
420                   t3 = 2.0 * ( vad(k,j) - vad(k,j+1) ) + &
421                        ( r(k,j) + r(k,j+1) ) * dy
422                   vad_in_out(k,j,i) = vad(k,j) + r(k,j) * t1 * dy - &
423                                       t2 * t1**2 + t3 * t1**3
424                   IF ( vad(k,j+1) == vad(k,j) )  THEN
425                      vad_in_out(k,j,i) = vad(k,j)
426                   ENDIF
427                ENDIF
428
429             ENDIF
430          ENDDO
431       ENDDO
432
433!
434!--    Limit values in order to prevent overshooting
435       IF ( cut_spline_overshoot )  THEN
436
437          DO  j = 0, ny
438             DO  k = nzb+1, nzt
439                IF ( ad_v(k,j,i) > 0.0 ) THEN
440                   IF ( vad(k,j) > vad(k,j-1) )  THEN
441                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
442                                               vad(k,j)   + overshoot_limit )
443                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
444                                               vad(k,j-1) - overshoot_limit )
445                   ELSE
446                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
447                                               vad(k,j)   - overshoot_limit )
448                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
449                                               vad(k,j-1) + overshoot_limit )
450                   ENDIF
451                ELSE
452                   IF ( vad(k,j) > vad(k,j+1) )  THEN
453                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
454                                               vad(k,j)   + overshoot_limit )
455                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
456                                               vad(k,j+1) - overshoot_limit )
457                   ELSE
458                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
459                                               vad(k,j)   - overshoot_limit )
460                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
461                                               vad(k,j+1) + overshoot_limit )
462                   ENDIF
463                ENDIF
464             ENDDO
465          ENDDO
466
467       ENDIF
468
469!
470!--    Long filter (acting on tendency only)
471       IF ( long_filter_factor /= 0.0 )  THEN
472
473!
474!--       Compute tendency
475          DO  j = nys, nyn
476             DO  k = nzb+1, nzt
477                tf(k,j) = vad_in_out(k,j,i) - vad(k,j)
478             ENDDO
479          ENDDO
480
481!
482!--       Apply the filter
483          wrk_long(:,0,1) = 2.0 * ( 1.0 + long_filter_factor )
484          wrk_long(:,0,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,0,1)
485          wrk_long(:,0,3) = ( long_filter_factor * tf(:,ny) + &
486                                   2.0 * tf(:,0) + tf(:,1) ) / wrk_long(:,0,1)
487
488          DO  j = 1, ny-1
489             DO  k = nzb+1, nzt
490                wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
491                                  ( 1.0 - long_filter_factor ) *       &
492                                  wrk_long(k,j-1,2)
493                wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) / wrk_long(k,j,1)
494                wrk_long(k,j,3) = ( tf(k,j-1) + 2.0 * tf(k,j) +                &
495                                    tf(k,j+1) - ( 1.0 - long_filter_factor ) * &
496                                    wrk_long(k,j-1,3) ) / wrk_long(k,j,1)
497             ENDDO
498             wrk_long(:,ny,1) = 2.0 * ( 1.0 + long_filter_factor ) - &
499                                      ( 1.0 - long_filter_factor ) * &
500                                        wrk_long(:,ny-1,2)
501             wrk_long(:,ny,2) = ( 1.0 - long_filter_factor ) / wrk_long(:,ny,1)
502             wrk_long(:,ny,3) = ( tf(:,ny-1) + 2.0 * tf(:,ny) +    &
503                                    long_filter_factor * tf(:,0) - &
504                                    ( 1.0 - long_filter_factor ) * &
505                                    wrk_long(:,ny-1,3) ) / wrk_long(:,ny,1)
506             r(:,ny) = wrk_long(:,ny,3)
507          ENDDO
508          DO  j = ny-1, 0, -1
509             DO  k = nzb+1, nzt
510                r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k,j+1)
511             ENDDO
512          ENDDO
513          DO  j = 0, ny
514             DO  k = nzb+1, nzt
515                vad_in_out(k,j,i) = vad(k,j) + r(k,j)
516             ENDDO
517          ENDDO
518
519       ENDIF  ! Long filter
520
521    ENDDO
522#endif
523
524    IF ( long_filter_factor /= 0.0 )  DEALLOCATE( tf, wrk_long )
525    DEALLOCATE( r, vad, wrk_spline )
526
527 END SUBROUTINE spline_y
Note: See TracBrowser for help on using the repository browser.