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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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