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

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

added _mod string to several filenames to meet the naming convection for modules

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