Home

Context Navigation

source: palm/trunk/SOURCE/advec_s_bc.f90 @ 16

Last change on this file since 16 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: 42.6 KB

Line
1	SUBROUTINE advec_s_bc( sk, sk_char )
2
3	!-------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: advec_s_bc.f90 4 2007-02-13 11:33:16Z raasch $
11	! RCS Log replace by Id keyword, revision history cleaned up
12	!
13	! Revision 1.22 2006/02/23 09:42:08 raasch
14	! anz renamed ngp
15	!
16	! Revision 1.1 1997/08/29 08:53:46 raasch
17	! Initial revision
18	!
19	!
20	! Description:
21	! ------------
22	! Advection term for scalar quantities using the Bott-Chlond scheme.
23	! Computation in individual steps for each of the three dimensions.
24	! Limiting assumptions:
25	! So far the scheme has been assuming equidistant grid spacing. As this is not
26	! the case in the stretched portion of the z-direction, there dzw(k) is used as
27	! a substitute for a constant grid length. This certainly causes incorrect
28	! results; however, it is hoped that they are not too apparent for weakly
29	! stretched grids.
30	! NOTE: This is a provisional, non-optimised version!
31	!-------------------------------------------------------------------------------!
32
33	USE advection
34	USE arrays_3d
35	USE cpulog
36	USE grid_variables
37	USE indices
38	USE interfaces
39	USE pegrid
40	USE statistics
41	USE control_parameters
42
43	IMPLICIT NONE
44
45	CHARACTER (LEN=*) :: sk_char
46
47	INTEGER :: i, ix, j, k, ngp, sr, type_xz_2
48
49	REAL :: cim, cimf, cip, cipf, d_new, ffmax, fminus, fplus, f2, f4, f8, &
50	f12, f24, f48, f1920, im, ip, m2, m3, nenner, snenn, sterm, &
51	tendenz, t1, t2, zaehler
52	REAL :: fmax(2), fmax_l(2)
53	REAL, DIMENSION(:,:,:), POINTER :: sk
54
55	REAL, DIMENSION(:,:), ALLOCATABLE :: a0, a1, a12, a2, a22, immb, imme, &
56	impb, impe, ipmb, ipme, ippb, ippe
57	REAL, DIMENSION(:,:,:), ALLOCATABLE :: sk_p
58
59	#if defined( __t3eb ) \|\| defined( __t3eh ) \|\| defined( __t3ej2 ) \|\| defined( __t3ej5 )
60	REAL (kind=4) :: m1n, m1z !Wichtig: Division
61	REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw
62	#else
63	REAL :: m1n, m1z
64	REAL, DIMENSION(:,:), ALLOCATABLE :: m1, sw
65	#endif
66
67
68	!
69	!-- Array sk_p requires 2 extra elements for each dimension
70	ALLOCATE( sk_p(nzb-2:nzt+2,nys-3:nyn+3,nxl-3:nxr+3) )
71	sk_p = 0.0
72
73	!
74	!-- Assign reciprocal values in order to avoid divisions later
75	f2 = 0.5
76	f4 = 0.25
77	f8 = 0.125
78	f12 = 0.8333333333333333E-01
79	f24 = 0.4166666666666666E-01
80	f48 = 0.2083333333333333E-01
81	f1920 = 0.5208333333333333E-03
82
83	!
84	!-- Advection in x-direction:
85
86	!
87	!-- Save the quantity to be advected in a local array
88	!-- add an enlarged boundary in x-direction
89	DO i = nxl-1, nxr+1
90	DO j = nys, nyn
91	DO k = nzb+1, nzt-1
92	sk_p(k,j,i) = sk(k,j,i)
93	ENDDO
94	ENDDO
95	ENDDO
96	#if defined( __parallel )
97	ngp = 2 * ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
98	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
99	!
100	!-- Send left boundary, receive right boundary
101	CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, &
102	sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, &
103	comm2d, status, ierr )
104	!
105	!-- Send right boundary, receive left boundary
106	CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, &
107	sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, &
108	comm2d, status, ierr )
109	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
110	#else
111
112	!
113	!-- Cyclic boundary conditions
114	sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)
115	sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)
116	sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)
117	sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)
118	#endif
119
120	!
121	!-- In case of a sloping surface, the additional gridpoints in x-direction
122	!-- of the temperature field at the left and right boundary of the total
123	!-- domain must be adjusted by the temperature difference between this distance
124	IF ( sloping_surface .AND. sk_char == 'pt' ) THEN
125	IF ( nxl == 0 ) THEN
126	sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset
127	sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset
128	ENDIF
129	IF ( nxr == nx ) THEN
130	sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset
131	sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset
132	ENDIF
133	ENDIF
134
135	!
136	!-- Initialise control density
137	d = 0.0
138
139	!
140	!-- Determine maxima of the first and second derivative in x-direction
141	fmax_l = 0.0
142	DO i = nxl, nxr
143	DO j = nys, nyn
144	DO k = nzb+1, nzt-1
145	zaehler = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
146	nenner = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
147	fmax_l(1) = MAX( fmax_l(1) , zaehler )
148	fmax_l(2) = MAX( fmax_l(2) , nenner )
149	ENDDO
150	ENDDO
151	ENDDO
152	#if defined( __parallel )
153	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
154	#else
155	fmax = fmax_l
156	#endif
157
158	fmax = 0.04 * fmax
159
160	!
161	!-- Allocate temporary arrays
162	ALLOCATE( a0(nzb+1:nzt-1,nxl-1:nxr+1), a1(nzb+1:nzt-1,nxl-1:nxr+1), &
163	a2(nzb+1:nzt-1,nxl-1:nxr+1), a12(nzb+1:nzt-1,nxl-1:nxr+1), &
164	a22(nzb+1:nzt-1,nxl-1:nxr+1), immb(nzb+1:nzt-1,nxl-1:nxr+1), &
165	imme(nzb+1:nzt-1,nxl-1:nxr+1), impb(nzb+1:nzt-1,nxl-1:nxr+1), &
166	impe(nzb+1:nzt-1,nxl-1:nxr+1), ipmb(nzb+1:nzt-1,nxl-1:nxr+1), &
167	ipme(nzb+1:nzt-1,nxl-1:nxr+1), ippb(nzb+1:nzt-1,nxl-1:nxr+1), &
168	ippe(nzb+1:nzt-1,nxl-1:nxr+1), m1(nzb+1:nzt-1,nxl-2:nxr+2), &
169	sw(nzb+1:nzt-1,nxl-1:nxr+1) &
170	)
171	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
172
173	!
174	!-- Initialise point of time measuring of the exponential portion (this would
175	!-- not work if done locally within the loop)
176	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
177	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
178
179	!
180	!-- Outer loop of all j
181	DO j = nys, nyn
182
183	!
184	!-- Compute polynomial coefficients
185	DO i = nxl-1, nxr+1
186	DO k = nzb+1, nzt-1
187	a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
188	a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) &
189	+ sk_p(k,j,i-1) )
190	a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1) &
191	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1) &
192	+ 9.0 * sk_p(k,j,i-2) ) * f1920
193	a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1) &
194	- 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2) &
195	) * f48
196	a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1) &
197	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1) &
198	- 3.0 * sk_p(k,j,i-2) ) * f48
199	ENDDO
200	ENDDO
201
202	!
203	!-- Fluxes using the Bott scheme
204	!-- *VOCL LOOP,UNROLL(2)
205	DO i = nxl, nxr
206	DO k = nzb+1, nzt-1
207	cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
208	cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
209	cipf = 1.0 - 2.0 * cip
210	cimf = 1.0 - 2.0 * cim
211	ip = a0(k,i) * f2 * ( 1.0 - cipf ) &
212	+ a1(k,i) * f8 * ( 1.0 - cipf*cipf ) &
213	+ a2(k,i) * f24 * ( 1.0 - cipfcipfcipf )
214	im = a0(k,i+1) * f2 * ( 1.0 - cimf ) &
215	- a1(k,i+1) * f8 * ( 1.0 - cimf*cimf ) &
216	+ a2(k,i+1) * f24 * ( 1.0 - cimfcimfcimf )
217	ip = MAX( ip, 0.0 )
218	im = MAX( im, 0.0 )
219	ippb(k,i) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
220	impb(k,i) = im * MIN( 1.0, sk_p(k,j,i+1) / (ip+im+1E-15) )
221
222	cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
223	cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
224	cipf = 1.0 - 2.0 * cip
225	cimf = 1.0 - 2.0 * cim
226	ip = a0(k,i-1) * f2 * ( 1.0 - cipf ) &
227	+ a1(k,i-1) * f8 * ( 1.0 - cipf*cipf ) &
228	+ a2(k,i-1) * f24 * ( 1.0 - cipfcipfcipf )
229	im = a0(k,i) * f2 * ( 1.0 - cimf ) &
230	- a1(k,i) * f8 * ( 1.0 - cimf*cimf ) &
231	+ a2(k,i) * f24 * ( 1.0 - cimfcimfcimf )
232	ip = MAX( ip, 0.0 )
233	im = MAX( im, 0.0 )
234	ipmb(k,i) = ip * MIN( 1.0, sk_p(k,j,i-1) / (ip+im+1E-15) )
235	immb(k,i) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
236	ENDDO
237	ENDDO
238
239	!
240	!-- Compute monitor function m1
241	DO i = nxl-2, nxr+2
242	DO k = nzb+1, nzt-1
243	m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
244	m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
245	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
246	m1(k,i) = m1z / m1n
247	IF ( m1(k,i) /= 2.0 .AND. m1n < fmax(2) ) m1(k,i) = 0.0
248	ELSEIF ( m1n < m1z ) THEN
249	m1(k,i) = -1.0
250	ELSE
251	m1(k,i) = 0.0
252	ENDIF
253	ENDDO
254	ENDDO
255
256	!
257	!-- Compute switch sw
258	sw = 0.0
259	DO i = nxl-1, nxr+1
260	DO k = nzb+1, nzt-1
261	m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / &
262	MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 )
263	IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0
264
265	m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) / &
266	MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35 )
267	IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0
268
269	t1 = 0.35
270	t2 = 0.35
271	IF ( m1(k,i) == -1.0 ) t2 = 0.12
272
273	!-- *VOCL STMT,IF(10)
274	IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0 &
275	.OR. m2 > t2 .OR. m3 > T2 .OR. &
276	( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND. &
277	m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 ) &
278	) sw(k,i) = 1.0
279	ENDDO
280	ENDDO
281
282	!
283	!-- Fluxes using the exponential scheme
284	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
285	DO i = nxl, nxr
286	DO k = nzb+1, nzt-1
287
288	!-- *VOCL STMT,IF(10)
289	IF ( sw(k,i) == 1.0 ) THEN
290	snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
291	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
292	sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
293	sterm = MIN( sterm, 0.9999 )
294	sterm = MAX( sterm, 0.0001 )
295
296	ix = INT( sterm * 1000 ) + 1
297
298	cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
299
300	ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( &
301	aex(ix) * cip + bex(ix) / dex(ix) * ( &
302	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
303	) &
304	)
305	IF ( sterm == 0.0001 ) ippe(k,i) = sk_p(k,j,i) * cip
306	IF ( sterm == 0.9999 ) ippe(k,i) = sk_p(k,j,i) * cip
307
308	snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
309	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
310	sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
311	sterm = MIN( sterm, 0.9999 )
312	sterm = MAX( sterm, 0.0001 )
313
314	ix = INT( sterm * 1000 ) + 1
315
316	cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
317
318	imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( &
319	aex(ix) * cim + bex(ix) / dex(ix) * ( &
320	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
321	) &
322	)
323	IF ( sterm == 0.0001 ) imme(k,i) = sk_p(k,j,i) * cim
324	IF ( sterm == 0.9999 ) imme(k,i) = sk_p(k,j,i) * cim
325	ENDIF
326
327	!-- *VOCL STMT,IF(10)
328	IF ( sw(k,i+1) == 1.0 ) THEN
329	snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
330	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
331	sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
332	sterm = MIN( sterm, 0.9999 )
333	sterm = MAX( sterm, 0.0001 )
334
335	ix = INT( sterm * 1000 ) + 1
336
337	cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
338
339	impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( &
340	aex(ix) * cim + bex(ix) / dex(ix) * ( &
341	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
342	) &
343	)
344	IF ( sterm == 0.0001 ) impe(k,i) = sk_p(k,j,i+1) * cim
345	IF ( sterm == 0.9999 ) impe(k,i) = sk_p(k,j,i+1) * cim
346	ENDIF
347
348	!-- *VOCL STMT,IF(10)
349	IF ( sw(k,i-1) == 1.0 ) THEN
350	snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
351	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
352	sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
353	sterm = MIN( sterm, 0.9999 )
354	sterm = MAX( sterm, 0.0001 )
355
356	ix = INT( sterm * 1000 ) + 1
357
358	cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
359
360	ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( &
361	aex(ix) * cip + bex(ix) / dex(ix) * ( &
362	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
363	) &
364	)
365	IF ( sterm == 0.0001 ) ipme(k,i) = sk_p(k,j,i-1) * cip
366	IF ( sterm == 0.9999 ) ipme(k,i) = sk_p(k,j,i-1) * cip
367	ENDIF
368
369	ENDDO
370	ENDDO
371	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
372
373	!
374	!-- Prognostic equation
375	DO i = nxl, nxr
376	DO k = nzb+1, nzt-1
377	fplus = ( 1.0 - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) &
378	- ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
379	fminus = ( 1.0 - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) &
380	- ( 1.0 - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i)
381	tendenz = fplus - fminus
382	!
383	!-- Removed in order to optimize speed
384	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
385	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
386	!
387	!-- Density correction because of possible remaining divergences
388	d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
389	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
390	( 1.0 + d_new )
391	d(k,j,i) = d_new
392	ENDDO
393	ENDDO
394
395	ENDDO ! End of the advection in x-direction
396
397	!
398	!-- Deallocate temporary arrays
399	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
400	ippb, ippe, m1, sw )
401
402
403	!
404	!-- Enlarge boundary of local array cyclically in y-direction
405	#if defined( __parallel )
406	ngp = ( nzt - nzb + 5 ) * ( nyn - nys + 7 )
407	CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+5), ngp, MPI_REAL, &
408	type_xz_2, ierr )
409	CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
410	!
411	!-- Send front boundary, receive rear boundary
412	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
413	CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, &
414	sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, &
415	comm2d, status, ierr )
416	!
417	!-- Send rear boundary, receive front boundary
418	CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, &
419	sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, &
420	comm2d, status, ierr )
421	CALL MPI_TYPE_FREE( type_xz_2, ierr )
422	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
423	#else
424	DO i = nxl, nxr
425	DO k = nzb+1, nzt-1
426	sk_p(k,nys-1,i) = sk_p(k,nyn,i)
427	sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
428	sk_p(k,nys-3,i) = sk_p(k,nyn-2,i)
429	sk_p(k,nyn+1,i) = sk_p(k,nys,i)
430	sk_p(k,nyn+2,i) = sk_p(k,nys+1,i)
431	sk_p(k,nyn+3,i) = sk_p(k,nys+2,i)
432	ENDDO
433	ENDDO
434	#endif
435
436	!
437	!-- Determine the maxima of the first and second derivative in y-direction
438	fmax_l = 0.0
439	DO i = nxl, nxr
440	DO j = nys, nyn
441	DO k = nzb+1, nzt-1
442	zaehler = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
443	nenner = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
444	fmax_l(1) = MAX( fmax_l(1) , zaehler )
445	fmax_l(2) = MAX( fmax_l(2) , nenner )
446	ENDDO
447	ENDDO
448	ENDDO
449	#if defined( __parallel )
450	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
451	#else
452	fmax = fmax_l
453	#endif
454
455	fmax = 0.04 * fmax
456
457	!
458	!-- Allocate temporary arrays
459	ALLOCATE( a0(nzb+1:nzt-1,nys-1:nyn+1), a1(nzb+1:nzt-1,nys-1:nyn+1), &
460	a2(nzb+1:nzt-1,nys-1:nyn+1), a12(nzb+1:nzt-1,nys-1:nyn+1), &
461	a22(nzb+1:nzt-1,nys-1:nyn+1), immb(nzb+1:nzt-1,nys-1:nyn+1), &
462	imme(nzb+1:nzt-1,nys-1:nyn+1), impb(nzb+1:nzt-1,nys-1:nyn+1), &
463	impe(nzb+1:nzt-1,nys-1:nyn+1), ipmb(nzb+1:nzt-1,nys-1:nyn+1), &
464	ipme(nzb+1:nzt-1,nys-1:nyn+1), ippb(nzb+1:nzt-1,nys-1:nyn+1), &
465	ippe(nzb+1:nzt-1,nys-1:nyn+1), m1(nzb+1:nzt-1,nys-2:nyn+2), &
466	sw(nzb+1:nzt-1,nys-1:nyn+1) &
467	)
468	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
469
470	!
471	!-- Outer loop of all i
472	DO i = nxl, nxr
473
474	!
475	!-- Compute polynomial coefficients
476	DO j = nys-1, nyn+1
477	DO k = nzb+1, nzt-1
478	a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
479	a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) &
480	+ sk_p(k,j-1,i) )
481	a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i) &
482	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i) &
483	+ 9.0 * sk_p(k,j-2,i) ) * f1920
484	a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i) &
485	- 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i) &
486	) * f48
487	a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i) &
488	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i) &
489	- 3.0 * sk_p(k,j-2,i) ) * f48
490	ENDDO
491	ENDDO
492
493	!
494	!-- Fluxes using the Bott scheme
495	!-- *VOCL LOOP,UNROLL(2)
496	DO j = nys, nyn
497	DO k = nzb+1, nzt-1
498	cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
499	cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
500	cipf = 1.0 - 2.0 * cip
501	cimf = 1.0 - 2.0 * cim
502	ip = a0(k,j) * f2 * ( 1.0 - cipf ) &
503	+ a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &
504	+ a2(k,j) * f24 * ( 1.0 - cipfcipfcipf )
505	im = a0(k,j+1) * f2 * ( 1.0 - cimf ) &
506	- a1(k,j+1) * f8 * ( 1.0 - cimf*cimf ) &
507	+ a2(k,j+1) * f24 * ( 1.0 - cimfcimfcimf )
508	ip = MAX( ip, 0.0 )
509	im = MAX( im, 0.0 )
510	ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
511	impb(k,j) = im * MIN( 1.0, sk_p(k,j+1,i) / (ip+im+1E-15) )
512
513	cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
514	cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
515	cipf = 1.0 - 2.0 * cip
516	cimf = 1.0 - 2.0 * cim
517	ip = a0(k,j-1) * f2 * ( 1.0 - cipf ) &
518	+ a1(k,j-1) * f8 * ( 1.0 - cipf*cipf ) &
519	+ a2(k,j-1) * f24 * ( 1.0 - cipfcipfcipf )
520	im = a0(k,j) * f2 * ( 1.0 - cimf ) &
521	- a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &
522	+ a2(k,j) * f24 * ( 1.0 - cimfcimfcimf )
523	ip = MAX( ip, 0.0 )
524	im = MAX( im, 0.0 )
525	ipmb(k,j) = ip * MIN( 1.0, sk_p(k,j-1,i) / (ip+im+1E-15) )
526	immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
527	ENDDO
528	ENDDO
529
530	!
531	!-- Compute monitor function m1
532	DO j = nys-2, nyn+2
533	DO k = nzb+1, nzt-1
534	m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
535	m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
536	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
537	m1(k,j) = m1z / m1n
538	IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0
539	ELSEIF ( m1n < m1z ) THEN
540	m1(k,j) = -1.0
541	ELSE
542	m1(k,j) = 0.0
543	ENDIF
544	ENDDO
545	ENDDO
546
547	!
548	!-- Compute switch sw
549	sw = 0.0
550	DO j = nys-1, nyn+1
551	DO k = nzb+1, nzt-1
552	m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
553	MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
554	IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0
555
556	m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
557	MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
558	IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0
559
560	t1 = 0.35
561	t2 = 0.35
562	IF ( m1(k,j) == -1.0 ) t2 = 0.12
563
564	!-- *VOCL STMT,IF(10)
565	IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0 &
566	.OR. m2 > t2 .OR. m3 > T2 .OR. &
567	( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND. &
568	m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 ) &
569	) sw(k,j) = 1.0
570	ENDDO
571	ENDDO
572
573	!
574	!-- Fluxes using exponential scheme
575	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
576	DO j = nys, nyn
577	DO k = nzb+1, nzt-1
578
579	!-- *VOCL STMT,IF(10)
580	IF ( sw(k,j) == 1.0 ) THEN
581	snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
582	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
583	sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
584	sterm = MIN( sterm, 0.9999 )
585	sterm = MAX( sterm, 0.0001 )
586
587	ix = INT( sterm * 1000 ) + 1
588
589	cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
590
591	ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( &
592	aex(ix) * cip + bex(ix) / dex(ix) * ( &
593	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
594	) &
595	)
596	IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip
597	IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip
598
599	snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
600	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
601	sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
602	sterm = MIN( sterm, 0.9999 )
603	sterm = MAX( sterm, 0.0001 )
604
605	ix = INT( sterm * 1000 ) + 1
606
607	cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
608
609	imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( &
610	aex(ix) * cim + bex(ix) / dex(ix) * ( &
611	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
612	) &
613	)
614	IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim
615	IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim
616	ENDIF
617
618	!-- *VOCL STMT,IF(10)
619	IF ( sw(k,j+1) == 1.0 ) THEN
620	snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
621	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
622	sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
623	sterm = MIN( sterm, 0.9999 )
624	sterm = MAX( sterm, 0.0001 )
625
626	ix = INT( sterm * 1000 ) + 1
627
628	cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
629
630	impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( &
631	aex(ix) * cim + bex(ix) / dex(ix) * ( &
632	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
633	) &
634	)
635	IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k,j+1,i) * cim
636	IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k,j+1,i) * cim
637	ENDIF
638
639	!-- *VOCL STMT,IF(10)
640	IF ( sw(k,j-1) == 1.0 ) THEN
641	snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
642	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
643	sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
644	sterm = MIN( sterm, 0.9999 )
645	sterm = MAX( sterm, 0.0001 )
646
647	ix = INT( sterm * 1000 ) + 1
648
649	cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
650
651	ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( &
652	aex(ix) * cip + bex(ix) / dex(ix) * ( &
653	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
654	) &
655	)
656	IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k,j-1,i) * cip
657	IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k,j-1,i) * cip
658	ENDIF
659
660	ENDDO
661	ENDDO
662	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
663
664	!
665	!-- Prognostic equation
666	DO j = nys, nyn
667	DO k = nzb+1, nzt-1
668	fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &
669	- ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
670	fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
671	- ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)
672	tendenz = fplus - fminus
673	!
674	!-- Removed in order to optimise speed
675	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
676	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
677	!
678	!-- Density correction because of possible remaining divergences
679	d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
680	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
681	( 1.0 + d_new )
682	d(k,j,i) = d_new
683	ENDDO
684	ENDDO
685
686	ENDDO ! End of the advection in y-direction
687	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
688	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
689
690	!
691	!-- Deallocate temporary arrays
692	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
693	ippb, ippe, m1, sw )
694
695
696	!
697	!-- Initialise for the computation of heat fluxes (see below; required in
698	!-- UP flow_statistics)
699	IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0
700
701	!
702	!-- Add top and bottom boundaries according to the relevant boundary conditions
703	IF ( sk_char == 'pt' ) THEN
704
705	!
706	!-- Temperature boundary condition at the bottom boundary
707	IF ( ibc_pt_b == 0 ) THEN
708	!
709	!-- Dirichlet (fixed surface temperature)
710	DO i = nxl, nxr
711	DO j = nys, nyn
712	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
713	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
714	ENDDO
715	ENDDO
716
717	ELSE
718	!
719	!-- Neumann (i.e. here zero gradient)
720	DO i = nxl, nxr
721	DO j = nys, nyn
722	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
723	sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
724	sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
725	ENDDO
726	ENDDO
727
728	ENDIF
729
730	!
731	!-- Temperature boundary condition at the top boundary
732	IF ( ibc_pt_t == 0 ) THEN
733	!
734	!-- Dirichlet
735	DO i = nxl, nxr
736	DO j = nys, nyn
737	sk_p(nzt+1,j,i) = sk_p(nzt,j,i)
738	sk_p(nzt+2,j,i) = sk_p(nzt,j,i)
739	ENDDO
740	ENDDO
741
742	ELSE
743	!
744	!-- Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
745	DO i = nxl, nxr
746	DO j = nys, nyn
747	sk_p(nzt,j,i) = sk_p(nzt-1,j,i) + bc_pt_t_val * dzu(nzt)
748	sk_p(nzt+1,j,i) = sk_p(nzt,j,i) + bc_pt_t_val * dzu(nzt+1)
749	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
750	ENDDO
751	ENDDO
752
753	ENDIF
754
755	ELSEIF ( sk_char == 'q' ) THEN
756
757	!
758	!-- Temperature boundary condition at the bottom boundary
759	IF ( ibc_q_b == 0 ) THEN
760	!
761	!-- Dirichlet (fixed surface temperature)
762	DO i = nxl, nxr
763	DO j = nys, nyn
764	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
765	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
766	ENDDO
767	ENDDO
768
769	ELSE
770	!
771	!-- Neumann (i.e. here zero gradient)
772	DO i = nxl, nxr
773	DO j = nys, nyn
774	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
775	sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
776	sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
777	ENDDO
778	ENDDO
779
780	ENDIF
781
782	!
783	!-- Temperature boundary condition at the top boundary
784	IF ( ibc_q_t == 0 ) THEN
785	!
786	!-- Dirichlet
787	DO i = nxl, nxr
788	DO j = nys, nyn
789	sk_p(nzt+1,j,i) = sk_p(nzt,j,i)
790	sk_p(nzt+2,j,i) = sk_p(nzt,j,i)
791	ENDDO
792	ENDDO
793
794	ELSE
795	!
796	!-- Neumann: dzu(nzt+2) is not defined, thus instead dzu(nzt+1) is used
797	DO i = nxl, nxr
798	DO j = nys, nyn
799	sk_p(nzt,j,i) = sk_p(nzt-1,j,i) + bc_q_t_val * dzu(nzt)
800	sk_p(nzt+1,j,i) = sk_p(nzt,j,i) + bc_q_t_val * dzu(nzt+1)
801	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
802	ENDDO
803	ENDDO
804
805	ENDIF
806
807	ELSEIF ( sk_char == 'e' ) THEN
808
809	!
810	!-- TKE boundary condition at bottom and top boundary (generally Neumann)
811	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
812	sk_p(nzb-1,j,i) = sk_p(nzb+1,j,i)
813	sk_p(nzb-2,j,i) = sk_p(nzb+1,j,i)
814	sk_p(nzt,j,i) = sk_p(nzt-1,j,i)
815	sk_p(nzt+1,j,i) = sk_p(nzt-1,j,i)
816	sk_p(nzt+2,j,i) = sk_p(nzt-1,j,i)
817
818	ELSE
819
820	IF ( myid == 0 ) PRINT*,'+++ advec_s_bc: no vertical boundary condi', &
821	'tion for variable "', sk_char, '"'
822	CALL local_stop
823
824	ENDIF
825
826	!
827	!-- Determine the maxima of the first and second derivative in z-direction
828	fmax_l = 0.0
829	DO i = nxl, nxr
830	DO j = nys, nyn
831	DO k = nzb, nzt
832	zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
833	nenner = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
834	fmax_l(1) = MAX( fmax_l(1) , zaehler )
835	fmax_l(2) = MAX( fmax_l(2) , nenner )
836	ENDDO
837	ENDDO
838	ENDDO
839	#if defined( __parallel )
840	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
841	#else
842	fmax = fmax_l
843	#endif
844
845	fmax = 0.04 * fmax
846
847	!
848	!-- Allocate temporary arrays
849	ALLOCATE( a0(nzb:nzt,nys:nyn), a1(nzb:nzt,nys:nyn), &
850	a2(nzb:nzt,nys:nyn), a12(nzb:nzt,nys:nyn), &
851	a22(nzb:nzt,nys:nyn), immb(nzb+1:nzt-1,nys:nyn), &
852	imme(nzb+1:nzt-1,nys:nyn), impb(nzb+1:nzt-1,nys:nyn), &
853	impe(nzb+1:nzt-1,nys:nyn), ipmb(nzb+1:nzt-1,nys:nyn), &
854	ipme(nzb+1:nzt-1,nys:nyn), ippb(nzb+1:nzt-1,nys:nyn), &
855	ippe(nzb+1:nzt-1,nys:nyn), m1(nzb-1:nzt+1,nys:nyn), &
856	sw(nzb:nzt,nys:nyn) &
857	)
858	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
859
860	!
861	!-- Outer loop of all i
862	DO i = nxl, nxr
863
864	!
865	!-- Compute polynomial coefficients
866	DO j = nys, nyn
867	DO k = nzb, nzt
868	a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
869	a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
870	+ sk_p(k-1,j,i) )
871	a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) &
872	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) &
873	+ 9.0 * sk_p(k-2,j,i) ) * f1920
874	a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) &
875	- 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) &
876	) * f48
877	a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) &
878	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) &
879	- 3.0 * sk_p(k-2,j,i) ) * f48
880	ENDDO
881	ENDDO
882
883	!
884	!-- Fluxes using the Bott scheme
885	!-- *VOCL LOOP,UNROLL(2)
886	DO j = nys, nyn
887	DO k = nzb+1, nzt-1
888	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
889	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
890	cipf = 1.0 - 2.0 * cip
891	cimf = 1.0 - 2.0 * cim
892	ip = a0(k,j) * f2 * ( 1.0 - cipf ) &
893	+ a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &
894	+ a2(k,j) * f24 * ( 1.0 - cipfcipfcipf )
895	im = a0(k+1,j) * f2 * ( 1.0 - cimf ) &
896	- a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) &
897	+ a2(k+1,j) * f24 * ( 1.0 - cimfcimfcimf )
898	ip = MAX( ip, 0.0 )
899	im = MAX( im, 0.0 )
900	ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
901	impb(k,j) = im * MIN( 1.0, sk_p(k+1,j,i) / (ip+im+1E-15) )
902
903	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
904	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
905	cipf = 1.0 - 2.0 * cip
906	cimf = 1.0 - 2.0 * cim
907	ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) &
908	+ a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) &
909	+ a2(k-1,j) * f24 * ( 1.0 - cipfcipfcipf )
910	im = a0(k,j) * f2 * ( 1.0 - cimf ) &
911	- a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &
912	+ a2(k,j) * f24 * ( 1.0 - cimfcimfcimf )
913	ip = MAX( ip, 0.0 )
914	im = MAX( im, 0.0 )
915	ipmb(k,j) = ip * MIN( 1.0, sk_p(k-1,j,i) / (ip+im+1E-15) )
916	immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
917	ENDDO
918	ENDDO
919
920	!
921	!-- Compute monitor function m1
922	DO j = nys, nyn
923	DO k = nzb-1, nzt+1
924	m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
925	m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
926	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
927	m1(k,j) = m1z / m1n
928	IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0
929	ELSEIF ( m1n < m1z ) THEN
930	m1(k,j) = -1.0
931	ELSE
932	m1(k,j) = 0.0
933	ENDIF
934	ENDDO
935	ENDDO
936
937	!
938	!-- Compute switch sw
939	sw = 0.0
940	DO j = nys, nyn
941	DO k = nzb, nzt
942	m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
943	MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
944	IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0
945
946	m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
947	MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
948	IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0
949
950	t1 = 0.35
951	t2 = 0.35
952	IF ( m1(k,j) == -1.0 ) t2 = 0.12
953
954	!-- *VOCL STMT,IF(10)
955	IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 &
956	.OR. m2 > t2 .OR. m3 > T2 .OR. &
957	( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. &
958	m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) &
959	) sw(k,j) = 1.0
960	ENDDO
961	ENDDO
962
963	!
964	!-- Fluxes using exponential scheme
965	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
966	DO j = nys, nyn
967	DO k = nzb+1, nzt-1
968
969	!-- *VOCL STMT,IF(10)
970	IF ( sw(k,j) == 1.0 ) THEN
971	snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
972	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
973	sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
974	sterm = MIN( sterm, 0.9999 )
975	sterm = MAX( sterm, 0.0001 )
976
977	ix = INT( sterm * 1000 ) + 1
978
979	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
980
981	ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( &
982	aex(ix) * cip + bex(ix) / dex(ix) * ( &
983	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
984	) &
985	)
986	IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip
987	IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip
988
989	snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
990	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
991	sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
992	sterm = MIN( sterm, 0.9999 )
993	sterm = MAX( sterm, 0.0001 )
994
995	ix = INT( sterm * 1000 ) + 1
996
997	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
998
999	imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( &
1000	aex(ix) * cim + bex(ix) / dex(ix) * ( &
1001	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
1002	) &
1003	)
1004	IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim
1005	IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim
1006	ENDIF
1007
1008	!-- *VOCL STMT,IF(10)
1009	IF ( sw(k+1,j) == 1.0 ) THEN
1010	snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
1011	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
1012	sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
1013	sterm = MIN( sterm, 0.9999 )
1014	sterm = MAX( sterm, 0.0001 )
1015
1016	ix = INT( sterm * 1000 ) + 1
1017
1018	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
1019
1020	impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( &
1021	aex(ix) * cim + bex(ix) / dex(ix) * ( &
1022	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
1023	) &
1024	)
1025	IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k+1,j,i) * cim
1026	IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k+1,j,i) * cim
1027	ENDIF
1028
1029	!-- *VOCL STMT,IF(10)
1030	IF ( sw(k-1,j) == 1.0 ) THEN
1031	snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
1032	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
1033	sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
1034	sterm = MIN( sterm, 0.9999 )
1035	sterm = MAX( sterm, 0.0001 )
1036
1037	ix = INT( sterm * 1000 ) + 1
1038
1039	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
1040
1041	ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( &
1042	aex(ix) * cip + bex(ix) / dex(ix) * ( &
1043	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
1044	) &
1045	)
1046	IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1047	IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1048	ENDIF
1049
1050	ENDDO
1051	ENDDO
1052	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
1053
1054	!
1055	!-- Prognostic equation
1056	DO j = nys, nyn
1057	DO k = nzb+1, nzt-1
1058	fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &
1059	- ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
1060	fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
1061	- ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)
1062	tendenz = fplus - fminus
1063	!
1064	!-- Removed in order to optimise speed
1065	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
1066	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
1067	!
1068	!-- Density correction because of possible remaining divergences
1069	d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
1070	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
1071	( 1.0 + d_new )
1072	!
1073	!-- Store heat flux for subsequent statistics output.
1074	!-- array m1 is here used as temporary storage
1075	m1(k,j) = fplus / dt_3d * dzw(k)
1076	ENDDO
1077	ENDDO
1078
1079	!
1080	!-- Sum up heat flux in order to order to obtain horizontal averages
1081	IF ( sk_char == 'pt' ) THEN
1082	DO sr = 0, statistic_regions
1083	DO j = nys, nyn
1084	DO k = nzb+1, nzt-1
1085	sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
1086	m1(k,j) * rmask(j,i,sr)
1087	ENDDO
1088	ENDDO
1089	ENDDO
1090	ENDIF
1091
1092	ENDDO ! End of the advection in z-direction
1093	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1094	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
1095
1096	!
1097	!-- Deallocate temporary arrays
1098	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
1099	ippb, ippe, m1, sw )
1100
1101	!
1102	!-- Store results as tendency and deallocate local array
1103	DO i = nxl, nxr
1104	DO j = nys, nyn
1105	DO k = nzb+1, nzt-1
1106	tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
1107	ENDDO
1108	ENDDO
1109	ENDDO
1110
1111	DEALLOCATE( sk_p )
1112
1113	END SUBROUTINE advec_s_bc

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |