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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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