Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |