Home

Context Navigation

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

Last change on this file since 226 was 226, checked in by raasch, 16 years ago
preparations for the next release
Property svn:keywords set to `Id`
File size: 43.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |