Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |