Home

Context Navigation

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

Last change on this file since 251 was 247, checked in by heinze, 16 years ago
Output of messages replaced by message handling routin
Property svn:keywords set to `Id`
File size: 43.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |