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

Last change on this file since 2201 was 2119, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 28.5 KB
Line 
1!> @file wall_fluxes.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: wall_fluxes.f90 2119 2017-01-17 16:51:50Z suehring $
27!
28! 2118 2017-01-17 16:38:49Z raasch
29! OpenACC versions of subroutines removed
30!
31! 2000 2016-08-20 18:09:15Z knoop
32! Forced header and separation lines into 80 columns
33!
34! 1873 2016-04-18 14:50:06Z maronga
35! Module renamed (removed _mod)
36!
37!
38! 1850 2016-04-08 13:29:27Z maronga
39! Module renamed
40!
41!
42! 1691 2015-10-26 16:17:44Z maronga
43! Renamed rif_min and rif_max with zeta_min and zeta_max, respectively.
44!
45! 1682 2015-10-07 23:56:08Z knoop
46! Code annotations made doxygen readable
47!
48! 1374 2014-04-25 12:55:07Z raasch
49! pt removed from acc-present-list
50!
51! 1353 2014-04-08 15:21:23Z heinze
52! REAL constants provided with KIND-attribute
53!
54! 1320 2014-03-20 08:40:49Z raasch
55! ONLY-attribute added to USE-statements,
56! kind-parameters added to all INTEGER and REAL declaration statements,
57! kinds are defined in new module kinds,
58! old module precision_kind is removed,
59! revision history before 2012 removed,
60! comment fields (!:) to be used for variable explanations added to
61! all variable declaration statements
62!
63! 1257 2013-11-08 15:18:40Z raasch
64! openacc loop and loop vector clauses removed
65!
66! 1153 2013-05-10 14:33:08Z raasch
67! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
68!
69! 1128 2013-04-12 06:19:32Z raasch
70! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
71! j_north
72!
73! 1036 2012-10-22 13:43:42Z raasch
74! code put under GPL (PALM 3.9)
75!
76! 1015 2012-09-27 09:23:24Z raasch
77! accelerator version (*_acc) added
78!
79! Initial version (2007/03/07)
80!
81! Description:
82! ------------
83!> Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
84!> similarity.
85!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
86!> The all-gridpoint version of wall_fluxes_e is not used so far, because
87!> it gives slightly different results from the ij-version for some unknown
88!> reason.
89!>
90!> @todo Rename rif to zeta throughout the routine
91!------------------------------------------------------------------------------!
92 MODULE wall_fluxes_mod
93 
94    PRIVATE
95    PUBLIC wall_fluxes, wall_fluxes_e
96   
97    INTERFACE wall_fluxes
98       MODULE PROCEDURE wall_fluxes
99       MODULE PROCEDURE wall_fluxes_ij
100    END INTERFACE wall_fluxes
101   
102    INTERFACE wall_fluxes_e
103       MODULE PROCEDURE wall_fluxes_e
104       MODULE PROCEDURE wall_fluxes_e_ij
105    END INTERFACE wall_fluxes_e
106 
107 CONTAINS
108
109!------------------------------------------------------------------------------!
110! Description:
111! ------------
112!> Call for all grid points
113!------------------------------------------------------------------------------!
114    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner,            &
115                            nzb_uvw_outer, wall )
116
117       USE arrays_3d,                                                          &
118           ONLY:  rif_wall, u, v, w, z0, pt
119       
120       USE control_parameters,                                                 &
121           ONLY:  g, kappa, zeta_max, zeta_min
122       
123       USE grid_variables,                                                     &
124           ONLY:  dx, dy
125       
126       USE indices,                                                            &
127           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
128       
129       USE kinds
130       
131       USE statistics,                                                         &
132           ONLY:  hom
133
134       IMPLICIT NONE
135
136       INTEGER(iwp) ::  i            !<
137       INTEGER(iwp) ::  j            !<
138       INTEGER(iwp) ::  k            !<
139       INTEGER(iwp) ::  wall_index   !<
140
141       INTEGER(iwp),                                                           &
142          DIMENSION(nysg:nyng,nxlg:nxrg) ::                                    &
143             nzb_uvw_inner   !<
144       INTEGER(iwp),                                                           &
145          DIMENSION(nysg:nyng,nxlg:nxrg) ::                                    &
146             nzb_uvw_outer   !<
147       
148       REAL(wp) ::  a           !<
149       REAL(wp) ::  b           !<
150       REAL(wp) ::  c1          !<
151       REAL(wp) ::  c2          !<
152       REAL(wp) ::  h1          !<
153       REAL(wp) ::  h2          !<
154       REAL(wp) ::  zp          !<
155       REAL(wp) ::  pts         !<
156       REAL(wp) ::  pt_i        !<
157       REAL(wp) ::  rifs        !<
158       REAL(wp) ::  u_i         !<
159       REAL(wp) ::  v_i         !<
160       REAL(wp) ::  us_wall     !<
161       REAL(wp) ::  vel_total   !<
162       REAL(wp) ::  ws          !<
163       REAL(wp) ::  wspts       !<
164
165       REAL(wp),                                                               &
166          DIMENSION(nysg:nyng,nxlg:nxrg) ::                                    &
167             wall   !<
168       
169       REAL(wp),                                                               &
170          DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::                              &
171             wall_flux   !<
172
173
174       zp         = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
175       wall_flux  = 0.0_wp
176       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
177
178       DO  i = nxl, nxr
179          DO  j = nys, nyn
180
181             IF ( wall(j,i) /= 0.0_wp )  THEN
182!
183!--             All subsequent variables are computed for the respective
184!--             location where the respective flux is defined.
185                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
186
187!
188!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
189                   rifs  = rif_wall(k,j,i,wall_index)
190
191                   u_i   = a * u(k,j,i) + c1 * 0.25_wp *                       &
192                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
193
194                   v_i   = b * v(k,j,i) + c2 * 0.25_wp *                       &
195                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
196
197                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * (                &
198                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
199                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
200                                                              )
201                   pt_i  = 0.5_wp * ( pt(k,j,i) + a *  pt(k,j,i-1) +           &
202                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
203
204                   pts   = pt_i - hom(k,1,4,0)
205                   wspts = ws * pts
206
207!
208!--                (2) Compute wall-parallel absolute velocity vel_total
209                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
210
211!
212!--                (3) Compute wall friction velocity us_wall
213                   IF ( rifs >= 0.0_wp )  THEN
214
215!
216!--                   Stable stratification (and neutral)
217                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
218                                         5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
219                                                    )
220                   ELSE
221
222!
223!--                   Unstable stratification
224                      h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
225                      h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
226
227                      us_wall = kappa * vel_total / (                          &
228                           LOG( zp / z0(j,i) ) -                               &
229                           LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (    &
230                                ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +&
231                                2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )           &
232                                                    )
233                   ENDIF
234
235!
236!--                (4) Compute zp/L (corresponds to neutral Richardson flux
237!--                    number rifs)
238                   rifs = -1.0_wp * zp * kappa * g * wspts /                   &
239                          ( pt_i * ( us_wall**3 + 1E-30 ) )
240
241!
242!--                Limit the value range of the Richardson numbers.
243!--                This is necessary for very small velocities (u,w --> 0),
244!--                because the absolute value of rif can then become very
245!--                large, which in consequence would result in very large
246!--                shear stresses and very small momentum fluxes (both are
247!--                generally unrealistic).
248                   IF ( rifs < zeta_min )  rifs = zeta_min
249                   IF ( rifs > zeta_max )  rifs = zeta_max
250
251!
252!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
253                   IF ( rifs >= 0.0_wp )  THEN
254
255!
256!--                   Stable stratification (and neutral)
257                      wall_flux(k,j,i) = kappa *                               &
258                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
259                              (  LOG( zp / z0(j,i) ) +                         &
260                                 5.0_wp * rifs * ( zp - z0(j,i) ) / zp         &
261                              )
262                   ELSE
263
264!
265!--                   Unstable stratification
266                      h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
267                      h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
268
269                      wall_flux(k,j,i) = kappa *                               &
270                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
271                           LOG( zp / z0(j,i) ) -                               &
272                           LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (    &
273                                ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +&
274                                2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )           &
275                                                                            )
276                   ENDIF
277                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
278
279!
280!--                store rifs for next time step
281                   rif_wall(k,j,i,wall_index) = rifs
282
283                ENDDO
284
285             ENDIF
286
287          ENDDO
288       ENDDO
289
290    END SUBROUTINE wall_fluxes
291
292
293!------------------------------------------------------------------------------!
294! Description:
295! ------------
296!> Call for all grid point i,j
297!------------------------------------------------------------------------------!
298    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
299
300       USE arrays_3d,                                                          &
301           ONLY:  rif_wall, pt, u, v, w, z0
302       
303       USE control_parameters,                                                 &
304           ONLY:  g, kappa, zeta_max, zeta_min
305       
306       USE grid_variables,                                                     &
307           ONLY:  dx, dy
308       
309       USE indices,                                                            &
310           ONLY:  nzb, nzt
311       
312       USE kinds
313       
314       USE statistics,                                                         &
315           ONLY:  hom
316
317       IMPLICIT NONE
318
319       INTEGER(iwp) ::  i            !<
320       INTEGER(iwp) ::  j            !<
321       INTEGER(iwp) ::  k            !<
322       INTEGER(iwp) ::  nzb_w        !<
323       INTEGER(iwp) ::  nzt_w        !<
324       INTEGER(iwp) ::  wall_index   !<
325       
326       REAL(wp) ::  a           !<
327       REAL(wp) ::  b           !<
328       REAL(wp) ::  c1          !<
329       REAL(wp) ::  c2          !<
330       REAL(wp) ::  h1          !<
331       REAL(wp) ::  h2          !<
332       REAL(wp) ::  zp          !<
333       REAL(wp) ::  pts         !<
334       REAL(wp) ::  pt_i        !<
335       REAL(wp) ::  rifs        !<
336       REAL(wp) ::  u_i         !<
337       REAL(wp) ::  v_i         !<
338       REAL(wp) ::  us_wall     !<
339       REAL(wp) ::  vel_total   !<
340       REAL(wp) ::  ws          !<
341       REAL(wp) ::  wspts       !<
342
343       REAL(wp), DIMENSION(nzb:nzt+1) ::  wall_flux   !<
344
345
346       zp         = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
347       wall_flux  = 0.0_wp
348       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
349
350!
351!--    All subsequent variables are computed for the respective location where
352!--    the respective flux is defined.
353       DO  k = nzb_w, nzt_w
354
355!
356!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
357          rifs  = rif_wall(k,j,i,wall_index)
358
359          u_i   = a * u(k,j,i) + c1 * 0.25_wp *                                &
360                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
361
362          v_i   = b * v(k,j,i) + c2 * 0.25_wp *                                &
363                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
364
365          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25_wp * (                         &
366                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
367                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
368                                                     )
369          pt_i  = 0.5_wp * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)    &
370                          + ( c1 + c2 ) * pt(k+1,j,i) )
371
372          pts   = pt_i - hom(k,1,4,0)
373          wspts = ws * pts
374
375!
376!--       (2) Compute wall-parallel absolute velocity vel_total
377          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
378
379!
380!--       (3) Compute wall friction velocity us_wall
381          IF ( rifs >= 0.0_wp )  THEN
382
383!
384!--          Stable stratification (and neutral)
385             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
386                                         5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
387                                           )
388          ELSE
389
390!
391!--          Unstable stratification
392             h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
393             h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
394
395             us_wall = kappa * vel_total / (                                   &
396                  LOG( zp / z0(j,i) ) -                                        &
397                  LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (             &
398                       ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +         &
399                       2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )                    &
400                                           )
401          ENDIF
402
403!
404!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
405!--           rifs)
406          rifs = -1.0_wp * zp * kappa * g * wspts /                            &
407                  ( pt_i * (us_wall**3 + 1E-30) )
408
409!
410!--       Limit the value range of the Richardson numbers.
411!--       This is necessary for very small velocities (u,w --> 0), because
412!--       the absolute value of rif can then become very large, which in
413!--       consequence would result in very large shear stresses and very
414!--       small momentum fluxes (both are generally unrealistic).
415          IF ( rifs < zeta_min )  rifs = zeta_min
416          IF ( rifs > zeta_max )  rifs = zeta_max
417
418!
419!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
420          IF ( rifs >= 0.0_wp )  THEN
421
422!
423!--          Stable stratification (and neutral)
424             wall_flux(k) = kappa *                                            &
425                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
426                            (  LOG( zp / z0(j,i) ) +                           &
427                               5.0_wp * rifs * ( zp - z0(j,i) ) / zp           &
428                            )
429          ELSE
430
431!
432!--          Unstable stratification
433             h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
434             h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
435
436             wall_flux(k) = kappa *                                            &
437                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (           &
438                  LOG( zp / z0(j,i) ) -                                        &
439                  LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (             &
440                       ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +         &
441                       2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )                    &
442                                                                   )
443          ENDIF
444          wall_flux(k) = -wall_flux(k) * us_wall
445
446!
447!--       store rifs for next time step
448          rif_wall(k,j,i,wall_index) = rifs
449
450       ENDDO
451
452    END SUBROUTINE wall_fluxes_ij
453
454
455
456!------------------------------------------------------------------------------!
457! Description:
458! ------------
459!> Call for all grid points
460!> Calculates momentum fluxes at vertical walls for routine production_e
461!> assuming Monin-Obukhov similarity.
462!> Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
463!------------------------------------------------------------------------------!
464
465    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
466
467
468       USE arrays_3d,                                                          &
469           ONLY:  rif_wall, u, v, w, z0
470       
471       USE control_parameters,                                                 &
472           ONLY:  kappa
473       
474       USE grid_variables,                                                     &
475           ONLY:  dx, dy
476       
477       USE indices,                                                            &
478           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb,             &
479                  nzb_diff_s_inner, nzb_diff_s_outer, nzt
480       
481       USE kinds
482
483       IMPLICIT NONE
484
485       INTEGER(iwp) ::  i            !<
486       INTEGER(iwp) ::  j            !<
487       INTEGER(iwp) ::  k            !<
488       INTEGER(iwp) ::  kk           !<
489       INTEGER(iwp) ::  wall_index   !<
490       
491       REAL(wp) ::  a           !<
492       REAL(wp) ::  b           !<
493       REAL(wp) ::  c1          !<
494       REAL(wp) ::  c2          !<
495       REAL(wp) ::  h1          !<
496       REAL(wp) ::  h2          !<
497       REAL(wp) ::  u_i         !<
498       REAL(wp) ::  v_i         !<
499       REAL(wp) ::  us_wall     !<
500       REAL(wp) ::  vel_total   !<
501       REAL(wp) ::  vel_zp      !<
502       REAL(wp) ::  ws          !<
503       REAL(wp) ::  zp          !<
504       REAL(wp) ::  rifs        !<
505
506       REAL(wp),                                                               &
507          DIMENSION(nysg:nyng,nxlg:nxrg) ::                                    &
508             wall   !<
509       
510       REAL(wp),                                                               &
511          DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::                              &
512             wall_flux   !<
513
514
515       zp         = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
516       wall_flux  = 0.0_wp
517       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
518
519       DO  i = nxl, nxr
520          DO  j = nys, nyn
521
522             IF ( wall(j,i) /= 0.0_wp )  THEN
523!
524!--             All subsequent variables are computed for scalar locations.
525                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
526!
527!--                (1) Compute rifs, u_i, v_i, and ws
528                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
529                      kk = nzb_diff_s_inner(j,i)-1
530                   ELSE
531                      kk = k-1
532                   ENDIF
533                   rifs  = 0.5_wp * (      rif_wall(k,j,i,wall_index) +        &
534                                      a  * rif_wall(k,j,i+1,1)        +        &
535                                      b  * rif_wall(k,j+1,i,2)        +        &
536                                      c1 * rif_wall(kk,j,i,3)         +        &
537                                      c2 * rif_wall(kk,j,i,4)                  &
538                                    )
539
540                   u_i   = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
541                   v_i   = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
542                   ws    = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
543!
544!--                (2) Compute wall-parallel absolute velocity vel_total and
545!--                interpolate appropriate velocity component vel_zp.
546                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
547                   vel_zp    = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
548!
549!--                (3) Compute wall friction velocity us_wall
550                   IF ( rifs >= 0.0_wp )  THEN
551
552!
553!--                   Stable stratification (and neutral)
554                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
555                                         5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
556                                                    )
557                   ELSE
558
559!
560!--                   Unstable stratification
561                      h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
562                      h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
563
564                      us_wall = kappa * vel_total / (                          &
565                           LOG( zp / z0(j,i) ) -                               &
566                           LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (    &
567                                ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +&
568                                2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )           &
569                                                    )
570                   ENDIF
571
572!
573!--                Skip step (4) of wall_fluxes, because here rifs is already
574!--                available from (1)
575!
576!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
577
578                   IF ( rifs >= 0.0_wp )  THEN
579
580!
581!--                   Stable stratification (and neutral)
582                      wall_flux(k,j,i) = kappa * vel_zp / ( LOG( zp/z0(j,i) ) +&
583                                         5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
584                   ELSE
585
586!
587!--                   Unstable stratification
588                      h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
589                      h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
590
591                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
592                           LOG( zp / z0(j,i) ) -                               &
593                           LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (    &
594                                ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +&
595                                2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )           &
596                                                          )
597                   ENDIF
598                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
599
600                ENDDO
601
602             ENDIF
603
604          ENDDO
605       ENDDO
606
607    END SUBROUTINE wall_fluxes_e
608
609
610!------------------------------------------------------------------------------!
611! Description:
612! ------------
613!> Call for grid point i,j
614!------------------------------------------------------------------------------!
615    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
616
617       USE arrays_3d,                                                          &
618           ONLY:  rif_wall, u, v, w, z0
619       
620       USE control_parameters,                                                 &
621           ONLY:  kappa
622       
623       USE grid_variables,                                                     &
624           ONLY:  dx, dy
625       
626       USE indices,                                                            &
627           ONLY:  nzb, nzt
628       
629       USE kinds
630
631       IMPLICIT NONE
632
633       INTEGER(iwp) ::  i            !<
634       INTEGER(iwp) ::  j            !<
635       INTEGER(iwp) ::  k            !<
636       INTEGER(iwp) ::  kk           !<
637       INTEGER(iwp) ::  nzb_w        !<
638       INTEGER(iwp) ::  nzt_w        !<
639       INTEGER(iwp) ::  wall_index   !<
640       
641       REAL(wp) ::  a           !<
642       REAL(wp) ::  b           !<
643       REAL(wp) ::  c1          !<
644       REAL(wp) ::  c2          !<
645       REAL(wp) ::  h1          !<
646       REAL(wp) ::  h2          !<
647       REAL(wp) ::  u_i         !<
648       REAL(wp) ::  v_i         !<
649       REAL(wp) ::  us_wall     !<
650       REAL(wp) ::  vel_total   !<
651       REAL(wp) ::  vel_zp      !<
652       REAL(wp) ::  ws          !<
653       REAL(wp) ::  zp          !<
654       REAL(wp) ::  rifs        !<
655
656       REAL(wp), DIMENSION(nzb:nzt+1) ::  wall_flux   !<
657
658
659       zp         = 0.5_wp * ( (a+c1) * dy + (b+c2) * dx )
660       wall_flux  = 0.0_wp
661       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
662
663!
664!--    All subsequent variables are computed for scalar locations.
665       DO  k = nzb_w, nzt_w
666
667!
668!--       (1) Compute rifs, u_i, v_i, and ws
669          IF ( k == nzb_w )  THEN
670             kk = nzb_w
671          ELSE
672             kk = k-1
673          ENDIF
674          rifs  = 0.5_wp * (      rif_wall(k,j,i,wall_index) +                 &
675                             a  * rif_wall(k,j,i+1,1)        +                 &
676                             b  * rif_wall(k,j+1,i,2)        +                 &
677                             c1 * rif_wall(kk,j,i,3)         +                 &
678                             c2 * rif_wall(kk,j,i,4)                           &
679                           )
680
681          u_i   = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
682          v_i   = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
683          ws    = 0.5_wp * ( w(k,j,i) + w(k-1,j,i) )
684!
685!--       (2) Compute wall-parallel absolute velocity vel_total and
686!--       interpolate appropriate velocity component vel_zp.
687          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
688          vel_zp    = 0.5_wp * ( a * u_i + b * v_i + (c1+c2) * ws )
689!
690!--       (3) Compute wall friction velocity us_wall
691          IF ( rifs >= 0.0_wp )  THEN
692
693!
694!--          Stable stratification (and neutral)
695             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
696                                         5.0_wp * rifs * ( zp - z0(j,i) ) / zp &
697                                           )
698          ELSE
699
700!
701!--          Unstable stratification
702             h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
703             h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
704
705             us_wall = kappa * vel_total / (                                   &
706                  LOG( zp / z0(j,i) ) -                                        &
707                  LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (             &
708                       ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +         &
709                       2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )                    &
710                                           )
711          ENDIF
712
713!
714!--       Skip step (4) of wall_fluxes, because here rifs is already
715!--       available from (1)
716!
717!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
718!--       First interpolate the velocity (this is different from
719!--       subroutine wall_fluxes because fluxes in subroutine
720!--       wall_fluxes_e are defined at scalar locations).
721          vel_zp = 0.5_wp * (       a * ( u(k,j,i) + u(k,j,i+1) ) +            &
722                                    b * ( v(k,j,i) + v(k,j+1,i) ) +            &
723                              (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )              &
724                            )
725
726          IF ( rifs >= 0.0_wp )  THEN
727
728!
729!--          Stable stratification (and neutral)
730             wall_flux(k) = kappa *  vel_zp /                                  &
731                     ( LOG( zp/z0(j,i) ) + 5.0_wp * rifs * ( zp-z0(j,i) ) / zp )
732          ELSE
733
734!
735!--          Unstable stratification
736             h1 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs ) )
737             h2 = SQRT( SQRT( 1.0_wp - 16.0_wp * rifs * z0(j,i) / zp ) )
738
739             wall_flux(k) = kappa * vel_zp / (                                 &
740                  LOG( zp / z0(j,i) ) -                                        &
741                  LOG( ( 1.0_wp + h1 )**2 * ( 1.0_wp + h1**2 ) / (             &
742                       ( 1.0_wp + h2 )**2 * ( 1.0_wp + h2**2 )   ) ) +         &
743                       2.0_wp * ( ATAN( h1 ) - ATAN( h2 ) )                    &
744                                             )
745          ENDIF
746          wall_flux(k) = - wall_flux(k) * us_wall
747
748       ENDDO
749
750    END SUBROUTINE wall_fluxes_e_ij
751
752 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.