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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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