Home

Context Navigation

advec_s_bc.f90 @ 97

Last change on this file since 97 was 97, checked in by raasch, 17 years ago

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:

inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:

Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

Property svn:keywords set to Id

File size: 42.7 KB

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

Download in other formats: