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

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

revised renaming of modules

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