Home

Context Navigation

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

Last change on this file since 70 was 63, checked in by raasch, 18 years ago
preliminary changes concerning update of BC-scheme
Property svn:keywords set to `Id`
File size: 42.1 KB

Line
1	SUBROUTINE advec_s_bc( sk, sk_char )
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! Calculation extended for gridpoint nzt
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: advec_s_bc.f90 63 2007-03-13 03:52:49Z 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 control_parameters
36	USE cpulog
37	USE grid_variables
38	USE indices
39	USE interfaces
40	USE pegrid
41	USE statistics
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( __nec )
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+3,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, 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 + 6 ) * ( 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
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,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), &
163	a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), &
164	a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), &
165	imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &
166	impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &
167	ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &
168	ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), &
169	sw(nzb+1:nzt,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
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
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
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
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
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
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 + 6 ) * ( nyn - nys + 7 )
407	CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), 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
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
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,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), &
460	a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), &
461	a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), &
462	imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), &
463	impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), &
464	ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), &
465	ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), &
466	sw(nzb+1:nzt,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
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
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
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
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
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
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-1,j,i) = sk_p(nzb,j,i)
723	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
724	ENDDO
725	ENDDO
726
727	ENDIF
728
729	!
730	!-- Temperature boundary condition at the top boundary
731	IF ( ibc_pt_t == 0 .OR. ibc_pt_t == 1 ) THEN
732	!
733	!-- Dirichlet or Neumann (zero gradient)
734	DO i = nxl, nxr
735	DO j = nys, nyn
736	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
737	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
738	ENDDO
739	ENDDO
740
741	ELSEIF ( ibc_pt_t == 2 ) THEN
742	!
743	!-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
744	DO i = nxl, nxr
745	DO j = nys, nyn
746	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
747	sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
748	ENDDO
749	ENDDO
750
751	ENDIF
752
753	ELSEIF ( sk_char == 'q' ) THEN
754
755	!
756	!-- Specific humidity boundary condition at the bottom boundary
757	IF ( ibc_q_b == 0 ) THEN
758	!
759	!-- Dirichlet (fixed surface humidity)
760	DO i = nxl, nxr
761	DO j = nys, nyn
762	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
763	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
764	ENDDO
765	ENDDO
766
767	ELSE
768	!
769	!-- Neumann (i.e. here zero gradient)
770	DO i = nxl, nxr
771	DO j = nys, nyn
772	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
773	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
774	ENDDO
775	ENDDO
776
777	ENDIF
778
779	!
780	!-- Specific humidity boundary condition at the top boundary
781	IF ( ibc_q_t == 0 ) THEN
782	!
783	!-- Dirichlet
784	DO i = nxl, nxr
785	DO j = nys, nyn
786	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
787	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
788	ENDDO
789	ENDDO
790
791	ELSE
792	!
793	!-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
794	DO i = nxl, nxr
795	DO j = nys, nyn
796	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
797	sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
798	ENDDO
799	ENDDO
800
801	ENDIF
802
803	ELSEIF ( sk_char == 'e' ) THEN
804
805	!
806	!-- TKE boundary condition at bottom and top boundary (generally Neumann)
807	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
808	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
809	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
810	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
811
812	ELSE
813
814	IF ( myid == 0 ) PRINT*,'+++ advec_s_bc: no vertical boundary condi', &
815	'tion for variable "', sk_char, '"'
816	CALL local_stop
817
818	ENDIF
819
820	!
821	!-- Determine the maxima of the first and second derivative in z-direction
822	fmax_l = 0.0
823	DO i = nxl, nxr
824	DO j = nys, nyn
825	DO k = nzb, nzt+1
826	zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
827	nenner = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
828	fmax_l(1) = MAX( fmax_l(1) , zaehler )
829	fmax_l(2) = MAX( fmax_l(2) , nenner )
830	ENDDO
831	ENDDO
832	ENDDO
833	#if defined( __parallel )
834	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
835	#else
836	fmax = fmax_l
837	#endif
838
839	fmax = 0.04 * fmax
840
841	!
842	!-- Allocate temporary arrays
843	ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), &
844	a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), &
845	a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), &
846	imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), &
847	impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), &
848	ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), &
849	ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), &
850	sw(nzb:nzt+1,nys:nyn) &
851	)
852	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
853
854	!
855	!-- Outer loop of all i
856	DO i = nxl, nxr
857
858	!
859	!-- Compute polynomial coefficients
860	DO j = nys, nyn
861	DO k = nzb, nzt+1
862	a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
863	a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
864	+ sk_p(k-1,j,i) )
865	a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) &
866	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) &
867	+ 9.0 * sk_p(k-2,j,i) ) * f1920
868	a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) &
869	- 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) &
870	) * f48
871	a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) &
872	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) &
873	- 3.0 * sk_p(k-2,j,i) ) * f48
874	ENDDO
875	ENDDO
876
877	!
878	!-- Fluxes using the Bott scheme
879	!-- *VOCL LOOP,UNROLL(2)
880	DO j = nys, nyn
881	DO k = nzb+1, nzt
882	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
883	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
884	cipf = 1.0 - 2.0 * cip
885	cimf = 1.0 - 2.0 * cim
886	ip = a0(k,j) * f2 * ( 1.0 - cipf ) &
887	+ a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &
888	+ a2(k,j) * f24 * ( 1.0 - cipfcipfcipf )
889	im = a0(k+1,j) * f2 * ( 1.0 - cimf ) &
890	- a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) &
891	+ a2(k+1,j) * f24 * ( 1.0 - cimfcimfcimf )
892	ip = MAX( ip, 0.0 )
893	im = MAX( im, 0.0 )
894	ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
895	impb(k,j) = im * MIN( 1.0, sk_p(k+1,j,i) / (ip+im+1E-15) )
896
897	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
898	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
899	cipf = 1.0 - 2.0 * cip
900	cimf = 1.0 - 2.0 * cim
901	ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) &
902	+ a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) &
903	+ a2(k-1,j) * f24 * ( 1.0 - cipfcipfcipf )
904	im = a0(k,j) * f2 * ( 1.0 - cimf ) &
905	- a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &
906	+ a2(k,j) * f24 * ( 1.0 - cimfcimfcimf )
907	ip = MAX( ip, 0.0 )
908	im = MAX( im, 0.0 )
909	ipmb(k,j) = ip * MIN( 1.0, sk_p(k-1,j,i) / (ip+im+1E-15) )
910	immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
911	ENDDO
912	ENDDO
913
914	!
915	!-- Compute monitor function m1
916	DO j = nys, nyn
917	DO k = nzb-1, nzt+2
918	m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
919	m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
920	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
921	m1(k,j) = m1z / m1n
922	IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0
923	ELSEIF ( m1n < m1z ) THEN
924	m1(k,j) = -1.0
925	ELSE
926	m1(k,j) = 0.0
927	ENDIF
928	ENDDO
929	ENDDO
930
931	!
932	!-- Compute switch sw
933	sw = 0.0
934	DO j = nys, nyn
935	DO k = nzb, nzt+1
936	m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
937	MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
938	IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0
939
940	m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
941	MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
942	IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0
943
944	t1 = 0.35
945	t2 = 0.35
946	IF ( m1(k,j) == -1.0 ) t2 = 0.12
947
948	!-- *VOCL STMT,IF(10)
949	IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 &
950	.OR. m2 > t2 .OR. m3 > T2 .OR. &
951	( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. &
952	m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) &
953	) sw(k,j) = 1.0
954	ENDDO
955	ENDDO
956
957	!
958	!-- Fluxes using exponential scheme
959	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
960	DO j = nys, nyn
961	DO k = nzb+1, nzt
962
963	!-- *VOCL STMT,IF(10)
964	IF ( sw(k,j) == 1.0 ) THEN
965	snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
966	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
967	sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
968	sterm = MIN( sterm, 0.9999 )
969	sterm = MAX( sterm, 0.0001 )
970
971	ix = INT( sterm * 1000 ) + 1
972
973	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
974
975	ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( &
976	aex(ix) * cip + bex(ix) / dex(ix) * ( &
977	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
978	) &
979	)
980	IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip
981	IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip
982
983	snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
984	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
985	sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
986	sterm = MIN( sterm, 0.9999 )
987	sterm = MAX( sterm, 0.0001 )
988
989	ix = INT( sterm * 1000 ) + 1
990
991	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
992
993	imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( &
994	aex(ix) * cim + bex(ix) / dex(ix) * ( &
995	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
996	) &
997	)
998	IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim
999	IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim
1000	ENDIF
1001
1002	!-- *VOCL STMT,IF(10)
1003	IF ( sw(k+1,j) == 1.0 ) THEN
1004	snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
1005	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
1006	sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
1007	sterm = MIN( sterm, 0.9999 )
1008	sterm = MAX( sterm, 0.0001 )
1009
1010	ix = INT( sterm * 1000 ) + 1
1011
1012	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
1013
1014	impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( &
1015	aex(ix) * cim + bex(ix) / dex(ix) * ( &
1016	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
1017	) &
1018	)
1019	IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k+1,j,i) * cim
1020	IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k+1,j,i) * cim
1021	ENDIF
1022
1023	!-- *VOCL STMT,IF(10)
1024	IF ( sw(k-1,j) == 1.0 ) THEN
1025	snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
1026	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
1027	sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
1028	sterm = MIN( sterm, 0.9999 )
1029	sterm = MAX( sterm, 0.0001 )
1030
1031	ix = INT( sterm * 1000 ) + 1
1032
1033	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
1034
1035	ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( &
1036	aex(ix) * cip + bex(ix) / dex(ix) * ( &
1037	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
1038	) &
1039	)
1040	IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1041	IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1042	ENDIF
1043
1044	ENDDO
1045	ENDDO
1046	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
1047
1048	!
1049	!-- Prognostic equation
1050	DO j = nys, nyn
1051	DO k = nzb+1, nzt
1052	fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &
1053	- ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
1054	fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
1055	- ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)
1056	tendenz = fplus - fminus
1057	!
1058	!-- Removed in order to optimise speed
1059	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
1060	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
1061	!
1062	!-- Density correction because of possible remaining divergences
1063	d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
1064	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
1065	( 1.0 + d_new )
1066	!
1067	!-- Store heat flux for subsequent statistics output.
1068	!-- array m1 is here used as temporary storage
1069	m1(k,j) = fplus / dt_3d * dzw(k)
1070	ENDDO
1071	ENDDO
1072
1073	!
1074	!-- Sum up heat flux in order to order to obtain horizontal averages
1075	IF ( sk_char == 'pt' ) THEN
1076	DO sr = 0, statistic_regions
1077	DO j = nys, nyn
1078	DO k = nzb+1, nzt
1079	sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
1080	m1(k,j) * rmask(j,i,sr)
1081	ENDDO
1082	ENDDO
1083	ENDDO
1084	ENDIF
1085
1086	ENDDO ! End of the advection in z-direction
1087	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1088	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
1089
1090	!
1091	!-- Deallocate temporary arrays
1092	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
1093	ippb, ippe, m1, sw )
1094
1095	!
1096	!-- Store results as tendency and deallocate local array
1097	DO i = nxl, nxr
1098	DO j = nys, nyn
1099	DO k = nzb+1, nzt
1100	tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
1101	ENDDO
1102	ENDDO
1103	ENDDO
1104
1105	DEALLOCATE( sk_p )
1106
1107	END SUBROUTINE advec_s_bc

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |