Home

Context Navigation

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

Last change on this file since 1010 was 1010, checked in by raasch, 12 years ago
pointer free version can be generated with cpp switch nopointer
Property svn:keywords set to `Id`
File size: 43.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |