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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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