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

Last change on this file since 4351 was 4346, checked in by motisi, 5 years ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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