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

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

last commit documented

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