Home

Context Navigation

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

Last change on this file since 4181 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 55.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |