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

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

last commit documented

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