Home

Context Navigation

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

Last change on this file since 260 was 247, checked in by heinze, 16 years ago
Output of messages replaced by message handling routin
Property svn:keywords set to `Id`
File size: 43.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |