source: palm/trunk/SOURCE/diffusion_v.f90

Last change on this file was 4828, checked in by Giersch, 7 months ago

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

  • Property svn:keywords set to Id
File size: 25.3 KB
RevLine 
[1873]1!> @file diffusion_v.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_v.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:24:01  raasch
48! Initial revision
49!
50!
[1]51! Description:
52! ------------
[1682]53!> Diffusion term of the v-component
[4583]54!--------------------------------------------------------------------------------------------------!
[1682]55 MODULE diffusion_v_mod
[1]56
[4583]57
[1]58    PRIVATE
[2118]59    PUBLIC diffusion_v
[1]60
61    INTERFACE diffusion_v
62       MODULE PROCEDURE diffusion_v
63       MODULE PROCEDURE diffusion_v_ij
64    END INTERFACE diffusion_v
65
66 CONTAINS
67
68
[4583]69!--------------------------------------------------------------------------------------------------!
[1682]70! Description:
71! ------------
72!> Call for all grid points
[4583]73!--------------------------------------------------------------------------------------------------!
[1001]74    SUBROUTINE diffusion_v
[1]75
[4583]76       USE arrays_3d,                                                                              &
77           ONLY:  ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
78
79       USE control_parameters,                                                                     &
80           ONLY:  constant_top_momentumflux, use_surface_fluxes,  use_top_fluxes
81
82       USE grid_variables,                                                                         &
[2232]83           ONLY:  ddx, ddy, ddy2
[4583]84
85       USE indices,                                                                                &
[4346]86           ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
[4583]87
[1320]88       USE kinds
[1]89
[4583]90       USE surface_mod,                                                                            &
91           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]92
[1]93       IMPLICIT NONE
94
[2232]95       INTEGER(iwp) ::  i             !< running index x direction
96       INTEGER(iwp) ::  j             !< running index y direction
97       INTEGER(iwp) ::  k             !< running index z direction
98       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
99       INTEGER(iwp) ::  m             !< running index surface elements
100       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
101       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1001]102
[2232]103       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]104       REAL(wp)     ::  kmxm          !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
105       REAL(wp)     ::  kmxp          !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
106       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
107       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
[4583]108       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
109       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
110       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
111       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
[1]112
[3634]113       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
114       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmzm, kmzp) &
115       !$ACC PRIVATE(mask_bottom, mask_east, mask_west, mask_top) &
[4346]116       !$ACC PRESENT(wall_flags_total_0, km) &
[3634]117       !$ACC PRESENT(u, v, w) &
118       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
119       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(2:3)) &
[4674]120       !$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(2:3)) &
121       !$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(2:3)) &
[3634]122       !$ACC PRESENT(tend)
[1]123       DO  i = nxl, nxr
[106]124          DO  j = nysv, nyn
[1]125!
126!--          Compute horizontal diffusion
[2232]127             DO  k = nzb+1, nzt
128
[1]129!
[4583]130!--             Predetermine flag to mask topography and wall-bounded grid points.
131!--             It is sufficient to masked only east- and west-facing surfaces, which need special
132!--             treatment for the v-component.
133                flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
[4346]134                mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
135                mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
[2232]136!
[1]137!--             Interpolate eddy diffusivities on staggered gridpoints
[2232]138                kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
139                kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]140
[4583]141                tend(k,j,i) = tend(k,j,i) +          (                                             &
142                                mask_east * kmxp * (                                               &
143                                       ( v(k,j,i+1) - v(k,j,i)     ) * ddx                         &
144                                     + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                         &
145                                                   )                                               &
146                              - mask_west * kmxm * (                                               &
147                                       ( v(k,j,i) - v(k,j,i-1) ) * ddx                             &
148                                     + ( u(k,j,i) - u(k,j-1,i) ) * ddy                             &
149                                                   )                                               &
150                                                     ) * ddx  * flag                               &
151                                          + 2.0_wp * (                                             &
152                                        km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                  &
153                                      - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                  &
154                                                     ) * ddy2 * flag
[2232]155
[1]156             ENDDO
157
158!
[4583]159!--          Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note,
160!--          in the the flat case, loops won't be entered as start_index > end_index. Furtermore,
161!--          note, no vertical natural surfaces so far.
[2232]162!--          Default-type surfaces
163             DO  l = 2, 3
164                surf_s = surf_def_v(l)%start_index(j,i)
165                surf_e = surf_def_v(l)%end_index(j,i)
166                DO  m = surf_s, surf_e
167                   k           = surf_def_v(l)%k(m)
[4583]168                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
169                ENDDO
[2232]170             ENDDO
[1]171!
[2232]172!--          Natural-type surfaces
173             DO  l = 2, 3
174                surf_s = surf_lsm_v(l)%start_index(j,i)
175                surf_e = surf_lsm_v(l)%end_index(j,i)
176                DO  m = surf_s, surf_e
177                   k           = surf_lsm_v(l)%k(m)
[4583]178                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
179                ENDDO
[2232]180             ENDDO
[1]181!
[2232]182!--          Urban-type surfaces
183             DO  l = 2, 3
184                surf_s = surf_usm_v(l)%start_index(j,i)
185                surf_e = surf_usm_v(l)%end_index(j,i)
186                DO  m = surf_s, surf_e
187                   k           = surf_usm_v(l)%k(m)
[4583]188                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
189                ENDDO
[2232]190             ENDDO
191!
[4583]192!--          Compute vertical diffusion. In case of simulating a surface layer, respective grid
193!--          diffusive fluxes are masked (flag 10) within this loop, and added further below, else,
194!--          simple gradient approach is applied. Model top is also mask if top-momentum flux is
195!--          given.
[2232]196             DO  k = nzb+1, nzt
197!
[4583]198!--             Determine flags to mask topography below and above. Flag 2 is used to mask
199!--             topography in general, while flag 8 implies also information about
200!--             use_surface_fluxes. Flag 9 is used to control momentum flux at model top.
201                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
202                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
203                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
204                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2232]205!
[1]206!--             Interpolate eddy diffusivities on staggered gridpoints
[4583]207                kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
208                kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]209
[4583]210                tend(k,j,i) = tend(k,j,i)                                                          &
211                              & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                 &
212                              &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                       &
213                              &            ) * rho_air_zw(k)   * mask_top                          &
214                              &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)               &
215                              &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                   &
216                              &            ) * rho_air_zw(k-1) * mask_bottom                       &
217                              &   ) * ddzw(k) * drho_air(k) * flag
[1]218             ENDDO
219
220!
[4583]221!--          Vertical diffusion at the first grid point above the surface, if the momentum flux at
222!--          the bottom is given by the Prandtl law or if it is prescribed by the user.
223!--          Difference quotient of the momentum flux is not formed over half of the grid spacing
224!--          (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
225!--          values of the momentum flux becomes too large in this case.
[1]226             IF ( use_surface_fluxes )  THEN
227!
[2232]228!--             Default-type surfaces, upward-facing
229                surf_s = surf_def_h(0)%start_index(j,i)
230                surf_e = surf_def_h(0)%end_index(j,i)
231                DO  m = surf_s, surf_e
232                   k   = surf_def_h(0)%k(m)
[1]233
[4583]234                   tend(k,j,i) = tend(k,j,i)                                                       &
235                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]236                ENDDO
237!
238!--             Default-type surfaces, dowward-facing
239                surf_s = surf_def_h(1)%start_index(j,i)
240                surf_e = surf_def_h(1)%end_index(j,i)
241                DO  m = surf_s, surf_e
242                   k   = surf_def_h(1)%k(m)
[1]243
[4583]244                   tend(k,j,i) = tend(k,j,i)                                                       &
245                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]246                ENDDO
[102]247!
[2232]248!--             Natural-type surfaces, upward-facing
[4671]249                surf_s = surf_lsm_h(0)%start_index(j,i)
250                surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]251                DO  m = surf_s, surf_e
[4671]252                   k   = surf_lsm_h(0)%k(m)
[2232]253
[4583]254                   tend(k,j,i) = tend(k,j,i)                                                       &
[4671]255                                 + ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]256
257                ENDDO
[102]258!
[4671]259!--             Natural-type surfaces, downward-facing
260                surf_s = surf_lsm_h(1)%start_index(j,i)
261                surf_e = surf_lsm_h(1)%end_index(j,i)
262                DO  m = surf_s, surf_e
263                   k   = surf_lsm_h(1)%k(m)
264
265                   tend(k,j,i) = tend(k,j,i)                                                       &
266                                 +  ( - surf_lsm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
267
268                ENDDO
269!
[2232]270!--             Urban-type surfaces, upward-facing
[4671]271                surf_s = surf_usm_h(0)%start_index(j,i)
272                surf_e = surf_usm_h(0)%end_index(j,i)
[2232]273                DO  m = surf_s, surf_e
[4671]274                   k   = surf_usm_h(0)%k(m)
[102]275
[4583]276                   tend(k,j,i) = tend(k,j,i)                                                       &
[4671]277                                 + ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]278
279                ENDDO
[4671]280!
281!--             Urban-type surfaces, downward-facing
282                surf_s = surf_usm_h(1)%start_index(j,i)
283                surf_e = surf_usm_h(1)%end_index(j,i)
284                DO  m = surf_s, surf_e
285                   k   = surf_usm_h(1)%k(m)
286
287                   tend(k,j,i) = tend(k,j,i)                                                       &
288                                 + ( - surf_usm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
289
290                ENDDO
[102]291             ENDIF
[2232]292!
293!--          Add momentum flux at model top
[2638]294             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]295                surf_s = surf_def_h(2)%start_index(j,i)
296                surf_e = surf_def_h(2)%end_index(j,i)
297                DO  m = surf_s, surf_e
[102]298
[2232]299                   k   = surf_def_h(2)%k(m)
300
[4583]301                   tend(k,j,i) = tend(k,j,i)                                                       &
302                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]303                ENDDO
304             ENDIF
305
[1]306          ENDDO
307       ENDDO
308
309    END SUBROUTINE diffusion_v
310
311
[4583]312!--------------------------------------------------------------------------------------------------!
[1682]313! Description:
314! ------------
315!> Call for grid point i,j
[4583]316!--------------------------------------------------------------------------------------------------!
[1001]317    SUBROUTINE diffusion_v_ij( i, j )
[1]318
[4583]319       USE arrays_3d,                                                                              &
320           ONLY:  ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
321
322       USE control_parameters,                                                                     &
323           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
324
325       USE grid_variables,                                                                         &
[2232]326           ONLY:  ddx, ddy, ddy2
[4583]327
328       USE indices,                                                                                &
[4346]329           ONLY:  nzb, nzt, wall_flags_total_0
[4583]330
[1320]331       USE kinds
[1]332
[4583]333       USE surface_mod,                                                                            &
334           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]335
[1]336       IMPLICIT NONE
337
338
[2232]339       INTEGER(iwp) ::  i             !< running index x direction
340       INTEGER(iwp) ::  j             !< running index y direction
341       INTEGER(iwp) ::  k             !< running index z direction
342       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
343       INTEGER(iwp) ::  m             !< running index surface elements
344       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
345       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1001]346
[2232]347       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]348       REAL(wp)     ::  kmxm          !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
349       REAL(wp)     ::  kmxp          !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
350       REAL(wp)     ::  kmzm          !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
351       REAL(wp)     ::  kmzp          !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]352       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface
353       REAL(wp)     ::  mask_east     !< flag to mask vertical surface south of the grid point
354       REAL(wp)     ::  mask_west     !< flag to mask vertical surface north of the grid point
[2232]355       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
356
[1]357!
358!--    Compute horizontal diffusion
[2232]359       DO  k = nzb+1, nzt
[1]360!
[4583]361!--       Predetermine flag to mask topography and wall-bounded grid points.
362!--       It is sufficient to masked only east- and west-facing surfaces, which need special
363!--       treatment for the v-component.
364          flag      = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   2 ) )
[4346]365          mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
366          mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
[2232]367!
[1]368!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]369          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
370          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]371
[4583]372          tend(k,j,i) = tend(k,j,i) +          (                                                   &
373                          mask_east * kmxp * (                                                     &
374                                 ( v(k,j,i+1) - v(k,j,i)     ) * ddx                               &
375                               + ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy                               &
376                                             )                                                     &
377                        - mask_west * kmxm * (                                                     &
378                                 ( v(k,j,i) - v(k,j,i-1) ) * ddx                                   &
379                               + ( u(k,j,i) - u(k,j-1,i) ) * ddy                                   &
380                                             )                                                     &
381                                               ) * ddx  * flag                                     &
382                                    + 2.0_wp * (                                                   &
383                                  km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i)   )                        &
384                                - km(k,j-1,i) * ( v(k,j,i)   - v(k,j-1,i) )                        &
[2232]385                                               ) * ddy2 * flag
[1]386       ENDDO
387
388!
[4583]389!--    Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note, in the
390!--    the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
391!--    vertical natural surfaces so far.
[2232]392!--    Default-type surfaces
393       DO  l = 2, 3
394          surf_s = surf_def_v(l)%start_index(j,i)
395          surf_e = surf_def_v(l)%end_index(j,i)
396          DO  m = surf_s, surf_e
397             k           = surf_def_v(l)%k(m)
398             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
[4583]399          ENDDO
[2232]400       ENDDO
[51]401!
[2232]402!--    Natural-type surfaces
403       DO  l = 2, 3
404          surf_s = surf_lsm_v(l)%start_index(j,i)
405          surf_e = surf_lsm_v(l)%end_index(j,i)
406          DO  m = surf_s, surf_e
407             k           = surf_lsm_v(l)%k(m)
408             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
[4583]409          ENDDO
[2232]410       ENDDO
[1]411!
[2232]412!--    Urban-type surfaces
413       DO  l = 2, 3
414          surf_s = surf_usm_v(l)%start_index(j,i)
415          surf_e = surf_usm_v(l)%end_index(j,i)
416          DO  m = surf_s, surf_e
417             k           = surf_usm_v(l)%k(m)
418             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
[4583]419          ENDDO
[2232]420       ENDDO
[1]421!
[4583]422!--    Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
423!--    fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
424!--    approach is applied. Model top is also mask if top-momentum flux is given.
[2232]425       DO  k = nzb+1, nzt
426!
[4583]427!--       Determine flags to mask topography below and above. Flag 2 is used to mask topography in
428!--       general, while flag 10 implies also information about use_surface_fluxes. Flag 9 is used
429!--       to control momentum flux at model top.
430          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
431          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *         &
432                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
433          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2232]434!
[1]435!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]436          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
437          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]438
[4583]439          tend(k,j,i) = tend(k,j,i)                                                                &
440                        & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)                       &
441                        &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy                             &
442                        &            ) * rho_air_zw(k)   * mask_top                                &
443                        &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)                     &
444                        &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy                         &
445                        &            ) * rho_air_zw(k-1) * mask_bottom                             &
446                        &   ) * ddzw(k) * drho_air(k) * flag
[1]447       ENDDO
448
449!
[4583]450!--    Vertical diffusion at the first grid point above the surface, if the momentum flux at the
451!--    bottom is given by the Prandtl law or if it is prescribed by the user.
452!--    Difference quotient of the momentum flux is not formed over half of the grid spacing
453!--    (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
454!--    of the momentum flux becomes too large in this case.
[1]455       IF ( use_surface_fluxes )  THEN
456!
[2232]457!--       Default-type surfaces, upward-facing
458          surf_s = surf_def_h(0)%start_index(j,i)
459          surf_e = surf_def_h(0)%end_index(j,i)
460          DO  m = surf_s, surf_e
461             k   = surf_def_h(0)%k(m)
[1]462
[4583]463             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]464          ENDDO
465!
466!--       Default-type surfaces, dowward-facing
467          surf_s = surf_def_h(1)%start_index(j,i)
468          surf_e = surf_def_h(1)%end_index(j,i)
469          DO  m = surf_s, surf_e
470             k   = surf_def_h(1)%k(m)
[1]471
[4583]472             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]473          ENDDO
[102]474!
[2232]475!--       Natural-type surfaces, upward-facing
[4671]476          surf_s = surf_lsm_h(0)%start_index(j,i)
477          surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]478          DO  m = surf_s, surf_e
[4671]479             k   = surf_lsm_h(0)%k(m)
[2232]480
[4671]481             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]482
483          ENDDO
[102]484!
[4671]485!--       Natural-type surfaces, downward-facing
486          surf_s = surf_lsm_h(1)%start_index(j,i)
487          surf_e = surf_lsm_h(1)%end_index(j,i)
488          DO  m = surf_s, surf_e
489             k   = surf_lsm_h(1)%k(m)
490
491             tend(k,j,i) = tend(k,j,i) + ( - surf_lsm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
492
493          ENDDO
494!
[2232]495!--       Urban-type surfaces, upward-facing
[4671]496          surf_s = surf_usm_h(0)%start_index(j,i)
497          surf_e = surf_usm_h(0)%end_index(j,i)
[2232]498          DO  m = surf_s, surf_e
[4671]499             k   = surf_usm_h(0)%k(m)
[102]500
[4671]501             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]502
503          ENDDO
[4671]504!
505!--       Urban-type surfaces, downward-facing
506          surf_s = surf_usm_h(1)%start_index(j,i)
507          surf_e = surf_usm_h(1)%end_index(j,i)
508          DO  m = surf_s, surf_e
509             k   = surf_usm_h(1)%k(m)
510
511             tend(k,j,i) = tend(k,j,i) + ( - surf_usm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
512
513          ENDDO
[102]514       ENDIF
[2232]515!
516!--    Add momentum flux at model top
[2638]517       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]518          surf_s = surf_def_h(2)%start_index(j,i)
519          surf_e = surf_def_h(2)%end_index(j,i)
520          DO  m = surf_s, surf_e
[102]521
[2232]522             k   = surf_def_h(2)%k(m)
523
[4583]524             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]525          ENDDO
526       ENDIF
527
[1]528    END SUBROUTINE diffusion_v_ij
529
[1321]530 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.