Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |