Home

Context Navigation

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

Last change on this file since 4495 was 4488, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 57.2 KB

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