Home

Context Navigation

advec_s_bc.f90 @ 4879

Last change on this file since 4879 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 55.6 KB

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

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

Context Navigation

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

Download in other formats: