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

Last change on this file since 1374 was 1374, checked in by raasch, 7 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

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