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

Last change on this file since 3208 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

  • Property svn:keywords set to Id
File size: 24.8 KB
RevLine 
[1873]1!> @file diffusion_u.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!
[2718]17! Copyright 1997-2018 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_u.f90 2718 2018-01-02 08:49:38Z sward $
[2716]27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
31!
32! 2638 2017-11-23 12:44:23Z raasch
[2638]33! bugfix for constant top momentumflux
34!
35! 2233 2017-05-30 18:08:54Z suehring
[1321]36!
[2233]37! 2232 2017-05-30 17:47:52Z suehring
38! Adjustments to new topography and surface concept
39!
[2119]40! 2118 2017-01-17 16:38:49Z raasch
41! OpenACC version of subroutine removed
42!
[2038]43! 2037 2016-10-26 11:15:40Z knoop
44! Anelastic approximation implemented
45!
[2001]46! 2000 2016-08-20 18:09:15Z knoop
47! Forced header and separation lines into 80 columns
48!
[1874]49! 1873 2016-04-18 14:50:06Z maronga
50! Module renamed (removed _mod)
51!
[1851]52! 1850 2016-04-08 13:29:27Z maronga
53! Module renamed
54!
[1741]55! 1740 2016-01-13 08:19:40Z raasch
56! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
57!
[1692]58! 1691 2015-10-26 16:17:44Z maronga
59! Formatting corrections.
60!
[1683]61! 1682 2015-10-07 23:56:08Z knoop
62! Code annotations made doxygen readable
63!
[1341]64! 1340 2014-03-25 19:45:13Z kanani
65! REAL constants defined as wp-kind
66!
[1321]67! 1320 2014-03-20 08:40:49Z raasch
[1320]68! ONLY-attribute added to USE-statements,
69! kind-parameters added to all INTEGER and REAL declaration statements,
70! kinds are defined in new module kinds,
71! revision history before 2012 removed,
72! comment fields (!:) to be used for variable explanations added to
73! all variable declaration statements
[1321]74!
[1258]75! 1257 2013-11-08 15:18:40Z raasch
76! openacc loop and loop vector clauses removed, declare create moved after
77! the FORTRAN declaration statement
78!
[1132]79! 1128 2013-04-12 06:19:32Z raasch
80! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
81! j_north
82!
[1037]83! 1036 2012-10-22 13:43:42Z raasch
84! code put under GPL (PALM 3.9)
85!
[1017]86! 1015 2012-09-27 09:23:24Z raasch
87! accelerator version (*_acc) added
88!
[1002]89! 1001 2012-09-13 14:08:46Z raasch
90! arrays comunicated by module instead of parameter list
91!
[979]92! 978 2012-08-09 08:28:32Z fricke
93! outflow damping layer removed
94! kmym_x/_y and kmyp_x/_y change to kmym and kmyp
95!
[1]96! Revision 1.1  1997/09/12 06:23:51  raasch
97! Initial revision
98!
99!
100! Description:
101! ------------
[1682]102!> Diffusion term of the u-component
103!> @todo additional damping (needed for non-cyclic bc) causes bad vectorization
104!>       and slows down the speed on NEC about 5-10%
[1]105!------------------------------------------------------------------------------!
[1682]106 MODULE diffusion_u_mod
107 
[1]108
109    PRIVATE
[2118]110    PUBLIC diffusion_u
[1]111
112    INTERFACE diffusion_u
113       MODULE PROCEDURE diffusion_u
114       MODULE PROCEDURE diffusion_u_ij
115    END INTERFACE diffusion_u
116
117 CONTAINS
118
119
120!------------------------------------------------------------------------------!
[1682]121! Description:
122! ------------
123!> Call for all grid points
[1]124!------------------------------------------------------------------------------!
[1001]125    SUBROUTINE diffusion_u
[1]126
[1320]127       USE arrays_3d,                                                          &
[2232]128           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
[1320]129       
130       USE control_parameters,                                                 &
[2232]131           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
[1320]132                  use_top_fluxes
133       
134       USE grid_variables,                                                     &
[2232]135           ONLY:  ddx, ddx2, ddy
[1320]136       
137       USE indices,                                                            &
[2232]138           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzt, wall_flags_0
139     
[1320]140       USE kinds
[1]141
[2232]142       USE surface_mod,                                                        &
143           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
144                   surf_usm_v
145
[1]146       IMPLICIT NONE
147
[2232]148       INTEGER(iwp) ::  i             !< running index x direction
149       INTEGER(iwp) ::  j             !< running index y direction
150       INTEGER(iwp) ::  k             !< running index z direction
151       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
152       INTEGER(iwp) ::  m             !< running index surface elements
153       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
154       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1001]155
[2232]156       REAL(wp)     ::  flag          !< flag to mask topography grid points
157       REAL(wp)     ::  kmym          !<
158       REAL(wp)     ::  kmyp          !<
159       REAL(wp)     ::  kmzm          !<
160       REAL(wp)     ::  kmzp          !<
161       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
162       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
163       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
164       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
[1]165
[56]166
[2232]167
[106]168       DO  i = nxlu, nxr
[1001]169          DO  j = nys, nyn
[1]170!
171!--          Compute horizontal diffusion
[2232]172             DO  k = nzb+1, nzt
[1]173!
[2232]174!--             Predetermine flag to mask topography and wall-bounded grid points.
175!--             It is sufficient to masked only north- and south-facing surfaces, which
176!--             need special treatment for the u-component.
177                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   1 ) ) 
178                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 1 ) )
179                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 1 ) )
180!
[1]181!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]182                kmyp = 0.25_wp *                                               &
[978]183                       ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
[1340]184                kmym = 0.25_wp *                                               &
[978]185                       ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]186
[1320]187                tend(k,j,i) = tend(k,j,i)                                      &
[2232]188                        + 2.0_wp * (                                           &
189                                  km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )    &
190                                - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )    &
191                                   ) * ddx2 * flag                             &
192                        +          ( mask_north * (                            &
193                            kmyp * ( u(k,j+1,i) - u(k,j,i)     ) * ddy         &
194                          + kmyp * ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx         &
195                                                  )                            &
196                                   - mask_south * (                            &
197                            kmym * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
198                          + kmym * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
199                                                  )                            &
200                                   ) * ddy  * flag                             
[1]201             ENDDO
202!
[2232]203!--          Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
204!--          surfaces. Note, in the the flat case, loops won't be entered as
205!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
206!--          so far.           
207!--          Default-type surfaces
208             DO  l = 0, 1
209                surf_s = surf_def_v(l)%start_index(j,i)
210                surf_e = surf_def_v(l)%end_index(j,i)
211                DO  m = surf_s, surf_e
212                   k           = surf_def_v(l)%k(m)
213                   tend(k,j,i) = tend(k,j,i) +                                 &                   
214                                    surf_def_v(l)%mom_flux_uv(m) * ddy
215                ENDDO   
216             ENDDO
217!
218!--          Natural-type surfaces
219             DO  l = 0, 1
220                surf_s = surf_lsm_v(l)%start_index(j,i)
221                surf_e = surf_lsm_v(l)%end_index(j,i)
222                DO  m = surf_s, surf_e
223                   k           = surf_lsm_v(l)%k(m)
224                   tend(k,j,i) = tend(k,j,i) +                                 &                   
225                                    surf_lsm_v(l)%mom_flux_uv(m) * ddy
226                ENDDO   
227             ENDDO
228!
229!--          Urban-type surfaces
230             DO  l = 0, 1
231                surf_s = surf_usm_v(l)%start_index(j,i)
232                surf_e = surf_usm_v(l)%end_index(j,i)
233                DO  m = surf_s, surf_e
234                   k           = surf_usm_v(l)%k(m)
235                   tend(k,j,i) = tend(k,j,i) +                                 &                   
236                                    surf_usm_v(l)%mom_flux_uv(m) * ddy
237                ENDDO   
238             ENDDO
[51]239
[1]240!
[2232]241!--          Compute vertical diffusion. In case of simulating a surface layer,
242!--          respective grid diffusive fluxes are masked (flag 8) within this
243!--          loop, and added further below, else, simple gradient approach is
244!--          applied. Model top is also mask if top-momentum flux is given.
245             DO  k = nzb+1, nzt
[1]246!
[2232]247!--             Determine flags to mask topography below and above. Flag 1 is
248!--             used to mask topography in general, and flag 8 implies
249!--             information about use_surface_fluxes. Flag 9 is used to control
250!--             momentum flux at model top. 
251                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
252                                     BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
253                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
254                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
255                              MERGE( 1.0_wp, 0.0_wp,                           &
256                                     BTEST( wall_flags_0(k+1,j,i), 9 ) ) 
257                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
258                                     BTEST( wall_flags_0(k,j,i), 1 ) ) 
259!
[1]260!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]261                kmzp = 0.25_wp *                                               &
[1]262                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
[1340]263                kmzm = 0.25_wp *                                               &
[1]264                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
265
[1320]266                tend(k,j,i) = tend(k,j,i)                                      &
[2232]267                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
268                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
269                                   ) * rho_air_zw(k)   * mask_top              &
270                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
271                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
272                                   ) * rho_air_zw(k-1) * mask_bottom           &
273                          ) * ddzw(k) * drho_air(k) * flag
[1]274             ENDDO
275
276!
277!--          Vertical diffusion at the first grid point above the surface,
278!--          if the momentum flux at the bottom is given by the Prandtl law or
279!--          if it is prescribed by the user.
280!--          Difference quotient of the momentum flux is not formed over half
281!--          of the grid spacing (2.0*ddzw(k)) any more, since the comparison
[1320]282!--          with other (LES) models showed that the values of the momentum
[1]283!--          flux becomes too large in this case.
284!--          The term containing w(k-1,..) (see above equation) is removed here
285!--          because the vertical velocity is assumed to be zero at the surface.
286             IF ( use_surface_fluxes )  THEN
287!
[2232]288!--             Default-type surfaces, upward-facing
289                surf_s = surf_def_h(0)%start_index(j,i)
290                surf_e = surf_def_h(0)%end_index(j,i)
291                DO  m = surf_s, surf_e
[1]292
[2232]293                   k   = surf_def_h(0)%k(m)
[1]294
[2232]295                   tend(k,j,i) = tend(k,j,i)                                   &
296                        + ( - ( - surf_def_h(0)%usws(m) )                      &
297                          ) * ddzw(k) * drho_air(k)
298                ENDDO
[102]299!
[2232]300!--             Default-type surfaces, dowward-facing
301                surf_s = surf_def_h(1)%start_index(j,i)
302                surf_e = surf_def_h(1)%end_index(j,i)
303                DO  m = surf_s, surf_e
304
305                   k   = surf_def_h(1)%k(m)
306
307                   tend(k,j,i) = tend(k,j,i)                                   &
308                        + ( - surf_def_h(1)%usws(m)                            &
309                          ) * ddzw(k) * drho_air(k)
310                ENDDO
[102]311!
[2232]312!--             Natural-type surfaces, upward-facing
313                surf_s = surf_lsm_h%start_index(j,i)
314                surf_e = surf_lsm_h%end_index(j,i)
315                DO  m = surf_s, surf_e
[102]316
[2232]317                   k   = surf_lsm_h%k(m)
318
319                   tend(k,j,i) = tend(k,j,i)                                   &
320                        + ( - ( - surf_lsm_h%usws(m) )                         &
321                          ) * ddzw(k) * drho_air(k)
322                ENDDO
323!
324!--             Urban-type surfaces, upward-facing
325                surf_s = surf_usm_h%start_index(j,i)
326                surf_e = surf_usm_h%end_index(j,i)
327                DO  m = surf_s, surf_e
328
329                   k   = surf_usm_h%k(m)
330
331                   tend(k,j,i) = tend(k,j,i)                                   &
332                       + ( - ( - surf_usm_h%usws(m) )                          &
333                         ) * ddzw(k) * drho_air(k)
334                ENDDO
335
[102]336             ENDIF
[2232]337!
338!--          Add momentum flux at model top
[2638]339             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]340                surf_s = surf_def_h(2)%start_index(j,i)
341                surf_e = surf_def_h(2)%end_index(j,i)
342                DO  m = surf_s, surf_e
[102]343
[2232]344                   k   = surf_def_h(2)%k(m)
345
346                   tend(k,j,i) = tend(k,j,i)                                   &
347                        + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
348                ENDDO
349             ENDIF
350
[1]351          ENDDO
352       ENDDO
353
354    END SUBROUTINE diffusion_u
355
356
357!------------------------------------------------------------------------------!
[1682]358! Description:
359! ------------
360!> Call for grid point i,j
[1]361!------------------------------------------------------------------------------!
[1001]362    SUBROUTINE diffusion_u_ij( i, j )
[1]363
[1320]364       USE arrays_3d,                                                          &
[2232]365           ONLY:  ddzu, ddzw, km, tend, u, v, w, drho_air, rho_air_zw
[1320]366       
367       USE control_parameters,                                                 &
[2232]368           ONLY:  constant_top_momentumflux, use_surface_fluxes,               &
369                  use_top_fluxes
[1320]370       
371       USE grid_variables,                                                     &
[2232]372           ONLY:  ddx, ddx2, ddy
[1320]373       
374       USE indices,                                                            &
[2232]375           ONLY:  nzb, nzt, wall_flags_0
376     
[1320]377       USE kinds
[1]378
[2232]379       USE surface_mod,                                                        &
380           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
381                   surf_usm_v
382
[1]383       IMPLICIT NONE
384
[2232]385       INTEGER(iwp) ::  i             !< running index x direction
386       INTEGER(iwp) ::  j             !< running index y direction
387       INTEGER(iwp) ::  k             !< running index z direction
388       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
389       INTEGER(iwp) ::  m             !< running index surface elements
390       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
391       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
[1]392
[2232]393       REAL(wp)     ::  flag          !< flag to mask topography grid points
394       REAL(wp)     ::  kmym          !<
395       REAL(wp)     ::  kmyp          !<
396       REAL(wp)     ::  kmzm          !<
397       REAL(wp)     ::  kmzp          !<
398       REAL(wp)     ::  mask_bottom   !< flag to mask vertical upward-facing surface       
399       REAL(wp)     ::  mask_north    !< flag to mask vertical surface north of the grid point
400       REAL(wp)     ::  mask_south    !< flag to mask vertical surface south of the grid point
401       REAL(wp)     ::  mask_top      !< flag to mask vertical downward-facing surface
402!
[1]403!--    Compute horizontal diffusion
[2232]404       DO  k = nzb+1, nzt
[1]405!
[2232]406!--       Predetermine flag to mask topography and wall-bounded grid points.
407!--       It is sufficient to masked only north- and south-facing surfaces, which
408!--       need special treatment for the u-component.
409          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   1 ) ) 
410          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 1 ) )
411          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 1 ) )
412!
[1]413!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]414          kmyp = 0.25_wp * ( km(k,j,i)+km(k,j+1,i)+km(k,j,i-1)+km(k,j+1,i-1) )
415          kmym = 0.25_wp * ( km(k,j,i)+km(k,j-1,i)+km(k,j,i-1)+km(k,j-1,i-1) )
[1]416
[1320]417          tend(k,j,i) = tend(k,j,i)                                            &
[2232]418                       + 2.0_wp * (                                            &
419                                 km(k,j,i)   * ( u(k,j,i+1) - u(k,j,i)   )     &
420                               - km(k,j,i-1) * ( u(k,j,i)   - u(k,j,i-1) )     &
421                                   ) * ddx2 * flag                             &
422                                 + (                                           &
423                  mask_north * kmyp * ( ( u(k,j+1,i) - u(k,j,i)     ) * ddy    &
424                                      + ( v(k,j+1,i) - v(k,j+1,i-1) ) * ddx    &
425                                      )                                        &
426                - mask_south * kmym * ( ( u(k,j,i)   - u(k,j-1,i)   ) * ddy    &
427                                      + ( v(k,j,i)   - v(k,j,i-1)   ) * ddx    &
428                                      )                                        &
429                                   ) * ddy  * flag
[1]430       ENDDO
431
432!
[2232]433!--    Add horizontal momentum flux u'v' at north- (l=0) and south-facing (l=1)
434!--    surfaces. Note, in the the flat case, loops won't be entered as
435!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
436!--    so far.           
437!--    Default-type surfaces
438       DO  l = 0, 1
439          surf_s = surf_def_v(l)%start_index(j,i)
440          surf_e = surf_def_v(l)%end_index(j,i)
441          DO  m = surf_s, surf_e
442             k           = surf_def_v(l)%k(m)
443             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddy
444          ENDDO   
445       ENDDO
[51]446!
[2232]447!--    Natural-type surfaces
448       DO  l = 0, 1
449          surf_s = surf_lsm_v(l)%start_index(j,i)
450          surf_e = surf_lsm_v(l)%end_index(j,i)
451          DO  m = surf_s, surf_e
452             k           = surf_lsm_v(l)%k(m)
453             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddy
454          ENDDO   
455       ENDDO
[1]456!
[2232]457!--    Urban-type surfaces
458       DO  l = 0, 1
459          surf_s = surf_usm_v(l)%start_index(j,i)
460          surf_e = surf_usm_v(l)%end_index(j,i)
461          DO  m = surf_s, surf_e
462             k           = surf_usm_v(l)%k(m)
463             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddy
464          ENDDO   
465       ENDDO
[1]466!
[2232]467!--    Compute vertical diffusion. In case of simulating a surface layer,
468!--    respective grid diffusive fluxes are masked (flag 8) within this
469!--    loop, and added further below, else, simple gradient approach is
470!--    applied. Model top is also mask if top-momentum flux is given.
471       DO  k = nzb+1, nzt
472!
473!--       Determine flags to mask topography below and above. Flag 1 is
474!--       used to mask topography in general, and flag 8 implies
475!--       information about use_surface_fluxes. Flag 9 is used to control
476!--       momentum flux at model top.
477          mask_bottom = MERGE( 1.0_wp, 0.0_wp,                                 &
478                               BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
479          mask_top    = MERGE( 1.0_wp, 0.0_wp,                                 &
480                               BTEST( wall_flags_0(k+1,j,i), 8 ) ) *           &
481                        MERGE( 1.0_wp, 0.0_wp,                                 &
482                               BTEST( wall_flags_0(k+1,j,i), 9 ) )
483          flag        = MERGE( 1.0_wp, 0.0_wp,                                 &
484                               BTEST( wall_flags_0(k,j,i), 1 ) ) 
485!
[1]486!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]487          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
488          kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
[1]489
[1320]490          tend(k,j,i) = tend(k,j,i)                                            &
[2232]491                        + ( kmzp * ( ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
492                                   + ( w(k,j,i)   - w(k,j,i-1) ) * ddx         &
493                                   ) * rho_air_zw(k)   * mask_top              &
494                          - kmzm * ( ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
495                                   + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
496                                   ) * rho_air_zw(k-1) * mask_bottom           &
497                          ) * ddzw(k) * drho_air(k) * flag
[1]498       ENDDO
499
500!
[2232]501!--    Vertical diffusion at the first surface grid points, if the
[1]502!--    momentum flux at the bottom is given by the Prandtl law or if it is
503!--    prescribed by the user.
504!--    Difference quotient of the momentum flux is not formed over half of
505!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]506!--    other (LES) models showed that the values of the momentum flux becomes
[1]507!--    too large in this case.
508       IF ( use_surface_fluxes )  THEN
509!
[2232]510!--       Default-type surfaces, upward-facing
511          surf_s = surf_def_h(0)%start_index(j,i)
512          surf_e = surf_def_h(0)%end_index(j,i)
513          DO  m = surf_s, surf_e
[1]514
[2232]515             k   = surf_def_h(0)%k(m)
[1]516
[2232]517             tend(k,j,i) = tend(k,j,i)                                         &
518                        + ( - ( - surf_def_h(0)%usws(m) )                      &
519                          ) * ddzw(k) * drho_air(k)
520          ENDDO
[102]521!
[2232]522!--       Default-type surfaces, dowward-facing (except for model-top fluxes)
523          surf_s = surf_def_h(1)%start_index(j,i)
524          surf_e = surf_def_h(1)%end_index(j,i)
525          DO  m = surf_s, surf_e
526
527             k   = surf_def_h(1)%k(m)
528
529             tend(k,j,i) = tend(k,j,i)                                         &
530                        + ( - surf_def_h(1)%usws(m)                            &
531                          ) * ddzw(k) * drho_air(k)
532          ENDDO
[102]533!
[2232]534!--       Natural-type surfaces, upward-facing
535          surf_s = surf_lsm_h%start_index(j,i)
536          surf_e = surf_lsm_h%end_index(j,i)
537          DO  m = surf_s, surf_e
[102]538
[2232]539             k   = surf_lsm_h%k(m)
540
541             tend(k,j,i) = tend(k,j,i)                                         &
542                        + ( - ( - surf_lsm_h%usws(m) )                         &
543                          ) * ddzw(k) * drho_air(k)
544          ENDDO
545!
546!--       Urban-type surfaces, upward-facing
547          surf_s = surf_usm_h%start_index(j,i)
548          surf_e = surf_usm_h%end_index(j,i)
549          DO  m = surf_s, surf_e
550
551             k   = surf_usm_h%k(m)
552
553             tend(k,j,i) = tend(k,j,i)                                         &
554                       + ( - ( - surf_usm_h%usws(m) )                          &
555                         ) * ddzw(k) * drho_air(k)
556          ENDDO
557
[102]558       ENDIF
[2232]559!
560!--    Add momentum flux at model top
[2638]561       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[2232]562          surf_s = surf_def_h(2)%start_index(j,i)
563          surf_e = surf_def_h(2)%end_index(j,i)
564          DO  m = surf_s, surf_e
[102]565
[2232]566             k   = surf_def_h(2)%k(m)
567
568             tend(k,j,i) = tend(k,j,i)                                         &
569                           + ( - surf_def_h(2)%usws(m) ) * ddzw(k) * drho_air(k)
570          ENDDO
571       ENDIF
572
573
[1]574    END SUBROUTINE diffusion_u_ij
575
576 END MODULE diffusion_u_mod
Note: See TracBrowser for help on using the repository browser.