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

Last change on this file since 1818 was 1818, checked in by maronga, 6 years ago

last commit documented / copyright update

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