Home

Context Navigation

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

Last change on this file since 1053 was 1037, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 44.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |