Home

Context Navigation

advec_s_bc_mod.f90 @ 1859

Last change on this file since 1859 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/SOURCE/advec_s_bc_mod.f90 @ 1859

Download in other formats: