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

Last change on this file since 4656 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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