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

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

Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added.
Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed.
mask array is set zero for ghost boundaries

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