Home

Context Navigation

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

Last change on this file since 2321 was 2300, checked in by raasch, 7 years ago
NEC related code partly removed, host variable partly removed, host specific code completely removed, default values for host, loop_optimization and termination time_needed changed
Property svn:keywords set to `Id`
File size: 57.4 KB

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