source: palm/trunk/SOURCE/diffusion_u.f90 @ 4661

Last change on this file since 4661 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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