source: palm/trunk/SOURCE/wall_fluxes_mod.f90 @ 1858

Last change on this file since 1858 was 1851, checked in by maronga, 9 years ago

last commit documented

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