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

Last change on this file since 4341 was 4329, checked in by motisi, 4 years ago

Renamed wall_flags_0 to wall_flags_static_0

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