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

Last change on this file since 1755 was 1692, checked in by maronga, 9 years ago

last commit documented

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