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

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

last commit documented

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