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

Last change on this file since 4356 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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