Home

Context Navigation

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

Last change on this file since 203 was 98, checked in by raasch, 17 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 42.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |