Home

Context Navigation

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

Last change on this file since 1685 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.2 KB

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