Home

Context Navigation

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

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