Home

Context Navigation

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

Last change on this file since 29 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 42.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |