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

Last change on this file since 1558 was 1341, checked in by kanani, 11 years ago

last commit documented

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