Home

Context Navigation

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

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