Home

Context Navigation

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

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