Home

Context Navigation

source: palm/trunk/SOURCE/spline_y.f90 @ 806

Last change on this file since 806 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 19.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |