Home

Context Navigation

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

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