Home

Context Navigation

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

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