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

Last change on this file since 4555 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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