Home

Context Navigation

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

Last change on this file since 703 was 623, checked in by raasch, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |