Home

Context Navigation

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

Last change on this file since 1354 was 1354, checked in by heinze, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 49.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |