Home

Context Navigation

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

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