Home

Context Navigation

advec_s_bc.f90 @ 622

Last change on this file since 622 was 622, checked in by raasch, 13 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:

Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:

Property svn:keywords set to Id

File size: 43.4 KB

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

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

Context Navigation

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

Download in other formats: