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

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.1 KB
Line 
1!> @file diffusion_w.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
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.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_w.f90 2101 2017-01-05 16:42:31Z suehring $
27!
28! 2037 2016-10-26 11:15:40Z knoop
29! Anelastic approximation implemented
30!
31! 2000 2016-08-20 18:09:15Z knoop
32! Forced header and separation lines into 80 columns
33!
34! 1873 2016-04-18 14:50:06Z maronga
35! Module renamed (removed _mod)
36!
37!
38! 1850 2016-04-08 13:29:27Z maronga
39! Module renamed
40!
41!
42! 1682 2015-10-07 23:56:08Z knoop
43! Code annotations made doxygen readable
44!
45! 1374 2014-04-25 12:55:07Z raasch
46! vsws + vswst removed from acc-present-list
47!
48! 1353 2014-04-08 15:21:23Z heinze
49! REAL constants provided with KIND-attribute
50!
51! 1340 2014-03-25 19:45:13Z kanani
52! REAL constants defined as wp-kind
53!
54! 1322 2014-03-20 16:38:49Z raasch
55! REAL constants defined as wp-kind
56!
57! 1320 2014-03-20 08:40:49Z raasch
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
64!
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!
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!
73! 1036 2012-10-22 13:43:42Z raasch
74! code put under GPL (PALM 3.9)
75!
76! 1015 2012-09-27 09:23:24Z raasch
77! accelerator version (*_acc) added
78!
79! 1001 2012-09-13 14:08:46Z raasch
80! arrays comunicated by module instead of parameter list
81!
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!
87! Revision 1.1  1997/09/12 06:24:11  raasch
88! Initial revision
89!
90!
91! Description:
92! ------------
93!> Diffusion term of the w-component
94!------------------------------------------------------------------------------!
95 MODULE diffusion_w_mod
96 
97
98    USE wall_fluxes_mod,                                                       &
99        ONLY :  wall_fluxes, wall_fluxes_acc
100
101    PRIVATE
102    PUBLIC diffusion_w, diffusion_w_acc
103
104    INTERFACE diffusion_w
105       MODULE PROCEDURE diffusion_w
106       MODULE PROCEDURE diffusion_w_ij
107    END INTERFACE diffusion_w
108
109    INTERFACE diffusion_w_acc
110       MODULE PROCEDURE diffusion_w_acc
111    END INTERFACE diffusion_w_acc
112
113 CONTAINS
114
115
116!------------------------------------------------------------------------------!
117! Description:
118! ------------
119!> Call for all grid points
120!------------------------------------------------------------------------------!
121    SUBROUTINE diffusion_w
122
123       USE arrays_3d,                                                          &         
124           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
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
136
137       IMPLICIT NONE
138
139       INTEGER(iwp) ::  i     !<
140       INTEGER(iwp) ::  j     !<
141       INTEGER(iwp) ::  k     !<
142       
143       REAL(wp) ::  kmxm  !<
144       REAL(wp) ::  kmxp  !<
145       REAL(wp) ::  kmym  !<
146       REAL(wp) ::  kmyp  !<
147
148       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
149       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
150
151
152!
153!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
154!--    walls, if neccessary
155       IF ( topography /= 'flat' )  THEN
156          CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner,             &
157                            nzb_w_outer, wall_w_x )
158          CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner,             &
159                            nzb_w_outer, wall_w_y )
160       ENDIF
161
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
167                kmxp = 0.25_wp *                                               &
168                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
169                kmxm = 0.25_wp *                                               &
170                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
171                kmyp = 0.25_wp *                                               &
172                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
173                kmym = 0.25_wp *                                               &
174                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
175
176                tend(k,j,i) = tend(k,j,i)                                      &
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)     &
181                      &   ) * ddx                                              &
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)     &
186                      &   ) * ddy                                              &
187                      & + 2.0_wp * (                                           &
188                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
189                      &               * rho_air(k+1)                           &
190                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
191                      &               * rho_air(k)                             &
192                      &            ) * ddzu(k+1) * drho_air_zw(k)
193             ENDDO
194
195!
196!--          Wall functions at all vertical walls, where necessary
197             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
198
199                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
200!
201!--                Interpolate eddy diffusivities on staggered gridpoints
202                   kmxp = 0.25_wp *                                            &
203                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
204                   kmxm = 0.25_wp *                                            &
205                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
206                   kmyp = 0.25_wp *                                            &
207                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
208                   kmym = 0.25_wp *                                            &
209                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
210
211                   tend(k,j,i) = tend(k,j,i)                                   &
212                                 + (   fwxp(j,i) * (                           &
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)   &
215                                                   )                           &
216                                     - fwxm(j,i) * (                           &
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)   &
219                                                   )                           &
220                                     + wall_w_x(j,i) * wsus(k,j,i)             &
221                                   ) * ddx                                     &
222                                 + (   fwyp(j,i) * (                           &
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)   &
225                                                   )                           &
226                                     - fwym(j,i) * (                           &
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)   &
229                                                   )                           &
230                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
231                                   ) * ddy                                     &
232                                 + 2.0_wp * (                                  &
233                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
234                                       * rho_air(k+1)                          &
235                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
236                                       * rho_air(k)                            &
237                                            ) * ddzu(k+1) * drho_air_zw(k)
238                ENDDO
239             ENDIF
240
241          ENDDO
242       ENDDO
243
244    END SUBROUTINE diffusion_w
245
246
247!------------------------------------------------------------------------------!
248! Description:
249! ------------
250!> Call for all grid points - accelerator version
251!------------------------------------------------------------------------------!
252    SUBROUTINE diffusion_w_acc
253
254       USE arrays_3d,                                                          &
255           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
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
268
269       IMPLICIT NONE
270
271       INTEGER(iwp) ::  i     !<
272       INTEGER(iwp) ::  j     !<
273       INTEGER(iwp) ::  k     !<
274       
275       REAL(wp) ::  kmxm  !<
276       REAL(wp) ::  kmxp  !<
277       REAL(wp) ::  kmym  !<
278       REAL(wp) ::  kmyp  !<
279
280       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !<
281       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !<
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
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 )
292       ENDIF
293
294       !$acc kernels present ( u, v, w, km, tend )                             &
295       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
296       !$acc         present ( nzb_w_inner, nzb_w_outer )
297       DO  i = i_left, i_right
298          DO  j = j_south, j_north
299             DO  k = 1, nzt
300                IF ( k > nzb_w_outer(j,i) )  THEN
301!
302!--                Interpolate eddy diffusivities on staggered gridpoints
303                   kmxp = 0.25_wp *                                            &
304                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
305                   kmxm = 0.25_wp *                                            &
306                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
307                   kmyp = 0.25_wp *                                            &
308                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
309                   kmym = 0.25_wp *                                            &
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                                             &
323                         & + 2.0_wp * (                                          &
324                         &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
325                         &               * rho_air(k+1)                          &
326                         & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
327                         &               * rho_air(k)                            &
328                         &            ) * ddzu(k+1) * drho_air_zw(k)
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. &
337                     wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
338!
339!--                Interpolate eddy diffusivities on staggered gridpoints
340                   kmxp = 0.25_wp *                                            &
341                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
342                   kmxm = 0.25_wp *                                            &
343                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
344                   kmyp = 0.25_wp *                                            &
345                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
346                   kmym = 0.25_wp *                                            &
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                                     &
370                                 + 2.0_wp * (                                  &
371                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
372                                       * rho_air(k+1)                          &
373                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
374                                       * rho_air(k)                            &
375                                            ) * ddzu(k+1) * drho_air_zw(k)
376                ENDIF
377             ENDDO
378
379          ENDDO
380       ENDDO
381       !$acc end kernels
382
383    END SUBROUTINE diffusion_w_acc
384
385
386!------------------------------------------------------------------------------!
387! Description:
388! ------------
389!> Call for grid point i,j
390!------------------------------------------------------------------------------!
391    SUBROUTINE diffusion_w_ij( i, j )
392
393       USE arrays_3d,                                                          &         
394           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
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
406
407       IMPLICIT NONE
408
409       INTEGER(iwp) ::  i     !<
410       INTEGER(iwp) ::  j     !<
411       INTEGER(iwp) ::  k     !<
412       
413       REAL(wp) ::  kmxm  !<
414       REAL(wp) ::  kmxp  !<
415       REAL(wp) ::  kmym  !<
416       REAL(wp) ::  kmyp  !<
417
418       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus
419       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs
420
421
422       DO  k = nzb_w_outer(j,i)+1, nzt-1
423!
424!--       Interpolate eddy diffusivities on staggered gridpoints
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) )
429
430          tend(k,j,i) = tend(k,j,i)                                            &
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)     &
435                      &   ) * ddx                                              &
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)     &
440                      &   ) * ddy                                              &
441                      & + 2.0_wp * (                                           &
442                      &   km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)  &
443                      &               * rho_air(k+1)                           &
444                      & - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)  &
445                      &               * rho_air(k)                             &
446                      &            ) * ddzu(k+1) * drho_air_zw(k)
447       ENDDO
448
449!
450!--    Wall functions at all vertical walls, where necessary
451       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
452
453!
454!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
455          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
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 )
458          ELSE
459             wsus = 0.0_wp
460          ENDIF
461
462          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
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 )
465          ELSE
466             wsvs = 0.0_wp
467          ENDIF
468
469          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
470!
471!--          Interpolate eddy diffusivities on staggered gridpoints
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) )
476
477             tend(k,j,i) = tend(k,j,i)                                         &
478                                 + (   fwxp(j,i) * (                           &
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)   &
481                                                   )                           &
482                                     - fwxm(j,i) * (                           &
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)   &
485                                                   )                           &
486                                     + wall_w_x(j,i) * wsus(k)                 &
487                                   ) * ddx                                     &
488                                 + (   fwyp(j,i) * (                           &
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)   &
491                                                   )                           &
492                                     - fwym(j,i) * (                           &
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)   &
495                                                   )                           &
496                                     + wall_w_y(j,i) * wsvs(k)                 &
497                                   ) * ddy                                     &
498                                 + 2.0_wp * (                                  &
499                           km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1) &
500                                       * rho_air(k+1)                          &
501                         - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k) &
502                                       * rho_air(k)                            &
503                                            ) * ddzu(k+1) * drho_air_zw(k)
504          ENDDO
505       ENDIF
506
507    END SUBROUTINE diffusion_w_ij
508
509 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.