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

Last change on this file since 4753 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
File size: 19.3 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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_w.f90 4583 2020-06-29 12:36:47Z raasch $
26! file re-formatted to follow the PALM coding standard
27!
28! 4360 2020-01-07 11:25:50Z suehring
29! Introduction of wall_flags_total_0, which currently sets bits based on static topography
30! information used in wall_flags_static_0
31!
32! 4329 2019-12-10 15:46:36Z motisi
33! Renamed wall_flags_0 to wall_flags_static_0
34!
35! 4182 2019-08-22 15:20:23Z scharf
36! Corrected "Former revisions" section
37!
38! 3655 2019-01-07 16:51:22Z knoop
39! OpenACC port for SPEC
40!
41! Revision 1.1  1997/09/12 06:24:11  raasch
42! Initial revision
43!
44!
45! Description:
46! ------------
47!> Diffusion term of the w-component
48!--------------------------------------------------------------------------------------------------!
49 MODULE diffusion_w_mod
50
51
52    PRIVATE
53    PUBLIC diffusion_w
54
55    INTERFACE diffusion_w
56       MODULE PROCEDURE diffusion_w
57       MODULE PROCEDURE diffusion_w_ij
58    END INTERFACE diffusion_w
59
60 CONTAINS
61
62
63!--------------------------------------------------------------------------------------------------!
64! Description:
65! ------------
66!> Call for all grid points
67!--------------------------------------------------------------------------------------------------!
68    SUBROUTINE diffusion_w
69
70       USE arrays_3d,                                                                              &
71           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
72
73       USE grid_variables,                                                                         &
74           ONLY :  ddx, ddy
75
76       USE indices,                                                                                &
77           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
78
79       USE kinds
80
81       USE surface_mod,                                                                            &
82           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
83
84       IMPLICIT NONE
85
86       INTEGER(iwp) ::  i             !< running index x direction
87       INTEGER(iwp) ::  j             !< running index y direction
88       INTEGER(iwp) ::  k             !< running index z direction
89       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
90       INTEGER(iwp) ::  m             !< running index surface elements
91       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
92       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
93
94       REAL(wp) ::  flag              !< flag to mask topography grid points
95       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
96       REAL(wp) ::  kmxp              !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
97       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
98       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
99       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
100       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
101       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
102       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
103
104
105
106       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
107       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
108       !$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
109       !$ACC PRESENT(wall_flags_total_0, km) &
110       !$ACC PRESENT(u, v, w) &
111       !$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
112       !$ACC PRESENT(surf_def_v(0:3)) &
113       !$ACC PRESENT(surf_lsm_v(0:3)) &
114       !$ACC PRESENT(surf_usm_v(0:3)) &
115       !$ACC PRESENT(tend)
116       DO  i = nxl, nxr
117          DO  j = nys, nyn
118             DO  k = nzb+1, nzt-1
119!
120!--             Predetermine flag to mask topography and wall-bounded grid points.
121                flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
122                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
123                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
124                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
125                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
126!
127!--             Interpolate eddy diffusivities on staggered gridpoints
128                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +                                   &
129                                   km(k+1,j,i) + km(k+1,j,i+1) )
130                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +                                   &
131                                   km(k+1,j,i) + km(k+1,j,i-1) )
132                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
133                                   km(k,j+1,i) + km(k+1,j+1,i) )
134                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +                                   &
135                                   km(k,j-1,i) + km(k+1,j-1,i) )
136
137                tend(k,j,i) = tend(k,j,i)                                                          &
138                              + ( mask_east *  kmxp * (                                            &
139                                          ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                      &
140                                        + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                &
141                                                      )                                            &
142                                - mask_west * kmxm *  (                                            &
143                                          ( w(k,j,i)     - w(k,j,i-1) ) * ddx                      &
144                                        + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                &
145                                                      )                                            &
146                                ) * ddx                                 * flag                     &
147                              + ( mask_north * kmyp * (                                            &
148                                          ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                      &
149                                        + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                &
150                                                      )                                            &
151                                - mask_south * kmym * (                                            &
152                                          ( w(k,j,i)     - w(k,j-1,i) ) * ddy                      &
153                                        + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                &
154                                                      )                                            &
155                                ) * ddy                                 * flag                     &
156                              + 2.0_wp * (                                                         &
157                                km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1)              &
158                                            * rho_air(k+1)                                         &
159                              - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                &
160                                            * rho_air(k)                                           &
161                                         ) * ddzu(k+1) * drho_air_zw(k) * flag
162             ENDDO
163
164!
165!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces.
166!--          Note, in the the flat case, loops won't be entered as start_index > end_index.
167!--          Furtermore, note, no vertical natural surfaces so far.
168!--          Default-type surfaces
169             DO  l = 0, 1
170                surf_s = surf_def_v(l)%start_index(j,i)
171                surf_e = surf_def_v(l)%end_index(j,i)
172                DO  m = surf_s, surf_e
173                   k           = surf_def_v(l)%k(m)
174                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
175                ENDDO
176             ENDDO
177!
178!--          Natural-type surfaces
179             DO  l = 0, 1
180                surf_s = surf_lsm_v(l)%start_index(j,i)
181                surf_e = surf_lsm_v(l)%end_index(j,i)
182                DO  m = surf_s, surf_e
183                   k           = surf_lsm_v(l)%k(m)
184                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
185                ENDDO
186             ENDDO
187!
188!--          Urban-type surfaces
189             DO  l = 0, 1
190                surf_s = surf_usm_v(l)%start_index(j,i)
191                surf_e = surf_usm_v(l)%end_index(j,i)
192                DO  m = surf_s, surf_e
193                   k           = surf_usm_v(l)%k(m)
194                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
195                ENDDO
196             ENDDO
197!
198!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surface.
199!--          Default-type surfaces
200             DO  l = 2, 3
201                surf_s = surf_def_v(l)%start_index(j,i)
202                surf_e = surf_def_v(l)%end_index(j,i)
203                DO  m = surf_s, surf_e
204                   k           = surf_def_v(l)%k(m)
205                   tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
206                ENDDO
207             ENDDO
208!
209!--          Natural-type surfaces
210             DO  l = 2, 3
211                surf_s = surf_lsm_v(l)%start_index(j,i)
212                surf_e = surf_lsm_v(l)%end_index(j,i)
213                DO  m = surf_s, surf_e
214                   k           = surf_lsm_v(l)%k(m)
215                   tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
216                ENDDO
217             ENDDO
218!
219!--          Urban-type surfaces
220             DO  l = 2, 3
221                surf_s = surf_usm_v(l)%start_index(j,i)
222                surf_e = surf_usm_v(l)%end_index(j,i)
223                DO  m = surf_s, surf_e
224                   k           = surf_usm_v(l)%k(m)
225                   tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
226                ENDDO
227             ENDDO
228
229          ENDDO
230       ENDDO
231
232    END SUBROUTINE diffusion_w
233
234
235!--------------------------------------------------------------------------------------------------!
236! Description:
237! ------------
238!> Call for grid point i,j
239!--------------------------------------------------------------------------------------------------!
240    SUBROUTINE diffusion_w_ij( i, j )
241
242       USE arrays_3d,                                                                              &
243           ONLY :  ddzu, ddzw, drho_air_zw, km, rho_air, tend, u, v, w
244
245       USE grid_variables,                                                                         &
246           ONLY :  ddx, ddy
247
248       USE indices,                                                                                &
249           ONLY :  nzb, nzt, wall_flags_total_0
250
251       USE kinds
252
253       USE surface_mod,                                                                            &
254           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
255
256       IMPLICIT NONE
257
258
259       INTEGER(iwp) ::  i             !< running index x direction
260       INTEGER(iwp) ::  j             !< running index y direction
261       INTEGER(iwp) ::  k             !< running index z direction
262       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
263       INTEGER(iwp) ::  m             !< running index surface elements
264       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
265       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
266
267       REAL(wp) ::  flag              !< flag to mask topography grid points
268       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
269       REAL(wp) ::  kmxp              !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
270       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
271       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
272       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
273       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
274       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
275       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
276
277
278       DO  k = nzb+1, nzt-1
279!
280!--       Predetermine flag to mask topography and wall-bounded grid points.
281          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i),   3 ) )
282          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
283          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
284          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
285          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
286!
287!--       Interpolate eddy diffusivities on staggered gridpoints
288          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
289          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
290          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
291          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
292
293          tend(k,j,i) = tend(k,j,i)                                                                &
294                        + ( mask_east *  kmxp * (                                                  &
295                                    ( w(k,j,i+1)   - w(k,j,i)   ) * ddx                            &
296                                  + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)                      &
297                                                )                                                  &
298                          - mask_west * kmxm *  (                                                  &
299                                    ( w(k,j,i)     - w(k,j,i-1) ) * ddx                            &
300                                  + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)                      &
301                                                )                                                  &
302                          ) * ddx                                 * flag                           &
303                        + ( mask_north * kmyp * (                                                  &
304                                    ( w(k,j+1,i)   - w(k,j,i)   ) * ddy                            &
305                                  + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)                      &
306                                                )                                                  &
307                          - mask_south * kmym * (                                                  &
308                                    ( w(k,j,i)     - w(k,j-1,i) ) * ddy                            &
309                                  + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)                      &
310                                                )                                                  &
311                          ) * ddy                                 * flag                           &
312                        + 2.0_wp * (                                                               &
313                          km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)                      &
314                                      * rho_air(k+1)                                               &
315                        - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)                      &
316                                      * rho_air(k)                                                 &
317                                   ) * ddzu(k+1) * drho_air_zw(k) * flag
318       ENDDO
319!
320!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1) surfaces. Note, in
321!--    the the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
322!--    vertical natural surfaces so far.
323!--    Default-type surfaces
324       DO  l = 0, 1
325          surf_s = surf_def_v(l)%start_index(j,i)
326          surf_e = surf_def_v(l)%end_index(j,i)
327          DO  m = surf_s, surf_e
328             k           = surf_def_v(l)%k(m)
329             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddy
330          ENDDO
331       ENDDO
332!
333!--    Natural-type surfaces
334       DO  l = 0, 1
335          surf_s = surf_lsm_v(l)%start_index(j,i)
336          surf_e = surf_lsm_v(l)%end_index(j,i)
337          DO  m = surf_s, surf_e
338             k           = surf_lsm_v(l)%k(m)
339             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddy
340          ENDDO
341       ENDDO
342!
343!--    Urban-type surfaces
344       DO  l = 0, 1
345          surf_s = surf_usm_v(l)%start_index(j,i)
346          surf_e = surf_usm_v(l)%end_index(j,i)
347          DO  m = surf_s, surf_e
348             k           = surf_usm_v(l)%k(m)
349             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddy
350          ENDDO
351       ENDDO
352!
353!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3) surfaces.
354!--    Default-type surfaces
355       DO  l = 2, 3
356          surf_s = surf_def_v(l)%start_index(j,i)
357          surf_e = surf_def_v(l)%end_index(j,i)
358          DO  m = surf_s, surf_e
359             k           = surf_def_v(l)%k(m)
360             tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_w(m) * ddx
361          ENDDO
362       ENDDO
363!
364!--    Natural-type surfaces
365       DO  l = 2, 3
366          surf_s = surf_lsm_v(l)%start_index(j,i)
367          surf_e = surf_lsm_v(l)%end_index(j,i)
368          DO  m = surf_s, surf_e
369             k           = surf_lsm_v(l)%k(m)
370             tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_w(m) * ddx
371          ENDDO
372       ENDDO
373!
374!--    Urban-type surfaces
375       DO  l = 2, 3
376          surf_s = surf_usm_v(l)%start_index(j,i)
377          surf_e = surf_usm_v(l)%end_index(j,i)
378          DO  m = surf_s, surf_e
379             k           = surf_usm_v(l)%k(m)
380             tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_w(m) * ddx
381          ENDDO
382       ENDDO
383
384
385    END SUBROUTINE diffusion_w_ij
386
387 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.