Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |