Home

Context Navigation

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

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