Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |