Home

Context Navigation

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

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