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

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

last commit documented

  • Property svn:keywords set to Id
File size: 22.5 KB
RevLine 
[1]1 MODULE diffusion_w_mod
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
[1375]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_w.f90 1375 2014-04-25 13:07:08Z suehring $
27!
[1375]28! 1374 2014-04-25 12:55:07Z raasch
29! vsws + vswst removed from acc-present-list
30!
[1354]31! 1353 2014-04-08 15:21:23Z heinze
32! REAL constants provided with KIND-attribute
33!
[1341]34! 1340 2014-03-25 19:45:13Z kanani
35! REAL constants defined as wp-kind
36!
[1323]37! 1322 2014-03-20 16:38:49Z raasch
38! REAL constants defined as wp-kind
39!
[1321]40! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1]47!
[1258]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!
[1132]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!
[1037]56! 1036 2012-10-22 13:43:42Z raasch
57! code put under GPL (PALM 3.9)
58!
[1017]59! 1015 2012-09-27 09:23:24Z raasch
60! accelerator version (*_acc) added
61!
[1002]62! 1001 2012-09-13 14:08:46Z raasch
63! arrays comunicated by module instead of parameter list
64!
[979]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!
[1]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
[1320]79    USE wall_fluxes_mod,                                                       &
80        ONLY :  wall_fluxes, wall_fluxes_acc
[56]81
[1]82    PRIVATE
[1015]83    PUBLIC diffusion_w, diffusion_w_acc
[1]84
85    INTERFACE diffusion_w
86       MODULE PROCEDURE diffusion_w
87       MODULE PROCEDURE diffusion_w_ij
88    END INTERFACE diffusion_w
89
[1015]90    INTERFACE diffusion_w_acc
91       MODULE PROCEDURE diffusion_w_acc
92    END INTERFACE diffusion_w_acc
93
[1]94 CONTAINS
95
96
97!------------------------------------------------------------------------------!
98! Call for all grid points
99!------------------------------------------------------------------------------!
[1001]100    SUBROUTINE diffusion_w
[1]101
[1320]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
[1]115
116       IMPLICIT NONE
117
[1320]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  !:
[1001]126
[1320]127       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !:
128       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !:
[1]129
130
[56]131!
132!--    First calculate horizontal momentum flux w'u' and/or w'v' at vertical
133!--    walls, if neccessary
134       IF ( topography /= 'flat' )  THEN
[1320]135          CALL wall_fluxes( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, nzb_w_inner,             &
[56]136                            nzb_w_outer, wall_w_x )
[1320]137          CALL wall_fluxes( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, nzb_w_inner,             &
[56]138                            nzb_w_outer, wall_w_y )
139       ENDIF
140
[1]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
[1322]146                kmxp = 0.25_wp *                                               &
[978]147                       ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]148                kmxm = 0.25_wp *                                               &
[978]149                       ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]150                kmyp = 0.25_wp *                                               &
[978]151                       ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]152                kmym = 0.25_wp *                                               &
[978]153                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]154
155                tend(k,j,i) = tend(k,j,i)                                      &
[978]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)     &
[1]160                      &   ) * ddx                                              &
[978]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)     &
[1]165                      &   ) * ddy                                              &
[1322]166                      & + 2.0_wp * (                                           &
[1]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)  &
[1340]169                      &            ) * ddzu(k+1)
[1]170             ENDDO
171
172!
173!--          Wall functions at all vertical walls, where necessary
[1340]174             IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
[51]175
[1]176                DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
177!
178!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]179                   kmxp = 0.25_wp *                                            &
[978]180                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]181                   kmxm = 0.25_wp *                                            &
[978]182                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]183                   kmyp = 0.25_wp *                                            &
[978]184                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]185                   kmym = 0.25_wp *                                            &
[978]186                          ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
[1]187
188                   tend(k,j,i) = tend(k,j,i)                                   &
189                                 + (   fwxp(j,i) * (                           &
[978]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)   &
[1]192                                                   )                           &
193                                     - fwxm(j,i) * (                           &
[978]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)   &
[1]196                                                   )                           &
[56]197                                     + wall_w_x(j,i) * wsus(k,j,i)             &
[1]198                                   ) * ddx                                     &
199                                 + (   fwyp(j,i) * (                           &
[978]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)   &
[1]202                                                   )                           &
203                                     - fwym(j,i) * (                           &
[978]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)   &
[1]206                                                   )                           &
[56]207                                     + wall_w_y(j,i) * wsvs(k,j,i)             &
[1]208                                   ) * ddy                                     &
[1340]209                                 + 2.0_wp * (                                  &
[1]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) &
[1340]212                                            ) * ddzu(k+1)
[1]213                ENDDO
214             ENDIF
215
216          ENDDO
217       ENDDO
218
219    END SUBROUTINE diffusion_w
220
221
222!------------------------------------------------------------------------------!
[1015]223! Call for all grid points - accelerator version
224!------------------------------------------------------------------------------!
225    SUBROUTINE diffusion_w_acc
226
[1320]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
[1015]241
242       IMPLICIT NONE
243
[1320]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  !:
[1015]252
[1320]253       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsus  !:
254       REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wsvs  !:
[1015]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
[1320]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 )
[1015]265       ENDIF
266
[1374]267       !$acc kernels present ( u, v, w, km, tend )                             &
[1015]268       !$acc         present ( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, wall_w_x, wall_w_y )           &
269       !$acc         present ( nzb_w_inner, nzb_w_outer )
[1128]270       DO  i = i_left, i_right
271          DO  j = j_south, j_north
[1015]272             DO  k = 1, nzt
273                IF ( k > nzb_w_outer(j,i) )  THEN
274!
275!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]276                   kmxp = 0.25_wp *                                            &
[1015]277                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]278                   kmxm = 0.25_wp *                                            &
[1015]279                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]280                   kmyp = 0.25_wp *                                            &
[1015]281                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]282                   kmym = 0.25_wp *                                            &
[1015]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                                             &
[1322]296                         & + 2.0_wp * (                                          &
[1015]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) &
[1340]299                         &            ) * ddzu(k+1)
[1015]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. &
[1340]308                     wall_w_x(j,i) /= 0.0_wp  .AND.  wall_w_y(j,i) /= 0.0_wp )  THEN
[1015]309!
310!--                Interpolate eddy diffusivities on staggered gridpoints
[1322]311                   kmxp = 0.25_wp *                                            &
[1015]312                          ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
[1322]313                   kmxm = 0.25_wp *                                            &
[1015]314                          ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
[1322]315                   kmyp = 0.25_wp *                                            &
[1015]316                          ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
[1322]317                   kmym = 0.25_wp *                                            &
[1015]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                                     &
[1340]341                                 + 2.0_wp * (                                  &
[1015]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) &
[1340]344                                            ) * ddzu(k+1)
[1015]345                ENDIF
346             ENDDO
347
348          ENDDO
349       ENDDO
350       !$acc end kernels
351
352    END SUBROUTINE diffusion_w_acc
353
354
355!------------------------------------------------------------------------------!
[1]356! Call for grid point i,j
357!------------------------------------------------------------------------------!
[1001]358    SUBROUTINE diffusion_w_ij( i, j )
[1]359
[1320]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
[1]373
374       IMPLICIT NONE
375
[1320]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  !:
[1]384
[1320]385       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsus
386       REAL(wp), DIMENSION(nzb:nzt+1) ::  wsvs
[1]387
[1001]388
[1]389       DO  k = nzb_w_outer(j,i)+1, nzt-1
390!
391!--       Interpolate eddy diffusivities on staggered gridpoints
[1322]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) )
[1]396
397          tend(k,j,i) = tend(k,j,i)                                            &
[978]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)     &
[1]402                      &   ) * ddx                                              &
[978]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)     &
[1]407                      &   ) * ddy                                              &
[1322]408                      & + 2.0_wp * (                                           &
[1]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)  &
[1340]411                      &            ) * ddzu(k+1)
[1]412       ENDDO
413
414!
415!--    Wall functions at all vertical walls, where necessary
[1340]416       IF ( wall_w_x(j,i) /= 0.0_wp  .OR.  wall_w_y(j,i) /= 0.0_wp )  THEN
[51]417
418!
419!--       Calculate the horizontal momentum fluxes w'u' and/or w'v'
[1340]420          IF ( wall_w_x(j,i) /= 0.0_wp )  THEN
[1320]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 )
[51]423          ELSE
[1340]424             wsus = 0.0_wp
[51]425          ENDIF
426
[1340]427          IF ( wall_w_y(j,i) /= 0.0_wp )  THEN
[1320]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 )
[51]430          ELSE
[1340]431             wsvs = 0.0_wp
[51]432          ENDIF
433
[1]434          DO  k = nzb_w_inner(j,i)+1, nzb_w_outer(j,i)
435!
436!--          Interpolate eddy diffusivities on staggered gridpoints
[1322]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) )
[1]441
442             tend(k,j,i) = tend(k,j,i)                                         &
443                                 + (   fwxp(j,i) * (                           &
[978]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)   &
[1]446                                                   )                           &
447                                     - fwxm(j,i) * (                           &
[978]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)   &
[1]450                                                   )                           &
[51]451                                     + wall_w_x(j,i) * wsus(k)                 &
[1]452                                   ) * ddx                                     &
453                                 + (   fwyp(j,i) * (                           &
[978]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)   &
[1]456                                                   )                           &
457                                     - fwym(j,i) * (                           &
[978]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)   &
[1]460                                                   )                           &
[51]461                                     + wall_w_y(j,i) * wsvs(k)                 &
[1]462                                   ) * ddy                                     &
[1340]463                                 + 2.0_wp * (                                  &
[1]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) &
[1340]466                                            ) * ddzu(k+1)
[1]467          ENDDO
468       ENDIF
469
470    END SUBROUTINE diffusion_w_ij
471
[1323]472 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.