Home

Context Navigation

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

Last change on this file since 1330 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 47.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |