source: palm/trunk/SOURCE/spline_z.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: 15.9 KB
Line 
1 SUBROUTINE spline_z( vad_in_out, ad_v, dz_spline, spline_tri, var_char )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: spline_z.f90,v $
11! Revision 1.9  2005/06/29 08:22:56  steinfeld
12! Dependency of ug and vg on height considered in the determination of the
13! upper boundary condition for vad
14!
15! Revision 1.8  2004/04/30 12:55:23  raasch
16! Enlarged transposition arrays introduced
17!
18! Revision 1.7  2003/03/16 09:49:33  raasch
19! Two underscores (_) are placed in front of all define-strings
20!
21! Revision 1.6  2001/03/30 07:54:07  raasch
22! Arrays r vad, and wrk_spline changed from 3D to 2D and removed from argument
23! list. Several loops over i combined to one loop.
24! Application of long filter moved to this routine. All comments and identifiers
25! translated into English.
26!
27! Revision 1.5  2001/01/22 08:10:21  raasch
28! Module test_variables removed
29!
30! Revision 1.4  1999/03/25 07:34:55  raasch
31! Filterung der Ueberschwinger geschieht optional, ups_limit_e eingefuehrt,
32! Ueberschwinger werden in gewissen Grenzen erlaubt
33!
34! Revision 1.3  1999/02/26 17:55:41  schroeter
35! - Gradientenkontrolle fuer den nicht-parallelen Teil
36! - statistische Auswertung ueber den prozentualen Anteil des
37!   Upstream-Verfahrens an der Gesamtadvektion fuer nicht-
38!   parallelen Teil
39!
40! Revision 1.2  1999/02/17 09:32:55  raasch
41! Wertebegrenzung zur Verhinderung von Ueberschwingern
42!
43! Revision 1.1  1999/02/05 09:17:16  raasch
44! Initial revision
45!
46!
47! Description:
48! ------------
49! Upstream-spline advection along x
50!
51! Input/output parameters:
52! ad_v       = advecting wind speed component
53! dz_spline  = vertical grid spacing (dzu or dzw, depending on quantity to be
54!              advected)
55! spline_tri = grid spacing factors (spl_tri_zu or spl_tri_zw, depending on
56!              quantity to be advected)
57! vad_in_out = quantity to be advected, excluding ghost- or cyclic boundaries
58!              result is given to the calling routine in this array
59! var_char   = string which defines the quantity to be advected
60!
61! Internal arrays:
62! r          = 2D-working array (right hand side of linear equation, buffer for
63!              Long filter)
64! tf         = tendency field (2D), used for long filter
65! vad        = quantity to be advected (2D), including ghost- or cyclic
66!              boundarys along the direction of advection
67! wrk_long   = working array (long coefficients)
68! wrk_spline = working array (spline coefficients)
69!------------------------------------------------------------------------------!
70
71    USE arrays_3d
72    USE grid_variables
73    USE indices
74    USE statistics
75    USE control_parameters
76    USE transpose_indices
77
78    IMPLICIT NONE
79
80    CHARACTER (LEN=*) ::  var_char
81
82    INTEGER ::  component, i, j, k, sr
83    REAL ::     dzwd, dzwu, overshoot_limit, t1, t2, t3, ups_limit
84    REAL ::     dz_spline(1:nzt+1)
85    REAL ::     spline_tri(5,nzb:nzt+1)
86    REAL ::     ad_v(nzb+1:nzta,nys:nyna,nxl:nxra)
87
88    REAL, DIMENSION(:,:), ALLOCATABLE   ::  r, tf, vad, wrk_spline
89    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  wrk_long
90
91#if defined( __parallel )
92    REAL ::     vad_in_out(nzb+1:nzta,nys:nyna,nxl:nxra)
93#else
94    REAL ::     vad_in_out(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
95#endif
96
97!
98!-- Set criteria for switching between upstream- and upstream-spline-method
99    IF ( var_char == 'u' )  THEN
100       overshoot_limit = overshoot_limit_u
101       ups_limit       = ups_limit_u
102       component       = 1
103    ELSEIF (  var_char == 'v' )  THEN
104       overshoot_limit = overshoot_limit_v
105       ups_limit       = ups_limit_v
106       component       = 2
107    ELSEIF (  var_char == 'w' )  THEN
108       overshoot_limit = overshoot_limit_w
109       ups_limit       = ups_limit_w
110       component       = 3
111    ELSEIF (  var_char == 'pt' )  THEN
112       overshoot_limit = overshoot_limit_pt
113       ups_limit       = ups_limit_pt
114       component       = 4
115    ELSEIF (  var_char == 'e' )  THEN
116       overshoot_limit = overshoot_limit_e
117       ups_limit       = ups_limit_e
118       component       = 5
119    ENDIF
120
121!
122!-- Allocate working arrays
123    ALLOCATE( r(nzb:nzt+1,nys:nyn), vad(nzb:nzt+1,nys:nyn), &
124              wrk_spline(nzb:nzt+1,nys:nyn) )
125    IF ( long_filter_factor /= 0.0 )  THEN
126       ALLOCATE( tf(nzb:nzt+1,nys:nyn), wrk_long(nzb+1:nzt,nys:nyn,1:3) )
127    ENDIF   
128
129!
130!-- Initialize calculation of relative upstream fraction
131    sums_up_fraction_l(component,3,:) = 0.0
132
133!
134!-- Loop over all gridpoints along x
135    DO  i = nxl, nxr
136
137!
138!--    Store array to be advected on work array
139       vad(nzb+1:nzt,:) = vad_in_out(nzb+1:nzt,nys:nyn,i)
140!
141!--    Add boundary conditions along z
142       IF ( var_char == 'u'  .OR.  var_char == 'v' )  THEN
143!
144!--       Bottom boundary
145!--       u- and v-component
146          IF ( ibc_uv_b == 0 )  THEN
147             vad(nzb,:) = -vad(nzb+1,:)
148          ELSE
149             vad(nzb,:) = vad(nzb+1,:)
150          ENDIF
151!
152!--       Top boundary
153!--       Dirichlet condition
154          IF ( ibc_uv_t == 0  .AND.  var_char == 'u' )  THEN
155!
156!--          u-component
157             vad(nzt+1,:) = ug(nzt+1)
158          ELSEIF ( ibc_uv_t == 0  .AND.  var_char == 'v' )  THEN
159!
160!--          v-component
161             vad(nzt+1,:) = vg(nzt+1)
162          ELSE
163!
164!--          Neumann condition
165             vad(nzt+1,:) = vad(nzt,:)
166          ENDIF
167
168       ELSEIF ( var_char == 'w' )  THEN
169!
170!--       Bottom and top boundary for w-component
171          vad(nzb,:)   = 0.0
172          vad(nzt+1,:) = 0.0
173
174       ELSEIF ( var_char == 'pt' )  THEN
175!
176!--       Bottom boundary for temperature
177          IF ( ibc_pt_b == 1 )  THEN
178             vad(nzb,:) = vad(nzb+1,:)
179          ELSE
180             vad(nzb,:) = pt(nzb,:,i)
181          ENDIF
182!
183!--       Top boundary for temperature
184          IF ( ibc_pt_t == 1 )  THEN
185             vad(nzt,:)   = vad(nzt-1,:) + bc_pt_t_val * dz_spline(nzt)
186             vad(nzt+1,:) = vad(nzt,:)   + bc_pt_t_val * dz_spline(nzt+1)
187          ELSE
188             vad(nzt,:)   = pt(nzt,nys:nyn,i)
189             vad(nzt+1,:) = pt(nzt+1,nys:nyn,i)
190          ENDIF
191
192       ELSEIF ( var_char == 'e' )  THEN
193!
194!--       Boundary conditions for TKE (Neumann in any case)
195          vad(nzb,:)   = vad(nzb+1,:)
196          vad(nzt,:)   = vad(nzt-1,:)
197          vad(nzt+1,:) = vad(nzt,:)
198
199       ENDIF
200
201!
202!--    Calculate right hand side
203       DO  j = nys, nyn
204          r(nzb,j)   = 3.0 * ( vad(nzb+1,j)-vad(nzb,j) ) / dz_spline(1)
205          r(nzt+1,j) = 3.0 * ( vad(nzt+1,j)-vad(nzt,j) ) / dz_spline(nzt+1)
206          DO  k = nzb+1, nzt
207             r(k,j) = 3.0 * (                                                   &
208                     spline_tri(2,k) * ( vad(k,j)-vad(k-1,j) ) / dz_spline(k)   &
209                   + spline_tri(3,k) * ( vad(k+1,j)-vad(k,j) ) / dz_spline(k+1) &
210                            )
211          ENDDO
212       ENDDO
213   
214!
215!--    Forward substitution
216       DO  j = nys, nyn
217          wrk_spline(nzb,j) = r(nzb,j)
218          DO  k = nzb+1, nzt+1
219             wrk_spline(k,j) = r(k,j) - spline_tri(5,k) * r(k-1,j)
220          ENDDO
221       ENDDO
222
223!
224!--    Backward substitution
225       DO  j = nys, nyn
226          r(nzt+1,j) = wrk_spline(nzt+1,j) / spline_tri(4,nzt+1)
227          DO  k = nzt, nzb, -1
228             r(k,j) = ( wrk_spline(k,j) - spline_tri(3,k) * r(k+1,j) ) / &
229                      spline_tri(4,k)
230          ENDDO
231       ENDDO
232   
233!
234!--    Calculate advection along z
235       DO  j = nys, nyn
236          DO  k = nzb+1, nzt
237
238             IF ( ad_v(k,j,i) == 0.0 ) THEN
239
240                vad_in_out(k,j,i) = vad(k,j)
241
242             ELSEIF ( ad_v(k,j,i) > 0.0 ) THEN
243
244                IF ( ABS( vad(k,j) - vad(k-1,j) ) <= ups_limit )  THEN
245                   vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
246                                       ( vad(k,j) - vad(k-1,j) ) * ddzu(k)
247!
248!--                Calculate upstream fraction in % (s. flow_statistics)
249                   DO  sr = 0, statistic_regions
250                      sums_up_fraction_l(component,3,sr) = &
251                               sums_up_fraction_l(component,3,sr) + 1.0
252                   ENDDO
253                ELSE
254                   t1 = ad_v(k,j,i) * dt_3d / dz_spline(k)
255                   t2 = 3.0 * ( vad(k-1,j) - vad(k,j) ) + &
256                        ( 2.0 * r(k,j) + r(k-1,j) ) * dz_spline(k)
257                   t3 = 2.0 * ( vad(k-1,j) - vad(k,j) ) + &
258                        ( r(k,j) + r(k-1,j) ) * dz_spline(k)
259                   vad_in_out(k,j,i) = vad(k,j) - r(k,j) * t1* dz_spline(k) + &
260                                       t2 * t1**2 - t3 * t1**3
261                   IF ( vad(k-1,j) == vad(k,j) )  THEN
262                      vad_in_out(k,j,i) = vad(k,j)
263                   ENDIF
264                ENDIF
265
266             ELSE
267
268                IF( ABS( vad(k,j) - vad(k+1,j) ) <= ups_limit )  THEN
269                   vad_in_out(k,j,i) = vad(k,j) - dt_3d * ad_v(k,j,i) * &
270                                       ( vad(k+1,j) - vad(k,j) ) * ddzu(k+1)
271!
272!--                Calculate upstream fraction in % (s. flow_statistics)
273                   DO  sr = 0, statistic_regions
274                      sums_up_fraction_l(component,3,sr) = &
275                               sums_up_fraction_l(component,3,sr) + 1.0
276                   ENDDO
277                ELSE
278                   t1 = -ad_v(k,j,i) * dt_3d / dz_spline(k+1)
279                   t2 = 3.0 * ( vad(k,j) - vad(k+1,j) ) + &
280                        ( 2.0 * r(k,j) + r(k+1,j) ) * dz_spline(k+1)
281                   t3 = 2.0 * ( vad(k,j) - vad(k+1,j) ) + &
282                        ( r(k,j) + r(k+1,j) ) * dz_spline(k+1)
283                   vad_in_out(k,j,i) = vad(k,j) + r(k,j)*t1*dz_spline(k+1) - &
284                                       t2 * t1**2 + t3 * t1**3
285                   IF ( vad(k+1,j) == vad(k,j) )  THEN
286                      vad_in_out(k,j,i) = vad(k,j)
287                   ENDIF
288                ENDIF
289
290             ENDIF
291          ENDDO
292       ENDDO
293
294!
295!--    Limit values in order to prevent overshooting
296       IF ( cut_spline_overshoot )  THEN
297
298          DO  j = nys, nyn
299             DO  k = nzb+1, nzt
300                IF ( ad_v(k,j,i) > 0.0 ) THEN
301                   IF ( vad(k,j) > vad(k-1,j) )  THEN
302                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
303                                               vad(k,j)   + overshoot_limit )
304                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
305                                               vad(k-1,j) - overshoot_limit )
306                   ELSE
307                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
308                                               vad(k,j)   - overshoot_limit )
309                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
310                                               vad(k-1,j) + overshoot_limit )
311                   ENDIF
312                ELSE
313                   IF ( vad(k,j) > vad(k+1,j) )  THEN
314                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
315                                               vad(k,j)   + overshoot_limit )
316                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
317                                               vad(k+1,j) - overshoot_limit )
318                   ELSE
319                      vad_in_out(k,j,i) = MAX( vad_in_out(k,j,i), &
320                                               vad(k,j)   - overshoot_limit )
321                      vad_in_out(k,j,i) = MIN( vad_in_out(k,j,i), &
322                                               vad(k+1,j) + overshoot_limit )
323                   ENDIF
324                ENDIF
325             ENDDO
326          ENDDO
327
328       ENDIF
329
330!
331!--    Long-filter (acting on tendency only)
332       IF ( long_filter_factor /= 0.0 )  THEN
333
334!
335!--       Compute tendency
336          DO  j = nys, nyn
337
338!
339!--          Depending on the quantity to be advected, the respective vertical
340!--          boundary conditions must be applied.
341             IF ( var_char == 'u'  .OR.  var_char == 'v' ) THEN
342
343                IF ( ibc_uv_b == 0 )  THEN
344                   tf(nzb,j) = - ( vad_in_out(nzb+1,j,i) - vad(nzb+1,j) )
345                ELSE
346                   tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j)
347                ENDIF
348
349                IF ( ibc_uv_t == 0 )  THEN
350                   tf(nzt+1,j) = 0.0
351                ELSE
352                   tf(nzt+1,j) = vad_in_out(nzt,j,i) - vad(nzt,j)
353                ENDIF
354
355             ELSEIF ( var_char == 'w' )  THEN
356
357                tf(nzb,j)   = 0.0
358                tf(nzt+1,j) = 0.0
359
360             ELSEIF ( var_char == 'pt' )  THEN
361
362                IF ( ibc_pt_b == 1 )  THEN
363                   tf(nzb,j) = vad_in_out(nzb+1,j,i) - vad(nzb+1,j)
364                ELSE
365                   tf(nzb,j) = 0.0
366                ENDIF
367
368                IF ( ibc_pt_t == 1 )  THEN
369                   vad_in_out(nzt,j,i) = vad_in_out(nzt-1,j,i) + bc_pt_t_val * &
370                                         dz_spline(nzt)
371                   tf(nzt+1,j) = vad_in_out(nzt,j,i) + bc_pt_t_val * &
372                                         dz_spline(nzt+1) - vad(nzt+1,j)
373                ELSE
374                   vad_in_out(nzt,j,i) = pt(nzt,j,i)
375                   tf(nzt+1,j) = 0.0
376                ENDIF
377
378             ENDIF
379
380             DO  k = nzb+1, nzt
381                tf(k,j) = vad_in_out(k,j,i) - vad(k,j)
382             ENDDO
383
384          ENDDO
385
386!
387!--       Apply the filter.
388          DO  j = nys, nyn
389
390             dzwd = dz_spline(1) / ( dz_spline(1) + dz_spline(2) )
391             dzwu = dz_spline(2) / ( dz_spline(1) + dz_spline(2) )
392
393             wrk_long(nzb+1,j,1) = 2.0 * ( 1.0 + long_filter_factor )
394             wrk_long(nzb+1,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
395                                     wrk_long(nzb+1,j,1)
396             wrk_long(nzb+1,j,3) = ( long_filter_factor * dzwu * tf(nzb,j) +  &
397                                     2.0 * tf(nzb+1,j) + dzwd * tf(nzb+2,j) &
398                                   ) / wrk_long(nzb+1,j,1)
399
400             DO  k = nzb+2, nzt-1
401
402                dzwd = dz_spline(k)   / ( dz_spline(k) + dz_spline(k+1) )
403                dzwu = dz_spline(k+1) / ( dz_spline(k) + dz_spline(k+1) )
404
405                wrk_long(k,j,1) = 2.0 * ( 1.0 + long_filter_factor ) -  &
406                                  ( 1.0 - long_filter_factor ) * dzwu * &
407                                  wrk_long(k-1,j,2)
408                wrk_long(k,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
409                                  wrk_long(k,j,1) 
410                wrk_long(k,j,3) = ( dzwu * tf(k-1,j) + 2.0 * tf(k,j) +   &
411                                    dzwd * tf(k+1,j) -                   &
412                                   ( 1.0 - long_filter_factor ) * dzwu * &
413                                  wrk_long(k-1,j,3)                      &
414                                  ) / wrk_long(k,j,1)
415             ENDDO
416
417             dzwd = dz_spline(nzt)   / ( dz_spline(nzt) + dz_spline(nzt+1) )
418             dzwu = dz_spline(nzt+1) / ( dz_spline(nzt) + dz_spline(nzt+1) )
419
420             wrk_long(nzt,j,1) = 2.0 * ( 1.0 + long_filter_factor ) -  &
421                                 ( 1.0 - long_filter_factor ) * dzwu * &
422                                 wrk_long(nzt-1,j,2)
423             wrk_long(nzt,j,2) = ( 1.0 - long_filter_factor ) * dzwd / &
424                                 wrk_long(nzt,j,1)
425             wrk_long(nzt,j,3) = ( dzwu * tf(nzt-1,j) + 2.0 * tf(nzt,j) +    & 
426                                   dzwd * long_filter_factor * tf(nzt+1,j) - &
427                                   ( 1.0 - long_filter_factor ) * dzwu *     &
428                                   wrk_long(nzt-1,j,3)                       &
429                                 ) / wrk_long(nzt,j,1)
430             r(nzt,j)     = wrk_long(nzt,j,3)
431
432          ENDDO
433
434          DO  j = nys, nyn
435             DO  k = nzt-1, nzb+1, -1
436                r(k,j) = wrk_long(k,j,3) - wrk_long(k,j,2) * r(k+1,j)
437             ENDDO
438          ENDDO
439
440          DO  j = nys, nyn
441             DO  k = nzb+1, nzt
442                vad_in_out(k,j,i) = vad(k,j) + r(k,j)
443             ENDDO
444          ENDDO
445
446       ENDIF  ! Long filter
447
448    ENDDO
449
450    DEALLOCATE( r, vad, wrk_spline )
451    IF ( long_filter_factor /= 0.0 )  DEALLOCATE( tf, wrk_long )   
452
453 END SUBROUTINE spline_z
Note: See TracBrowser for help on using the repository browser.