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

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

last commit documented

  • Property svn:keywords set to Id
File size: 27.0 KB
RevLine 
[1]1 MODULE diffusion_v_mod
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_v.f90 1341 2014-03-25 19:48:09Z heinze $
27!
[1341]28! 1340 2014-03-25 19:45:13Z kanani
29! REAL constants defined as wp-kind
30!
[1321]31! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1321]38!
[1258]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!
[1132]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!
[1037]47! 1036 2012-10-22 13:43:42Z raasch
48! code put under GPL (PALM 3.9)
49!
[1017]50! 1015 2012-09-27 09:23:24Z raasch
51! accelerator version (*_acc) added
52!
[1002]53! 1001 2012-09-13 14:08:46Z raasch
54! arrays comunicated by module instead of parameter list
55!
[979]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!
[1]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
[56]69    USE wall_fluxes_mod
70
[1]71    PRIVATE
[1015]72    PUBLIC diffusion_v, diffusion_v_acc
[1]73
74    INTERFACE diffusion_v
75       MODULE PROCEDURE diffusion_v
76       MODULE PROCEDURE diffusion_v_ij
77    END INTERFACE diffusion_v
78
[1015]79    INTERFACE diffusion_v_acc
80       MODULE PROCEDURE diffusion_v_acc
81    END INTERFACE diffusion_v_acc
82
[1]83 CONTAINS
84
85
86!------------------------------------------------------------------------------!
87! Call for all grid points
88!------------------------------------------------------------------------------!
[1001]89    SUBROUTINE diffusion_v
[1]90
[1320]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
[1]106
107       IMPLICIT NONE
108
[1320]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  !:
[1001]116
[1320]117       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !:
[1]118
[56]119!
120!--    First calculate horizontal momentum flux v'u' at vertical walls,
121!--    if neccessary
122       IF ( topography /= 'flat' )  THEN
[1320]123          CALL wall_fluxes( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, nzb_v_inner, &
[56]124                            nzb_v_outer, wall_v )
125       ENDIF
126
[1]127       DO  i = nxl, nxr
[106]128          DO  j = nysv, nyn
[1]129!
130!--          Compute horizontal diffusion
131             DO  k = nzb_v_outer(j,i)+1, nzt
132!
133!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]134                kmxp = 0.25_wp * &
[978]135                       ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]136                kmxm = 0.25_wp * &
[978]137                       ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]138
[1320]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                                              &
[1340]145                      & + 2.0_wp * (                                           &
[1320]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) )      &
[1340]148                      &            ) * ddy2
[1]149             ENDDO
150
151!
152!--          Wall functions at the left and right walls, respectively
[1340]153             IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]154
[1]155                DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]156                   kmxp = 0.25_wp *                                            &
[978]157                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]158                   kmxm = 0.25_wp *                                            &
[978]159                          ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
160                   
[1]161                   tend(k,j,i) = tend(k,j,i)                                   &
[1340]162                                 + 2.0_wp * (                                  &
[1]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) ) &
[1340]165                                            ) * ddy2                           &
[1]166                                 + (   fxp(j,i) * (                            &
[978]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   &
[1]169                                                  )                            &
170                                     - fxm(j,i) * (                            &
[978]171                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
172                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]173                                                  )                            &
[56]174                                     + wall_v(j,i) * vsus(k,j,i)               &
[1]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.
[102]182             DO  k = nzb_diff_v(j,i), nzt_diff
[1]183!
184!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]185                kmzp = 0.25_wp * &
[1]186                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]187                kmzm = 0.25_wp * &
[1]188                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
189
[1320]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                      &            )                                           &
[1]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
[1320]206!--          comparison with other (LES) models showed that the values of
[1]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
[1340]214                kmzp = 0.25_wp *                                               &
[1]215                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]216                kmzm = 0.25_wp *                                               &
[1]217                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
218
[1320]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)                                          &
[1]224                      &   ) * ddzw(k)
225             ENDIF
226
[102]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
[103]230             IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[102]231                k = nzt
232!
233!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]234                kmzp = 0.25_wp *                                               &
[102]235                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]236                kmzm = 0.25_wp *                                               &
[102]237                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
238
[1320]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)    &
[102]244                      &   ) * ddzw(k)
245             ENDIF
246
[1]247          ENDDO
248       ENDDO
249
250    END SUBROUTINE diffusion_v
251
252
253!------------------------------------------------------------------------------!
[1015]254! Call for all grid points - accelerator version
255!------------------------------------------------------------------------------!
256    SUBROUTINE diffusion_v_acc
257
[1320]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
[1015]273
274       IMPLICIT NONE
275
[1320]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  !:
[1015]283
[1320]284       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus  !:
[1015]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
[1320]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 )
[1015]293       ENDIF
294
[1320]295       !$acc kernels present ( u, v, w, km, tend, vsws, vswst )                &
296       !$acc         present ( ddzu, ddzw, fxm, fxp, wall_v )                  &
[1015]297       !$acc         present ( nzb_v_inner, nzb_v_outer, nzb_diff_v )
[1128]298       DO  i = i_left, i_right
299          DO  j = j_south, j_north
[1015]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
[1340]306                   kmxp = 0.25_wp *                                            &
[1015]307                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]308                   kmxm = 0.25_wp *                                            &
[1015]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                                           &
[1340]317                         & + 2.0_wp * (                                        &
[1015]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) )   &
[1340]320                         &            ) * ddy2
[1015]321                ENDIF
322             ENDDO
323
324!
325!--          Wall functions at the left and right walls, respectively
326             DO  k = 1, nzt
[1320]327                IF( k > nzb_v_inner(j,i)  .AND.  k <= nzb_v_outer(j,i)  .AND.  &
[1340]328                    wall_v(j,i) /= 0.0_wp )  THEN
[1015]329
[1340]330                   kmxp = 0.25_wp *                                            &
[1015]331                          ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]332                   kmxm = 0.25_wp *                                            &
[1015]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)                                   &
[1340]336                                 + 2.0_wp * (                                  &
[1015]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) ) &
[1340]339                                            ) * ddy2                           &
[1015]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
[1340]360                   kmzp = 0.25_wp *                                            &
[1015]361                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]362                   kmzm = 0.25_wp *                                            &
[1015]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
[1320]385!--    comparison with other (LES) models showed that the values of
[1015]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
[1128]391          DO  i = i_left, i_right
392             DO  j = j_south, j_north
[1015]393         
394                k = nzb_v_inner(j,i)+1
395!
396!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]397                kmzp = 0.25_wp *                                               &
[1015]398                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]399                kmzm = 0.25_wp *                                               &
[1015]400                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
401
[1320]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)                                          &
[1015]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
[1128]420          DO  i = i_left, i_right
421             DO  j = j_south, j_north
[1015]422
423!
424!--             Interpolate eddy diffusivities on staggered gridpoints
[1340]425                kmzp = 0.25_wp *                                               &
[1015]426                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]427                kmzm = 0.25_wp *                                               &
[1015]428                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
429
[1320]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)    &
[1015]435                      &   ) * ddzw(k)
436             ENDDO
437          ENDDO
438
439       ENDIF
440       !$acc end kernels
441
442    END SUBROUTINE diffusion_v_acc
443
444
445!------------------------------------------------------------------------------!
[1]446! Call for grid point i,j
447!------------------------------------------------------------------------------!
[1001]448    SUBROUTINE diffusion_v_ij( i, j )
[1]449
[1320]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
[1]463
464       IMPLICIT NONE
465
[1320]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  !:
[1]473
[1320]474       REAL(wp), DIMENSION(nzb:nzt+1) ::  vsus  !:
[1001]475
[1]476!
477!--    Compute horizontal diffusion
478       DO  k = nzb_v_outer(j,i)+1, nzt
479!
480!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]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) )
[1]483
[1320]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                                              &
[1340]490                      & + 2.0_wp * (                                           &
[1320]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) )      &
[1340]493                      &            ) * ddy2
[1]494       ENDDO
495
496!
497!--    Wall functions at the left and right walls, respectively
[1340]498       IF ( wall_v(j,i) /= 0.0_wp )  THEN
[51]499
500!
501!--       Calculate the horizontal momentum flux v'u'
[1320]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 )
[51]504
[1]505          DO  k = nzb_v_inner(j,i)+1, nzb_v_outer(j,i)
[1340]506             kmxp = 0.25_wp *                                                  &
[978]507                    ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
[1340]508             kmxm = 0.25_wp *                                                  &
[978]509                    ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]510
511             tend(k,j,i) = tend(k,j,i)                                         &
[1340]512                                 + 2.0_wp * (                                  &
[1]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) ) &
[1340]515                                            ) * ddy2                           &
[1]516                                 + (   fxp(j,i) * (                            &
[978]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   &
[1]519                                                  )                            &
520                                     - fxm(j,i) * (                            &
[978]521                                  kmxm * ( v(k,j,i) - v(k,j,i-1) ) * ddx       &
522                                + kmxm * ( u(k,j,i) - u(k,j-1,i) ) * ddy       &
[1]523                                                  )                            &
[51]524                                     + wall_v(j,i) * vsus(k)                   &
[1]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.
[102]532       DO  k = nzb_diff_v(j,i), nzt_diff
[1]533!
534!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]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) )
[1]537
[1320]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                      &            )                                           &
[1]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
[1320]554!--    other (LES) models showed that the values of the momentum flux becomes
[1]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
[1340]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) )
[1]564
[1320]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)                                          &
[1]570                      &   ) * ddzw(k)
571       ENDIF
572
[102]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
[103]576       IF ( use_top_fluxes  .AND.  constant_top_momentumflux )  THEN
[102]577          k = nzt
578!
579!--       Interpolate eddy diffusivities on staggered gridpoints
[1340]580          kmzp = 0.25_wp * &
[102]581                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1340]582          kmzm = 0.25_wp * &
[102]583                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
584
[1320]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)    &
[102]590                      &   ) * ddzw(k)
591       ENDIF
592
[1]593    END SUBROUTINE diffusion_v_ij
594
[1321]595 END MODULE diffusion_v_mod
Note: See TracBrowser for help on using the repository browser.