Home

Context Navigation

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

Last change on this file since 1684 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |