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

Last change on this file since 2696 was 2696, checked in by kanani, 4 years ago

Merge of branch palm4u into trunk

  • Property svn:keywords set to Id
File size: 19.9 KB
Line 
1!> @file diffusion_w.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_w.f90 2696 2017-12-14 17:12:51Z kanani $
27!
28! 2232 2017-05-30 17:47:52Z suehring
29! Adjustments to new topography and surface concept
30!
31! 2118 2017-01-17 16:38:49Z raasch
32! OpenACC version of subroutine removed
33!
34! 2037 2016-10-26 11:15:40Z knoop
35! Anelastic approximation implemented
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1873 2016-04-18 14:50:06Z maronga
41! Module renamed (removed _mod)
42!
43! 1850 2016-04-08 13:29:27Z maronga
44! Module renamed
45!
46! 1682 2015-10-07 23:56:08Z knoop
47! Code annotations made doxygen readable
48!
49! 1374 2014-04-25 12:55:07Z raasch
50! vsws + vswst removed from acc-present-list
51!
52! 1353 2014-04-08 15:21:23Z heinze
53! REAL constants provided with KIND-attribute
54!
55! 1340 2014-03-25 19:45:13Z kanani
56! REAL constants defined as wp-kind
57!
58! 1322 2014-03-20 16:38:49Z raasch
59! REAL constants defined as wp-kind
60!
61! 1320 2014-03-20 08:40:49Z raasch
62! ONLY-attribute added to USE-statements,
63! kind-parameters added to all INTEGER and REAL declaration statements,
64! kinds are defined in new module kinds,
65! revision history before 2012 removed,
66! comment fields (!:) to be used for variable explanations added to
67! all variable declaration statements
68!
69! 1257 2013-11-08 15:18:40Z raasch
70! openacc loop and loop vector clauses removed, declare create moved after
71! the FORTRAN declaration statement
72!
73! 1128 2013-04-12 06:19:32Z raasch
74! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
75! j_north
76!
77! 1036 2012-10-22 13:43:42Z raasch
78! code put under GPL (PALM 3.9)
79!
80! 1015 2012-09-27 09:23:24Z raasch
81! accelerator version (*_acc) added
82!
83! 1001 2012-09-13 14:08:46Z raasch
84! arrays comunicated by module instead of parameter list
85!
86! 978 2012-08-09 08:28:32Z fricke
87! outflow damping layer removed
88! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
89! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
90!
91! Revision 1.1  1997/09/12 06:24:11  raasch
92! Initial revision
93!
94!
95! Description:
96! ------------
97!> Diffusion term of the w-component
98!------------------------------------------------------------------------------!
99 MODULE diffusion_w_mod
100 
101
102    PRIVATE
103    PUBLIC diffusion_w
104
105    INTERFACE diffusion_w
106       MODULE PROCEDURE diffusion_w
107       MODULE PROCEDURE diffusion_w_ij
108    END INTERFACE diffusion_w
109
110 CONTAINS
111
112
113!------------------------------------------------------------------------------!
114! Description:
115! ------------
116!> Call for all grid points
117!------------------------------------------------------------------------------!
118    SUBROUTINE diffusion_w
119
120       USE arrays_3d,                                                          &         
121           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
122           
123       USE control_parameters,                                                 & 
124           ONLY :  topography
125           
126       USE grid_variables,                                                     &     
127           ONLY :  ddx, ddy
128           
129       USE indices,                                                            &           
130           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
131           
132       USE kinds
133
134       USE surface_mod,                                                        &
135           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
136
137       IMPLICIT NONE
138
139       INTEGER(iwp) ::  i             !< running index x direction
140       INTEGER(iwp) ::  j             !< running index y direction
141       INTEGER(iwp) ::  k             !< running index z direction
142       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
143       INTEGER(iwp) ::  m             !< running index surface elements
144       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
145       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
146       
147       REAL(wp) ::  flag              !< flag to mask topography grid points
148       REAL(wp) ::  kmxm              !<
149       REAL(wp) ::  kmxp              !<
150       REAL(wp) ::  kmym              !<
151       REAL(wp) ::  kmyp              !<
152       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
153       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
154       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
155       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
156
157
158
159       DO  i = nxl, nxr
160          DO  j = nys, nyn
161             DO  k = nzb+1, nzt-1
162!
163!--             Predetermine flag to mask topography and wall-bounded grid points.
164                flag       = MERGE( 1.0_wp, 0.0_wp,                            &
165                                    BTEST( wall_flags_0(k,j,i),   3 ) ) 
166                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
167                                    BTEST( wall_flags_0(k,j,i+1), 3 ) )
168                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
169                                    BTEST( wall_flags_0(k,j,i-1), 3 ) )
170                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
171                                    BTEST( wall_flags_0(k,j-1,i), 3 ) )
172                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
173                                    BTEST( wall_flags_0(k,j+1,i), 3 ) )
174!
175!--             Interpolate eddy diffusivities on staggered gridpoints
176                kmxp = 0.25_wp * ( km(k,j,i)   +   km(k,j,i+1) +               &
177                                   km(k+1,j,i) + km(k+1,j,i+1) )
178                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +               &
179                                   km(k+1,j,i) + km(k+1,j,i-1) )
180                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
181                                   km(k,j+1,i) + km(k+1,j+1,i) )
182                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
183                                   km(k,j-1,i) + km(k+1,j-1,i) )
184
185                tend(k,j,i) = tend(k,j,i)                                      &
186                       + ( mask_east *  kmxp * (                               &
187                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
188                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
189                                               )                               &
190                         - mask_west * kmxm *  (                               &
191                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
192                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
193                                               )                               &
194                         ) * ddx                                 * flag        &
195                       + ( mask_north * kmyp * (                               &
196                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
197                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
198                                               )                               &
199                         - mask_south * kmym * (                               &
200                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
201                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
202                                               )                               &
203                         ) * ddy                                 * flag        &
204                       + 2.0_wp * (                                            &
205                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
206                                     * rho_air(k+1)                            &
207                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
208                                     * rho_air(k)                              &
209                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
210             ENDDO
211
212!
213!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
214!--          surfaces. Note, in the the flat case, loops won't be entered as
215!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
216!--          so far.           
217!--          Default-type surfaces
218             DO  l = 0, 1
219                surf_s = surf_def_v(l)%start_index(j,i)
220                surf_e = surf_def_v(l)%end_index(j,i)
221                DO  m = surf_s, surf_e
222                   k           = surf_def_v(l)%k(m)
223                   tend(k,j,i) = tend(k,j,i) +                                 &
224                                     surf_def_v(l)%mom_flux_w(m) * ddy
225                ENDDO   
226             ENDDO
227!
228!--          Natural-type surfaces
229             DO  l = 0, 1
230                surf_s = surf_lsm_v(l)%start_index(j,i)
231                surf_e = surf_lsm_v(l)%end_index(j,i)
232                DO  m = surf_s, surf_e
233                   k           = surf_lsm_v(l)%k(m)
234                   tend(k,j,i) = tend(k,j,i) +                                 &
235                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
236                ENDDO   
237             ENDDO
238!
239!--          Urban-type surfaces
240             DO  l = 0, 1
241                surf_s = surf_usm_v(l)%start_index(j,i)
242                surf_e = surf_usm_v(l)%end_index(j,i)
243                DO  m = surf_s, surf_e
244                   k           = surf_usm_v(l)%k(m)
245                   tend(k,j,i) = tend(k,j,i) +                                 &
246                                     surf_usm_v(l)%mom_flux_w(m) * ddy
247                ENDDO   
248             ENDDO
249!
250!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
251!--          surface.
252!--          Default-type surfaces
253             DO  l = 2, 3
254                surf_s = surf_def_v(l)%start_index(j,i)
255                surf_e = surf_def_v(l)%end_index(j,i)
256                DO  m = surf_s, surf_e
257                   k           = surf_def_v(l)%k(m)
258                   tend(k,j,i) = tend(k,j,i) +                                 &
259                                     surf_def_v(l)%mom_flux_w(m) * ddx
260                ENDDO   
261             ENDDO
262!
263!--          Natural-type surfaces
264             DO  l = 2, 3
265                surf_s = surf_lsm_v(l)%start_index(j,i)
266                surf_e = surf_lsm_v(l)%end_index(j,i)
267                DO  m = surf_s, surf_e
268                   k           = surf_lsm_v(l)%k(m)
269                   tend(k,j,i) = tend(k,j,i) +                                 &
270                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
271                ENDDO   
272             ENDDO
273!
274!--          Urban-type surfaces
275             DO  l = 2, 3
276                surf_s = surf_usm_v(l)%start_index(j,i)
277                surf_e = surf_usm_v(l)%end_index(j,i)
278                DO  m = surf_s, surf_e
279                   k           = surf_usm_v(l)%k(m)
280                   tend(k,j,i) = tend(k,j,i) +                                 &
281                                     surf_usm_v(l)%mom_flux_w(m) * ddx
282                ENDDO   
283             ENDDO
284
285          ENDDO
286       ENDDO
287
288    END SUBROUTINE diffusion_w
289
290
291!------------------------------------------------------------------------------!
292! Description:
293! ------------
294!> Call for grid point i,j
295!------------------------------------------------------------------------------!
296    SUBROUTINE diffusion_w_ij( i, j )
297
298       USE arrays_3d,                                                          &         
299           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
300           
301       USE control_parameters,                                                 & 
302           ONLY :  topography
303           
304       USE grid_variables,                                                     &     
305           ONLY :  ddx, ddy
306           
307       USE indices,                                                            &           
308           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
309           
310       USE kinds
311
312       USE surface_mod,                                                        &
313           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
314
315       IMPLICIT NONE
316
317
318       INTEGER(iwp) ::  i             !< running index x direction
319       INTEGER(iwp) ::  j             !< running index y direction
320       INTEGER(iwp) ::  k             !< running index z direction
321       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
322       INTEGER(iwp) ::  m             !< running index surface elements
323       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
324       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
325       
326       REAL(wp) ::  flag              !< flag to mask topography grid points
327       REAL(wp) ::  kmxm              !<
328       REAL(wp) ::  kmxp              !<
329       REAL(wp) ::  kmym              !<
330       REAL(wp) ::  kmyp              !<
331       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
332       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
333       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
334       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
335
336
337       DO  k = nzb+1, nzt-1
338!
339!--       Predetermine flag to mask topography and wall-bounded grid points.
340          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   3 ) ) 
341          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
342          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
343          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
344          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
345!
346!--       Interpolate eddy diffusivities on staggered gridpoints
347          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
348          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
349          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
350          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
351
352          tend(k,j,i) = tend(k,j,i)                                            &
353                       + ( mask_east *  kmxp * (                               &
354                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
355                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
356                                               )                               &
357                         - mask_west * kmxm *  (                               &
358                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
359                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
360                                               )                               &
361                         ) * ddx                                 * flag        &
362                       + ( mask_north * kmyp * (                               &
363                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
364                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
365                                               )                               &
366                         - mask_south * kmym * (                               &
367                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
368                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
369                                               )                               &
370                         ) * ddy                                 * flag        &
371                       + 2.0_wp * (                                            &
372                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
373                                     * rho_air(k+1)                            &
374                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
375                                     * rho_air(k)                              &
376                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
377       ENDDO
378!
379!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
380!--    surfaces. Note, in the the flat case, loops won't be entered as
381!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
382!--    so far.           
383!--    Default-type surfaces
384       DO  l = 0, 1
385          surf_s = surf_def_v(l)%start_index(j,i)
386          surf_e = surf_def_v(l)%end_index(j,i)
387          DO  m = surf_s, surf_e
388             k           = surf_def_v(l)%k(m)
389             tend(k,j,i) = tend(k,j,i) +                                       &
390                                     surf_def_v(l)%mom_flux_w(m) * ddy
391          ENDDO   
392       ENDDO
393!
394!--    Natural-type surfaces
395       DO  l = 0, 1
396          surf_s = surf_lsm_v(l)%start_index(j,i)
397          surf_e = surf_lsm_v(l)%end_index(j,i)
398          DO  m = surf_s, surf_e
399             k           = surf_lsm_v(l)%k(m)
400             tend(k,j,i) = tend(k,j,i) +                                       &
401                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
402          ENDDO   
403       ENDDO
404!
405!--    Urban-type surfaces
406       DO  l = 0, 1
407          surf_s = surf_usm_v(l)%start_index(j,i)
408          surf_e = surf_usm_v(l)%end_index(j,i)
409          DO  m = surf_s, surf_e
410             k           = surf_usm_v(l)%k(m)
411             tend(k,j,i) = tend(k,j,i) +                                       &
412                                     surf_usm_v(l)%mom_flux_w(m) * ddy
413          ENDDO   
414       ENDDO
415!
416!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
417!--    surfaces.
418!--    Default-type surfaces
419       DO  l = 2, 3
420          surf_s = surf_def_v(l)%start_index(j,i)
421          surf_e = surf_def_v(l)%end_index(j,i)
422          DO  m = surf_s, surf_e
423             k           = surf_def_v(l)%k(m)
424             tend(k,j,i) = tend(k,j,i) +                                       &
425                                     surf_def_v(l)%mom_flux_w(m) * ddx
426          ENDDO   
427       ENDDO
428!
429!--    Natural-type surfaces
430       DO  l = 2, 3
431          surf_s = surf_lsm_v(l)%start_index(j,i)
432          surf_e = surf_lsm_v(l)%end_index(j,i)
433          DO  m = surf_s, surf_e
434             k           = surf_lsm_v(l)%k(m)
435             tend(k,j,i) = tend(k,j,i) +                                       &
436                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
437          ENDDO   
438       ENDDO
439!
440!--    Urban-type surfaces
441       DO  l = 2, 3
442          surf_s = surf_usm_v(l)%start_index(j,i)
443          surf_e = surf_usm_v(l)%end_index(j,i)
444          DO  m = surf_s, surf_e
445             k           = surf_usm_v(l)%k(m)
446             tend(k,j,i) = tend(k,j,i) +                                       &
447                                     surf_usm_v(l)%mom_flux_w(m) * ddx
448          ENDDO   
449       ENDDO
450
451
452    END SUBROUTINE diffusion_w_ij
453
454 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.