source: palm/trunk/SOURCE/diffusion_w.f90 @ 2039

Last change on this file since 2039 was 2038, checked in by knoop, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.1 KB
RevLine 
[1873]1!> @file diffusion_w.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]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!
[1818]17! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
[2038]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_w.f90 2038 2016-10-26 11:16:56Z gronemeier $
27!
[2038]28! 2037 2016-10-26 11:15:40Z knoop
29! Anelastic approximation implemented
30!
[2001]31! 2000 2016-08-20 18:09:15Z knoop
32! Forced header and separation lines into 80 columns
33!
[1874]34! 1873 2016-04-18 14:50:06Z maronga
35! Module renamed (removed _mod)
36!
37!
[1851]38! 1850 2016-04-08 13:29:27Z maronga
39! Module renamed
40!
41!
[1683]42! 1682 2015-10-07 23:56:08Z knoop
43! Code annotations made doxygen readable
44!
[1375]45! 1374 2014-04-25 12:55:07Z raasch
46! vsws + vswst removed from acc-present-list
47!
[1354]48! 1353 2014-04-08 15:21:23Z heinze
49! REAL constants provided with KIND-attribute
50!
[1341]51! 1340 2014-03-25 19:45:13Z kanani
52! REAL constants defined as wp-kind
53!
[1323]54! 1322 2014-03-20 16:38:49Z raasch
55! REAL constants defined as wp-kind
56!
[1321]57! 1320 2014-03-20 08:40:49Z raasch
[1320]58! ONLY-attribute added to USE-statements,
59! kind-parameters added to all INTEGER and REAL declaration statements,
60! kinds are defined in new module kinds,
61! revision history before 2012 removed,
62! comment fields (!:) to be used for variable explanations added to
63! all variable declaration statements
[1]64!
[1258]65! 1257 2013-11-08 15:18:40Z raasch
66! openacc loop and loop vector clauses removed, declare create moved after
67! the FORTRAN declaration statement
68!
[1132]69! 1128 2013-04-12 06:19:32Z raasch
70! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
71! j_north
72!
[1037]73! 1036 2012-10-22 13:43:42Z raasch
74! code put under GPL (PALM 3.9)
75!
[1017]76! 1015 2012-09-27 09:23:24Z raasch
77! accelerator version (*_acc) added
78!
[1002]79! 1001 2012-09-13 14:08:46Z raasch
80! arrays comunicated by module instead of parameter list
81!
[979]82! 978 2012-08-09 08:28:32Z fricke
83! outflow damping layer removed
84! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
85! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
86!
[1]87! Revision 1.1  1997/09/12 06:24:11  raasch
88! Initial revision
89!
90!
91! Description:
92! ------------
[1682]93!> Diffusion term of the w-component
[1]94!------------------------------------------------------------------------------!
[1682]95 MODULE diffusion_w_mod
96 
[1]97
[1320]98    USE wall_fluxes_mod,                                                       &
99        ONLY :  wall_fluxes, wall_fluxes_acc
[56]100
[1]101    PRIVATE
[1015]102    PUBLIC diffusion_w, diffusion_w_acc
[1]103
104    INTERFACE diffusion_w
105       MODULE PROCEDURE diffusion_w
106       MODULE PROCEDURE diffusion_w_ij
107    END INTERFACE diffusion_w
108
[1015]109    INTERFACE diffusion_w_acc
110       MODULE PROCEDURE diffusion_w_acc
111    END INTERFACE diffusion_w_acc
112
[1]113 CONTAINS
114
115
116!------------------------------------------------------------------------------!
[1682]117! Description:
118! ------------
119!> Call for all grid points
[1]120!------------------------------------------------------------------------------!
[1001]121    SUBROUTINE diffusion_w
[1]122
[1320]123       USE arrays_3d,                                                          &         
[2037]124           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]125           
126       USE control_parameters,                                                 & 
127           ONLY :  topography
128           
129       USE grid_variables,                                                     &     
130           ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
131           
132       USE indices,                                                            &           
133           ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
134           
135       USE kinds
[1]136
137       IMPLICIT NONE
138
[1682]139       INTEGER(iwp) ::  i     !<
140       INTEGER(iwp) ::  j     !<
141       INTEGER(iwp) ::  k     !<
[1320]142       
[1682]143       REAL(wp) ::  kmxm  !<
144       REAL(wp) ::  kmxp  !<
145       REAL(wp) ::  kmym  !<
146       REAL(wp) ::  kmyp  !<
[1001]147
[1682]148       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
149       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
[1]150
151
[56]152!
153!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
154!--    walls, if neccessary
155       IF ( topography /= 'flat' )  THEN
[1320]156          CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner,             &
[56]157                            nzb_w_outer, wall_w_x )
[1320]158          CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner,             &
[56]159                            nzb_w_outer, wall_w_y )
160       ENDIF
161
[1]162       DO  i = nxl, nxr
163          DO  j = nys, nyn
164             DO  k = nzb_w_outer(j,i)+1, nzt-1
165!
166!--             Interpolate eddy diffusivities on staggered gridpoints
[1322]167                kmxp = 0.25_wp *                                               &
[978]168                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]169                kmxm = 0.25_wp *                                               &
[978]170                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]171                kmyp = 0.25_wp *                                               &
[978]172                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]173                kmym = 0.25_wp *                                               &
[978]174                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]175
176                tend(k,j,i) = tend(k,j,i)                                      &
[978]177                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
178                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
179                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
180                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
[1]181                      &   ) * ddx                                              &
[978]182                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
183                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
184                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
185                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
[1]186                      &   ) * ddy                                              &
[1322]187                      & + 2.0_wp * (                                           &
[1]188                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
[2037]189                      &               * rho_air(k+1)                           &
[1]190                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
[2037]191                      &               * rho_air(k)                             &
192                      &            ) * ddzu(k+1) * drho_air_zw(k)
[1]193             ENDDO
194
195!
196!--          Wall functions at all vertical walls, where necessary
[1340]197             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
[51]198
[1]199                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
200!
201!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]202                   kmxp = 0.25_wp *                                            &
[978]203                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]204                   kmxm = 0.25_wp *                                            &
[978]205                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]206                   kmyp = 0.25_wp *                                            &
[978]207                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]208                   kmym = 0.25_wp *                                            &
[978]209                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]210
211                   tend(k,j,i) = tend(k,j,i)                                   &
212                                 + (   fwxp(j,i) * (                           &
[978]213                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
214                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
[1]215                                                   )                           &
216                                     - fwxm(j,i) * (                           &
[978]217                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
218                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
[1]219                                                   )                           &
[56]220                                     + wall_w_x(j,i) * wsus(k,j,i)             &
[1]221                                   ) * ddx                                     &
222                                 + (   fwyp(j,i) * (                           &
[978]223                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
224                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
[1]225                                                   )                           &
226                                     - fwym(j,i) * (                           &
[978]227                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
228                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
[1]229                                                   )                           &
[56]230                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
[1]231                                   ) * ddy                                     &
[1340]232                                 + 2.0_wp * (                                  &
[1]233                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
[2037]234                                       * rho_air(k+1)                          &
[1]235                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
[2037]236                                       * rho_air(k)                            &
237                                            ) * ddzu(k+1) * drho_air_zw(k)
[1]238                ENDDO
239             ENDIF
240
241          ENDDO
242       ENDDO
243
244    END SUBROUTINE diffusion_w
245
246
247!------------------------------------------------------------------------------!
[1682]248! Description:
249! ------------
250!> Call for all grid points - accelerator version
[1015]251!------------------------------------------------------------------------------!
252    SUBROUTINE diffusion_w_acc
253
[1320]254       USE arrays_3d,                                                          &
[2037]255           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]256           
257       USE control_parameters,                                                 &
258           ONLY :  topography
259           
260       USE grid_variables,                                                     &
261           ONLY : ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
262           
263       USE indices,                                                            &
264           ONLY :  i_left, i_right, j_north, j_south, nxl, nxr, nyn, nys, nzb, &
265                   nzb_w_inner, nzb_w_outer, nzt
266                   
267       USE kinds
[1015]268
269       IMPLICIT NONE
270
[1682]271       INTEGER(iwp) ::  i     !<
272       INTEGER(iwp) ::  j     !<
273       INTEGER(iwp) ::  k     !<
[1320]274       
[1682]275       REAL(wp) ::  kmxm  !<
276       REAL(wp) ::  kmxp  !<
277       REAL(wp) ::  kmym  !<
278       REAL(wp) ::  kmyp  !<
[1015]279
[1682]280       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
281       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
[1015]282       !$acc declare create ( wsus, wsvs )
283
284!
285!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
286!--    walls, if neccessary
287       IF ( topography /= 'flat' )  THEN
[1320]288          CALL wall_fluxes_acc( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp,          &
289                                nzb_w_inner, nzb_w_outer, wall_w_x )
290          CALL wall_fluxes_acc( wsvs, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp,          &
291                                nzb_w_inner, nzb_w_outer, wall_w_y )
[1015]292       ENDIF
293
[1374]294       !$acc kernels present ( u, v, w, km, tend )                             &
[1015]295       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
296       !$acc         present ( nzb_w_inner, nzb_w_outer )
[1128]297       DO  i = i_left, i_right
298          DO  j = j_south, j_north
[1015]299             DO  k = 1, nzt
300                IF ( k > nzb_w_outer(j,i) )  THEN
301!
302!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]303                   kmxp = 0.25_wp *                                            &
[1015]304                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]305                   kmxm = 0.25_wp *                                            &
[1015]306                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]307                   kmyp = 0.25_wp *                                            &
[1015]308                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]309                   kmym = 0.25_wp *                                            &
[1015]310                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
311
312                   tend(k,j,i) = tend(k,j,i)                                     &
313                         & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx        &
314                         &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)  &
315                         &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx          &
316                         &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)    &
317                         &   ) * ddx                                             &
318                         & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy        &
319                         &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)  &
320                         &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy          &
321                         &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)    &
322                         &   ) * ddy                                             &
[1322]323                         & + 2.0_wp * (                                          &
[1015]324                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
[2037]325                         &               * rho_air(k+1)                          &
[1015]326                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
[2037]327                         &               * rho_air(k)                            &
328                         &            ) * ddzu(k+1) * drho_air_zw(k)
[1015]329                ENDIF
330             ENDDO
331
332!
333!--          Wall functions at all vertical walls, where necessary
334             DO  k = 1,nzt
335
336                IF ( k > nzb_w_inner(j,i)  .AND.  k <= nzb_w_outer(j,i)  .AND. &
[1340]337                     wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
[1015]338!
339!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]340                   kmxp = 0.25_wp *                                            &
[1015]341                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]342                   kmxm = 0.25_wp *                                            &
[1015]343                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]344                   kmyp = 0.25_wp *                                            &
[1015]345                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]346                   kmym = 0.25_wp *                                            &
[1015]347                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
348
349                   tend(k,j,i) = tend(k,j,i)                                   &
350                                 + (   fwxp(j,i) * (                           &
351                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
352                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
353                                                   )                           &
354                                     - fwxm(j,i) * (                           &
355                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
356                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
357                                                   )                           &
358                                     + wall_w_x(j,i) * wsus(k,j,i)             &
359                                   ) * ddx                                     &
360                                 + (   fwyp(j,i) * (                           &
361                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
362                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
363                                                   )                           &
364                                     - fwym(j,i) * (                           &
365                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
366                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
367                                                   )                           &
368                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
369                                   ) * ddy                                     &
[1340]370                                 + 2.0_wp * (                                  &
[1015]371                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
[2037]372                                       * rho_air(k+1)                          &
[1015]373                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
[2037]374                                       * rho_air(k)                            &
375                                            ) * ddzu(k+1) * drho_air_zw(k)
[1015]376                ENDIF
377             ENDDO
378
379          ENDDO
380       ENDDO
381       !$acc end kernels
382
383    END SUBROUTINE diffusion_w_acc
384
385
386!------------------------------------------------------------------------------!
[1682]387! Description:
388! ------------
389!> Call for grid point i,j
[1]390!------------------------------------------------------------------------------!
[1001]391    SUBROUTINE diffusion_w_ij( i, j )
[1]392
[1320]393       USE arrays_3d,                                                          &         
[2037]394           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]395           
396       USE control_parameters,                                                 & 
397           ONLY :  topography
398           
399       USE grid_variables,                                                     &     
400           ONLY :  ddx, ddy, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y
401           
402       USE indices,                                                            &           
403           ONLY :  nxl, nxr, nyn, nys, nzb, nzb_w_inner, nzb_w_outer, nzt
404           
405       USE kinds
[1]406
407       IMPLICIT NONE
408
[1682]409       INTEGER(iwp) ::  i     !<
410       INTEGER(iwp) ::  j     !<
411       INTEGER(iwp) ::  k     !<
[1320]412       
[1682]413       REAL(wp) ::  kmxm  !<
414       REAL(wp) ::  kmxp  !<
415       REAL(wp) ::  kmym  !<
416       REAL(wp) ::  kmyp  !<
[1]417
[1320]418       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus
419       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs
[1]420
[1001]421
[1]422       DO  k = nzb_w_outer(j,i)+1, nzt-1
423!
424!--       Interpolate eddy diffusivities on staggered gridpoints
[1322]425          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
426          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
427          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
428          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]429
430          tend(k,j,i) = tend(k,j,i)                                            &
[978]431                      & + ( kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
432                      &   + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
433                      &   - kmxm * ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
434                      &   - kmxm * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
[1]435                      &   ) * ddx                                              &
[978]436                      & + ( kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
437                      &   + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
438                      &   - kmym * ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
439                      &   - kmym * ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
[1]440                      &   ) * ddy                                              &
[1322]441                      & + 2.0_wp * (                                           &
[1]442                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
[2037]443                      &               * rho_air(k+1)                           &
[1]444                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
[2037]445                      &               * rho_air(k)                             &
446                      &            ) * ddzu(k+1) * drho_air_zw(k)
[1]447       ENDDO
448
449!
450!--    Wall functions at all vertical walls, where necessary
[1340]451       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
[51]452
453!
454!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
[1340]455          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
[1320]456             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
457                               wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
[51]458          ELSE
[1340]459             wsus = 0.0_wp
[51]460          ENDIF
461
[1340]462          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
[1320]463             CALL wall_fluxes( i, j, nzb_w_inner(j,i)+1, nzb_w_outer(j,i),     &
464                               wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
[51]465          ELSE
[1340]466             wsvs = 0.0_wp
[51]467          ENDIF
468
[1]469          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
470!
471!--          Interpolate eddy diffusivities on staggered gridpoints
[1322]472             kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
473             kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
474             kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
475             kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]476
477             tend(k,j,i) = tend(k,j,i)                                         &
478                                 + (   fwxp(j,i) * (                           &
[978]479                            kmxp * ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
480                          + kmxp * ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
[1]481                                                   )                           &
482                                     - fwxm(j,i) * (                           &
[978]483                            kmxm * ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
484                          + kmxm * ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
[1]485                                                   )                           &
[51]486                                     + wall_w_x(j,i) * wsus(k)                 &
[1]487                                   ) * ddx                                     &
488                                 + (   fwyp(j,i) * (                           &
[978]489                            kmyp * ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
490                          + kmyp * ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
[1]491                                                   )                           &
492                                     - fwym(j,i) * (                           &
[978]493                            kmym * ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
494                          + kmym * ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
[1]495                                                   )                           &
[51]496                                     + wall_w_y(j,i) * wsvs(k)                 &
[1]497                                   ) * ddy                                     &
[1340]498                                 + 2.0_wp * (                                  &
[1]499                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
[2037]500                                       * rho_air(k+1)                          &
[1]501                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
[2037]502                                       * rho_air(k)                            &
503                                            ) * ddzu(k+1) * drho_air_zw(k)
[1]504          ENDDO
505       ENDIF
506
507    END SUBROUTINE diffusion_w_ij
508
[1323]509 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.