Home

Context Navigation

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

Last change on this file since 870 was 623, checked in by raasch, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |