Home

Context Navigation

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

Last change on this file since 1887 was 1874, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.5 KB

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