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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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