Home

Context Navigation

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

Last change on this file since 221 was 216, checked in by raasch, 16 years ago
reading mechanism for restart files completely revised
Property svn:keywords set to `Id`
File size: 43.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |