Home

Context Navigation

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

Last change on this file since 1233 was 1093, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 44.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |