source: palm/trunk/SOURCE/spline_z.f90 @ 29

Last change on this file since 29 was 19, checked in by raasch, 18 years ago

preliminary version of modified boundary conditions at top

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