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

Last change on this file since 1279 was 1258, checked in by raasch, 11 years ago

last commit documented

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