Home

Context Navigation

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

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