Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |