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

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