Home

Context Navigation

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

Last change on this file since 1691 was 1691, checked in by maronga, 8 years ago
various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings
Property svn:keywords set to `Id`
File size: 43.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |