Home

Context Navigation

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

Last change on this file since 1036 was 1036, checked in by raasch, 12 years ago
code has been put under the GNU General Public License (v3)
Property svn:keywords set to `Id`
File size: 44.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |