Home

Context Navigation

source: palm/trunk/SOURCE/wall_fluxes.f90 @ 1016

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 32.5 KB

Line
1	MODULE wall_fluxes_mod
2	!------------------------------------------------------------------------------!
3	! Current revisions:
4	! -----------------
5	! accelerator version (*_acc) added
6	!
7	! Former revisions:
8	! -----------------
9	! $Id: wall_fluxes.f90 1015 2012-09-27 09:23:24Z raasch $
10	!
11	! 187 2008-08-06 16:25:09Z letzel
12	! Bugfix: Modification of the evaluation of the vertical turbulent momentum
13	! fluxes u'w' and v'w (see prandtl_fluxes), this requires the calculation of
14	! us_wall (and vel_total, u_i, v_i, ws) also in wall_fluxes_e.
15	! Bugfix: change definition of us_wall from 1D to 2D
16	! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
17	! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
18	! consistency with subroutine wall_fluxes
19	! Change: Modification of the integrated version of the profile function for
20	! momentum for unstable stratification
21	!
22	! Initial version (2007/03/07)
23	!
24	! Description:
25	! ------------
26	! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
27	! similarity.
28	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
29	! The all-gridpoint version of wall_fluxes_e is not used so far, because
30	! it gives slightly different results from the ij-version for some unknown
31	! reason.
32	!------------------------------------------------------------------------------!
33	PRIVATE
34	PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
35
36	INTERFACE wall_fluxes
37	MODULE PROCEDURE wall_fluxes
38	MODULE PROCEDURE wall_fluxes_ij
39	END INTERFACE wall_fluxes
40
41	INTERFACE wall_fluxes_acc
42	MODULE PROCEDURE wall_fluxes_acc
43	END INTERFACE wall_fluxes_acc
44
45	INTERFACE wall_fluxes_e
46	MODULE PROCEDURE wall_fluxes_e
47	MODULE PROCEDURE wall_fluxes_e_ij
48	END INTERFACE wall_fluxes_e
49
50	INTERFACE wall_fluxes_e_acc
51	MODULE PROCEDURE wall_fluxes_e_acc
52	END INTERFACE wall_fluxes_e_acc
53
54	CONTAINS
55
56	!------------------------------------------------------------------------------!
57	! Call for all grid points
58	!------------------------------------------------------------------------------!
59	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
60	nzb_uvw_outer, wall )
61
62	USE arrays_3d
63	USE control_parameters
64	USE grid_variables
65	USE indices
66	USE statistics
67
68	IMPLICIT NONE
69
70	INTEGER :: i, j, k, wall_index
71
72	INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) :: nzb_uvw_inner, &
73	nzb_uvw_outer
74	REAL :: a, b, c1, c2, h1, h2, zp
75	REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
76
77	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall
78	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux
79
80
81	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
82	wall_flux = 0.0
83	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
84
85	DO i = nxl, nxr
86	DO j = nys, nyn
87
88	IF ( wall(j,i) /= 0.0 ) THEN
89	!
90	!-- All subsequent variables are computed for the respective
91	!-- location where the respective flux is defined.
92	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
93
94	!
95	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
96	rifs = rif_wall(k,j,i,wall_index)
97
98	u_i = a * u(k,j,i) + c1 * 0.25 * &
99	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
100
101	v_i = b * v(k,j,i) + c2 * 0.25 * &
102	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
103
104	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
105	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
106	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
107	)
108	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + &
109	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
110
111	pts = pt_i - hom(k,1,4,0)
112	wspts = ws * pts
113
114	!
115	!-- (2) Compute wall-parallel absolute velocity vel_total
116	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
117
118	!
119	!-- (3) Compute wall friction velocity us_wall
120	IF ( rifs >= 0.0 ) THEN
121
122	!
123	!-- Stable stratification (and neutral)
124	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
125	5.0 * rifs * ( zp - z0(j,i) ) / zp &
126	)
127	ELSE
128
129	!
130	!-- Unstable stratification
131	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
132	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
133
134	us_wall = kappa * vel_total / ( &
135	LOG( zp / z0(j,i) ) - &
136	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
137	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
138	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
139	)
140	ENDIF
141
142	!
143	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
144	!-- number rifs)
145	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
146	( us_wall**3 + 1E-30 ) )
147
148	!
149	!-- Limit the value range of the Richardson numbers.
150	!-- This is necessary for very small velocities (u,w --> 0),
151	!-- because the absolute value of rif can then become very
152	!-- large, which in consequence would result in very large
153	!-- shear stresses and very small momentum fluxes (both are
154	!-- generally unrealistic).
155	IF ( rifs < rif_min ) rifs = rif_min
156	IF ( rifs > rif_max ) rifs = rif_max
157
158	!
159	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
160	IF ( rifs >= 0.0 ) THEN
161
162	!
163	!-- Stable stratification (and neutral)
164	wall_flux(k,j,i) = kappa * &
165	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
166	( LOG( zp / z0(j,i) ) + &
167	5.0 * rifs * ( zp - z0(j,i) ) / zp &
168	)
169	ELSE
170
171	!
172	!-- Unstable stratification
173	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
174	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
175
176	wall_flux(k,j,i) = kappa * &
177	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
178	LOG( zp / z0(j,i) ) - &
179	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
180	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
181	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
182	)
183	ENDIF
184	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
185
186	!
187	!-- store rifs for next time step
188	rif_wall(k,j,i,wall_index) = rifs
189
190	ENDDO
191
192	ENDIF
193
194	ENDDO
195	ENDDO
196
197	END SUBROUTINE wall_fluxes
198
199
200	!------------------------------------------------------------------------------!
201	! Call for all grid points - accelerator version
202	!------------------------------------------------------------------------------!
203	SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
204	nzb_uvw_outer, wall )
205
206	USE arrays_3d
207	USE control_parameters
208	USE grid_variables
209	USE indices
210	USE statistics
211
212	IMPLICIT NONE
213
214	INTEGER :: i, j, k, max_outer, min_inner, wall_index
215
216	INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) :: nzb_uvw_inner, &
217	nzb_uvw_outer
218	REAL :: a, b, c1, c2, h1, h2, zp
219	REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
220
221	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall
222	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux
223
224
225	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
226	wall_flux = 0.0
227	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
228
229	min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
230	max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
231
232	!$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
233	!$acc present( u, v, w, wall, wall_flux, z0 )
234	!$acc loop
235	DO i = nxl, nxr
236	DO j = nys, nyn
237	!$acc loop vector( 32 )
238	DO k = min_inner, max_outer
239	!
240	!-- All subsequent variables are computed for the respective
241	!-- location where the respective flux is defined.
242	IF ( k >= nzb_uvw_inner(j,i)+1 .AND. &
243	k <= nzb_uvw_outer(j,i) .AND. wall(j,i) /= 0.0 ) THEN
244	!
245	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
246	rifs = rif_wall(k,j,i,wall_index)
247
248	u_i = a * u(k,j,i) + c1 * 0.25 * &
249	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
250
251	v_i = b * v(k,j,i) + c2 * 0.25 * &
252	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
253
254	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
255	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
256	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
257	)
258	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + &
259	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
260
261	pts = pt_i - hom(k,1,4,0)
262	wspts = ws * pts
263
264	!
265	!-- (2) Compute wall-parallel absolute velocity vel_total
266	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
267
268	!
269	!-- (3) Compute wall friction velocity us_wall
270	IF ( rifs >= 0.0 ) THEN
271
272	!
273	!-- Stable stratification (and neutral)
274	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
275	5.0 * rifs * ( zp - z0(j,i) ) / zp &
276	)
277	ELSE
278
279	!
280	!-- Unstable stratification
281	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
282	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
283
284	us_wall = kappa * vel_total / ( &
285	LOG( zp / z0(j,i) ) - &
286	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
287	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
288	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
289	)
290	ENDIF
291
292	!
293	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
294	!-- number rifs)
295	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
296	( us_wall**3 + 1E-30 ) )
297
298	!
299	!-- Limit the value range of the Richardson numbers.
300	!-- This is necessary for very small velocities (u,w --> 0),
301	!-- because the absolute value of rif can then become very
302	!-- large, which in consequence would result in very large
303	!-- shear stresses and very small momentum fluxes (both are
304	!-- generally unrealistic).
305	IF ( rifs < rif_min ) rifs = rif_min
306	IF ( rifs > rif_max ) rifs = rif_max
307
308	!
309	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
310	IF ( rifs >= 0.0 ) THEN
311
312	!
313	!-- Stable stratification (and neutral)
314	wall_flux(k,j,i) = kappa * &
315	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
316	( LOG( zp / z0(j,i) ) + &
317	5.0 * rifs * ( zp - z0(j,i) ) / zp &
318	)
319	ELSE
320
321	!
322	!-- Unstable stratification
323	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
324	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
325
326	wall_flux(k,j,i) = kappa * &
327	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
328	LOG( zp / z0(j,i) ) - &
329	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
330	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
331	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
332	)
333	ENDIF
334	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
335
336	!
337	!-- store rifs for next time step
338	rif_wall(k,j,i,wall_index) = rifs
339
340	ENDIF
341
342	ENDDO
343	ENDDO
344	ENDDO
345	!$acc end kernels
346
347	END SUBROUTINE wall_fluxes_acc
348
349
350	!------------------------------------------------------------------------------!
351	! Call for all grid point i,j
352	!------------------------------------------------------------------------------!
353	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
354
355	USE arrays_3d
356	USE control_parameters
357	USE grid_variables
358	USE indices
359	USE statistics
360
361	IMPLICIT NONE
362
363	INTEGER :: i, j, k, nzb_w, nzt_w, wall_index
364	REAL :: a, b, c1, c2, h1, h2, zp
365
366	REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
367
368	REAL, DIMENSION(nzb:nzt+1) :: wall_flux
369
370
371	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
372	wall_flux = 0.0
373	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
374
375	!
376	!-- All subsequent variables are computed for the respective location where
377	!-- the respective flux is defined.
378	DO k = nzb_w, nzt_w
379
380	!
381	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
382	rifs = rif_wall(k,j,i,wall_index)
383
384	u_i = a * u(k,j,i) + c1 * 0.25 * &
385	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
386
387	v_i = b * v(k,j,i) + c2 * 0.25 * &
388	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
389
390	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
391	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
392	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
393	)
394	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
395	+ ( c1 + c2 ) * pt(k+1,j,i) )
396
397	pts = pt_i - hom(k,1,4,0)
398	wspts = ws * pts
399
400	!
401	!-- (2) Compute wall-parallel absolute velocity vel_total
402	vel_total = SQRT( ws*2 + ( a+c1 ) u_i*2 + ( b+c2 ) v_i**2 )
403
404	!
405	!-- (3) Compute wall friction velocity us_wall
406	IF ( rifs >= 0.0 ) THEN
407
408	!
409	!-- Stable stratification (and neutral)
410	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
411	5.0 * rifs * ( zp - z0(j,i) ) / zp &
412	)
413	ELSE
414
415	!
416	!-- Unstable stratification
417	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
418	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
419
420	us_wall = kappa * vel_total / ( &
421	LOG( zp / z0(j,i) ) - &
422	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
423	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
424	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
425	)
426	ENDIF
427
428	!
429	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
430	!-- rifs)
431	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
432
433	!
434	!-- Limit the value range of the Richardson numbers.
435	!-- This is necessary for very small velocities (u,w --> 0), because
436	!-- the absolute value of rif can then become very large, which in
437	!-- consequence would result in very large shear stresses and very
438	!-- small momentum fluxes (both are generally unrealistic).
439	IF ( rifs < rif_min ) rifs = rif_min
440	IF ( rifs > rif_max ) rifs = rif_max
441
442	!
443	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
444	IF ( rifs >= 0.0 ) THEN
445
446	!
447	!-- Stable stratification (and neutral)
448	wall_flux(k) = kappa * &
449	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
450	( LOG( zp / z0(j,i) ) + &
451	5.0 * rifs * ( zp - z0(j,i) ) / zp &
452	)
453	ELSE
454
455	!
456	!-- Unstable stratification
457	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
458	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
459
460	wall_flux(k) = kappa * &
461	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
462	LOG( zp / z0(j,i) ) - &
463	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
464	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
465	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
466	)
467	ENDIF
468	wall_flux(k) = -wall_flux(k) * us_wall
469
470	!
471	!-- store rifs for next time step
472	rif_wall(k,j,i,wall_index) = rifs
473
474	ENDDO
475
476	END SUBROUTINE wall_fluxes_ij
477
478
479
480	!------------------------------------------------------------------------------!
481	! Call for all grid points
482	!------------------------------------------------------------------------------!
483	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
484
485	!------------------------------------------------------------------------------!
486	! Description:
487	! ------------
488	! Calculates momentum fluxes at vertical walls for routine production_e
489	! assuming Monin-Obukhov similarity.
490	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
491	!------------------------------------------------------------------------------!
492
493	USE arrays_3d
494	USE control_parameters
495	USE grid_variables
496	USE indices
497	USE statistics
498
499	IMPLICIT NONE
500
501	INTEGER :: i, j, k, kk, wall_index
502	REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
503	ws, zp
504
505	REAL :: rifs
506
507	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall
508	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux
509
510
511	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
512	wall_flux = 0.0
513	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
514
515	DO i = nxl, nxr
516	DO j = nys, nyn
517
518	IF ( wall(j,i) /= 0.0 ) THEN
519	!
520	!-- All subsequent variables are computed for scalar locations.
521	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
522	!
523	!-- (1) Compute rifs, u_i, v_i, and ws
524	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
525	kk = nzb_diff_s_inner(j,i)-1
526	ELSE
527	kk = k-1
528	ENDIF
529	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
530	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
531	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
532	)
533
534	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
535	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
536	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
537	!
538	!-- (2) Compute wall-parallel absolute velocity vel_total and
539	!-- interpolate appropriate velocity component vel_zp.
540	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
541	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
542	!
543	!-- (3) Compute wall friction velocity us_wall
544	IF ( rifs >= 0.0 ) THEN
545
546	!
547	!-- Stable stratification (and neutral)
548	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
549	5.0 * rifs * ( zp - z0(j,i) ) / zp &
550	)
551	ELSE
552
553	!
554	!-- Unstable stratification
555	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
556	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
557
558	us_wall = kappa * vel_total / ( &
559	LOG( zp / z0(j,i) ) - &
560	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
561	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
562	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
563	)
564	ENDIF
565
566	!
567	!-- Skip step (4) of wall_fluxes, because here rifs is already
568	!-- available from (1)
569	!
570	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
571
572	IF ( rifs >= 0.0 ) THEN
573
574	!
575	!-- Stable stratification (and neutral)
576	wall_flux(k,j,i) = kappa * vel_zp / &
577	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
578	ELSE
579
580	!
581	!-- Unstable stratification
582	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
583	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
584
585	wall_flux(k,j,i) = kappa * vel_zp / ( &
586	LOG( zp / z0(j,i) ) - &
587	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
588	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
589	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
590	)
591	ENDIF
592	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
593
594	ENDDO
595
596	ENDIF
597
598	ENDDO
599	ENDDO
600
601	END SUBROUTINE wall_fluxes_e
602
603
604	!------------------------------------------------------------------------------!
605	! Call for all grid points - accelerator version
606	!------------------------------------------------------------------------------!
607	SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
608
609	!------------------------------------------------------------------------------!
610	! Description:
611	! ------------
612	! Calculates momentum fluxes at vertical walls for routine production_e
613	! assuming Monin-Obukhov similarity.
614	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
615	!------------------------------------------------------------------------------!
616
617	USE arrays_3d
618	USE control_parameters
619	USE grid_variables
620	USE indices
621	USE statistics
622
623	IMPLICIT NONE
624
625	INTEGER :: i, j, k, kk, max_outer, min_inner, wall_index
626	REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
627	ws, zp
628
629	REAL :: rifs
630
631	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: wall
632	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wall_flux
633
634
635	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
636	wall_flux = 0.0
637	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
638
639	min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
640	max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
641
642	!$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
643	!$acc present( u, v, w, wall, wall_flux, z0 )
644	!$acc loop
645	DO i = nxl, nxr
646	DO j = nys, nyn
647	!$acc loop vector(32)
648	DO k = min_inner, max_outer
649	!
650	!-- All subsequent variables are computed for scalar locations
651	IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &
652	k <= nzb_diff_s_outer(j,i)-2 .AND. wall(j,i) /= 0.0 ) THEN
653	!
654	!-- (1) Compute rifs, u_i, v_i, and ws
655	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
656	kk = nzb_diff_s_inner(j,i)-1
657	ELSE
658	kk = k-1
659	ENDIF
660	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
661	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
662	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
663	)
664
665	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
666	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
667	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
668	!
669	!-- (2) Compute wall-parallel absolute velocity vel_total and
670	!-- interpolate appropriate velocity component vel_zp.
671	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
672	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
673	!
674	!-- (3) Compute wall friction velocity us_wall
675	IF ( rifs >= 0.0 ) THEN
676
677	!
678	!-- Stable stratification (and neutral)
679	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
680	5.0 * rifs * ( zp - z0(j,i) ) / zp &
681	)
682	ELSE
683
684	!
685	!-- Unstable stratification
686	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
687	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
688
689	us_wall = kappa * vel_total / ( &
690	LOG( zp / z0(j,i) ) - &
691	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
692	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
693	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
694	)
695	ENDIF
696
697	!
698	!-- Skip step (4) of wall_fluxes, because here rifs is already
699	!-- available from (1)
700	!
701	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
702
703	IF ( rifs >= 0.0 ) THEN
704
705	!
706	!-- Stable stratification (and neutral)
707	wall_flux(k,j,i) = kappa * vel_zp / &
708	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
709	ELSE
710
711	!
712	!-- Unstable stratification
713	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
714	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
715
716	wall_flux(k,j,i) = kappa * vel_zp / ( &
717	LOG( zp / z0(j,i) ) - &
718	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
719	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
720	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
721	)
722	ENDIF
723	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
724
725	ENDIF
726
727	ENDDO
728	ENDDO
729	ENDDO
730	!$acc end kernels
731
732	END SUBROUTINE wall_fluxes_e_acc
733
734
735	!------------------------------------------------------------------------------!
736	! Call for grid point i,j
737	!------------------------------------------------------------------------------!
738	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
739
740	USE arrays_3d
741	USE control_parameters
742	USE grid_variables
743	USE indices
744	USE statistics
745
746	IMPLICIT NONE
747
748	INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index
749	REAL :: a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
750	ws, zp
751
752	REAL :: rifs
753
754	REAL, DIMENSION(nzb:nzt+1) :: wall_flux
755
756
757	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
758	wall_flux = 0.0
759	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
760
761	!
762	!-- All subsequent variables are computed for scalar locations.
763	DO k = nzb_w, nzt_w
764
765	!
766	!-- (1) Compute rifs, u_i, v_i, and ws
767	IF ( k == nzb_w ) THEN
768	kk = nzb_w
769	ELSE
770	kk = k-1
771	ENDIF
772	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
773	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
774	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
775	)
776
777	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
778	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
779	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
780	!
781	!-- (2) Compute wall-parallel absolute velocity vel_total and
782	!-- interpolate appropriate velocity component vel_zp.
783	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
784	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
785	!
786	!-- (3) Compute wall friction velocity us_wall
787	IF ( rifs >= 0.0 ) THEN
788
789	!
790	!-- Stable stratification (and neutral)
791	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
792	5.0 * rifs * ( zp - z0(j,i) ) / zp &
793	)
794	ELSE
795
796	!
797	!-- Unstable stratification
798	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
799	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
800
801	us_wall = kappa * vel_total / ( &
802	LOG( zp / z0(j,i) ) - &
803	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
804	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
805	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
806	)
807	ENDIF
808
809	!
810	!-- Skip step (4) of wall_fluxes, because here rifs is already
811	!-- available from (1)
812	!
813	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
814	!-- First interpolate the velocity (this is different from
815	!-- subroutine wall_fluxes because fluxes in subroutine
816	!-- wall_fluxes_e are defined at scalar locations).
817	vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &
818	b * ( v(k,j,i) + v(k,j+1,i) ) + &
819	(c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &
820	)
821
822	IF ( rifs >= 0.0 ) THEN
823
824	!
825	!-- Stable stratification (and neutral)
826	wall_flux(k) = kappa * vel_zp / &
827	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
828	ELSE
829
830	!
831	!-- Unstable stratification
832	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
833	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
834
835	wall_flux(k) = kappa * vel_zp / ( &
836	LOG( zp / z0(j,i) ) - &
837	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
838	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
839	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
840	)
841	ENDIF
842	wall_flux(k) = - wall_flux(k) * us_wall
843
844	ENDDO
845
846	END SUBROUTINE wall_fluxes_e_ij
847
848	END MODULE wall_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |