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

Last change on this file since 1375 was 1375, checked in by raasch, 8 years ago

last commit documented

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