Home

Context Navigation

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

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