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

Last change on this file since 1353 was 1353, checked in by heinze, 10 years ago

REAL constants provided with KIND-attribute

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