Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 55.8 KB

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