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

Last change on this file since 1037 was 1037, checked in by raasch, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 33.3 KB
Line 
1 MODULE wall_fluxes_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: wall_fluxes.f90 1037 2012-10-22 14:10:22Z raasch $
27!
28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
31! 1015 2012-09-27 09:23:24Z raasch
32! accelerator version (*_acc) added
33!
34! 187 2008-08-06 16:25:09Z letzel
35! Bugfix: Modification of the evaluation of the vertical turbulent momentum
36! fluxes u'w' and v'w (see prandtl_fluxes), this requires the calculation of
37! us_wall (and vel_total, u_i, v_i, ws) also in wall_fluxes_e.
38! Bugfix: change definition of us_wall from 1D to 2D
39! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
40! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
41! consistency with subroutine wall_fluxes
42! Change: Modification of the integrated version of the profile function for
43! momentum for unstable stratification
44!
45! Initial version (2007/03/07)
46!
47! Description:
48! ------------
49! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
50! similarity.
51! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
52! The all-gridpoint version of wall_fluxes_e is not used so far, because
53! it gives slightly different results from the ij-version for some unknown
54! reason.
55!------------------------------------------------------------------------------!
56    PRIVATE
57    PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
58   
59    INTERFACE wall_fluxes
60       MODULE PROCEDURE wall_fluxes
61       MODULE PROCEDURE wall_fluxes_ij
62    END INTERFACE wall_fluxes
63   
64    INTERFACE wall_fluxes_acc
65       MODULE PROCEDURE wall_fluxes_acc
66    END INTERFACE wall_fluxes_acc
67
68    INTERFACE wall_fluxes_e
69       MODULE PROCEDURE wall_fluxes_e
70       MODULE PROCEDURE wall_fluxes_e_ij
71    END INTERFACE wall_fluxes_e
72 
73    INTERFACE wall_fluxes_e_acc
74       MODULE PROCEDURE wall_fluxes_e_acc
75    END INTERFACE wall_fluxes_e_acc
76
77 CONTAINS
78
79!------------------------------------------------------------------------------!
80! Call for all grid points
81!------------------------------------------------------------------------------!
82    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
83                            nzb_uvw_outer, wall )
84
85       USE arrays_3d
86       USE control_parameters
87       USE grid_variables
88       USE indices
89       USE statistics
90
91       IMPLICIT NONE
92
93       INTEGER ::  i, j, k, wall_index
94
95       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
96                                                       nzb_uvw_outer
97       REAL ::  a, b, c1, c2, h1, h2, zp
98       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
99
100       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
101       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
102
103
104       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
105       wall_flux  = 0.0
106       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
107
108       DO  i = nxl, nxr
109          DO  j = nys, nyn
110
111             IF ( wall(j,i) /= 0.0 )  THEN
112!
113!--             All subsequent variables are computed for the respective
114!--             location where the respective flux is defined.
115                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
116
117!
118!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
119                   rifs  = rif_wall(k,j,i,wall_index)
120
121                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
122                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
123
124                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
125                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
126
127                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
128                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
129                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
130                                                           )
131                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
132                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
133
134                   pts   = pt_i - hom(k,1,4,0)
135                   wspts = ws * pts
136
137!
138!--                (2) Compute wall-parallel absolute velocity vel_total
139                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
140
141!
142!--                (3) Compute wall friction velocity us_wall
143                   IF ( rifs >= 0.0 )  THEN
144
145!
146!--                   Stable stratification (and neutral)
147                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
148                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
149                                                    )
150                   ELSE
151
152!
153!--                   Unstable stratification
154                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
155                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
156
157                      us_wall = kappa * vel_total / (                          &
158                           LOG( zp / z0(j,i) ) -                               &
159                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
160                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
161                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
162                                                    )
163                   ENDIF
164
165!
166!--                (4) Compute zp/L (corresponds to neutral Richardson flux
167!--                    number rifs)
168                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
169                                                        ( us_wall**3 + 1E-30 ) )
170
171!
172!--                Limit the value range of the Richardson numbers.
173!--                This is necessary for very small velocities (u,w --> 0),
174!--                because the absolute value of rif can then become very
175!--                large, which in consequence would result in very large
176!--                shear stresses and very small momentum fluxes (both are
177!--                generally unrealistic).
178                   IF ( rifs < rif_min )  rifs = rif_min
179                   IF ( rifs > rif_max )  rifs = rif_max
180
181!
182!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
183                   IF ( rifs >= 0.0 )  THEN
184
185!
186!--                   Stable stratification (and neutral)
187                      wall_flux(k,j,i) = kappa *                               &
188                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
189                              (  LOG( zp / z0(j,i) ) +                         &
190                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
191                              )
192                   ELSE
193
194!
195!--                   Unstable stratification
196                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
197                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
198
199                      wall_flux(k,j,i) = kappa *                               &
200                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
201                           LOG( zp / z0(j,i) ) -                               &
202                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
203                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
204                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
205                                                                            )
206                   ENDIF
207                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
208
209!
210!--                store rifs for next time step
211                   rif_wall(k,j,i,wall_index) = rifs
212
213                ENDDO
214
215             ENDIF
216
217          ENDDO
218       ENDDO
219
220    END SUBROUTINE wall_fluxes
221
222
223!------------------------------------------------------------------------------!
224! Call for all grid points - accelerator version
225!------------------------------------------------------------------------------!
226    SUBROUTINE wall_fluxes_acc( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
227                                nzb_uvw_outer, wall )
228
229       USE arrays_3d
230       USE control_parameters
231       USE grid_variables
232       USE indices
233       USE statistics
234
235       IMPLICIT NONE
236
237       INTEGER ::  i, j, k, max_outer, min_inner, wall_index
238
239       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
240                                                   nzb_uvw_outer
241       REAL ::  a, b, c1, c2, h1, h2, zp
242       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
243
244       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
245       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
246
247
248       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
249       wall_flux  = 0.0
250       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
251
252       min_inner = MINVAL( nzb_uvw_inner(nys:nyn,nxl:nxr) ) + 1
253       max_outer = MINVAL( nzb_uvw_outer(nys:nyn,nxl:nxr) )
254
255       !$acc kernels present( hom, nzb_uvw_inner, nzb_uvw_outer, pt, rif_wall ) &
256       !$acc         present( u, v, w, wall, wall_flux, z0 )
257       !$acc loop
258       DO  i = nxl, nxr
259          DO  j = nys, nyn
260             !$acc loop vector( 32 )
261             DO  k = min_inner, max_outer
262!
263!--             All subsequent variables are computed for the respective
264!--             location where the respective flux is defined.
265                IF ( k >= nzb_uvw_inner(j,i)+1  .AND. &
266                     k <= nzb_uvw_outer(j,i)    .AND.  wall(j,i) /= 0.0 )  THEN
267!
268!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
269                   rifs  = rif_wall(k,j,i,wall_index)
270
271                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
272                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
273
274                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
275                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
276
277                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
278                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
279                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
280                                                           )
281                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
282                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
283
284                   pts   = pt_i - hom(k,1,4,0)
285                   wspts = ws * pts
286
287!
288!--                (2) Compute wall-parallel absolute velocity vel_total
289                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
290
291!
292!--                (3) Compute wall friction velocity us_wall
293                   IF ( rifs >= 0.0 )  THEN
294
295!
296!--                   Stable stratification (and neutral)
297                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
298                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
299                                                    )
300                   ELSE
301
302!
303!--                   Unstable stratification
304                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
305                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
306
307                      us_wall = kappa * vel_total / (                          &
308                           LOG( zp / z0(j,i) ) -                               &
309                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
310                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
311                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
312                                                    )
313                   ENDIF
314
315!
316!--                (4) Compute zp/L (corresponds to neutral Richardson flux
317!--                    number rifs)
318                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
319                                                        ( us_wall**3 + 1E-30 ) )
320
321!
322!--                Limit the value range of the Richardson numbers.
323!--                This is necessary for very small velocities (u,w --> 0),
324!--                because the absolute value of rif can then become very
325!--                large, which in consequence would result in very large
326!--                shear stresses and very small momentum fluxes (both are
327!--                generally unrealistic).
328                   IF ( rifs < rif_min )  rifs = rif_min
329                   IF ( rifs > rif_max )  rifs = rif_max
330
331!
332!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
333                   IF ( rifs >= 0.0 )  THEN
334
335!
336!--                   Stable stratification (and neutral)
337                      wall_flux(k,j,i) = kappa *                               &
338                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
339                              (  LOG( zp / z0(j,i) ) +                         &
340                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
341                              )
342                   ELSE
343
344!
345!--                   Unstable stratification
346                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
347                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
348
349                      wall_flux(k,j,i) = kappa *                               &
350                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
351                           LOG( zp / z0(j,i) ) -                               &
352                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
353                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
354                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
355                                                                            )
356                   ENDIF
357                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
358
359!
360!--                store rifs for next time step
361                   rif_wall(k,j,i,wall_index) = rifs
362
363                ENDIF
364
365             ENDDO
366          ENDDO
367       ENDDO
368       !$acc end kernels
369
370    END SUBROUTINE wall_fluxes_acc
371
372
373!------------------------------------------------------------------------------!
374! Call for all grid point i,j
375!------------------------------------------------------------------------------!
376    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
377
378       USE arrays_3d
379       USE control_parameters
380       USE grid_variables
381       USE indices
382       USE statistics
383
384       IMPLICIT NONE
385
386       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
387       REAL    ::  a, b, c1, c2, h1, h2, zp
388
389       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
390
391       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
392
393
394       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
395       wall_flux  = 0.0
396       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
397
398!
399!--    All subsequent variables are computed for the respective location where
400!--    the respective flux is defined.
401       DO  k = nzb_w, nzt_w
402
403!
404!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
405          rifs  = rif_wall(k,j,i,wall_index)
406
407          u_i   = a * u(k,j,i) + c1 * 0.25 * &
408                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
409
410          v_i   = b * v(k,j,i) + c2 * 0.25 * &
411                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
412
413          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
414                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
415                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
416                                                  )
417          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
418                          + ( c1 + c2 ) * pt(k+1,j,i) )
419
420          pts   = pt_i - hom(k,1,4,0)
421          wspts = ws * pts
422
423!
424!--       (2) Compute wall-parallel absolute velocity vel_total
425          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
426
427!
428!--       (3) Compute wall friction velocity us_wall
429          IF ( rifs >= 0.0 )  THEN
430
431!
432!--          Stable stratification (and neutral)
433             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
434                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
435                                           )
436          ELSE
437
438!
439!--          Unstable stratification
440             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
441             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
442
443             us_wall = kappa * vel_total / (                          &
444                  LOG( zp / z0(j,i) ) -                               &
445                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
446                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
447                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
448                                           )
449          ENDIF
450
451!
452!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
453!--           rifs)
454          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
455
456!
457!--       Limit the value range of the Richardson numbers.
458!--       This is necessary for very small velocities (u,w --> 0), because
459!--       the absolute value of rif can then become very large, which in
460!--       consequence would result in very large shear stresses and very
461!--       small momentum fluxes (both are generally unrealistic).
462          IF ( rifs < rif_min )  rifs = rif_min
463          IF ( rifs > rif_max )  rifs = rif_max
464
465!
466!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
467          IF ( rifs >= 0.0 )  THEN
468
469!
470!--          Stable stratification (and neutral)
471             wall_flux(k) = kappa *                                          &
472                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
473                            (  LOG( zp / z0(j,i) ) +                         &
474                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
475                            )
476          ELSE
477
478!
479!--          Unstable stratification
480             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
481             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
482
483             wall_flux(k) = kappa *                               &
484                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
485                  LOG( zp / z0(j,i) ) -                               &
486                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
487                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
488                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
489                                                                   )
490          ENDIF
491          wall_flux(k) = -wall_flux(k) * us_wall
492
493!
494!--       store rifs for next time step
495          rif_wall(k,j,i,wall_index) = rifs
496
497       ENDDO
498
499    END SUBROUTINE wall_fluxes_ij
500
501
502
503!------------------------------------------------------------------------------!
504! Call for all grid points
505!------------------------------------------------------------------------------!
506    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
507
508!------------------------------------------------------------------------------!
509! Description:
510! ------------
511! Calculates momentum fluxes at vertical walls for routine production_e
512! assuming Monin-Obukhov similarity.
513! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
514!------------------------------------------------------------------------------!
515
516       USE arrays_3d
517       USE control_parameters
518       USE grid_variables
519       USE indices
520       USE statistics
521
522       IMPLICIT NONE
523
524       INTEGER ::  i, j, k, kk, wall_index
525       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
526                   ws, zp
527
528       REAL ::  rifs
529
530       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
531       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
532
533
534       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
535       wall_flux  = 0.0
536       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
537
538       DO  i = nxl, nxr
539          DO  j = nys, nyn
540
541             IF ( wall(j,i) /= 0.0 )  THEN
542!
543!--             All subsequent variables are computed for scalar locations.
544                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
545!
546!--                (1) Compute rifs, u_i, v_i, and ws
547                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
548                      kk = nzb_diff_s_inner(j,i)-1
549                   ELSE
550                      kk = k-1
551                   ENDIF
552                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
553                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
554                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
555                                 )
556
557                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
558                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
559                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
560!
561!--                (2) Compute wall-parallel absolute velocity vel_total and
562!--                interpolate appropriate velocity component vel_zp.
563                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
564                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
565!
566!--                (3) Compute wall friction velocity us_wall
567                   IF ( rifs >= 0.0 )  THEN
568
569!
570!--                   Stable stratification (and neutral)
571                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
572                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
573                                                    )
574                   ELSE
575
576!
577!--                   Unstable stratification
578                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
579                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
580
581                      us_wall = kappa * vel_total / (                          &
582                           LOG( zp / z0(j,i) ) -                               &
583                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
584                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
585                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
586                                                    )
587                   ENDIF
588
589!
590!--                Skip step (4) of wall_fluxes, because here rifs is already
591!--                available from (1)
592!
593!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
594
595                   IF ( rifs >= 0.0 )  THEN
596
597!
598!--                   Stable stratification (and neutral)
599                      wall_flux(k,j,i) = kappa *  vel_zp / &
600                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
601                   ELSE
602
603!
604!--                   Unstable stratification
605                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
606                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
607
608                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
609                           LOG( zp / z0(j,i) ) -                               &
610                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
611                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
612                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
613                                                          )
614                   ENDIF
615                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
616
617                ENDDO
618
619             ENDIF
620
621          ENDDO
622       ENDDO
623
624    END SUBROUTINE wall_fluxes_e
625
626
627!------------------------------------------------------------------------------!
628! Call for all grid points - accelerator version
629!------------------------------------------------------------------------------!
630    SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
631
632!------------------------------------------------------------------------------!
633! Description:
634! ------------
635! Calculates momentum fluxes at vertical walls for routine production_e
636! assuming Monin-Obukhov similarity.
637! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
638!------------------------------------------------------------------------------!
639
640       USE arrays_3d
641       USE control_parameters
642       USE grid_variables
643       USE indices
644       USE statistics
645
646       IMPLICIT NONE
647
648       INTEGER ::  i, j, k, kk, max_outer, min_inner, wall_index
649       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
650                   ws, zp
651
652       REAL ::  rifs
653
654       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
655       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
656
657
658       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
659       wall_flux  = 0.0
660       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
661
662       min_inner = MINVAL( nzb_diff_s_inner(nys:nyn,nxl:nxr) ) - 1
663       max_outer = MAXVAL( nzb_diff_s_outer(nys:nyn,nxl:nxr) ) - 2
664
665       !$acc kernels present( nzb_diff_s_inner, nzb_diff_s_outer, pt, rif_wall ) &
666       !$acc         present( u, v, w, wall, wall_flux, z0 )
667       !$acc loop
668       DO  i = nxl, nxr
669          DO  j = nys, nyn
670             !$acc loop vector(32)
671             DO  k = min_inner, max_outer
672!
673!--             All subsequent variables are computed for scalar locations
674                IF ( k >= nzb_diff_s_inner(j,i)-1  .AND. &
675                     k <= nzb_diff_s_outer(j,i)-2  .AND.  wall(j,i) /= 0.0 )  THEN
676!
677!--                (1) Compute rifs, u_i, v_i, and ws
678                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
679                      kk = nzb_diff_s_inner(j,i)-1
680                   ELSE
681                      kk = k-1
682                   ENDIF
683                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
684                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
685                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
686                                 )
687
688                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
689                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
690                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
691!
692!--                (2) Compute wall-parallel absolute velocity vel_total and
693!--                interpolate appropriate velocity component vel_zp.
694                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
695                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
696!
697!--                (3) Compute wall friction velocity us_wall
698                   IF ( rifs >= 0.0 )  THEN
699
700!
701!--                   Stable stratification (and neutral)
702                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
703                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
704                                                    )
705                   ELSE
706
707!
708!--                   Unstable stratification
709                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
710                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
711
712                      us_wall = kappa * vel_total / (                          &
713                           LOG( zp / z0(j,i) ) -                               &
714                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
715                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
716                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
717                                                    )
718                   ENDIF
719
720!
721!--                Skip step (4) of wall_fluxes, because here rifs is already
722!--                available from (1)
723!
724!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
725
726                   IF ( rifs >= 0.0 )  THEN
727
728!
729!--                   Stable stratification (and neutral)
730                      wall_flux(k,j,i) = kappa *  vel_zp / &
731                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
732                   ELSE
733
734!
735!--                   Unstable stratification
736                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
737                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
738
739                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
740                           LOG( zp / z0(j,i) ) -                               &
741                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
742                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
743                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
744                                                          )
745                   ENDIF
746                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
747
748                ENDIF
749
750             ENDDO
751          ENDDO
752       ENDDO
753       !$acc end kernels
754
755    END SUBROUTINE wall_fluxes_e_acc
756
757
758!------------------------------------------------------------------------------!
759! Call for grid point i,j
760!------------------------------------------------------------------------------!
761    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
762
763       USE arrays_3d
764       USE control_parameters
765       USE grid_variables
766       USE indices
767       USE statistics
768
769       IMPLICIT NONE
770
771       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
772       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
773                   ws, zp
774
775       REAL ::  rifs
776
777       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
778
779
780       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
781       wall_flux  = 0.0
782       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
783
784!
785!--    All subsequent variables are computed for scalar locations.
786       DO  k = nzb_w, nzt_w
787
788!
789!--       (1) Compute rifs, u_i, v_i, and ws
790          IF ( k == nzb_w )  THEN
791             kk = nzb_w
792          ELSE
793             kk = k-1
794          ENDIF
795          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
796                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
797                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
798                        )
799
800          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
801          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
802          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
803!
804!--       (2) Compute wall-parallel absolute velocity vel_total and
805!--       interpolate appropriate velocity component vel_zp.
806          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
807          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
808!
809!--       (3) Compute wall friction velocity us_wall
810          IF ( rifs >= 0.0 )  THEN
811
812!
813!--          Stable stratification (and neutral)
814             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
815                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
816                                           )
817          ELSE
818
819!
820!--          Unstable stratification
821             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
822             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
823
824             us_wall = kappa * vel_total / (                          &
825                  LOG( zp / z0(j,i) ) -                               &
826                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
827                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
828                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
829                                           )
830          ENDIF
831
832!
833!--       Skip step (4) of wall_fluxes, because here rifs is already
834!--       available from (1)
835!
836!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
837!--       First interpolate the velocity (this is different from
838!--       subroutine wall_fluxes because fluxes in subroutine
839!--       wall_fluxes_e are defined at scalar locations).
840          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
841                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
842                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
843                         )
844
845          IF ( rifs >= 0.0 )  THEN
846
847!
848!--          Stable stratification (and neutral)
849             wall_flux(k) = kappa *  vel_zp / &
850                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
851          ELSE
852
853!
854!--          Unstable stratification
855             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
856             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
857
858             wall_flux(k) = kappa * vel_zp / (                        &
859                  LOG( zp / z0(j,i) ) -                               &
860                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
861                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
862                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
863                                                 )
864          ENDIF
865          wall_flux(k) = - wall_flux(k) * us_wall
866
867       ENDDO
868
869    END SUBROUTINE wall_fluxes_e_ij
870
871 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.