Home

Context Navigation

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

Last change on this file since 1067 was 1037, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 44.5 KB

Line
1	SUBROUTINE advec_s_bc( sk, sk_char )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later 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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: advec_s_bc.f90 1037 2012-10-22 14:10:22Z maronga $
27	!
28	! 1036 2012-10-22 13:43:42Z raasch
29	! code put under GPL (PALM 3.9)
30	!
31	! 1010 2012-09-20 07:59:54Z raasch
32	! cpp switch __nopointer added for pointer free version
33	!
34	! 622 2010-12-10 08:08:13Z raasch
35	! optional barriers included in order to speed up collective operations
36	!
37	! 247 2009-02-27 14:01:30Z heinze
38	! Output of messages replaced by message handling routine
39	!
40	! 216 2008-11-25 07:12:43Z raasch
41	! Neumann boundary condition at k=nzb is explicitly set for better reading,
42	! although this has been already done in boundary_conds
43	!
44	! 97 2007-06-21 08:23:15Z raasch
45	! Advection of salinity included
46	! Bugfix: Error in boundary condition for TKE removed
47	!
48	! 63 2007-03-13 03:52:49Z raasch
49	! Calculation extended for gridpoint nzt
50	!
51	! RCS Log replace by Id keyword, revision history cleaned up
52	!
53	! Revision 1.22 2006/02/23 09:42:08 raasch
54	! anz renamed ngp
55	!
56	! Revision 1.1 1997/08/29 08:53:46 raasch
57	! Initial revision
58	!
59	!
60	! Description:
61	! ------------
62	! Advection term for scalar quantities using the Bott-Chlond scheme.
63	! Computation in individual steps for each of the three dimensions.
64	! Limiting assumptions:
65	! So far the scheme has been assuming equidistant grid spacing. As this is not
66	! the case in the stretched portion of the z-direction, there dzw(k) is used as
67	! a substitute for a constant grid length. This certainly causes incorrect
68	! results; however, it is hoped that they are not too apparent for weakly
69	! stretched grids.
70	! NOTE: This is a provisional, non-optimised version!
71	!------------------------------------------------------------------------------!
72
73	USE advection
74	USE arrays_3d
75	USE control_parameters
76	USE cpulog
77	USE grid_variables
78	USE indices
79	USE interfaces
80	USE pegrid
81	USE statistics
82
83	IMPLICIT NONE
84
85	CHARACTER (LEN=*) :: sk_char
86
87	INTEGER :: i, ix, j, k, ngp, sr, type_xz_2
88
89	REAL :: cim, cimf, cip, cipf, d_new, ffmax, fminus, fplus, f2, f4, f8, &
90	f12, f24, f48, f1920, im, ip, m2, m3, nenner, snenn, sterm, &
91	tendenz, t1, t2, zaehler
92	REAL :: fmax(2), fmax_l(2)
93	#if defined( __nopointer )
94	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
95	#else
96	REAL, DIMENSION(:,:,:), POINTER :: sk
97	#endif
98
99	REAL, DIMENSION(:,:), ALLOCATABLE :: a0, a1, a12, a2, a22, immb, imme, &
100	impb, impe, ipmb, ipme, ippb, ippe
101	REAL, DIMENSION(:,:,:), ALLOCATABLE :: sk_p
102
103	#if defined( __nec )
104	REAL (kind=4) :: m1n, m1z !Wichtig: Division
105	REAL (kind=4), DIMENSION(:,:), ALLOCATABLE :: m1, sw
106	#else
107	REAL :: m1n, m1z
108	REAL, DIMENSION(:,:), ALLOCATABLE :: m1, sw
109	#endif
110
111
112	!
113	!-- Array sk_p requires 2 extra elements for each dimension
114	ALLOCATE( sk_p(nzb-2:nzt+3,nys-3:nyn+3,nxl-3:nxr+3) )
115	sk_p = 0.0
116
117	!
118	!-- Assign reciprocal values in order to avoid divisions later
119	f2 = 0.5
120	f4 = 0.25
121	f8 = 0.125
122	f12 = 0.8333333333333333E-01
123	f24 = 0.4166666666666666E-01
124	f48 = 0.2083333333333333E-01
125	f1920 = 0.5208333333333333E-03
126
127	!
128	!-- Advection in x-direction:
129
130	!
131	!-- Save the quantity to be advected in a local array
132	!-- add an enlarged boundary in x-direction
133	DO i = nxl-1, nxr+1
134	DO j = nys, nyn
135	DO k = nzb, nzt+1
136	sk_p(k,j,i) = sk(k,j,i)
137	ENDDO
138	ENDDO
139	ENDDO
140	#if defined( __parallel )
141	ngp = 2 * ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
142	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'start' )
143	!
144	!-- Send left boundary, receive right boundary
145	CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxl+1), ngp, MPI_REAL, pleft, 0, &
146	sk_p(nzb-2,nys-3,nxr+2), ngp, MPI_REAL, pright, 0, &
147	comm2d, status, ierr )
148	!
149	!-- Send right boundary, receive left boundary
150	CALL MPI_SENDRECV( sk_p(nzb-2,nys-3,nxr-2), ngp, MPI_REAL, pright, 1, &
151	sk_p(nzb-2,nys-3,nxl-3), ngp, MPI_REAL, pleft, 1, &
152	comm2d, status, ierr )
153	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
154	#else
155
156	!
157	!-- Cyclic boundary conditions
158	sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxr-2)
159	sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxr-1)
160	sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxl+1)
161	sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxl+2)
162	#endif
163
164	!
165	!-- In case of a sloping surface, the additional gridpoints in x-direction
166	!-- of the temperature field at the left and right boundary of the total
167	!-- domain must be adjusted by the temperature difference between this distance
168	IF ( sloping_surface .AND. sk_char == 'pt' ) THEN
169	IF ( nxl == 0 ) THEN
170	sk_p(:,nys:nyn,nxl-3) = sk_p(:,nys:nyn,nxl-3) - pt_slope_offset
171	sk_p(:,nys:nyn,nxl-2) = sk_p(:,nys:nyn,nxl-2) - pt_slope_offset
172	ENDIF
173	IF ( nxr == nx ) THEN
174	sk_p(:,nys:nyn,nxr+2) = sk_p(:,nys:nyn,nxr+2) + pt_slope_offset
175	sk_p(:,nys:nyn,nxr+3) = sk_p(:,nys:nyn,nxr+3) + pt_slope_offset
176	ENDIF
177	ENDIF
178
179	!
180	!-- Initialise control density
181	d = 0.0
182
183	!
184	!-- Determine maxima of the first and second derivative in x-direction
185	fmax_l = 0.0
186	DO i = nxl, nxr
187	DO j = nys, nyn
188	DO k = nzb+1, nzt
189	zaehler = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
190	nenner = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
191	fmax_l(1) = MAX( fmax_l(1) , zaehler )
192	fmax_l(2) = MAX( fmax_l(2) , nenner )
193	ENDDO
194	ENDDO
195	ENDDO
196	#if defined( __parallel )
197	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
198	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
199	#else
200	fmax = fmax_l
201	#endif
202
203	fmax = 0.04 * fmax
204
205	!
206	!-- Allocate temporary arrays
207	ALLOCATE( a0(nzb+1:nzt,nxl-1:nxr+1), a1(nzb+1:nzt,nxl-1:nxr+1), &
208	a2(nzb+1:nzt,nxl-1:nxr+1), a12(nzb+1:nzt,nxl-1:nxr+1), &
209	a22(nzb+1:nzt,nxl-1:nxr+1), immb(nzb+1:nzt,nxl-1:nxr+1), &
210	imme(nzb+1:nzt,nxl-1:nxr+1), impb(nzb+1:nzt,nxl-1:nxr+1), &
211	impe(nzb+1:nzt,nxl-1:nxr+1), ipmb(nzb+1:nzt,nxl-1:nxr+1), &
212	ipme(nzb+1:nzt,nxl-1:nxr+1), ippb(nzb+1:nzt,nxl-1:nxr+1), &
213	ippe(nzb+1:nzt,nxl-1:nxr+1), m1(nzb+1:nzt,nxl-2:nxr+2), &
214	sw(nzb+1:nzt,nxl-1:nxr+1) &
215	)
216	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
217
218	!
219	!-- Initialise point of time measuring of the exponential portion (this would
220	!-- not work if done locally within the loop)
221	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'start' )
222	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
223
224	!
225	!-- Outer loop of all j
226	DO j = nys, nyn
227
228	!
229	!-- Compute polynomial coefficients
230	DO i = nxl-1, nxr+1
231	DO k = nzb+1, nzt
232	a12(k,i) = 0.5 * ( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
233	a22(k,i) = 0.5 * ( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) &
234	+ sk_p(k,j,i-1) )
235	a0(k,i) = ( 9.0 * sk_p(k,j,i+2) - 116.0 * sk_p(k,j,i+1) &
236	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j,i-1) &
237	+ 9.0 * sk_p(k,j,i-2) ) * f1920
238	a1(k,i) = ( -5.0 * sk_p(k,j,i+2) + 34.0 * sk_p(k,j,i+1) &
239	- 34.0 * sk_p(k,j,i-1) + 5.0 * sk_p(k,j,i-2) &
240	) * f48
241	a2(k,i) = ( -3.0 * sk_p(k,j,i+2) + 36.0 * sk_p(k,j,i+1) &
242	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j,i-1) &
243	- 3.0 * sk_p(k,j,i-2) ) * f48
244	ENDDO
245	ENDDO
246
247	!
248	!-- Fluxes using the Bott scheme
249	!-- *VOCL LOOP,UNROLL(2)
250	DO i = nxl, nxr
251	DO k = nzb+1, nzt
252	cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
253	cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
254	cipf = 1.0 - 2.0 * cip
255	cimf = 1.0 - 2.0 * cim
256	ip = a0(k,i) * f2 * ( 1.0 - cipf ) &
257	+ a1(k,i) * f8 * ( 1.0 - cipf*cipf ) &
258	+ a2(k,i) * f24 * ( 1.0 - cipfcipfcipf )
259	im = a0(k,i+1) * f2 * ( 1.0 - cimf ) &
260	- a1(k,i+1) * f8 * ( 1.0 - cimf*cimf ) &
261	+ a2(k,i+1) * f24 * ( 1.0 - cimfcimfcimf )
262	ip = MAX( ip, 0.0 )
263	im = MAX( im, 0.0 )
264	ippb(k,i) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
265	impb(k,i) = im * MIN( 1.0, sk_p(k,j,i+1) / (ip+im+1E-15) )
266
267	cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
268	cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
269	cipf = 1.0 - 2.0 * cip
270	cimf = 1.0 - 2.0 * cim
271	ip = a0(k,i-1) * f2 * ( 1.0 - cipf ) &
272	+ a1(k,i-1) * f8 * ( 1.0 - cipf*cipf ) &
273	+ a2(k,i-1) * f24 * ( 1.0 - cipfcipfcipf )
274	im = a0(k,i) * f2 * ( 1.0 - cimf ) &
275	- a1(k,i) * f8 * ( 1.0 - cimf*cimf ) &
276	+ a2(k,i) * f24 * ( 1.0 - cimfcimfcimf )
277	ip = MAX( ip, 0.0 )
278	im = MAX( im, 0.0 )
279	ipmb(k,i) = ip * MIN( 1.0, sk_p(k,j,i-1) / (ip+im+1E-15) )
280	immb(k,i) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
281	ENDDO
282	ENDDO
283
284	!
285	!-- Compute monitor function m1
286	DO i = nxl-2, nxr+2
287	DO k = nzb+1, nzt
288	m1z = ABS( sk_p(k,j,i+1) - 2.0 * sk_p(k,j,i) + sk_p(k,j,i-1) )
289	m1n = ABS( sk_p(k,j,i+1) - sk_p(k,j,i-1) )
290	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
291	m1(k,i) = m1z / m1n
292	IF ( m1(k,i) /= 2.0 .AND. m1n < fmax(2) ) m1(k,i) = 0.0
293	ELSEIF ( m1n < m1z ) THEN
294	m1(k,i) = -1.0
295	ELSE
296	m1(k,i) = 0.0
297	ENDIF
298	ENDDO
299	ENDDO
300
301	!
302	!-- Compute switch sw
303	sw = 0.0
304	DO i = nxl-1, nxr+1
305	DO k = nzb+1, nzt
306	m2 = 2.0 * ABS( a1(k,i) - a12(k,i) ) / &
307	MAX( ABS( a1(k,i) + a12(k,i) ), 1E-35 )
308	IF ( ABS( a1(k,i) + a12(k,i) ) < fmax(2) ) m2 = 0.0
309
310	m3 = 2.0 * ABS( a2(k,i) - a22(k,i) ) / &
311	MAX( ABS( a2(k,i) + a22(k,i) ), 1E-35 )
312	IF ( ABS( a2(k,i) + a22(k,i) ) < fmax(1) ) m3 = 0.0
313
314	t1 = 0.35
315	t2 = 0.35
316	IF ( m1(k,i) == -1.0 ) t2 = 0.12
317
318	!-- *VOCL STMT,IF(10)
319	IF ( m1(k,i-1) == 1.0 .OR. m1(k,i) == 1.0 .OR. m1(k,i+1) == 1.0 &
320	.OR. m2 > t2 .OR. m3 > T2 .OR. &
321	( m1(k,i) > t1 .AND. m1(k,i-1) /= -1.0 .AND. &
322	m1(k,i) /= -1.0 .AND. m1(k,i+1) /= -1.0 ) &
323	) sw(k,i) = 1.0
324	ENDDO
325	ENDDO
326
327	!
328	!-- Fluxes using the exponential scheme
329	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
330	DO i = nxl, nxr
331	DO k = nzb+1, nzt
332
333	!-- *VOCL STMT,IF(10)
334	IF ( sw(k,i) == 1.0 ) THEN
335	snenn = sk_p(k,j,i+1) - sk_p(k,j,i-1)
336	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
337	sterm = ( sk_p(k,j,i) - sk_p(k,j,i-1) ) / snenn
338	sterm = MIN( sterm, 0.9999 )
339	sterm = MAX( sterm, 0.0001 )
340
341	ix = INT( sterm * 1000 ) + 1
342
343	cip = MAX( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
344
345	ippe(k,i) = sk_p(k,j,i-1) * cip + snenn * ( &
346	aex(ix) * cip + bex(ix) / dex(ix) * ( &
347	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
348	) &
349	)
350	IF ( sterm == 0.0001 ) ippe(k,i) = sk_p(k,j,i) * cip
351	IF ( sterm == 0.9999 ) ippe(k,i) = sk_p(k,j,i) * cip
352
353	snenn = sk_p(k,j,i-1) - sk_p(k,j,i+1)
354	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
355	sterm = ( sk_p(k,j,i) - sk_p(k,j,i+1) ) / snenn
356	sterm = MIN( sterm, 0.9999 )
357	sterm = MAX( sterm, 0.0001 )
358
359	ix = INT( sterm * 1000 ) + 1
360
361	cim = -MIN( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
362
363	imme(k,i) = sk_p(k,j,i+1) * cim + snenn * ( &
364	aex(ix) * cim + bex(ix) / dex(ix) * ( &
365	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
366	) &
367	)
368	IF ( sterm == 0.0001 ) imme(k,i) = sk_p(k,j,i) * cim
369	IF ( sterm == 0.9999 ) imme(k,i) = sk_p(k,j,i) * cim
370	ENDIF
371
372	!-- *VOCL STMT,IF(10)
373	IF ( sw(k,i+1) == 1.0 ) THEN
374	snenn = sk_p(k,j,i) - sk_p(k,j,i+2)
375	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
376	sterm = ( sk_p(k,j,i+1) - sk_p(k,j,i+2) ) / snenn
377	sterm = MIN( sterm, 0.9999 )
378	sterm = MAX( sterm, 0.0001 )
379
380	ix = INT( sterm * 1000 ) + 1
381
382	cim = -MIN( 0.0, ( u(k,j,i+1) - u_gtrans ) * dt_3d * ddx )
383
384	impe(k,i) = sk_p(k,j,i+2) * cim + snenn * ( &
385	aex(ix) * cim + bex(ix) / dex(ix) * ( &
386	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
387	) &
388	)
389	IF ( sterm == 0.0001 ) impe(k,i) = sk_p(k,j,i+1) * cim
390	IF ( sterm == 0.9999 ) impe(k,i) = sk_p(k,j,i+1) * cim
391	ENDIF
392
393	!-- *VOCL STMT,IF(10)
394	IF ( sw(k,i-1) == 1.0 ) THEN
395	snenn = sk_p(k,j,i) - sk_p(k,j,i-2)
396	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
397	sterm = ( sk_p(k,j,i-1) - sk_p(k,j,i-2) ) / snenn
398	sterm = MIN( sterm, 0.9999 )
399	sterm = MAX( sterm, 0.0001 )
400
401	ix = INT( sterm * 1000 ) + 1
402
403	cip = MAX( 0.0, ( u(k,j,i) - u_gtrans ) * dt_3d * ddx )
404
405	ipme(k,i) = sk_p(k,j,i-2) * cip + snenn * ( &
406	aex(ix) * cip + bex(ix) / dex(ix) * ( &
407	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
408	) &
409	)
410	IF ( sterm == 0.0001 ) ipme(k,i) = sk_p(k,j,i-1) * cip
411	IF ( sterm == 0.9999 ) ipme(k,i) = sk_p(k,j,i-1) * cip
412	ENDIF
413
414	ENDDO
415	ENDDO
416	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
417
418	!
419	!-- Prognostic equation
420	DO i = nxl, nxr
421	DO k = nzb+1, nzt
422	fplus = ( 1.0 - sw(k,i) ) * ippb(k,i) + sw(k,i) * ippe(k,i) &
423	- ( 1.0 - sw(k,i+1) ) * impb(k,i) - sw(k,i+1) * impe(k,i)
424	fminus = ( 1.0 - sw(k,i-1) ) * ipmb(k,i) + sw(k,i-1) * ipme(k,i) &
425	- ( 1.0 - sw(k,i) ) * immb(k,i) - sw(k,i) * imme(k,i)
426	tendenz = fplus - fminus
427	!
428	!-- Removed in order to optimize speed
429	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
430	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
431	!
432	!-- Density correction because of possible remaining divergences
433	d_new = d(k,j,i) - ( u(k,j,i+1) - u(k,j,i) ) * dt_3d * ddx
434	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
435	( 1.0 + d_new )
436	d(k,j,i) = d_new
437	ENDDO
438	ENDDO
439
440	ENDDO ! End of the advection in x-direction
441
442	!
443	!-- Deallocate temporary arrays
444	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
445	ippb, ippe, m1, sw )
446
447
448	!
449	!-- Enlarge boundary of local array cyclically in y-direction
450	#if defined( __parallel )
451	ngp = ( nzt - nzb + 6 ) * ( nyn - nys + 7 )
452	CALL MPI_TYPE_VECTOR( nxr-nxl+7, 3*(nzt-nzb+6), ngp, MPI_REAL, &
453	type_xz_2, ierr )
454	CALL MPI_TYPE_COMMIT( type_xz_2, ierr )
455	!
456	!-- Send front boundary, receive rear boundary
457	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
458	CALL MPI_SENDRECV( sk_p(nzb-2,nys,nxl-3), 1, type_xz_2, psouth, 0, &
459	sk_p(nzb-2,nyn+1,nxl-3), 1, type_xz_2, pnorth, 0, &
460	comm2d, status, ierr )
461	!
462	!-- Send rear boundary, receive front boundary
463	CALL MPI_SENDRECV( sk_p(nzb-2,nyn-2,nxl-3), 1, type_xz_2, pnorth, 1, &
464	sk_p(nzb-2,nys-3,nxl-3), 1, type_xz_2, psouth, 1, &
465	comm2d, status, ierr )
466	CALL MPI_TYPE_FREE( type_xz_2, ierr )
467	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'pause' )
468	#else
469	DO i = nxl, nxr
470	DO k = nzb+1, nzt
471	sk_p(k,nys-1,i) = sk_p(k,nyn,i)
472	sk_p(k,nys-2,i) = sk_p(k,nyn-1,i)
473	sk_p(k,nys-3,i) = sk_p(k,nyn-2,i)
474	sk_p(k,nyn+1,i) = sk_p(k,nys,i)
475	sk_p(k,nyn+2,i) = sk_p(k,nys+1,i)
476	sk_p(k,nyn+3,i) = sk_p(k,nys+2,i)
477	ENDDO
478	ENDDO
479	#endif
480
481	!
482	!-- Determine the maxima of the first and second derivative in y-direction
483	fmax_l = 0.0
484	DO i = nxl, nxr
485	DO j = nys, nyn
486	DO k = nzb+1, nzt
487	zaehler = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
488	nenner = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
489	fmax_l(1) = MAX( fmax_l(1) , zaehler )
490	fmax_l(2) = MAX( fmax_l(2) , nenner )
491	ENDDO
492	ENDDO
493	ENDDO
494	#if defined( __parallel )
495	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
496	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
497	#else
498	fmax = fmax_l
499	#endif
500
501	fmax = 0.04 * fmax
502
503	!
504	!-- Allocate temporary arrays
505	ALLOCATE( a0(nzb+1:nzt,nys-1:nyn+1), a1(nzb+1:nzt,nys-1:nyn+1), &
506	a2(nzb+1:nzt,nys-1:nyn+1), a12(nzb+1:nzt,nys-1:nyn+1), &
507	a22(nzb+1:nzt,nys-1:nyn+1), immb(nzb+1:nzt,nys-1:nyn+1), &
508	imme(nzb+1:nzt,nys-1:nyn+1), impb(nzb+1:nzt,nys-1:nyn+1), &
509	impe(nzb+1:nzt,nys-1:nyn+1), ipmb(nzb+1:nzt,nys-1:nyn+1), &
510	ipme(nzb+1:nzt,nys-1:nyn+1), ippb(nzb+1:nzt,nys-1:nyn+1), &
511	ippe(nzb+1:nzt,nys-1:nyn+1), m1(nzb+1:nzt,nys-2:nyn+2), &
512	sw(nzb+1:nzt,nys-1:nyn+1) &
513	)
514	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
515
516	!
517	!-- Outer loop of all i
518	DO i = nxl, nxr
519
520	!
521	!-- Compute polynomial coefficients
522	DO j = nys-1, nyn+1
523	DO k = nzb+1, nzt
524	a12(k,j) = 0.5 * ( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
525	a22(k,j) = 0.5 * ( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) &
526	+ sk_p(k,j-1,i) )
527	a0(k,j) = ( 9.0 * sk_p(k,j+2,i) - 116.0 * sk_p(k,j+1,i) &
528	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k,j-1,i) &
529	+ 9.0 * sk_p(k,j-2,i) ) * f1920
530	a1(k,j) = ( -5.0 * sk_p(k,j+2,i) + 34.0 * sk_p(k,j+1,i) &
531	- 34.0 * sk_p(k,j-1,i) + 5.0 * sk_p(k,j-2,i) &
532	) * f48
533	a2(k,j) = ( -3.0 * sk_p(k,j+2,i) + 36.0 * sk_p(k,j+1,i) &
534	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k,j-1,i) &
535	- 3.0 * sk_p(k,j-2,i) ) * f48
536	ENDDO
537	ENDDO
538
539	!
540	!-- Fluxes using the Bott scheme
541	!-- *VOCL LOOP,UNROLL(2)
542	DO j = nys, nyn
543	DO k = nzb+1, nzt
544	cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
545	cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
546	cipf = 1.0 - 2.0 * cip
547	cimf = 1.0 - 2.0 * cim
548	ip = a0(k,j) * f2 * ( 1.0 - cipf ) &
549	+ a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &
550	+ a2(k,j) * f24 * ( 1.0 - cipfcipfcipf )
551	im = a0(k,j+1) * f2 * ( 1.0 - cimf ) &
552	- a1(k,j+1) * f8 * ( 1.0 - cimf*cimf ) &
553	+ a2(k,j+1) * f24 * ( 1.0 - cimfcimfcimf )
554	ip = MAX( ip, 0.0 )
555	im = MAX( im, 0.0 )
556	ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
557	impb(k,j) = im * MIN( 1.0, sk_p(k,j+1,i) / (ip+im+1E-15) )
558
559	cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
560	cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
561	cipf = 1.0 - 2.0 * cip
562	cimf = 1.0 - 2.0 * cim
563	ip = a0(k,j-1) * f2 * ( 1.0 - cipf ) &
564	+ a1(k,j-1) * f8 * ( 1.0 - cipf*cipf ) &
565	+ a2(k,j-1) * f24 * ( 1.0 - cipfcipfcipf )
566	im = a0(k,j) * f2 * ( 1.0 - cimf ) &
567	- a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &
568	+ a2(k,j) * f24 * ( 1.0 - cimfcimfcimf )
569	ip = MAX( ip, 0.0 )
570	im = MAX( im, 0.0 )
571	ipmb(k,j) = ip * MIN( 1.0, sk_p(k,j-1,i) / (ip+im+1E-15) )
572	immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
573	ENDDO
574	ENDDO
575
576	!
577	!-- Compute monitor function m1
578	DO j = nys-2, nyn+2
579	DO k = nzb+1, nzt
580	m1z = ABS( sk_p(k,j+1,i) - 2.0 * sk_p(k,j,i) + sk_p(k,j-1,i) )
581	m1n = ABS( sk_p(k,j+1,i) - sk_p(k,j-1,i) )
582	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
583	m1(k,j) = m1z / m1n
584	IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0
585	ELSEIF ( m1n < m1z ) THEN
586	m1(k,j) = -1.0
587	ELSE
588	m1(k,j) = 0.0
589	ENDIF
590	ENDDO
591	ENDDO
592
593	!
594	!-- Compute switch sw
595	sw = 0.0
596	DO j = nys-1, nyn+1
597	DO k = nzb+1, nzt
598	m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
599	MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
600	IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0
601
602	m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
603	MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
604	IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0
605
606	t1 = 0.35
607	t2 = 0.35
608	IF ( m1(k,j) == -1.0 ) t2 = 0.12
609
610	!-- *VOCL STMT,IF(10)
611	IF ( m1(k,j-1) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k,j+1) == 1.0 &
612	.OR. m2 > t2 .OR. m3 > T2 .OR. &
613	( m1(k,j) > t1 .AND. m1(k,j-1) /= -1.0 .AND. &
614	m1(k,j) /= -1.0 .AND. m1(k,j+1) /= -1.0 ) &
615	) sw(k,j) = 1.0
616	ENDDO
617	ENDDO
618
619	!
620	!-- Fluxes using exponential scheme
621	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
622	DO j = nys, nyn
623	DO k = nzb+1, nzt
624
625	!-- *VOCL STMT,IF(10)
626	IF ( sw(k,j) == 1.0 ) THEN
627	snenn = sk_p(k,j+1,i) - sk_p(k,j-1,i)
628	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
629	sterm = ( sk_p(k,j,i) - sk_p(k,j-1,i) ) / snenn
630	sterm = MIN( sterm, 0.9999 )
631	sterm = MAX( sterm, 0.0001 )
632
633	ix = INT( sterm * 1000 ) + 1
634
635	cip = MAX( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
636
637	ippe(k,j) = sk_p(k,j-1,i) * cip + snenn * ( &
638	aex(ix) * cip + bex(ix) / dex(ix) * ( &
639	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
640	) &
641	)
642	IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip
643	IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip
644
645	snenn = sk_p(k,j-1,i) - sk_p(k,j+1,i)
646	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
647	sterm = ( sk_p(k,j,i) - sk_p(k,j+1,i) ) / snenn
648	sterm = MIN( sterm, 0.9999 )
649	sterm = MAX( sterm, 0.0001 )
650
651	ix = INT( sterm * 1000 ) + 1
652
653	cim = -MIN( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
654
655	imme(k,j) = sk_p(k,j+1,i) * cim + snenn * ( &
656	aex(ix) * cim + bex(ix) / dex(ix) * ( &
657	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
658	) &
659	)
660	IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim
661	IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim
662	ENDIF
663
664	!-- *VOCL STMT,IF(10)
665	IF ( sw(k,j+1) == 1.0 ) THEN
666	snenn = sk_p(k,j,i) - sk_p(k,j+2,i)
667	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
668	sterm = ( sk_p(k,j+1,i) - sk_p(k,j+2,i) ) / snenn
669	sterm = MIN( sterm, 0.9999 )
670	sterm = MAX( sterm, 0.0001 )
671
672	ix = INT( sterm * 1000 ) + 1
673
674	cim = -MIN( 0.0, ( v(k,j+1,i) - v_gtrans ) * dt_3d * ddy )
675
676	impe(k,j) = sk_p(k,j+2,i) * cim + snenn * ( &
677	aex(ix) * cim + bex(ix) / dex(ix) * ( &
678	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
679	) &
680	)
681	IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k,j+1,i) * cim
682	IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k,j+1,i) * cim
683	ENDIF
684
685	!-- *VOCL STMT,IF(10)
686	IF ( sw(k,j-1) == 1.0 ) THEN
687	snenn = sk_p(k,j,i) - sk_p(k,j-2,i)
688	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
689	sterm = ( sk_p(k,j-1,i) - sk_p(k,j-2,i) ) / snenn
690	sterm = MIN( sterm, 0.9999 )
691	sterm = MAX( sterm, 0.0001 )
692
693	ix = INT( sterm * 1000 ) + 1
694
695	cip = MAX( 0.0, ( v(k,j,i) - v_gtrans ) * dt_3d * ddy )
696
697	ipme(k,j) = sk_p(k,j-2,i) * cip + snenn * ( &
698	aex(ix) * cip + bex(ix) / dex(ix) * ( &
699	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
700	) &
701	)
702	IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k,j-1,i) * cip
703	IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k,j-1,i) * cip
704	ENDIF
705
706	ENDDO
707	ENDDO
708	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
709
710	!
711	!-- Prognostic equation
712	DO j = nys, nyn
713	DO k = nzb+1, nzt
714	fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &
715	- ( 1.0 - sw(k,j+1) ) * impb(k,j) - sw(k,j+1) * impe(k,j)
716	fminus = ( 1.0 - sw(k,j-1) ) * ipmb(k,j) + sw(k,j-1) * ipme(k,j) &
717	- ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)
718	tendenz = fplus - fminus
719	!
720	!-- Removed in order to optimise speed
721	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
722	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
723	!
724	!-- Density correction because of possible remaining divergences
725	d_new = d(k,j,i) - ( v(k,j+1,i) - v(k,j,i) ) * dt_3d * ddy
726	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
727	( 1.0 + d_new )
728	d(k,j,i) = d_new
729	ENDDO
730	ENDDO
731
732	ENDDO ! End of the advection in y-direction
733	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'continue' )
734	CALL cpu_log( log_point_s(11), 'advec_s_bc:sendrecv', 'stop' )
735
736	!
737	!-- Deallocate temporary arrays
738	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
739	ippb, ippe, m1, sw )
740
741
742	!
743	!-- Initialise for the computation of heat fluxes (see below; required in
744	!-- UP flow_statistics)
745	IF ( sk_char == 'pt' ) sums_wsts_bc_l = 0.0
746
747	!
748	!-- Add top and bottom boundaries according to the relevant boundary conditions
749	IF ( sk_char == 'pt' ) THEN
750
751	!
752	!-- Temperature boundary condition at the bottom boundary
753	IF ( ibc_pt_b == 0 ) THEN
754	!
755	!-- Dirichlet (fixed surface temperature)
756	DO i = nxl, nxr
757	DO j = nys, nyn
758	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
759	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
760	ENDDO
761	ENDDO
762
763	ELSE
764	!
765	!-- Neumann (i.e. here zero gradient)
766	DO i = nxl, nxr
767	DO j = nys, nyn
768	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
769	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
770	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
771	ENDDO
772	ENDDO
773
774	ENDIF
775
776	!
777	!-- Temperature boundary condition at the top boundary
778	IF ( ibc_pt_t == 0 .OR. ibc_pt_t == 1 ) THEN
779	!
780	!-- Dirichlet or Neumann (zero gradient)
781	DO i = nxl, nxr
782	DO j = nys, nyn
783	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
784	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
785	ENDDO
786	ENDDO
787
788	ELSEIF ( ibc_pt_t == 2 ) THEN
789	!
790	!-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
791	DO i = nxl, nxr
792	DO j = nys, nyn
793	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_pt_t_val * dzu(nzt+1)
794	sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_pt_t_val * dzu(nzt+1)
795	ENDDO
796	ENDDO
797
798	ENDIF
799
800	ELSEIF ( sk_char == 'sa' ) THEN
801
802	!
803	!-- Salinity boundary condition at the bottom boundary.
804	!-- So far, always Neumann (i.e. here zero gradient) is used
805	DO i = nxl, nxr
806	DO j = nys, nyn
807	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
808	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
809	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
810	ENDDO
811	ENDDO
812
813	!
814	!-- Salinity boundary condition at the top boundary.
815	!-- Dirichlet or Neumann (zero gradient)
816	DO i = nxl, nxr
817	DO j = nys, nyn
818	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
819	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
820	ENDDO
821	ENDDO
822
823	ELSEIF ( sk_char == 'q' ) THEN
824
825	!
826	!-- Specific humidity boundary condition at the bottom boundary.
827	!-- Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
828	DO i = nxl, nxr
829	DO j = nys, nyn
830	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
831	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
832	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
833	ENDDO
834	ENDDO
835
836	!
837	!-- Specific humidity boundary condition at the top boundary
838	IF ( ibc_q_t == 0 ) THEN
839	!
840	!-- Dirichlet
841	DO i = nxl, nxr
842	DO j = nys, nyn
843	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
844	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
845	ENDDO
846	ENDDO
847
848	ELSE
849	!
850	!-- Neumann: dzu(nzt+2:3) are not defined, dzu(nzt+1) is used instead
851	DO i = nxl, nxr
852	DO j = nys, nyn
853	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) + bc_q_t_val * dzu(nzt+1)
854	sk_p(nzt+3,j,i) = sk_p(nzt+2,j,i) + bc_q_t_val * dzu(nzt+1)
855	ENDDO
856	ENDDO
857
858	ENDIF
859
860	ELSEIF ( sk_char == 'e' ) THEN
861
862	!
863	!-- TKE boundary condition at bottom and top boundary (generally Neumann)
864	DO i = nxl, nxr
865	DO j = nys, nyn
866	sk_p(nzb,j,i) = sk_p(nzb+1,j,i)
867	sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
868	sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
869	sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
870	sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
871	ENDDO
872	ENDDO
873
874	ELSE
875
876	WRITE( message_string, * ) 'no vertical boundary condi', &
877	'tion for variable "', sk_char, '"'
878	CALL message( 'advec_s_bc', 'PA0158', 1, 2, 0, 6, 0 )
879
880	ENDIF
881
882	!
883	!-- Determine the maxima of the first and second derivative in z-direction
884	fmax_l = 0.0
885	DO i = nxl, nxr
886	DO j = nys, nyn
887	DO k = nzb, nzt+1
888	zaehler = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
889	nenner = ABS( sk_p(k+1,j,i+1) - sk_p(k-1,j,i) )
890	fmax_l(1) = MAX( fmax_l(1) , zaehler )
891	fmax_l(2) = MAX( fmax_l(2) , nenner )
892	ENDDO
893	ENDDO
894	ENDDO
895	#if defined( __parallel )
896	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
897	CALL MPI_ALLREDUCE( fmax_l, fmax, 2, MPI_REAL, MPI_MAX, comm2d, ierr )
898	#else
899	fmax = fmax_l
900	#endif
901
902	fmax = 0.04 * fmax
903
904	!
905	!-- Allocate temporary arrays
906	ALLOCATE( a0(nzb:nzt+1,nys:nyn), a1(nzb:nzt+1,nys:nyn), &
907	a2(nzb:nzt+1,nys:nyn), a12(nzb:nzt+1,nys:nyn), &
908	a22(nzb:nzt+1,nys:nyn), immb(nzb+1:nzt,nys:nyn), &
909	imme(nzb+1:nzt,nys:nyn), impb(nzb+1:nzt,nys:nyn), &
910	impe(nzb+1:nzt,nys:nyn), ipmb(nzb+1:nzt,nys:nyn), &
911	ipme(nzb+1:nzt,nys:nyn), ippb(nzb+1:nzt,nys:nyn), &
912	ippe(nzb+1:nzt,nys:nyn), m1(nzb-1:nzt+2,nys:nyn), &
913	sw(nzb:nzt+1,nys:nyn) &
914	)
915	imme = 0.0; impe = 0.0; ipme = 0.0; ippe = 0.0
916
917	!
918	!-- Outer loop of all i
919	DO i = nxl, nxr
920
921	!
922	!-- Compute polynomial coefficients
923	DO j = nys, nyn
924	DO k = nzb, nzt+1
925	a12(k,j) = 0.5 * ( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
926	a22(k,j) = 0.5 * ( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) &
927	+ sk_p(k-1,j,i) )
928	a0(k,j) = ( 9.0 * sk_p(k+2,j,i) - 116.0 * sk_p(k+1,j,i) &
929	+ 2134.0 * sk_p(k,j,i) - 116.0 * sk_p(k-1,j,i) &
930	+ 9.0 * sk_p(k-2,j,i) ) * f1920
931	a1(k,j) = ( -5.0 * sk_p(k+2,j,i) + 34.0 * sk_p(k+1,j,i) &
932	- 34.0 * sk_p(k-1,j,i) + 5.0 * sk_p(k-2,j,i) &
933	) * f48
934	a2(k,j) = ( -3.0 * sk_p(k+2,j,i) + 36.0 * sk_p(k+1,j,i) &
935	- 66.0 * sk_p(k,j,i) + 36.0 * sk_p(k-1,j,i) &
936	- 3.0 * sk_p(k-2,j,i) ) * f48
937	ENDDO
938	ENDDO
939
940	!
941	!-- Fluxes using the Bott scheme
942	!-- *VOCL LOOP,UNROLL(2)
943	DO j = nys, nyn
944	DO k = nzb+1, nzt
945	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
946	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
947	cipf = 1.0 - 2.0 * cip
948	cimf = 1.0 - 2.0 * cim
949	ip = a0(k,j) * f2 * ( 1.0 - cipf ) &
950	+ a1(k,j) * f8 * ( 1.0 - cipf*cipf ) &
951	+ a2(k,j) * f24 * ( 1.0 - cipfcipfcipf )
952	im = a0(k+1,j) * f2 * ( 1.0 - cimf ) &
953	- a1(k+1,j) * f8 * ( 1.0 - cimf*cimf ) &
954	+ a2(k+1,j) * f24 * ( 1.0 - cimfcimfcimf )
955	ip = MAX( ip, 0.0 )
956	im = MAX( im, 0.0 )
957	ippb(k,j) = ip * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
958	impb(k,j) = im * MIN( 1.0, sk_p(k+1,j,i) / (ip+im+1E-15) )
959
960	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
961	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
962	cipf = 1.0 - 2.0 * cip
963	cimf = 1.0 - 2.0 * cim
964	ip = a0(k-1,j) * f2 * ( 1.0 - cipf ) &
965	+ a1(k-1,j) * f8 * ( 1.0 - cipf*cipf ) &
966	+ a2(k-1,j) * f24 * ( 1.0 - cipfcipfcipf )
967	im = a0(k,j) * f2 * ( 1.0 - cimf ) &
968	- a1(k,j) * f8 * ( 1.0 - cimf*cimf ) &
969	+ a2(k,j) * f24 * ( 1.0 - cimfcimfcimf )
970	ip = MAX( ip, 0.0 )
971	im = MAX( im, 0.0 )
972	ipmb(k,j) = ip * MIN( 1.0, sk_p(k-1,j,i) / (ip+im+1E-15) )
973	immb(k,j) = im * MIN( 1.0, sk_p(k,j,i) / (ip+im+1E-15) )
974	ENDDO
975	ENDDO
976
977	!
978	!-- Compute monitor function m1
979	DO j = nys, nyn
980	DO k = nzb-1, nzt+2
981	m1z = ABS( sk_p(k+1,j,i) - 2.0 * sk_p(k,j,i) + sk_p(k-1,j,i) )
982	m1n = ABS( sk_p(k+1,j,i) - sk_p(k-1,j,i) )
983	IF ( m1n /= 0.0 .AND. m1n >= m1z ) THEN
984	m1(k,j) = m1z / m1n
985	IF ( m1(k,j) /= 2.0 .AND. m1n < fmax(2) ) m1(k,j) = 0.0
986	ELSEIF ( m1n < m1z ) THEN
987	m1(k,j) = -1.0
988	ELSE
989	m1(k,j) = 0.0
990	ENDIF
991	ENDDO
992	ENDDO
993
994	!
995	!-- Compute switch sw
996	sw = 0.0
997	DO j = nys, nyn
998	DO k = nzb, nzt+1
999	m2 = 2.0 * ABS( a1(k,j) - a12(k,j) ) / &
1000	MAX( ABS( a1(k,j) + a12(k,j) ), 1E-35 )
1001	IF ( ABS( a1(k,j) + a12(k,j) ) < fmax(2) ) m2 = 0.0
1002
1003	m3 = 2.0 * ABS( a2(k,j) - a22(k,j) ) / &
1004	MAX( ABS( a2(k,j) + a22(k,j) ), 1E-35 )
1005	IF ( ABS( a2(k,j) + a22(k,j) ) < fmax(1) ) m3 = 0.0
1006
1007	t1 = 0.35
1008	t2 = 0.35
1009	IF ( m1(k,j) == -1.0 ) t2 = 0.12
1010
1011	!-- *VOCL STMT,IF(10)
1012	IF ( m1(k-1,j) == 1.0 .OR. m1(k,j) == 1.0 .OR. m1(k+1,j) == 1.0 &
1013	.OR. m2 > t2 .OR. m3 > T2 .OR. &
1014	( m1(k,j) > t1 .AND. m1(k-1,j) /= -1.0 .AND. &
1015	m1(k,j) /= -1.0 .AND. m1(k+1,j) /= -1.0 ) &
1016	) sw(k,j) = 1.0
1017	ENDDO
1018	ENDDO
1019
1020	!
1021	!-- Fluxes using exponential scheme
1022	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1023	DO j = nys, nyn
1024	DO k = nzb+1, nzt
1025
1026	!-- *VOCL STMT,IF(10)
1027	IF ( sw(k,j) == 1.0 ) THEN
1028	snenn = sk_p(k+1,j,i) - sk_p(k-1,j,i)
1029	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
1030	sterm = ( sk_p(k,j,i) - sk_p(k-1,j,i) ) / snenn
1031	sterm = MIN( sterm, 0.9999 )
1032	sterm = MAX( sterm, 0.0001 )
1033
1034	ix = INT( sterm * 1000 ) + 1
1035
1036	cip = MAX( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
1037
1038	ippe(k,j) = sk_p(k-1,j,i) * cip + snenn * ( &
1039	aex(ix) * cip + bex(ix) / dex(ix) * ( &
1040	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
1041	) &
1042	)
1043	IF ( sterm == 0.0001 ) ippe(k,j) = sk_p(k,j,i) * cip
1044	IF ( sterm == 0.9999 ) ippe(k,j) = sk_p(k,j,i) * cip
1045
1046	snenn = sk_p(k-1,j,i) - sk_p(k+1,j,i)
1047	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
1048	sterm = ( sk_p(k,j,i) - sk_p(k+1,j,i) ) / snenn
1049	sterm = MIN( sterm, 0.9999 )
1050	sterm = MAX( sterm, 0.0001 )
1051
1052	ix = INT( sterm * 1000 ) + 1
1053
1054	cim = -MIN( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
1055
1056	imme(k,j) = sk_p(k+1,j,i) * cim + snenn * ( &
1057	aex(ix) * cim + bex(ix) / dex(ix) * ( &
1058	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
1059	) &
1060	)
1061	IF ( sterm == 0.0001 ) imme(k,j) = sk_p(k,j,i) * cim
1062	IF ( sterm == 0.9999 ) imme(k,j) = sk_p(k,j,i) * cim
1063	ENDIF
1064
1065	!-- *VOCL STMT,IF(10)
1066	IF ( sw(k+1,j) == 1.0 ) THEN
1067	snenn = sk_p(k,j,i) - sk_p(k+2,j,i)
1068	IF ( ABS( snenn ) .LT. 1E-9 ) snenn = 1E-9
1069	sterm = ( sk_p(k+1,j,i) - sk_p(k+2,j,i) ) / snenn
1070	sterm = MIN( sterm, 0.9999 )
1071	sterm = MAX( sterm, 0.0001 )
1072
1073	ix = INT( sterm * 1000 ) + 1
1074
1075	cim = -MIN( 0.0, w(k,j,i) * dt_3d * ddzw(k) )
1076
1077	impe(k,j) = sk_p(k+2,j,i) * cim + snenn * ( &
1078	aex(ix) * cim + bex(ix) / dex(ix) * ( &
1079	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cim ) ) &
1080	) &
1081	)
1082	IF ( sterm == 0.0001 ) impe(k,j) = sk_p(k+1,j,i) * cim
1083	IF ( sterm == 0.9999 ) impe(k,j) = sk_p(k+1,j,i) * cim
1084	ENDIF
1085
1086	!-- *VOCL STMT,IF(10)
1087	IF ( sw(k-1,j) == 1.0 ) THEN
1088	snenn = sk_p(k,j,i) - sk_p(k-2,j,i)
1089	IF ( ABS( snenn ) < 1E-9 ) snenn = 1E-9
1090	sterm = ( sk_p(k-1,j,i) - sk_p(k-2,j,i) ) / snenn
1091	sterm = MIN( sterm, 0.9999 )
1092	sterm = MAX( sterm, 0.0001 )
1093
1094	ix = INT( sterm * 1000 ) + 1
1095
1096	cip = MAX( 0.0, w(k-1,j,i) * dt_3d * ddzw(k) )
1097
1098	ipme(k,j) = sk_p(k-2,j,i) * cip + snenn * ( &
1099	aex(ix) * cip + bex(ix) / dex(ix) * ( &
1100	eex(ix) - EXP( dex(ix)0.5 ( 1.0 - 2.0 * cip ) ) &
1101	) &
1102	)
1103	IF ( sterm == 0.0001 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1104	IF ( sterm == 0.9999 ) ipme(k,j) = sk_p(k-1,j,i) * cip
1105	ENDIF
1106
1107	ENDDO
1108	ENDDO
1109	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'pause' )
1110
1111	!
1112	!-- Prognostic equation
1113	DO j = nys, nyn
1114	DO k = nzb+1, nzt
1115	fplus = ( 1.0 - sw(k,j) ) * ippb(k,j) + sw(k,j) * ippe(k,j) &
1116	- ( 1.0 - sw(k+1,j) ) * impb(k,j) - sw(k+1,j) * impe(k,j)
1117	fminus = ( 1.0 - sw(k-1,j) ) * ipmb(k,j) + sw(k-1,j) * ipme(k,j) &
1118	- ( 1.0 - sw(k,j) ) * immb(k,j) - sw(k,j) * imme(k,j)
1119	tendenz = fplus - fminus
1120	!
1121	!-- Removed in order to optimise speed
1122	! ffmax = MAX( ABS( fplus ), ABS( fminus ), 1E-35 )
1123	! IF ( ( ABS( tendenz ) / ffmax ) < 1E-7 ) tendenz = 0.0
1124	!
1125	!-- Density correction because of possible remaining divergences
1126	d_new = d(k,j,i) - ( w(k,j,i) - w(k-1,j,i) ) * dt_3d * ddzw(k)
1127	sk_p(k,j,i) = ( ( 1.0 + d(k,j,i) ) * sk_p(k,j,i) - tendenz ) / &
1128	( 1.0 + d_new )
1129	!
1130	!-- Store heat flux for subsequent statistics output.
1131	!-- array m1 is here used as temporary storage
1132	m1(k,j) = fplus / dt_3d * dzw(k)
1133	ENDDO
1134	ENDDO
1135
1136	!
1137	!-- Sum up heat flux in order to order to obtain horizontal averages
1138	IF ( sk_char == 'pt' ) THEN
1139	DO sr = 0, statistic_regions
1140	DO j = nys, nyn
1141	DO k = nzb+1, nzt
1142	sums_wsts_bc_l(k,sr) = sums_wsts_bc_l(k,sr) + &
1143	m1(k,j) * rmask(j,i,sr)
1144	ENDDO
1145	ENDDO
1146	ENDDO
1147	ENDIF
1148
1149	ENDDO ! End of the advection in z-direction
1150	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'continue' )
1151	CALL cpu_log( log_point_s(12), 'advec_s_bc:exp', 'stop' )
1152
1153	!
1154	!-- Deallocate temporary arrays
1155	DEALLOCATE( a0, a1, a2, a12, a22, immb, imme, impb, impe, ipmb, ipme, &
1156	ippb, ippe, m1, sw )
1157
1158	!
1159	!-- Store results as tendency and deallocate local array
1160	DO i = nxl, nxr
1161	DO j = nys, nyn
1162	DO k = nzb+1, nzt
1163	tend(k,j,i) = tend(k,j,i) + ( sk_p(k,j,i) - sk(k,j,i) ) / dt_3d
1164	ENDDO
1165	ENDDO
1166	ENDDO
1167
1168	DEALLOCATE( sk_p )
1169
1170	END SUBROUTINE advec_s_bc

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |