source: palm/trunk/SOURCE/diffusion_v_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 5 years ago

added _mod string to several filenames to meet the naming convection for modules

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