Home

Context Navigation

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

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