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

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