Home

Context Navigation

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

Last change on this file since 240 was 226, checked in by raasch, 16 years ago
preparations for the next release
Property svn:keywords set to `Id`
File size: 43.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |