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

Last change on this file since 1314 was 1310, checked in by raasch, 11 years ago

update of GPL copyright

  • 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: wall_fluxes.f90 1310 2014-03-14 08:01:56Z suehring $
27!
28! 1257 2013-11-08 15:18:40Z raasch
29! openacc loop and loop vector clauses removed
30!
31! 1153 2013-05-10 14:33:08Z raasch
32! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
33!
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!
38! 1036 2012-10-22 13:43:42Z raasch
39! code put under GPL (PALM 3.9)
40!
41! 1015 2012-09-27 09:23:24Z raasch
42! accelerator version (*_acc) added
43!
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
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
51! consistency with subroutine wall_fluxes
52! Change: Modification of the integrated version of the profile function for
53! momentum for unstable stratification
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).
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.
65!------------------------------------------------------------------------------!
66    PRIVATE
67    PUBLIC wall_fluxes, wall_fluxes_acc, wall_fluxes_e, wall_fluxes_e_acc
68   
69    INTERFACE wall_fluxes
70       MODULE PROCEDURE wall_fluxes
71       MODULE PROCEDURE wall_fluxes_ij
72    END INTERFACE wall_fluxes
73   
74    INTERFACE wall_fluxes_acc
75       MODULE PROCEDURE wall_fluxes_acc
76    END INTERFACE wall_fluxes_acc
77
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 
83    INTERFACE wall_fluxes_e_acc
84       MODULE PROCEDURE wall_fluxes_e_acc
85    END INTERFACE wall_fluxes_e_acc
86
87 CONTAINS
88
89!------------------------------------------------------------------------------!
90! Call for all grid points
91!------------------------------------------------------------------------------!
92    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
93                            nzb_uvw_outer, wall )
94
95       USE arrays_3d
96       USE control_parameters
97       USE grid_variables
98       USE indices
99       USE statistics
100
101       IMPLICIT NONE
102
103       INTEGER ::  i, j, k, wall_index
104
105       INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  nzb_uvw_inner, &
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
109
110       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
111       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
112
113
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
118       DO  i = nxl, nxr
119          DO  j = nys, nyn
120
121             IF ( wall(j,i) /= 0.0 )  THEN
122!
123!--             All subsequent variables are computed for the respective
124!--             location where the respective flux is defined.
125                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
126
127!
128!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
129                   rifs  = rif_wall(k,j,i,wall_index)
130
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) )
133
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) )
136
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) )
143
144                   pts   = pt_i - hom(k,1,4,0)
145                   wspts = ws * pts
146
147!
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 )
150
151!
152!--                (3) Compute wall friction velocity us_wall
153                   IF ( rifs >= 0.0 )  THEN
154
155!
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
161
162!
163!--                   Unstable stratification
164                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
165                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
166
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                                                    )
173                   ENDIF
174
175!
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 ) )
180
181!
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
190
191!
192!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
193                   IF ( rifs >= 0.0 )  THEN
194
195!
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
203
204!
205!--                   Unstable stratification
206                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
207                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
208
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                                                                            )
216                   ENDIF
217                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
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
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 )
238
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 )
267       !$acc loop independent
268       DO  i = i_left, i_right
269          DO  j = j_south, j_north
270
271             IF ( wall(j,i) /= 0.0 )  THEN
272!
273!--             All subsequent variables are computed for the respective
274!--             location where the respective flux is defined.
275                !$acc loop independent
276                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
277
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
374!                ENDIF
375
376                ENDDO
377
378             ENDIF
379
380          ENDDO
381       ENDDO
382       !$acc end kernels
383
384    END SUBROUTINE wall_fluxes_acc
385
386
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
414!--    the respective flux is defined.
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
454             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
455             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
456
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                                           )
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)
485             wall_flux(k) = kappa *                                          &
486                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
487                            (  LOG( zp / z0(j,i) ) +                         &
488                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
489                            )
490          ELSE
491
492!
493!--          Unstable stratification
494             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
495             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
496
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                                                                   )
504          ENDIF
505          wall_flux(k) = -wall_flux(k) * us_wall
506
507!
508!--       store rifs for next time step
509          rif_wall(k,j,i,wall_index) = rifs
510
511       ENDDO
512
513    END SUBROUTINE wall_fluxes_ij
514
515
516
517!------------------------------------------------------------------------------!
518! Call for all grid points
519!------------------------------------------------------------------------------!
520    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
521
522!------------------------------------------------------------------------------!
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
530       USE arrays_3d
531       USE control_parameters
532       USE grid_variables
533       USE indices
534       USE statistics
535
536       IMPLICIT NONE
537
538       INTEGER ::  i, j, k, kk, wall_index
539       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
540                   ws, zp
541
542       REAL ::  rifs
543
544       REAL, DIMENSION(nysg:nyng,nxlg:nxrg)   ::  wall
545       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
546
547
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 )
551
552       DO  i = nxl, nxr
553          DO  j = nys, nyn
554
555             IF ( wall(j,i) /= 0.0 )  THEN
556!
557!--             All subsequent variables are computed for scalar locations.
558                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
559!
560!--                (1) Compute rifs, u_i, v_i, and ws
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                                 )
570
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) )
574!
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
582
583!
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!
607!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
608
609                   IF ( rifs >= 0.0 )  THEN
610
611!
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
616
617!
618!--                   Unstable stratification
619                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
620                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
621
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                                                          )
628                   ENDIF
629                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
630
631                ENDDO
632
633             ENDIF
634
635          ENDDO
636       ENDDO
637
638    END SUBROUTINE wall_fluxes_e
639
640
641!------------------------------------------------------------------------------!
642! Call for all grid points - accelerator version
643!------------------------------------------------------------------------------!
644    SUBROUTINE wall_fluxes_e_acc( wall_flux, a, b, c1, c2, wall )
645
646!------------------------------------------------------------------------------!
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 )
681       DO  i = i_left, i_right
682          DO  j = j_south, j_north
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!------------------------------------------------------------------------------!
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
784       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
785                   ws, zp
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!
797!--    All subsequent variables are computed for scalar locations.
798       DO  k = nzb_w, nzt_w
799
800!
801!--       (1) Compute rifs, u_i, v_i, and ws
802          IF ( k == nzb_w )  THEN
803             kk = nzb_w
804          ELSE
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
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) )
815!
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
823
824!
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!
848!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
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).
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
867             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
868             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
869
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                                                 )
876          ENDIF
877          wall_flux(k) = - wall_flux(k) * us_wall
878
879       ENDDO
880
881    END SUBROUTINE wall_fluxes_e_ij
882
883 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.