Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |