Home

Context Navigation

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

Last change on this file since 1832 was 1818, checked in by maronga, 9 years ago
last commit documented / copyright update
Property svn:keywords set to `Id`
File size: 54.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |