Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |