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

Last change on this file since 4671 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

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