source: palm/trunk/SOURCE/diffusion_u.f90

Last change on this file was 4828, checked in by Giersch, 3 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 25.5 KB
RevLine 
[1873]1!> @file diffusion_u.f90
[4583]2!--------------------------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[4583]5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
[1036]8!
[4583]9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
[1036]12!
[4583]13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
[1036]15!
[4828]16! Copyright 1997-2021 Leibniz Universitaet Hannover
[4583]17!--------------------------------------------------------------------------------------------------!
[1036]18!
[484]19! Current revisions:
[1]20! -----------------
[4671]21!
22!
[1321]23! Former revisions:
24! -----------------
25! $Id: diffusion_u.f90 4828 2021-01-05 11:21:41Z banzhafs $
[4674]26! Update ACC directives for downward facing USM and LSM surfaces
27!
28! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]29! Implementation of downward facing USM and LSM surfaces
30!
31! 4583 2020-06-29 12:36:47Z raasch
[4583]32! file re-formatted to follow the PALM coding standard
33!
34! 4360 2020-01-07 11:25:50Z suehring
35! Introduction of wall_flags_total_0, which currently sets bits based on static topography
36! information used in wall_flags_static_0
37!
[4346]38! 4329 2019-12-10 15:46:36Z motisi
[4329]39! Renamed wall_flags_0 to wall_flags_static_0
[4583]40!
[4329]41! 4182 2019-08-22 15:20:23Z scharf
[4182]42! Corrected "Former revisions" section
[4583]43!
[4182]44! 3655 2019-01-07 16:51:22Z knoop
[3634]45! OpenACC port for SPEC
[2716]46!
[4182]47! Revision 1.1  1997/09/12 06:23:51  raasch
48! Initial revision
49!
50!
[1]51! Description:
52! ------------
[1682]53!> Diffusion term of the u-component
[4583]54!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization and slows down the
55!        speed on NEC about 5-10%
56!--------------------------------------------------------------------------------------------------!
[1682]57 MODULE diffusion_u_mod
[1]58
[4583]59
[1]60    PRIVATE
[2118]61    PUBLIC diffusion_u
[1]62
63    INTERFACE diffusion_u
64       MODULE PROCEDURE diffusion_u
65       MODULE PROCEDURE diffusion_u_ij
66    END INTERFACE diffusion_u
67
68 CONTAINS
69
70
[4583]71!--------------------------------------------------------------------------------------------------!
[1682]72! Description:
73! ------------
74!> Call for all grid points
[4583]75!--------------------------------------------------------------------------------------------------!
[1001]76    SUBROUTINE diffusion_u
[1]77
[4583]78       USE arrays_3d,                                                                              &
79           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
80
81       USE control_parameters,                                                                     &
82           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
83
84       USE grid_variables,                                                                         &
[2232]85           ONLY:  ddx, ddx2, ddy
[4583]86
87       USE indices,                                                                                &
[4346]88           ONLY:  nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
[4583]89
[1320]90       USE kinds
[1]91
[4583]92       USE surface_mod,                                                                            &
93           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]94
[1]95       IMPLICIT NONE
96
[2232]97       INTEGER(iwp) ::  i             !< running index x direction
98       INTEGER(iwp) ::  j             !< running index y direction
99       INTEGER(iwp) ::  k             !< running index z direction
100       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
101       INTEGER(iwp) ::  m             !< running index surface elements
[3547]102       INTEGER(iwp) ::  surf_e        !< end index of surface elements at (j,i)-gridpoint
103       INTEGER(iwp) ::  surf_s        !< start index of surface elements at (j,i)-gridpoint
[1001]104
[2232]105       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]106       REAL(wp)     ::  kmym          !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
107       REAL(wp)     ::  kmyp          !< diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
108       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
109       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]110       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
111       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
112       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
113       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
[1]114
[56]115
[2232]116
[3634]117       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
118       !$ACC PRIVATE(surf_e, surf_s, flag, kmym, kmyp, kmzm, kmzp) &
119       !$ACC PRIVATE(mask_bottom, mask_north, mask_south, mask_top) &
[4346]120       !$ACC PRESENT(wall_flags_total_0, km) &
[3634]121       !$ACC PRESENT(u, v, w) &
122       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
123       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:1)) &
[4674]124       !$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(0:1)) &
125       !$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(0:1)) &
[3634]126       !$ACC PRESENT(tend)
[106]127       DO  i = nxlu, nxr
[1001]128          DO  j = nys, nyn
[1]129!
130!--          Compute horizontal diffusion
[2232]131             DO  k = nzb+1, nzt
[1]132!
[4583]133!--             Predetermine flag to mask topography and wall-bounded grid points.
134!--             It is sufficient to masked only north- and south-facing surfaces, which need special
135!--             treatment for the u-component.
136                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
[4346]137                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
138                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
[2232]139!
[1]140!--             Interpolate eddy diffusivities on staggered gridpoints
[4583]141                kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
142                kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]143
[4583]144                tend(k,j,i) = tend(k,j,i)                                                          &
145                              + 2.0_wp * (                                                         &
146                                        km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                  &
147                                      - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                  &
148                                         ) * ddx2 * flag                                           &
149                              +          ( mask_north * (                                          &
150                                  kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy                       &
151                                + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                       &
152                                                        )                                          &
153                                         - mask_south * (                                          &
154                                  kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy                           &
155                                + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx                           &
156                                                        )                                          &
157                                         ) * ddy  * flag
[1]158             ENDDO
159!
[4583]160!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces.
161!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
162!--          Furtermore, note, no vertical natural surfaces so far.
[2232]163!--          Default-type surfaces
164             DO  l = 0, 1
165                surf_s = surf_def_v(l)%start_index(j,i)
166                surf_e = surf_def_v(l)%end_index(j,i)
167                DO  m = surf_s, surf_e
168                   k           = surf_def_v(l)%k(m)
[4583]169                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
170                ENDDO
[2232]171             ENDDO
172!
173!--          Natural-type surfaces
174             DO  l = 0, 1
175                surf_s = surf_lsm_v(l)%start_index(j,i)
176                surf_e = surf_lsm_v(l)%end_index(j,i)
177                DO  m = surf_s, surf_e
178                   k           = surf_lsm_v(l)%k(m)
[4583]179                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
180                ENDDO
[2232]181             ENDDO
182!
183!--          Urban-type surfaces
184             DO  l = 0, 1
185                surf_s = surf_usm_v(l)%start_index(j,i)
186                surf_e = surf_usm_v(l)%end_index(j,i)
187                DO  m = surf_s, surf_e
188                   k           = surf_usm_v(l)%k(m)
[4583]189                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
190                ENDDO
[2232]191             ENDDO
[51]192
[1]193!
[4583]194!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
195!--          diffusive fluxes are masked (flag 8) within this loop, and added further below, else,
196!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
197!--          given.
[2232]198             DO  k = nzb+1, nzt
[1]199!
[4583]200!--             Determine flags to mask topography below and above. Flag 1 is used to mask
201!--             topography in general, and flag 8 implies information about use_surface_fluxes.
202!--             Flag 9 is used to control momentum flux at model top.
203                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
204                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
205                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
206                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]207!
[1]208!--             Interpolate eddy diffusivities on staggered gridpoints
[4583]209                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
210                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]211
[4583]212                tend(k,j,i) = tend(k,j,i)                                                          &
213                              + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                 &
214                                         + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                       &
215                                         ) * rho_air_zw(k)   * mask_top                            &
216                                - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                 &
217                                         + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                     &
218                                         ) * rho_air_zw(k-1) * mask_bottom                         &
219                                ) * ddzw(k) * drho_air(k) * flag
[1]220             ENDDO
221
222!
[4583]223!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
224!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
225!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
226!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
227!--          values of the momentum flux becomes too large in this case.
228!--          The term containing w(k-1,..) (see above equation) is removed here because the vertical
229!--          velocity is assumed to be zero at the surface.
[1]230             IF ( use_surface_fluxes )  THEN
231!
[2232]232!--             Default-type surfaces, upward-facing
233                surf_s = surf_def_h(0)%start_index(j,i)
234                surf_e = surf_def_h(0)%end_index(j,i)
235                DO  m = surf_s, surf_e
[1]236
[2232]237                   k   = surf_def_h(0)%k(m)
[1]238
[4583]239                   tend(k,j,i) = tend(k,j,i)                                                       &
240                                 + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]241                ENDDO
[102]242!
[2232]243!--             Default-type surfaces, dowward-facing
244                surf_s = surf_def_h(1)%start_index(j,i)
245                surf_e = surf_def_h(1)%end_index(j,i)
246                DO  m = surf_s, surf_e
247
248                   k   = surf_def_h(1)%k(m)
249
[4583]250                   tend(k,j,i) = tend(k,j,i)                                                       &
251                                 + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]252                ENDDO
[102]253!
[2232]254!--             Natural-type surfaces, upward-facing
[4671]255                surf_s = surf_lsm_h(0)%start_index(j,i)
256                surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]257                DO  m = surf_s, surf_e
[102]258
[4671]259                   k   = surf_lsm_h(0)%k(m)
[2232]260
[4583]261                   tend(k,j,i) = tend(k,j,i)                                                       &
[4671]262                                 + ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]263                ENDDO
264!
[4671]265!--             Natural-type surfaces, downward-facing
266                surf_s = surf_lsm_h(1)%start_index(j,i)
267                surf_e = surf_lsm_h(1)%end_index(j,i)
268                DO  m = surf_s, surf_e
269
270                   k   = surf_lsm_h(1)%k(m)
271
272                   tend(k,j,i) = tend(k,j,i)                                                       &
273                                 + ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
274                ENDDO
275!
[2232]276!--             Urban-type surfaces, upward-facing
[4671]277                surf_s = surf_usm_h(0)%start_index(j,i)
278                surf_e = surf_usm_h(0)%end_index(j,i)
[2232]279                DO  m = surf_s, surf_e
280
[4671]281                   k   = surf_usm_h(0)%k(m)
[2232]282
[4583]283                   tend(k,j,i) = tend(k,j,i)                                                       &
[4671]284                                 + ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]285                ENDDO
[4671]286!
287!--             Urban-type surfaces, downward-facing
288                surf_s = surf_usm_h(1)%start_index(j,i)
289                surf_e = surf_usm_h(1)%end_index(j,i)
290                DO  m = surf_s, surf_e
[2232]291
[4671]292                   k   = surf_usm_h(1)%k(m)
293
294                   tend(k,j,i) = tend(k,j,i)                                                       &
295                                 + ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
296                ENDDO
297
[102]298             ENDIF
[2232]299!
300!--          Add momentum flux at model top
[2638]301             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]302                surf_s = surf_def_h(2)%start_index(j,i)
303                surf_e = surf_def_h(2)%end_index(j,i)
304                DO  m = surf_s, surf_e
[102]305
[2232]306                   k   = surf_def_h(2)%k(m)
307
[4583]308                   tend(k,j,i) = tend(k,j,i)                                                       &
309                                 + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]310                ENDDO
311             ENDIF
312
[1]313          ENDDO
314       ENDDO
315
316    END SUBROUTINE diffusion_u
317
318
[4583]319!--------------------------------------------------------------------------------------------------!
[1682]320! Description:
321! ------------
322!> Call for grid point i,j
[4583]323!--------------------------------------------------------------------------------------------------!
[1001]324    SUBROUTINE diffusion_u_ij( i, j )
[1]325
[4583]326       USE arrays_3d,                                                                              &
327           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
328
329       USE control_parameters,                                                                     &
330           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
331
332       USE grid_variables,                                                                         &
[2232]333           ONLY:  ddx, ddx2, ddy
[4583]334
335       USE indices,                                                                                &
[4346]336           ONLY:  nzb, nzt, wall_flags_total_0
[4583]337
[1320]338       USE kinds
[1]339
[4583]340       USE surface_mod,                                                                            &
341           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]342
[1]343       IMPLICIT NONE
344
[2232]345       INTEGER(iwp) ::  i             !< running index x direction
346       INTEGER(iwp) ::  j             !< running index y direction
347       INTEGER(iwp) ::  k             !< running index z direction
348       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
349       INTEGER(iwp) ::  m             !< running index surface elements
350       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
351       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1]352
[2232]353       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]354       REAL(wp)     ::  kmym          !< diffusion coefficient on southward side of the u-gridbox - interpolated onto xu-yv grid
355       REAL(wp)     ::  kmyp          !<diffusion coefficient on northward side of the u-gridbox - interpolated onto xu-yv grid
356       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
357       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]358       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
359       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
360       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
361       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
362!
[1]363!--    Compute horizontal diffusion
[2232]364       DO  k = nzb+1, nzt
[1]365!
[4583]366!--       Predetermine flag to mask topography and wall-bounded grid points.
367!--       It is sufficient to masked only north- and south-facing surfaces, which need special
368!--       treatment for the u-component.
369          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   1 ) )
[4346]370          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 1 ) )
371          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 1 ) )
[2232]372!
[1]373!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]374          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
375          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]376
[4583]377          tend(k,j,i) = tend(k,j,i)                                                                &
378                                + 2.0_wp * (                                                       &
379                                         km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )                 &
380                                       - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )                 &
381                                           ) * ddx2 * flag                                         &
382                                         + (                                                       &
383                          mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy                &
384                                              + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx                &
385                                              )                                                    &
386                        - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy                &
387                                              + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx                &
388                                              )                                                    &
389                                           ) * ddy  * flag
[1]390       ENDDO
391
392!
[4583]393!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1) surfaces. Note, in
394!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
395!--    vertical natural surfaces so far.
[2232]396!--    Default-type surfaces
397       DO  l = 0, 1
398          surf_s = surf_def_v(l)%start_index(j,i)
399          surf_e = surf_def_v(l)%end_index(j,i)
400          DO  m = surf_s, surf_e
401             k           = surf_def_v(l)%k(m)
402             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
[4583]403          ENDDO
[2232]404       ENDDO
[51]405!
[2232]406!--    Natural-type surfaces
407       DO  l = 0, 1
408          surf_s = surf_lsm_v(l)%start_index(j,i)
409          surf_e = surf_lsm_v(l)%end_index(j,i)
410          DO  m = surf_s, surf_e
411             k           = surf_lsm_v(l)%k(m)
412             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
[4583]413          ENDDO
[2232]414       ENDDO
[1]415!
[2232]416!--    Urban-type surfaces
417       DO  l = 0, 1
418          surf_s = surf_usm_v(l)%start_index(j,i)
419          surf_e = surf_usm_v(l)%end_index(j,i)
420          DO  m = surf_s, surf_e
421             k           = surf_usm_v(l)%k(m)
422             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
[4583]423          ENDDO
[2232]424       ENDDO
[1]425!
[4583]426!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
427!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
428!--    approach is applied. Model top is also mask if top-momentum flux is given.
[2232]429       DO  k = nzb+1, nzt
430!
[4583]431!--       Determine flags to mask topography below and above. Flag 1 is used to mask topography in
432!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
433!--       control momentum flux at model top.
434          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
435          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
436                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
437          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[2232]438!
[1]439!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]440          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
441          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]442
[4583]443          tend(k,j,i) = tend(k,j,i)                                                                &
444                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                       &
445                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                             &
446                                   ) * rho_air_zw(k)   * mask_top                                  &
447                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)                       &
448                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx                           &
449                                   ) * rho_air_zw(k-1) * mask_bottom                               &
[2232]450                          ) * ddzw(k) * drho_air(k) * flag
[1]451       ENDDO
452
453!
[4583]454!--    Vertical diffusion at the first surface grid points, if the momentum flux at the bottom is
455!--    given by the Prandtl law or if it is prescribed by the user.
456!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
457!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
458!--    of the momentum flux becomes too large in this case.
[1]459       IF ( use_surface_fluxes )  THEN
460!
[2232]461!--       Default-type surfaces, upward-facing
462          surf_s = surf_def_h(0)%start_index(j,i)
463          surf_e = surf_def_h(0)%end_index(j,i)
464          DO  m = surf_s, surf_e
[1]465
[2232]466             k   = surf_def_h(0)%k(m)
[1]467
[4583]468             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]469          ENDDO
[102]470!
[2232]471!--       Default-type surfaces, dowward-facing (except for model-top fluxes)
472          surf_s = surf_def_h(1)%start_index(j,i)
473          surf_e = surf_def_h(1)%end_index(j,i)
474          DO  m = surf_s, surf_e
475
476             k   = surf_def_h(1)%k(m)
477
[4583]478             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]479          ENDDO
[102]480!
[2232]481!--       Natural-type surfaces, upward-facing
[4671]482          surf_s = surf_lsm_h(0)%start_index(j,i)
483          surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]484          DO  m = surf_s, surf_e
[102]485
[4671]486             k   = surf_lsm_h(0)%k(m)
[2232]487
[4671]488             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]489          ENDDO
490!
[4671]491!--       Natural-type surfaces, downward-facing
492          surf_s = surf_lsm_h(1)%start_index(j,i)
493          surf_e = surf_lsm_h(1)%end_index(j,i)
494          DO  m = surf_s, surf_e
495
496             k   = surf_lsm_h(1)%k(m)
497
498             tend(k,j,i) = tend(k,j,i) + ( - surf_lsm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
499          ENDDO
500!
[2232]501!--       Urban-type surfaces, upward-facing
[4671]502          surf_s = surf_usm_h(0)%start_index(j,i)
503          surf_e = surf_usm_h(0)%end_index(j,i)
[2232]504          DO  m = surf_s, surf_e
505
[4671]506             k   = surf_usm_h(0)%k(m)
[2232]507
[4671]508             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%usws(m) ) ) * ddzw(k) * drho_air(k)
[2232]509          ENDDO
[4671]510!
511!--       Urban-type surfaces, downward-facing
512          surf_s = surf_usm_h(1)%start_index(j,i)
513          surf_e = surf_usm_h(1)%end_index(j,i)
514          DO  m = surf_s, surf_e
[2232]515
[4671]516             k   = surf_usm_h(1)%k(m)
517
518             tend(k,j,i) = tend(k,j,i) + ( - surf_usm_h(1)%usws(m) ) * ddzw(k) * drho_air(k)
519          ENDDO
520
[102]521       ENDIF
[2232]522!
523!--    Add momentum flux at model top
[2638]524       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]525          surf_s = surf_def_h(2)%start_index(j,i)
526          surf_e = surf_def_h(2)%end_index(j,i)
527          DO  m = surf_s, surf_e
[102]528
[2232]529             k   = surf_def_h(2)%k(m)
530
[4583]531             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
[2232]532          ENDDO
533       ENDIF
534
535
[1]536    END SUBROUTINE diffusion_u_ij
537
538 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.