Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |