Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 55.8 KB

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