Home

Context Navigation

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

Last change on this file since 3683 was 3655, checked in by knoop, 6 years ago
Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update
Property svn:keywords set to `Id`
File size: 57.5 KB

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