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

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

all OpenACC directives and related parts removed from the code

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