Home

Context Navigation

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

Last change on this file since 156 was 98, checked in by raasch, 17 years ago
updating comments and rc-file
Property svn:keywords set to `Id`
File size: 42.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |