Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |