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

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

last commit documented / copyright update

  • Property svn:keywords set to Id
File size: 26.5 KB
Line 
1!> @file diffusion_v.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!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_v.f90 1818 2016-04-06 15:53:27Z maronga $
26!
27! 1740 2016-01-13 08:19:40Z raasch
28! unnecessary calculations of kmzm and kmzp in wall bounded parts removed
29!
30! 1682 2015-10-07 23:56:08Z knoop
31! Code annotations made doxygen readable
32!
33! 1340 2014-03-25 19:45:13Z kanani
34! REAL constants defined as wp-kind
35!
36! 1320 2014-03-20 08:40:49Z raasch
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
43!
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!
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!
52! 1036 2012-10-22 13:43:42Z raasch
53! code put under GPL (PALM 3.9)
54!
55! 1015 2012-09-27 09:23:24Z raasch
56! accelerator version (*_acc) added
57!
58! 1001 2012-09-13 14:08:46Z raasch
59! arrays comunicated by module instead of parameter list
60!
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!
65! Revision 1.1  1997/09/12 06:24:01  raasch
66! Initial revision
67!
68!
69! Description:
70! ------------
71!> Diffusion term of the v-component
72!------------------------------------------------------------------------------!
73 MODULE diffusion_v_mod
74 
75
76    USE wall_fluxes_mod
77
78    PRIVATE
79    PUBLIC diffusion_v, diffusion_v_acc
80
81    INTERFACE diffusion_v
82       MODULE PROCEDURE diffusion_v
83       MODULE PROCEDURE diffusion_v_ij
84    END INTERFACE diffusion_v
85
86    INTERFACE diffusion_v_acc
87       MODULE PROCEDURE diffusion_v_acc
88    END INTERFACE diffusion_v_acc
89
90 CONTAINS
91
92
93!------------------------------------------------------------------------------!
94! Description:
95! ------------
96!> Call for all grid points
97!------------------------------------------------------------------------------!
98    SUBROUTINE diffusion_v
99
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
115
116       IMPLICIT NONE
117
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  !<
125
126       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
127
128!
129!--    First calculate horizontal momentum flux v'u' at vertical walls,
130!--    if neccessary
131       IF ( topography /= 'flat' )  THEN
132          CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
133                            nzb_v_outer, wall_v )
134       ENDIF
135
136       DO  i = nxl, nxr
137          DO  j = nysv, nyn
138!
139!--          Compute horizontal diffusion
140             DO  k = nzb_v_outer(j,i)+1, nzt
141!
142!--             Interpolate eddy diffusivities on staggered gridpoints
143                kmxp = 0.25_wp * &
144                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
145                kmxm = 0.25_wp * &
146                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
147
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                                              &
154                      & + 2.0_wp * (                                           &
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) )      &
157                      &            ) * ddy2
158             ENDDO
159
160!
161!--          Wall functions at the left and right walls, respectively
162             IF ( wall_v(j,i) /= 0.0_wp )  THEN
163
164                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
165                   kmxp = 0.25_wp *                                            &
166                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
167                   kmxm = 0.25_wp *                                            &
168                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
169                   
170                   tend(k,j,i) = tend(k,j,i)                                   &
171                                 + 2.0_wp * (                                  &
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) ) &
174                                            ) * ddy2                           &
175                                 + (   fxp(j,i) * (                            &
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   &
178                                                  )                            &
179                                     - fxm(j,i) * (                            &
180                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
181                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
182                                                  )                            &
183                                     + wall_v(j,i) * vsus(k,j,i)               &
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.
191             DO  k = nzb_diff_v(j,i), nzt_diff
192!
193!--             Interpolate eddy diffusivities on staggered gridpoints
194                kmzp = 0.25_wp * &
195                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
196                kmzm = 0.25_wp * &
197                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
198
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                      &            )                                           &
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
215!--          comparison with other (LES) models showed that the values of
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
223                kmzp = 0.25_wp *                                               &
224                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
225
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)                                          &
231                      &   ) * ddzw(k)
232             ENDIF
233
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
237             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
238                k = nzt
239!
240!--             Interpolate eddy diffusivities on staggered gridpoints
241                kmzm = 0.25_wp *                                               &
242                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
243
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)    &
249                      &   ) * ddzw(k)
250             ENDIF
251
252          ENDDO
253       ENDDO
254
255    END SUBROUTINE diffusion_v
256
257
258!------------------------------------------------------------------------------!
259! Description:
260! ------------
261!> Call for all grid points - accelerator version
262!------------------------------------------------------------------------------!
263    SUBROUTINE diffusion_v_acc
264
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
280
281       IMPLICIT NONE
282
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  !<
290
291       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !<
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
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 )
300       ENDIF
301
302       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
303       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
304       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
305       DO  i = i_left, i_right
306          DO  j = j_south, j_north
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
313                   kmxp = 0.25_wp *                                            &
314                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
315                   kmxm = 0.25_wp *                                            &
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                                           &
324                         & + 2.0_wp * (                                        &
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) )   &
327                         &            ) * ddy2
328                ENDIF
329             ENDDO
330
331!
332!--          Wall functions at the left and right walls, respectively
333             DO  k = 1, nzt
334                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
335                    wall_v(j,i) /= 0.0_wp )  THEN
336
337                   kmxp = 0.25_wp *                                            &
338                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
339                   kmxm = 0.25_wp *                                            &
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)                                   &
343                                 + 2.0_wp * (                                  &
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) ) &
346                                            ) * ddy2                           &
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
367                   kmzp = 0.25_wp *                                            &
368                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
369                   kmzm = 0.25_wp *                                            &
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
392!--    comparison with other (LES) models showed that the values of
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
398          DO  i = i_left, i_right
399             DO  j = j_south, j_north
400         
401                k = nzb_v_inner(j,i)+1
402!
403!--             Interpolate eddy diffusivities on staggered gridpoints
404                kmzp = 0.25_wp *                                               &
405                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
406
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)                                          &
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
425          DO  i = i_left, i_right
426             DO  j = j_south, j_north
427
428!
429!--             Interpolate eddy diffusivities on staggered gridpoints
430                kmzm = 0.25_wp *                                               &
431                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
432
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)    &
438                      &   ) * ddzw(k)
439             ENDDO
440          ENDDO
441
442       ENDIF
443       !$acc end kernels
444
445    END SUBROUTINE diffusion_v_acc
446
447
448!------------------------------------------------------------------------------!
449! Description:
450! ------------
451!> Call for grid point i,j
452!------------------------------------------------------------------------------!
453    SUBROUTINE diffusion_v_ij( i, j )
454
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
468
469       IMPLICIT NONE
470
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  !<
478
479       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !<
480
481!
482!--    Compute horizontal diffusion
483       DO  k = nzb_v_outer(j,i)+1, nzt
484!
485!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
488
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                                              &
495                      & + 2.0_wp * (                                           &
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) )      &
498                      &            ) * ddy2
499       ENDDO
500
501!
502!--    Wall functions at the left and right walls, respectively
503       IF ( wall_v(j,i) /= 0.0_wp )  THEN
504
505!
506!--       Calculate the horizontal momentum flux v'u'
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 )
509
510          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
511             kmxp = 0.25_wp *                                                  &
512                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
513             kmxm = 0.25_wp *                                                  &
514                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
515
516             tend(k,j,i) = tend(k,j,i)                                         &
517                                 + 2.0_wp * (                                  &
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) ) &
520                                            ) * ddy2                           &
521                                 + (   fxp(j,i) * (                            &
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   &
524                                                  )                            &
525                                     - fxm(j,i) * (                            &
526                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
527                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
528                                                  )                            &
529                                     + wall_v(j,i) * vsus(k)                   &
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.
537       DO  k = nzb_diff_v(j,i), nzt_diff
538!
539!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
542
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                      &            )                                           &
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
559!--    other (LES) models showed that the values of the momentum flux becomes
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
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) )
568
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)                                          &
574                      &   ) * ddzw(k)
575       ENDIF
576
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
580       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
581          k = nzt
582!
583!--       Interpolate eddy diffusivities on staggered gridpoints
584          kmzm = 0.25_wp * &
585                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
586
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)    &
592                      &   ) * ddzw(k)
593       ENDIF
594
595    END SUBROUTINE diffusion_v_ij
596
597 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.