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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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