Home

Context Navigation

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

Last change on this file since 219 was 216, checked in by raasch, 16 years ago
reading mechanism for restart files completely revised
Property svn:keywords set to `Id`
File size: 43.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |