Home

Context Navigation

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

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