Home

Context Navigation

advec_s_bc.f90 @ 96

Last change on this file since 96 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 42.2 KB

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

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

Context Navigation

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

Download in other formats: