Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |