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

Last change on this file since 4181 was 4180, checked in by scharf, 5 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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