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

Last change on this file since 1952 was 1874, checked in by maronga, 8 years ago

last commit documented

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