Home

Context Navigation

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

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