Home

Context Navigation

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

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