Home

Context Navigation

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

Last change on this file since 1 was 1, checked in by raasch, 17 years ago
Initial repository layout and content
File size: 44.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |