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

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

Update ACC directives for downward facing USM and LSM surfaces

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