Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |