Home

Context Navigation

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

Last change on this file since 1355 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 49.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |