Home

Context Navigation

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

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