Home

Context Navigation

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

Last change on this file since 1961 was 1961, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 55.7 KB

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