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

Last change on this file since 1741 was 1741, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 26.5 KB
RevLine 
[1682]1!> @file diffusion_v.f90
[1036]2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
[1310]16! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]17!--------------------------------------------------------------------------------!
18!
[484]19! Current revisions:
[1]20! -----------------
[1341]21!
[1741]22!
[1321]23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 1741 2016-01-13 08:32:52Z raasch $
26!
[1741]27! 1740 2016-01-13 08:19:40Z raasch
28! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
29!
[1683]30! 1682 2015-10-07 23:56:08Z knoop
31! Code annotations made doxygen readable
32!
[1341]33! 1340 2014-03-25 19:45:13Z kanani
34! REAL constants defined as wp-kind
35!
[1321]36! 1320 2014-03-20 08:40:49Z raasch
[1320]37! ONLY-attribute added to USE-statements,
38! kind-parameters added to all INTEGER and REAL declaration statements,
39! kinds are defined in new module kinds,
40! revision history before 2012 removed,
41! comment fields (!:) to be used for variable explanations added to
42! all variable declaration statements
[1321]43!
[1258]44! 1257 2013-11-08 15:18:40Z raasch
45! openacc loop and loop vector clauses removed, declare create moved after
46! the FORTRAN declaration statement
47!
[1132]48! 1128 2013-04-12 06:19:32Z raasch
49! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
50! j_north
51!
[1037]52! 1036 2012-10-22 13:43:42Z raasch
53! code put under GPL (PALM 3.9)
54!
[1017]55! 1015 2012-09-27 09:23:24Z raasch
56! accelerator version (*_acc) added
57!
[1002]58! 1001 2012-09-13 14:08:46Z raasch
59! arrays comunicated by module instead of parameter list
60!
[979]61! 978 2012-08-09 08:28:32Z fricke
62! outflow damping layer removed
63! kmxm_x/_y and kmxp_x/_y change to kmxm and kmxp
64!
[1]65! Revision 1.1  1997/09/12 06:24:01  raasch
66! Initial revision
67!
68!
69! Description:
70! ------------
[1682]71!> Diffusion term of the v-component
[1]72!------------------------------------------------------------------------------!
[1682]73 MODULE diffusion_v_mod
74 
[1]75
[56]76    USE wall_fluxes_mod
77
[1]78    PRIVATE
[1015]79    PUBLIC diffusion_v, diffusion_v_acc
[1]80
81    INTERFACE diffusion_v
82       MODULE PROCEDURE diffusion_v
83       MODULE PROCEDURE diffusion_v_ij
84    END INTERFACE diffusion_v
85
[1015]86    INTERFACE diffusion_v_acc
87       MODULE PROCEDURE diffusion_v_acc
88    END INTERFACE diffusion_v_acc
89
[1]90 CONTAINS
91
92
93!------------------------------------------------------------------------------!
[1682]94! Description:
95! ------------
96!> Call for all grid points
[1]97!------------------------------------------------------------------------------!
[1001]98    SUBROUTINE diffusion_v
[1]99
[1320]100       USE arrays_3d,                                                          &
101           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
102       
103       USE control_parameters,                                                 &
104           ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
105                  use_top_fluxes
106       
107       USE grid_variables,                                                     &
108           ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
109       
110       USE indices,                                                            &
111           ONLY:  nxl, nxr, nyn, nys, nysv, nzb, nzb_diff_v, nzb_v_inner,      &
112                  nzb_v_outer, nzt, nzt_diff
113       
114       USE kinds
[1]115
116       IMPLICIT NONE
117
[1682]118       INTEGER(iwp) ::  i     !<
119       INTEGER(iwp) ::  j     !<
120       INTEGER(iwp) ::  k     !<
121       REAL(wp)     ::  kmxm  !<
122       REAL(wp)     ::  kmxp  !<
123       REAL(wp)     ::  kmzm  !<
124       REAL(wp)     ::  kmzp  !<
[1001]125
[1682]126       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
[1]127
[56]128!
129!--    First calculate horizontal momentum flux v'u' at vertical walls,
130!--    if neccessary
131       IF ( topography /= 'flat' )  THEN
[1320]132          CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
[56]133                            nzb_v_outer, wall_v )
134       ENDIF
135
[1]136       DO  i = nxl, nxr
[106]137          DO  j = nysv, nyn
[1]138!
139!--          Compute horizontal diffusion
140             DO  k = nzb_v_outer(j,i)+1, nzt
141!
142!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]143                kmxp = 0.25_wp * &
[978]144                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]145                kmxm = 0.25_wp * &
[978]146                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]147
[1320]148                tend(k,j,i) = tend(k,j,i)                                      &
149                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
150                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
151                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
152                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
153                      &   ) * ddx                                              &
[1340]154                      & + 2.0_wp * (                                           &
[1320]155                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
156                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
[1340]157                      &            ) * ddy2
[1]158             ENDDO
159
160!
161!--          Wall functions at the left and right walls, respectively
[1340]162             IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]163
[1]164                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]165                   kmxp = 0.25_wp *                                            &
[978]166                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]167                   kmxm = 0.25_wp *                                            &
[978]168                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
169                   
[1]170                   tend(k,j,i) = tend(k,j,i)                                   &
[1340]171                                 + 2.0_wp * (                                  &
[1]172                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
173                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]174                                            ) * ddy2                           &
[1]175                                 + (   fxp(j,i) * (                            &
[978]176                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
177                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
[1]178                                                  )                            &
179                                     - fxm(j,i) * (                            &
[978]180                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
181                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]182                                                  )                            &
[56]183                                     + wall_v(j,i) * vsus(k,j,i)               &
[1]184                                   ) * ddx
185                ENDDO
186             ENDIF
187
188!
189!--          Compute vertical diffusion. In case of simulating a Prandtl
190!--          layer, index k starts at nzb_v_inner+2.
[102]191             DO  k = nzb_diff_v(j,i), nzt_diff
[1]192!
193!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]194                kmzp = 0.25_wp * &
[1]195                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]196                kmzm = 0.25_wp * &
[1]197                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
198
[1320]199                tend(k,j,i) = tend(k,j,i)                                      &
200                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
201                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
202                      &            )                                           &
203                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
204                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
205                      &            )                                           &
[1]206                      &   ) * ddzw(k)
207             ENDDO
208
209!
210!--          Vertical diffusion at the first grid point above the surface,
211!--          if the momentum flux at the bottom is given by the Prandtl law
212!--          or if it is prescribed by the user.
213!--          Difference quotient of the momentum flux is not formed over
214!--          half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]215!--          comparison with other (LES) models showed that the values of
[1]216!--          the momentum flux becomes too large in this case.
217!--          The term containing w(k-1,..) (see above equation) is removed here
218!--          because the vertical velocity is assumed to be zero at the surface.
219             IF ( use_surface_fluxes )  THEN
220                k = nzb_v_inner(j,i)+1
221!
222!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]223                kmzp = 0.25_wp *                                               &
[1]224                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
225
[1320]226                tend(k,j,i) = tend(k,j,i)                                      &
227                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
228                      &   ) * ddzw(k)                                          &
229                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
230                      &   + vsws(j,i)                                          &
[1]231                      &   ) * ddzw(k)
232             ENDIF
233
[102]234!
235!--          Vertical diffusion at the first gridpoint below the top boundary,
236!--          if the momentum flux at the top is prescribed by the user
[103]237             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[102]238                k = nzt
239!
240!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]241                kmzm = 0.25_wp *                                               &
[102]242                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
243
[1320]244                tend(k,j,i) = tend(k,j,i)                                      &
245                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
246                      &   ) * ddzw(k)                                          &
247                      & + ( -vswst(j,i)                                        &
248                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
[102]249                      &   ) * ddzw(k)
250             ENDIF
251
[1]252          ENDDO
253       ENDDO
254
255    END SUBROUTINE diffusion_v
256
257
258!------------------------------------------------------------------------------!
[1682]259! Description:
260! ------------
261!> Call for all grid points - accelerator version
[1015]262!------------------------------------------------------------------------------!
263    SUBROUTINE diffusion_v_acc
264
[1320]265       USE arrays_3d,                                                          &
266           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
267       
268       USE control_parameters,                                                 &
269           ONLY:  constant_top_momentumflux, topography, use_surface_fluxes,   &
270                  use_top_fluxes
271       
272       USE grid_variables,                                                     &
273           ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
274       
275       USE indices,                                                            &
276           ONLY:  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb,  &
277                  nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
278       
279       USE kinds
[1015]280
281       IMPLICIT NONE
282
[1682]283       INTEGER(iwp) ::  i     !<
284       INTEGER(iwp) ::  j     !<
285       INTEGER(iwp) ::  k     !<
286       REAL(wp)     ::  kmxm  !<
287       REAL(wp)     ::  kmxp  !<
288       REAL(wp)     ::  kmzm  !<
289       REAL(wp)     ::  kmzp  !<
[1015]290
[1682]291       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
[1015]292       !$acc declare create ( vsus )
293
294!
295!--    First calculate horizontal momentum flux v'u' at vertical walls,
296!--    if neccessary
297       IF ( topography /= 'flat' )  THEN
[1320]298          CALL wall_fluxes_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp,          &
299                                nzb_v_inner, nzb_v_outer, wall_v )
[1015]300       ENDIF
301
[1320]302       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
303       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
[1015]304       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
[1128]305       DO  i = i_left, i_right
306          DO  j = j_south, j_north
[1015]307!
308!--          Compute horizontal diffusion
309             DO  k = 1, nzt
310                IF ( k > nzb_v_outer(j,i) )  THEN
311!
312!--                Interpolate eddy diffusivities on staggered gridpoints
[1340]313                   kmxp = 0.25_wp *                                            &
[1015]314                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]315                   kmxm = 0.25_wp *                                            &
[1015]316                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
317
318                   tend(k,j,i) = tend(k,j,i)                                   &
319                         & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx      &
320                         &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy      &
321                         &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx          &
322                         &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy          &
323                         &   ) * ddx                                           &
[1340]324                         & + 2.0_wp * (                                        &
[1015]325                         &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )   &
326                         &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )   &
[1340]327                         &            ) * ddy2
[1015]328                ENDIF
329             ENDDO
330
331!
332!--          Wall functions at the left and right walls, respectively
333             DO  k = 1, nzt
[1320]334                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
[1340]335                    wall_v(j,i) /= 0.0_wp )  THEN
[1015]336
[1340]337                   kmxp = 0.25_wp *                                            &
[1015]338                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]339                   kmxm = 0.25_wp *                                            &
[1015]340                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
341                   
342                   tend(k,j,i) = tend(k,j,i)                                   &
[1340]343                                 + 2.0_wp * (                                  &
[1015]344                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
345                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]346                                            ) * ddy2                           &
[1015]347                                 + (   fxp(j,i) * (                            &
348                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
349                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
350                                                  )                            &
351                                     - fxm(j,i) * (                            &
352                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
353                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
354                                                  )                            &
355                                     + wall_v(j,i) * vsus(k,j,i)               &
356                                   ) * ddx
357                ENDIF
358             ENDDO
359
360!
361!--          Compute vertical diffusion. In case of simulating a Prandtl
362!--          layer, index k starts at nzb_v_inner+2.
363             DO  k = 1, nzt_diff
364                IF ( k >= nzb_diff_v(j,i) )  THEN
365!
366!--                Interpolate eddy diffusivities on staggered gridpoints
[1340]367                   kmzp = 0.25_wp *                                            &
[1015]368                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]369                   kmzm = 0.25_wp *                                            &
[1015]370                          ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
371
372                   tend(k,j,i) = tend(k,j,i)                                   &
373                         & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)&
374                         &            + ( w(k,j,i)   - w(k,j-1,i) ) * ddy      &
375                         &            )                                        &
376                         &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)&
377                         &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy    &
378                         &            )                                        &
379                         &   ) * ddzw(k)
380                ENDIF
381             ENDDO
382
383          ENDDO
384       ENDDO
385
386!
387!--    Vertical diffusion at the first grid point above the surface,
388!--    if the momentum flux at the bottom is given by the Prandtl law
389!--    or if it is prescribed by the user.
390!--    Difference quotient of the momentum flux is not formed over
391!--    half of the grid spacing (2.0*ddzw(k)) any more, since the
[1320]392!--    comparison with other (LES) models showed that the values of
[1015]393!--    the momentum flux becomes too large in this case.
394!--    The term containing w(k-1,..) (see above equation) is removed here
395!--    because the vertical velocity is assumed to be zero at the surface.
396       IF ( use_surface_fluxes )  THEN
397
[1128]398          DO  i = i_left, i_right
399             DO  j = j_south, j_north
[1015]400         
401                k = nzb_v_inner(j,i)+1
402!
403!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]404                kmzp = 0.25_wp *                                               &
[1015]405                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
406
[1320]407                tend(k,j,i) = tend(k,j,i)                                      &
408                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
409                      &   ) * ddzw(k)                                          &
410                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
411                      &   + vsws(j,i)                                          &
[1015]412                      &   ) * ddzw(k)
413             ENDDO
414          ENDDO
415
416       ENDIF
417
418!
419!--    Vertical diffusion at the first gridpoint below the top boundary,
420!--    if the momentum flux at the top is prescribed by the user
421       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
422
423          k = nzt
424
[1128]425          DO  i = i_left, i_right
426             DO  j = j_south, j_north
[1015]427
428!
429!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]430                kmzm = 0.25_wp *                                               &
[1015]431                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
432
[1320]433                tend(k,j,i) = tend(k,j,i)                                      &
434                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
435                      &   ) * ddzw(k)                                          &
436                      & + ( -vswst(j,i)                                        &
437                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
[1015]438                      &   ) * ddzw(k)
439             ENDDO
440          ENDDO
441
442       ENDIF
443       !$acc end kernels
444
445    END SUBROUTINE diffusion_v_acc
446
447
448!------------------------------------------------------------------------------!
[1682]449! Description:
450! ------------
451!> Call for grid point i,j
[1]452!------------------------------------------------------------------------------!
[1001]453    SUBROUTINE diffusion_v_ij( i, j )
[1]454
[1320]455       USE arrays_3d,                                                          &
456           ONLY:  ddzu, ddzw, km, tend, u, v, vsws, vswst, w
457       
458       USE control_parameters,                                                 &
459           ONLY:  constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
460       
461       USE grid_variables,                                                     &
462           ONLY:  ddx, ddy, ddy2, fxm, fxp, wall_v
463       
464       USE indices,                                                            &
465           ONLY:  nzb, nzb_diff_v, nzb_v_inner, nzb_v_outer, nzt, nzt_diff
466       
467       USE kinds
[1]468
469       IMPLICIT NONE
470
[1682]471       INTEGER(iwp) ::  i     !<
472       INTEGER(iwp) ::  j     !<
473       INTEGER(iwp) ::  k     !<
474       REAL(wp)     ::  kmxm  !<
475       REAL(wp)     ::  kmxp  !<
476       REAL(wp)     ::  kmzm  !<
477       REAL(wp)     ::  kmzp  !<
[1]478
[1682]479       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
[1001]480
[1]481!
482!--    Compute horizontal diffusion
483       DO  k = nzb_v_outer(j,i)+1, nzt
484!
485!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]486          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
487          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]488
[1320]489          tend(k,j,i) = tend(k,j,i)                                            &
490                      & + ( kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx         &
491                      &   + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy         &
492                      &   - kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx             &
493                      &   - kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy             &
494                      &   ) * ddx                                              &
[1340]495                      & + 2.0_wp * (                                           &
[1320]496                      &           km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) )      &
497                      &         - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) )      &
[1340]498                      &            ) * ddy2
[1]499       ENDDO
500
501!
502!--    Wall functions at the left and right walls, respectively
[1340]503       IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]504
505!
506!--       Calculate the horizontal momentum flux v'u'
[1320]507          CALL wall_fluxes( i, j, nzb_v_inner(j,i)+1, nzb_v_outer(j,i),        &
508                            vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[51]509
[1]510          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]511             kmxp = 0.25_wp *                                                  &
[978]512                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]513             kmxm = 0.25_wp *                                                  &
[978]514                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]515
516             tend(k,j,i) = tend(k,j,i)                                         &
[1340]517                                 + 2.0_wp * (                                  &
[1]518                                       km(k,j,i)   * ( v(k,j+1,i) - v(k,j,i) ) &
519                                     - km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[1340]520                                            ) * ddy2                           &
[1]521                                 + (   fxp(j,i) * (                            &
[978]522                                  kmxp * ( v(k,j,i+1) - v(k,j,i)     ) * ddx   &
523                                + kmxp * ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy   &
[1]524                                                  )                            &
525                                     - fxm(j,i) * (                            &
[978]526                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
527                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]528                                                  )                            &
[51]529                                     + wall_v(j,i) * vsus(k)                   &
[1]530                                   ) * ddx
531          ENDDO
532       ENDIF
533
534!
535!--    Compute vertical diffusion. In case of simulating a Prandtl layer,
536!--    index k starts at nzb_v_inner+2.
[102]537       DO  k = nzb_diff_v(j,i), nzt_diff
[1]538!
539!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]540          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
541          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]542
[1320]543          tend(k,j,i) = tend(k,j,i)                                            &
544                      & + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
545                      &            + ( w(k,j,i) - w(k,j-1,i) ) * ddy           &
546                      &            )                                           &
547                      &   - kmzm * ( ( v(k,j,i)   - v(k-1,j,i)   ) * ddzu(k)   &
548                      &            + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy       &
549                      &            )                                           &
[1]550                      &   ) * ddzw(k)
551       ENDDO
552
553!
554!--    Vertical diffusion at the first grid point above the surface, if the
555!--    momentum flux at the bottom is given by the Prandtl law or if it is
556!--    prescribed by the user.
557!--    Difference quotient of the momentum flux is not formed over half of
558!--    the grid spacing (2.0*ddzw(k)) any more, since the comparison with
[1320]559!--    other (LES) models showed that the values of the momentum flux becomes
[1]560!--    too large in this case.
561!--    The term containing w(k-1,..) (see above equation) is removed here
562!--    because the vertical velocity is assumed to be zero at the surface.
563       IF ( use_surface_fluxes )  THEN
564          k = nzb_v_inner(j,i)+1
565!
566!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]567          kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]568
[1320]569          tend(k,j,i) = tend(k,j,i)                                            &
570                      & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy         &
571                      &   ) * ddzw(k)                                          &
572                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1)   &
573                      &   + vsws(j,i)                                          &
[1]574                      &   ) * ddzw(k)
575       ENDIF
576
[102]577!
578!--    Vertical diffusion at the first gridpoint below the top boundary,
579!--    if the momentum flux at the top is prescribed by the user
[103]580       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[102]581          k = nzt
582!
583!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]584          kmzm = 0.25_wp * &
[102]585                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
586
[1320]587          tend(k,j,i) = tend(k,j,i)                                            &
588                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy        &
589                      &   ) * ddzw(k)                                          &
590                      & + ( -vswst(j,i)                                        &
591                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)    &
[102]592                      &   ) * ddzw(k)
593       ENDIF
594
[1]595    END SUBROUTINE diffusion_v_ij
596
[1321]597 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.