Home

Context Navigation

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

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