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

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

last commit documented

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