Home

Context Navigation

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

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