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

Last change on this file since 1433 was 1375, checked in by raasch, 11 years ago

last commit documented

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