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

Last change on this file since 3436 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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