Home

Context Navigation

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

Last change on this file since 1 was 1, checked in by raasch, 18 years ago
Initial repository layout and content
File size: 19.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |