Home

Context Navigation

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

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