source: palm/trunk/SOURCE/diffusion_v.f90 @ 4828

Last change on this file since 4828 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.3 KB
Line 
1!> @file diffusion_v.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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.
8!
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.
12!
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/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 4828 2021-01-05 11:21:41Z Giersch $
26! Update ACC directives for downward facing USM and LSM surfaces
27!
28! 4671 2020-09-09 20:27:58Z pavelkrc
29! Implementation of downward facing USM and LSM surfaces
30!
31! 4583 2020-06-29 12:36:47Z raasch
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!
38! 4329 2019-12-10 15:46:36Z motisi
39! Renamed wall_flags_0 to wall_flags_static_0
40!
41! 4182 2019-08-22 15:20:23Z scharf
42! Corrected "Former revisions" section
43!
44! 3655 2019-01-07 16:51:22Z knoop
45! OpenACC port for SPEC
46!
47! Revision 1.1  1997/09/12 06:24:01  raasch
48! Initial revision
49!
50!
51! Description:
52! ------------
53!> Diffusion term of the v-component
54!--------------------------------------------------------------------------------------------------!
55 MODULE diffusion_v_mod
56
57
58    PRIVATE
59    PUBLIC diffusion_v
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
69!--------------------------------------------------------------------------------------------------!
70! Description:
71! ------------
72!> Call for all grid points
73!--------------------------------------------------------------------------------------------------!
74    SUBROUTINE diffusion_v
75
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,                                                                         &
83           ONLY:  ddx, ddy, ddy2
84
85       USE indices,                                                                                &
86           ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
87
88       USE kinds
89
90       USE surface_mod,                                                                            &
91           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
92
93       IMPLICIT NONE
94
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
102
103       REAL(wp)     ::  flag          !< flag to mask topography grid points
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
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
112
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) &
116       !$ACC PRESENT(wall_flags_total_0, km) &
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)) &
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)) &
122       !$ACC PRESENT(tend)
123       DO  i = nxl, nxr
124          DO  j = nysv, nyn
125!
126!--          Compute horizontal diffusion
127             DO  k = nzb+1, nzt
128
129!
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 ) )
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 ) )
136!
137!--             Interpolate eddy diffusivities on staggered gridpoints
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) )
140
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
155
156             ENDDO
157
158!
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.
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)
168                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
169                ENDDO
170             ENDDO
171!
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)
178                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
179                ENDDO
180             ENDDO
181!
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)
188                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
189                ENDDO
190             ENDDO
191!
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.
196             DO  k = nzb+1, nzt
197!
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 ) )
205!
206!--             Interpolate eddy diffusivities on staggered gridpoints
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) )
209
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
218             ENDDO
219
220!
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.
226             IF ( use_surface_fluxes )  THEN
227!
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)
233
234                   tend(k,j,i) = tend(k,j,i)                                                       &
235                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
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)
243
244                   tend(k,j,i) = tend(k,j,i)                                                       &
245                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
246                ENDDO
247!
248!--             Natural-type surfaces, upward-facing
249                surf_s = surf_lsm_h(0)%start_index(j,i)
250                surf_e = surf_lsm_h(0)%end_index(j,i)
251                DO  m = surf_s, surf_e
252                   k   = surf_lsm_h(0)%k(m)
253
254                   tend(k,j,i) = tend(k,j,i)                                                       &
255                                 + ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
256
257                ENDDO
258!
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!
270!--             Urban-type surfaces, upward-facing
271                surf_s = surf_usm_h(0)%start_index(j,i)
272                surf_e = surf_usm_h(0)%end_index(j,i)
273                DO  m = surf_s, surf_e
274                   k   = surf_usm_h(0)%k(m)
275
276                   tend(k,j,i) = tend(k,j,i)                                                       &
277                                 + ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
278
279                ENDDO
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
291             ENDIF
292!
293!--          Add momentum flux at model top
294             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
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
298
299                   k   = surf_def_h(2)%k(m)
300
301                   tend(k,j,i) = tend(k,j,i)                                                       &
302                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
303                ENDDO
304             ENDIF
305
306          ENDDO
307       ENDDO
308
309    END SUBROUTINE diffusion_v
310
311
312!--------------------------------------------------------------------------------------------------!
313! Description:
314! ------------
315!> Call for grid point i,j
316!--------------------------------------------------------------------------------------------------!
317    SUBROUTINE diffusion_v_ij( i, j )
318
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,                                                                         &
326           ONLY:  ddx, ddy, ddy2
327
328       USE indices,                                                                                &
329           ONLY:  nzb, nzt, wall_flags_total_0
330
331       USE kinds
332
333       USE surface_mod,                                                                            &
334           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
335
336       IMPLICIT NONE
337
338
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
346
347       REAL(wp)     ::  flag          !< flag to mask topography grid points
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
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
355       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
356
357!
358!--    Compute horizontal diffusion
359       DO  k = nzb+1, nzt
360!
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 ) )
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 ) )
367!
368!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
371
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) )                        &
385                                               ) * ddy2 * flag
386       ENDDO
387
388!
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.
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
399          ENDDO
400       ENDDO
401!
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
409          ENDDO
410       ENDDO
411!
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
419          ENDDO
420       ENDDO
421!
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.
425       DO  k = nzb+1, nzt
426!
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 ) )
434!
435!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
438
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
447       ENDDO
448
449!
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.
455       IF ( use_surface_fluxes )  THEN
456!
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)
462
463             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
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)
471
472             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
473          ENDDO
474!
475!--       Natural-type surfaces, upward-facing
476          surf_s = surf_lsm_h(0)%start_index(j,i)
477          surf_e = surf_lsm_h(0)%end_index(j,i)
478          DO  m = surf_s, surf_e
479             k   = surf_lsm_h(0)%k(m)
480
481             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
482
483          ENDDO
484!
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!
495!--       Urban-type surfaces, upward-facing
496          surf_s = surf_usm_h(0)%start_index(j,i)
497          surf_e = surf_usm_h(0)%end_index(j,i)
498          DO  m = surf_s, surf_e
499             k   = surf_usm_h(0)%k(m)
500
501             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
502
503          ENDDO
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
514       ENDIF
515!
516!--    Add momentum flux at model top
517       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
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
521
522             k   = surf_def_h(2)%k(m)
523
524             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
525          ENDDO
526       ENDIF
527
528    END SUBROUTINE diffusion_v_ij
529
530 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.