Home

Context Navigation

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

Last change on this file since 76 was 63, checked in by raasch, 18 years ago
preliminary changes concerning update of BC-scheme
Property svn:keywords set to `Id`
File size: 42.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |