Home

Context Navigation

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

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

Line
1	MODULE wall_fluxes_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: wall_fluxes.f90 1321 2014-03-20 09:40:40Z kanani $
27	!
28	! 1320 2014-03-20 08:40:49Z raasch
29	! ONLY-attribute added to USE-statements,
30	! kind-parameters added to all INTEGER and REAL declaration statements,
31	! kinds are defined in new module kinds,
32	! old module precision_kind is removed,
33	! revision history before 2012 removed,
34	! comment fields (!:) to be used for variable explanations added to
35	! all variable declaration statements
36	!
37	! 1257 2013-11-08 15:18:40Z raasch
38	! openacc loop and loop vector clauses removed
39	!
40	! 1153 2013-05-10 14:33:08Z raasch
41	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
42	!
43	! 1128 2013-04-12 06:19:32Z raasch
44	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
45	! j_north
46	!
47	! 1036 2012-10-22 13:43:42Z raasch
48	! code put under GPL (PALM 3.9)
49	!
50	! 1015 2012-09-27 09:23:24Z raasch
51	! accelerator version (*_acc) added
52	!
53	! Initial version (2007/03/07)
54	!
55	! Description:
56	! ------------
57	! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
58	! similarity.
59	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
60	! The all-gridpoint version of wall_fluxes_e is not used so far, because
61	! it gives slightly different results from the ij-version for some unknown
62	! reason.
63	!------------------------------------------------------------------------------!
64	PRIVATE
65	PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
66
67	INTERFACE wall_fluxes
68	MODULE PROCEDURE wall_fluxes
69	MODULE PROCEDURE wall_fluxes_ij
70	END INTERFACE wall_fluxes
71
72	INTERFACE wall_fluxes_acc
73	MODULE PROCEDURE wall_fluxes_acc
74	END INTERFACE wall_fluxes_acc
75
76	INTERFACE wall_fluxes_e
77	MODULE PROCEDURE wall_fluxes_e
78	MODULE PROCEDURE wall_fluxes_e_ij
79	END INTERFACE wall_fluxes_e
80
81	INTERFACE wall_fluxes_e_acc
82	MODULE PROCEDURE wall_fluxes_e_acc
83	END INTERFACE wall_fluxes_e_acc
84
85	CONTAINS
86
87	!------------------------------------------------------------------------------!
88	! Call for all grid points
89	!------------------------------------------------------------------------------!
90	SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
91	nzb_uvw_outer, wall )
92
93	USE arrays_3d, &
94	ONLY: rif_wall, u, v, w, z0, pt
95
96	USE control_parameters, &
97	ONLY: g, kappa, rif_max, rif_min
98
99	USE grid_variables, &
100	ONLY: dx, dy
101
102	USE indices, &
103	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
104
105	USE kinds
106
107	USE statistics, &
108	ONLY: hom
109
110	IMPLICIT NONE
111
112	INTEGER(iwp) :: i !:
113	INTEGER(iwp) :: j !:
114	INTEGER(iwp) :: k !:
115	INTEGER(iwp) :: wall_index !:
116
117	INTEGER(iwp), &
118	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
119	nzb_uvw_inner !:
120	INTEGER(iwp), &
121	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
122	nzb_uvw_outer !:
123
124	REAL(wp) :: a !:
125	REAL(wp) :: b !:
126	REAL(wp) :: c1 !:
127	REAL(wp) :: c2 !:
128	REAL(wp) :: h1 !:
129	REAL(wp) :: h2 !:
130	REAL(wp) :: zp !:
131	REAL(wp) :: pts !:
132	REAL(wp) :: pt_i !:
133	REAL(wp) :: rifs !:
134	REAL(wp) :: u_i !:
135	REAL(wp) :: v_i !:
136	REAL(wp) :: us_wall !:
137	REAL(wp) :: vel_total !:
138	REAL(wp) :: ws !:
139	REAL(wp) :: wspts !:
140
141	REAL(wp), &
142	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
143	wall !:
144
145	REAL(wp), &
146	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
147	wall_flux !:
148
149
150	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
151	wall_flux = 0.0
152	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
153
154	DO i = nxl, nxr
155	DO j = nys, nyn
156
157	IF ( wall(j,i) /= 0.0 ) THEN
158	!
159	!-- All subsequent variables are computed for the respective
160	!-- location where the respective flux is defined.
161	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
162
163	!
164	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
165	rifs = rif_wall(k,j,i,wall_index)
166
167	u_i = a * u(k,j,i) + c1 * 0.25 * &
168	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
169
170	v_i = b * v(k,j,i) + c2 * 0.25 * &
171	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
172
173	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
174	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
175	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
176	)
177	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + &
178	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
179
180	pts = pt_i - hom(k,1,4,0)
181	wspts = ws * pts
182
183	!
184	!-- (2) Compute wall-parallel absolute velocity vel_total
185	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
186
187	!
188	!-- (3) Compute wall friction velocity us_wall
189	IF ( rifs >= 0.0 ) THEN
190
191	!
192	!-- Stable stratification (and neutral)
193	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
194	5.0 * rifs * ( zp - z0(j,i) ) / zp &
195	)
196	ELSE
197
198	!
199	!-- Unstable stratification
200	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
201	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
202
203	us_wall = kappa * vel_total / ( &
204	LOG( zp / z0(j,i) ) - &
205	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
206	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
207	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
208	)
209	ENDIF
210
211	!
212	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
213	!-- number rifs)
214	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
215	( us_wall**3 + 1E-30 ) )
216
217	!
218	!-- Limit the value range of the Richardson numbers.
219	!-- This is necessary for very small velocities (u,w --> 0),
220	!-- because the absolute value of rif can then become very
221	!-- large, which in consequence would result in very large
222	!-- shear stresses and very small momentum fluxes (both are
223	!-- generally unrealistic).
224	IF ( rifs < rif_min ) rifs = rif_min
225	IF ( rifs > rif_max ) rifs = rif_max
226
227	!
228	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
229	IF ( rifs >= 0.0 ) THEN
230
231	!
232	!-- Stable stratification (and neutral)
233	wall_flux(k,j,i) = kappa * &
234	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
235	( LOG( zp / z0(j,i) ) + &
236	5.0 * rifs * ( zp - z0(j,i) ) / zp &
237	)
238	ELSE
239
240	!
241	!-- Unstable stratification
242	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
243	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
244
245	wall_flux(k,j,i) = kappa * &
246	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
247	LOG( zp / z0(j,i) ) - &
248	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
249	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
250	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
251	)
252	ENDIF
253	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
254
255	!
256	!-- store rifs for next time step
257	rif_wall(k,j,i,wall_index) = rifs
258
259	ENDDO
260
261	ENDIF
262
263	ENDDO
264	ENDDO
265
266	END SUBROUTINE wall_fluxes
267
268
269	!------------------------------------------------------------------------------!
270	! Call for all grid points - accelerator version
271	!------------------------------------------------------------------------------!
272	SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
273	nzb_uvw_outer, wall )
274
275	USE arrays_3d, &
276	ONLY: rif_wall, pt, u, v, w, z0
277
278	USE control_parameters, &
279	ONLY: g, kappa, rif_max, rif_min
280
281	USE grid_variables, &
282	ONLY: dx, dy
283
284	USE indices, &
285	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, &
286	nyn, nyng, nys, nysg, nzb, nzt
287
288	USE kinds
289
290	USE statistics, &
291	ONLY: hom
292
293	IMPLICIT NONE
294
295	INTEGER(iwp) :: i !:
296	INTEGER(iwp) :: j !:
297	INTEGER(iwp) :: k !:
298	INTEGER(iwp) :: max_outer !:
299	INTEGER(iwp) :: min_inner !:
300	INTEGER(iwp) :: wall_index !:
301
302	INTEGER(iwp), &
303	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
304	nzb_uvw_inner !:
305	INTEGER(iwp), &
306	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
307	nzb_uvw_outer !:
308
309	REAL(wp) :: a !:
310	REAL(wp) :: b !:
311	REAL(wp) :: c1 !:
312	REAL(wp) :: c2 !:
313	REAL(wp) :: h1 !:
314	REAL(wp) :: h2 !:
315	REAL(wp) :: zp !:
316	REAL(wp) :: pts !:
317	REAL(wp) :: pt_i !:
318	REAL(wp) :: rifs !:
319	REAL(wp) :: u_i !:
320	REAL(wp) :: v_i !:
321	REAL(wp) :: us_wall !:
322	REAL(wp) :: vel_total !:
323	REAL(wp) :: ws !:
324	REAL(wp) :: wspts !:
325
326	REAL(wp), &
327	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
328	wall !:
329
330	REAL(wp), &
331	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
332	wall_flux !:
333
334
335	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
336	wall_flux = 0.0
337	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
338
339	min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
340	max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
341
342	!$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
343	!$acc present( u, v, w, wall, wall_flux, z0 )
344	!$acc loop independent
345	DO i = i_left, i_right
346	DO j = j_south, j_north
347
348	IF ( wall(j,i) /= 0.0 ) THEN
349	!
350	!-- All subsequent variables are computed for the respective
351	!-- location where the respective flux is defined.
352	!$acc loop independent
353	DO k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
354
355	!
356	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
357	rifs = rif_wall(k,j,i,wall_index)
358
359	u_i = a * u(k,j,i) + c1 * 0.25 * &
360	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
361
362	v_i = b * v(k,j,i) + c2 * 0.25 * &
363	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
364
365	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
366	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
367	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
368	)
369	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + &
370	b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
371
372	pts = pt_i - hom(k,1,4,0)
373	wspts = ws * pts
374
375	!
376	!-- (2) Compute wall-parallel absolute velocity vel_total
377	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
378
379	!
380	!-- (3) Compute wall friction velocity us_wall
381	IF ( rifs >= 0.0 ) THEN
382
383	!
384	!-- Stable stratification (and neutral)
385	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
386	5.0 * rifs * ( zp - z0(j,i) ) / zp &
387	)
388	ELSE
389
390	!
391	!-- Unstable stratification
392	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
393	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
394
395	us_wall = kappa * vel_total / ( &
396	LOG( zp / z0(j,i) ) - &
397	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
398	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
399	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
400	)
401	ENDIF
402
403	!
404	!-- (4) Compute zp/L (corresponds to neutral Richardson flux
405	!-- number rifs)
406	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
407	( us_wall**3 + 1E-30 ) )
408
409	!
410	!-- Limit the value range of the Richardson numbers.
411	!-- This is necessary for very small velocities (u,w --> 0),
412	!-- because the absolute value of rif can then become very
413	!-- large, which in consequence would result in very large
414	!-- shear stresses and very small momentum fluxes (both are
415	!-- generally unrealistic).
416	IF ( rifs < rif_min ) rifs = rif_min
417	IF ( rifs > rif_max ) rifs = rif_max
418
419	!
420	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
421	IF ( rifs >= 0.0 ) THEN
422
423	!
424	!-- Stable stratification (and neutral)
425	wall_flux(k,j,i) = kappa * &
426	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
427	( LOG( zp / z0(j,i) ) + &
428	5.0 * rifs * ( zp - z0(j,i) ) / zp &
429	)
430	ELSE
431
432	!
433	!-- Unstable stratification
434	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
435	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
436
437	wall_flux(k,j,i) = kappa * &
438	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
439	LOG( zp / z0(j,i) ) - &
440	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
441	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
442	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
443	)
444	ENDIF
445	wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
446
447	!
448	!-- store rifs for next time step
449	rif_wall(k,j,i,wall_index) = rifs
450
451	ENDDO
452
453	ENDIF
454
455	ENDDO
456	ENDDO
457	!$acc end kernels
458
459	END SUBROUTINE wall_fluxes_acc
460
461
462	!------------------------------------------------------------------------------!
463	! Call for all grid point i,j
464	!------------------------------------------------------------------------------!
465	SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
466
467	USE arrays_3d, &
468	ONLY: rif_wall, pt, u, v, w, z0
469
470	USE control_parameters, &
471	ONLY: g, kappa, rif_max, rif_min
472
473	USE grid_variables, &
474	ONLY: dx, dy
475
476	USE indices, &
477	ONLY: nzb, nzt
478
479	USE kinds
480
481	USE statistics, &
482	ONLY: hom
483
484	IMPLICIT NONE
485
486	INTEGER(iwp) :: i !:
487	INTEGER(iwp) :: j !:
488	INTEGER(iwp) :: k !:
489	INTEGER(iwp) :: nzb_w !:
490	INTEGER(iwp) :: nzt_w !:
491	INTEGER(iwp) :: wall_index !:
492
493	REAL(wp) :: a !:
494	REAL(wp) :: b !:
495	REAL(wp) :: c1 !:
496	REAL(wp) :: c2 !:
497	REAL(wp) :: h1 !:
498	REAL(wp) :: h2 !:
499	REAL(wp) :: zp !:
500	REAL(wp) :: pts !:
501	REAL(wp) :: pt_i !:
502	REAL(wp) :: rifs !:
503	REAL(wp) :: u_i !:
504	REAL(wp) :: v_i !:
505	REAL(wp) :: us_wall !:
506	REAL(wp) :: vel_total !:
507	REAL(wp) :: ws !:
508	REAL(wp) :: wspts !:
509
510	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !:
511
512
513	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
514	wall_flux = 0.0
515	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
516
517	!
518	!-- All subsequent variables are computed for the respective location where
519	!-- the respective flux is defined.
520	DO k = nzb_w, nzt_w
521
522	!
523	!-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
524	rifs = rif_wall(k,j,i,wall_index)
525
526	u_i = a * u(k,j,i) + c1 * 0.25 * &
527	( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
528
529	v_i = b * v(k,j,i) + c2 * 0.25 * &
530	( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
531
532	ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( &
533	a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
534	+ b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
535	)
536	pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) &
537	+ ( c1 + c2 ) * pt(k+1,j,i) )
538
539	pts = pt_i - hom(k,1,4,0)
540	wspts = ws * pts
541
542	!
543	!-- (2) Compute wall-parallel absolute velocity vel_total
544	vel_total = SQRT( ws*2 + ( a+c1 ) u_i*2 + ( b+c2 ) v_i**2 )
545
546	!
547	!-- (3) Compute wall friction velocity us_wall
548	IF ( rifs >= 0.0 ) THEN
549
550	!
551	!-- Stable stratification (and neutral)
552	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
553	5.0 * rifs * ( zp - z0(j,i) ) / zp &
554	)
555	ELSE
556
557	!
558	!-- Unstable stratification
559	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
560	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
561
562	us_wall = kappa * vel_total / ( &
563	LOG( zp / z0(j,i) ) - &
564	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
565	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
566	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
567	)
568	ENDIF
569
570	!
571	!-- (4) Compute zp/L (corresponds to neutral Richardson flux number
572	!-- rifs)
573	rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
574
575	!
576	!-- Limit the value range of the Richardson numbers.
577	!-- This is necessary for very small velocities (u,w --> 0), because
578	!-- the absolute value of rif can then become very large, which in
579	!-- consequence would result in very large shear stresses and very
580	!-- small momentum fluxes (both are generally unrealistic).
581	IF ( rifs < rif_min ) rifs = rif_min
582	IF ( rifs > rif_max ) rifs = rif_max
583
584	!
585	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
586	IF ( rifs >= 0.0 ) THEN
587
588	!
589	!-- Stable stratification (and neutral)
590	wall_flux(k) = kappa * &
591	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / &
592	( LOG( zp / z0(j,i) ) + &
593	5.0 * rifs * ( zp - z0(j,i) ) / zp &
594	)
595	ELSE
596
597	!
598	!-- Unstable stratification
599	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
600	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
601
602	wall_flux(k) = kappa * &
603	( au(k,j,i) + bv(k,j,i) + (c1+c2)*w(k,j,i) ) / ( &
604	LOG( zp / z0(j,i) ) - &
605	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
606	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
607	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
608	)
609	ENDIF
610	wall_flux(k) = -wall_flux(k) * us_wall
611
612	!
613	!-- store rifs for next time step
614	rif_wall(k,j,i,wall_index) = rifs
615
616	ENDDO
617
618	END SUBROUTINE wall_fluxes_ij
619
620
621
622	!------------------------------------------------------------------------------!
623	! Call for all grid points
624	!------------------------------------------------------------------------------!
625	SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
626
627	!------------------------------------------------------------------------------!
628	! Description:
629	! ------------
630	! Calculates momentum fluxes at vertical walls for routine production_e
631	! assuming Monin-Obukhov similarity.
632	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
633	!------------------------------------------------------------------------------!
634
635	USE arrays_3d, &
636	ONLY: rif_wall, u, v, w, z0
637
638	USE control_parameters, &
639	ONLY: kappa
640
641	USE grid_variables, &
642	ONLY: dx, dy
643
644	USE indices, &
645	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
646	nzb_diff_s_inner, nzb_diff_s_outer, nzt
647
648	USE kinds
649
650	IMPLICIT NONE
651
652	INTEGER(iwp) :: i !:
653	INTEGER(iwp) :: j !:
654	INTEGER(iwp) :: k !:
655	INTEGER(iwp) :: kk !:
656	INTEGER(iwp) :: wall_index !:
657
658	REAL(wp) :: a !:
659	REAL(wp) :: b !:
660	REAL(wp) :: c1 !:
661	REAL(wp) :: c2 !:
662	REAL(wp) :: h1 !:
663	REAL(wp) :: h2 !:
664	REAL(wp) :: u_i !:
665	REAL(wp) :: v_i !:
666	REAL(wp) :: us_wall !:
667	REAL(wp) :: vel_total !:
668	REAL(wp) :: vel_zp !:
669	REAL(wp) :: ws !:
670	REAL(wp) :: zp !:
671	REAL(wp) :: rifs !:
672
673	REAL(wp), &
674	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
675	wall !:
676
677	REAL(wp), &
678	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
679	wall_flux !:
680
681
682	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
683	wall_flux = 0.0
684	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
685
686	DO i = nxl, nxr
687	DO j = nys, nyn
688
689	IF ( wall(j,i) /= 0.0 ) THEN
690	!
691	!-- All subsequent variables are computed for scalar locations.
692	DO k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
693	!
694	!-- (1) Compute rifs, u_i, v_i, and ws
695	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
696	kk = nzb_diff_s_inner(j,i)-1
697	ELSE
698	kk = k-1
699	ENDIF
700	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
701	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
702	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
703	)
704
705	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
706	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
707	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
708	!
709	!-- (2) Compute wall-parallel absolute velocity vel_total and
710	!-- interpolate appropriate velocity component vel_zp.
711	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
712	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
713	!
714	!-- (3) Compute wall friction velocity us_wall
715	IF ( rifs >= 0.0 ) THEN
716
717	!
718	!-- Stable stratification (and neutral)
719	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
720	5.0 * rifs * ( zp - z0(j,i) ) / zp &
721	)
722	ELSE
723
724	!
725	!-- Unstable stratification
726	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
727	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
728
729	us_wall = kappa * vel_total / ( &
730	LOG( zp / z0(j,i) ) - &
731	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
732	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
733	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
734	)
735	ENDIF
736
737	!
738	!-- Skip step (4) of wall_fluxes, because here rifs is already
739	!-- available from (1)
740	!
741	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
742
743	IF ( rifs >= 0.0 ) THEN
744
745	!
746	!-- Stable stratification (and neutral)
747	wall_flux(k,j,i) = kappa * vel_zp / &
748	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
749	ELSE
750
751	!
752	!-- Unstable stratification
753	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
754	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
755
756	wall_flux(k,j,i) = kappa * vel_zp / ( &
757	LOG( zp / z0(j,i) ) - &
758	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
759	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
760	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
761	)
762	ENDIF
763	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
764
765	ENDDO
766
767	ENDIF
768
769	ENDDO
770	ENDDO
771
772	END SUBROUTINE wall_fluxes_e
773
774
775	!------------------------------------------------------------------------------!
776	! Call for all grid points - accelerator version
777	!------------------------------------------------------------------------------!
778	SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
779
780	!------------------------------------------------------------------------------!
781	! Description:
782	! ------------
783	! Calculates momentum fluxes at vertical walls for routine production_e
784	! assuming Monin-Obukhov similarity.
785	! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
786	!------------------------------------------------------------------------------!
787
788	USE arrays_3d, &
789	ONLY: rif_wall, u, v, w, z0
790
791	USE control_parameters, &
792	ONLY: kappa
793
794	USE grid_variables, &
795	ONLY: dx, dy
796
797	USE indices, &
798	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxr, nxrg, &
799	nyn, nyng, nys, nysg, nzb, nzb_diff_s_inner, &
800	nzb_diff_s_outer, nzt
801
802	USE kinds
803
804	IMPLICIT NONE
805
806	INTEGER(iwp) :: i !:
807	INTEGER(iwp) :: j !:
808	INTEGER(iwp) :: k !:
809	INTEGER(iwp) :: kk !:
810	INTEGER(iwp) :: max_outer !:
811	INTEGER(iwp) :: min_inner !:
812	INTEGER(iwp) :: wall_index !:
813
814	REAL(wp) :: a !:
815	REAL(wp) :: b !:
816	REAL(wp) :: c1 !:
817	REAL(wp) :: c2 !:
818	REAL(wp) :: h1 !:
819	REAL(wp) :: h2 !:
820	REAL(wp) :: u_i !:
821	REAL(wp) :: v_i !:
822	REAL(wp) :: us_wall !:
823	REAL(wp) :: vel_total !:
824	REAL(wp) :: vel_zp !:
825	REAL(wp) :: ws !:
826	REAL(wp) :: zp !:
827	REAL(wp) :: rifs !:
828
829	REAL(wp), &
830	DIMENSION(nysg:nyng,nxlg:nxrg) :: &
831	wall !:
832
833	REAL(wp), &
834	DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: &
835	wall_flux !:
836
837
838	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
839	wall_flux = 0.0
840	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
841
842	min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
843	max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
844
845	!$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
846	!$acc present( u, v, w, wall, wall_flux, z0 )
847	DO i = i_left, i_right
848	DO j = j_south, j_north
849	DO k = min_inner, max_outer
850	!
851	!-- All subsequent variables are computed for scalar locations
852	IF ( k >= nzb_diff_s_inner(j,i)-1 .AND. &
853	k <= nzb_diff_s_outer(j,i)-2 .AND. wall(j,i) /= 0.0 ) THEN
854	!
855	!-- (1) Compute rifs, u_i, v_i, and ws
856	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
857	kk = nzb_diff_s_inner(j,i)-1
858	ELSE
859	kk = k-1
860	ENDIF
861	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
862	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
863	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
864	)
865
866	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
867	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
868	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
869	!
870	!-- (2) Compute wall-parallel absolute velocity vel_total and
871	!-- interpolate appropriate velocity component vel_zp.
872	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
873	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
874	!
875	!-- (3) Compute wall friction velocity us_wall
876	IF ( rifs >= 0.0 ) THEN
877
878	!
879	!-- Stable stratification (and neutral)
880	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
881	5.0 * rifs * ( zp - z0(j,i) ) / zp &
882	)
883	ELSE
884
885	!
886	!-- Unstable stratification
887	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
888	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
889
890	us_wall = kappa * vel_total / ( &
891	LOG( zp / z0(j,i) ) - &
892	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
893	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
894	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
895	)
896	ENDIF
897
898	!
899	!-- Skip step (4) of wall_fluxes, because here rifs is already
900	!-- available from (1)
901	!
902	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
903
904	IF ( rifs >= 0.0 ) THEN
905
906	!
907	!-- Stable stratification (and neutral)
908	wall_flux(k,j,i) = kappa * vel_zp / &
909	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
910	ELSE
911
912	!
913	!-- Unstable stratification
914	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
915	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
916
917	wall_flux(k,j,i) = kappa * vel_zp / ( &
918	LOG( zp / z0(j,i) ) - &
919	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
920	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
921	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
922	)
923	ENDIF
924	wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
925
926	ENDIF
927
928	ENDDO
929	ENDDO
930	ENDDO
931	!$acc end kernels
932
933	END SUBROUTINE wall_fluxes_e_acc
934
935
936	!------------------------------------------------------------------------------!
937	! Call for grid point i,j
938	!------------------------------------------------------------------------------!
939	SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
940
941	USE arrays_3d, &
942	ONLY: rif_wall, u, v, w, z0
943
944	USE control_parameters, &
945	ONLY: kappa
946
947	USE grid_variables, &
948	ONLY: dx, dy
949
950	USE indices, &
951	ONLY: nzb, nzt
952
953	USE kinds
954
955	IMPLICIT NONE
956
957	INTEGER(iwp) :: i !:
958	INTEGER(iwp) :: j !:
959	INTEGER(iwp) :: k !:
960	INTEGER(iwp) :: kk !:
961	INTEGER(iwp) :: nzb_w !:
962	INTEGER(iwp) :: nzt_w !:
963	INTEGER(iwp) :: wall_index !:
964
965	REAL(wp) :: a !:
966	REAL(wp) :: b !:
967	REAL(wp) :: c1 !:
968	REAL(wp) :: c2 !:
969	REAL(wp) :: h1 !:
970	REAL(wp) :: h2 !:
971	REAL(wp) :: u_i !:
972	REAL(wp) :: v_i !:
973	REAL(wp) :: us_wall !:
974	REAL(wp) :: vel_total !:
975	REAL(wp) :: vel_zp !:
976	REAL(wp) :: ws !:
977	REAL(wp) :: zp !:
978	REAL(wp) :: rifs !:
979
980	REAL(wp), DIMENSION(nzb:nzt+1) :: wall_flux !:
981
982
983	zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
984	wall_flux = 0.0
985	wall_index = NINT( a+ 2b + 3c1 + 4*c2 )
986
987	!
988	!-- All subsequent variables are computed for scalar locations.
989	DO k = nzb_w, nzt_w
990
991	!
992	!-- (1) Compute rifs, u_i, v_i, and ws
993	IF ( k == nzb_w ) THEN
994	kk = nzb_w
995	ELSE
996	kk = k-1
997	ENDIF
998	rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + &
999	a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + &
1000	c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) &
1001	)
1002
1003	u_i = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
1004	v_i = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
1005	ws = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
1006	!
1007	!-- (2) Compute wall-parallel absolute velocity vel_total and
1008	!-- interpolate appropriate velocity component vel_zp.
1009	vel_total = SQRT( ws*2 + (a+c1) u_i*2 + (b+c2) v_i**2 )
1010	vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
1011	!
1012	!-- (3) Compute wall friction velocity us_wall
1013	IF ( rifs >= 0.0 ) THEN
1014
1015	!
1016	!-- Stable stratification (and neutral)
1017	us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + &
1018	5.0 * rifs * ( zp - z0(j,i) ) / zp &
1019	)
1020	ELSE
1021
1022	!
1023	!-- Unstable stratification
1024	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
1025	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
1026
1027	us_wall = kappa * vel_total / ( &
1028	LOG( zp / z0(j,i) ) - &
1029	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
1030	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
1031	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
1032	)
1033	ENDIF
1034
1035	!
1036	!-- Skip step (4) of wall_fluxes, because here rifs is already
1037	!-- available from (1)
1038	!
1039	!-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
1040	!-- First interpolate the velocity (this is different from
1041	!-- subroutine wall_fluxes because fluxes in subroutine
1042	!-- wall_fluxes_e are defined at scalar locations).
1043	vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + &
1044	b * ( v(k,j,i) + v(k,j+1,i) ) + &
1045	(c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) &
1046	)
1047
1048	IF ( rifs >= 0.0 ) THEN
1049
1050	!
1051	!-- Stable stratification (and neutral)
1052	wall_flux(k) = kappa * vel_zp / &
1053	( LOG( zp/z0(j,i) ) + 5.0rifs ( zp-z0(j,i) ) / zp )
1054	ELSE
1055
1056	!
1057	!-- Unstable stratification
1058	h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
1059	h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
1060
1061	wall_flux(k) = kappa * vel_zp / ( &
1062	LOG( zp / z0(j,i) ) - &
1063	LOG( ( 1.0 + h1 )*2 ( 1.0 + h1**2 ) / ( &
1064	( 1.0 + h2 )*2 ( 1.0 + h2**2 ) ) ) + &
1065	2.0 * ( ATAN( h1 ) - ATAN( h2 ) ) &
1066	)
1067	ENDIF
1068	wall_flux(k) = - wall_flux(k) * us_wall
1069
1070	ENDDO
1071
1072	END SUBROUTINE wall_fluxes_e_ij
1073
1074	END MODULE wall_fluxes_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |