Home

Context Navigation

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

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