Home

Context Navigation

advec_s_bc.f90 @ 4869

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

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

Context Navigation

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

Download in other formats: