Home

Context Navigation

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

Last change on this file since 3948 was 3761, checked in by raasch, 6 years ago
unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings
Property svn:keywords set to `Id`
File size: 57.6 KB

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