Home

Context Navigation

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

Last change on this file since 1595 was 1518, checked in by hoffmann, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 53.9 KB

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