Home

Context Navigation

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

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