Home

Context Navigation

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

Last change on this file since 4434 was 4429, checked in by raasch, 5 years ago
serial (non-MPI) test case added, several bugfixes for the serial mode
Property svn:keywords set to `Id`
File size: 55.7 KB

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