Home

Context Navigation

wall_fluxes_mod.f90 @ 1856

Last change on this file since 1856 was 1851, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 43.3 KB

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

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

Context Navigation

source: palm/trunk/SOURCE/wall_fluxes_mod.f90 @ 1856

Download in other formats: