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

Last change on this file since 4651 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
RevLine 
[1873]1!> @file diffusion_v.f90
[4583]2!--------------------------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[4583]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.
[1036]8!
[4583]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.
[1036]12!
[4583]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/>.
[1036]15!
[4360]16! Copyright 1997-2020 Leibniz Universitaet Hannover
[4583]17!--------------------------------------------------------------------------------------------------!
[1036]18!
[484]19! Current revisions:
[1]20! -----------------
[1341]21!
[2233]22!
[1321]23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 4583 2020-06-29 12:36:47Z raasch $
[4583]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!
[4346]32! 4329 2019-12-10 15:46:36Z motisi
[4329]33! Renamed wall_flags_0 to wall_flags_static_0
[4583]34!
[4329]35! 4182 2019-08-22 15:20:23Z scharf
[4182]36! Corrected "Former revisions" section
[4583]37!
[4182]38! 3655 2019-01-07 16:51:22Z knoop
[3634]39! OpenACC port for SPEC
[2716]40!
[4182]41! Revision 1.1  1997/09/12 06:24:01  raasch
42! Initial revision
43!
44!
[1]45! Description:
46! ------------
[1682]47!> Diffusion term of the v-component
[4583]48!--------------------------------------------------------------------------------------------------!
[1682]49 MODULE diffusion_v_mod
[1]50
[4583]51
[1]52    PRIVATE
[2118]53    PUBLIC diffusion_v
[1]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
[4583]63!--------------------------------------------------------------------------------------------------!
[1682]64! Description:
65! ------------
66!> Call for all grid points
[4583]67!--------------------------------------------------------------------------------------------------!
[1001]68    SUBROUTINE diffusion_v
[1]69
[4583]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,                                                                         &
[2232]77           ONLY:  ddx, ddy, ddy2
[4583]78
79       USE indices,                                                                                &
[4346]80           ONLY:  nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
[4583]81
[1320]82       USE kinds
[1]83
[4583]84       USE surface_mod,                                                                            &
85           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]86
[1]87       IMPLICIT NONE
88
[2232]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
[1001]96
[2232]97       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]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
[4583]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
[1]106
[3634]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) &
[4346]110       !$ACC PRESENT(wall_flags_total_0, km) &
[3634]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)
[1]117       DO  i = nxl, nxr
[106]118          DO  j = nysv, nyn
[1]119!
120!--          Compute horizontal diffusion
[2232]121             DO  k = nzb+1, nzt
122
[1]123!
[4583]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 ) )
[4346]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 ) )
[2232]130!
[1]131!--             Interpolate eddy diffusivities on staggered gridpoints
[2232]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) )
[1]134
[4583]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
[2232]149
[1]150             ENDDO
151
152!
[4583]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.
[2232]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)
[4583]162                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
163                ENDDO
[2232]164             ENDDO
[1]165!
[2232]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)
[4583]172                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
173                ENDDO
[2232]174             ENDDO
[1]175!
[2232]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)
[4583]182                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
183                ENDDO
[2232]184             ENDDO
185!
[4583]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.
[2232]190             DO  k = nzb+1, nzt
191!
[4583]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 ) )
[2232]199!
[1]200!--             Interpolate eddy diffusivities on staggered gridpoints
[4583]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) )
[1]203
[4583]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
[1]212             ENDDO
213
214!
[4583]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.
[1]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
[4583]228                   tend(k,j,i) = tend(k,j,i)                                                       &
229                                 + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]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)
[1]237
[4583]238                   tend(k,j,i) = tend(k,j,i)                                                       &
239                                 + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]240                ENDDO
[102]241!
[2232]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
[4583]248                   tend(k,j,i) = tend(k,j,i)                                                       &
249                                 + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]250
251                ENDDO
[102]252!
[2232]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)
[102]258
[4583]259                   tend(k,j,i) = tend(k,j,i)                                                       &
260                                 + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]261
262                ENDDO
[102]263             ENDIF
[2232]264!
265!--          Add momentum flux at model top
[2638]266             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]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
[102]270
[2232]271                   k   = surf_def_h(2)%k(m)
272
[4583]273                   tend(k,j,i) = tend(k,j,i)                                                       &
274                                 + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]275                ENDDO
276             ENDIF
277
[1]278          ENDDO
279       ENDDO
280
281    END SUBROUTINE diffusion_v
282
283
[4583]284!--------------------------------------------------------------------------------------------------!
[1682]285! Description:
286! ------------
287!> Call for grid point i,j
[4583]288!--------------------------------------------------------------------------------------------------!
[1001]289    SUBROUTINE diffusion_v_ij( i, j )
[1]290
[4583]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,                                                                         &
[2232]298           ONLY:  ddx, ddy, ddy2
[4583]299
300       USE indices,                                                                                &
[4346]301           ONLY:  nzb, nzt, wall_flags_total_0
[4583]302
[1320]303       USE kinds
[1]304
[4583]305       USE surface_mod,                                                                            &
306           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]307
[1]308       IMPLICIT NONE
309
310
[2232]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
[1001]318
[2232]319       REAL(wp)     ::  flag          !< flag to mask topography grid points
[3547]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
[4583]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
[2232]327       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
328
[1]329!
330!--    Compute horizontal diffusion
[2232]331       DO  k = nzb+1, nzt
[1]332!
[4583]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 ) )
[4346]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 ) )
[2232]339!
[1]340!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]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) )
[1]343
[4583]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) )                        &
[2232]357                                               ) * ddy2 * flag
[1]358       ENDDO
359
360!
[4583]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.
[2232]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
[4583]371          ENDDO
[2232]372       ENDDO
[51]373!
[2232]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
[4583]381          ENDDO
[2232]382       ENDDO
[1]383!
[2232]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
[4583]391          ENDDO
[2232]392       ENDDO
[1]393!
[4583]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.
[2232]397       DO  k = nzb+1, nzt
398!
[4583]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 ) )
[2232]406!
[1]407!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]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) )
[1]410
[4583]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
[1]419       ENDDO
420
421!
[4583]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.
[1]427       IF ( use_surface_fluxes )  THEN
428!
[2232]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)
[1]434
[4583]435             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]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)
[1]443
[4583]444             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]445          ENDDO
[102]446!
[2232]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
[4583]453             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]454
455          ENDDO
[102]456!
[2232]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)
[102]462
[4583]463             tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]464
465          ENDDO
[102]466       ENDIF
[2232]467!
468!--    Add momentum flux at model top
[2638]469       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]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
[102]473
[2232]474             k   = surf_def_h(2)%k(m)
475
[4583]476             tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]477          ENDDO
478       ENDIF
479
[1]480    END SUBROUTINE diffusion_v_ij
481
[1321]482 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.