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

Last change on this file since 4372 was 4360, checked in by suehring, 4 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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