Home

Context Navigation

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

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