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

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

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