Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |