Home

Context Navigation

advec_s_bc_mod.f90 @ 1857

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

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

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

Context Navigation

source: palm/trunk/SOURCE/advec_s_bc_mod.f90 @ 1857

Download in other formats: