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

Last change on this file since 1257 was 1257, checked in by raasch, 8 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:


openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:


bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

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