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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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