Home

Context Navigation

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

Last change on this file since 1873 was 1873, checked in by maronga, 8 years ago
revised renaming of modules
Property svn:keywords set to `Id`
File size: 54.4 KB

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