Home

Context Navigation

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

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