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

Last change on this file since 3768 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 21.1 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 3655 2019-01-07 16:51:22Z raasch $
27! OpenACC port for SPEC
28!
29! 3547 2018-11-21 13:21:24Z suehring
30! variables documented
31!
32! 3241 2018-09-12 15:02:00Z raasch
33! unused variables removed
34!
35! 2718 2018-01-02 08:49:38Z maronga
36! Corrected "Former revisions" section
37!
38! 2696 2017-12-14 17:12:51Z kanani
39! Change in file header (GPL part)
40!
41! 2233 2017-05-30 18:08:54Z suehring
42!
43! 2232 2017-05-30 17:47:52Z suehring
44! Adjustments to new topography and surface concept
45!
46! 2118 2017-01-17 16:38:49Z raasch
47! OpenACC version of subroutine removed
48!
49! 2037 2016-10-26 11:15:40Z knoop
50! Anelastic approximation implemented
51!
52! 2000 2016-08-20 18:09:15Z knoop
53! Forced header and separation lines into 80 columns
54!
55! 1873 2016-04-18 14:50:06Z maronga
56! Module renamed (removed _mod)
57!
58! 1850 2016-04-08 13:29:27Z maronga
59! Module renamed
60!
61! 1682 2015-10-07 23:56:08Z knoop
62! Code annotations made doxygen readable
63!
64! 1374 2014-04-25 12:55:07Z raasch
65! vsws + vswst removed from acc-present-list
66!
67! 1353 2014-04-08 15:21:23Z heinze
68! REAL constants provided with KIND-attribute
69!
70! 1340 2014-03-25 19:45:13Z kanani
71! REAL constants defined as wp-kind
72!
73! 1322 2014-03-20 16:38:49Z raasch
74! REAL constants defined as wp-kind
75!
76! 1320 2014-03-20 08:40:49Z raasch
77! ONLY-attribute added to USE-statements,
78! kind-parameters added to all INTEGER and REAL declaration statements,
79! kinds are defined in new module kinds,
80! revision history before 2012 removed,
81! comment fields (!:) to be used for variable explanations added to
82! all variable declaration statements
83!
84! 1257 2013-11-08 15:18:40Z raasch
85! openacc loop and loop vector clauses removed, declare create moved after
86! the FORTRAN declaration statement
87!
88! 1128 2013-04-12 06:19:32Z raasch
89! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
90! j_north
91!
92! 1036 2012-10-22 13:43:42Z raasch
93! code put under GPL (PALM 3.9)
94!
95! 1015 2012-09-27 09:23:24Z raasch
96! accelerator version (*_acc) added
97!
98! 1001 2012-09-13 14:08:46Z raasch
99! arrays comunicated by module instead of parameter list
100!
101! 978 2012-08-09 08:28:32Z fricke
102! outflow damping layer removed
103! kmxm_x/_z and kmxp_x/_z change to kmxm and kmxp
104! kmym_y/_z and kmyp_y/_z change to kmym and kmyp
105!
106! Revision 1.1  1997/09/12 06:24:11  raasch
107! Initial revision
108!
109!
110! Description:
111! ------------
112!> Diffusion term of the w-component
113!------------------------------------------------------------------------------!
114 MODULE diffusion_w_mod
115 
116
117    PRIVATE
118    PUBLIC diffusion_w
119
120    INTERFACE diffusion_w
121       MODULE PROCEDURE diffusion_w
122       MODULE PROCEDURE diffusion_w_ij
123    END INTERFACE diffusion_w
124
125 CONTAINS
126
127
128!------------------------------------------------------------------------------!
129! Description:
130! ------------
131!> Call for all grid points
132!------------------------------------------------------------------------------!
133    SUBROUTINE diffusion_w
134
135       USE arrays_3d,                                                          &         
136           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
137           
138       USE grid_variables,                                                     &     
139           ONLY :  ddx, ddy
140           
141       USE indices,                                                            &           
142           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
143           
144       USE kinds
145
146       USE surface_mod,                                                        &
147           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
148
149       IMPLICIT NONE
150
151       INTEGER(iwp) ::  i             !< running index x direction
152       INTEGER(iwp) ::  j             !< running index y direction
153       INTEGER(iwp) ::  k             !< running index z direction
154       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
155       INTEGER(iwp) ::  m             !< running index surface elements
156       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
157       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
158       
159       REAL(wp) ::  flag              !< flag to mask topography grid points
160       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
161       REAL(wp) ::  kmxp              !<diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
162       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
163       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
164       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
165       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
166       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
167       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
168
169
170
171       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
172       !$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmym, kmyp) &
173       !$ACC PRIVATE(mask_west, mask_east, mask_south, mask_north) &
174       !$ACC PRESENT(wall_flags_0, km) &
175       !$ACC PRESENT(u, v, w) &
176       !$ACC PRESENT(ddzu, ddzw, rho_air, drho_air_zw) &
177       !$ACC PRESENT(surf_def_v(0:3)) &
178       !$ACC PRESENT(surf_lsm_v(0:3)) &
179       !$ACC PRESENT(surf_usm_v(0:3)) &
180       !$ACC PRESENT(tend)
181       DO  i = nxl, nxr
182          DO  j = nys, nyn
183             DO  k = nzb+1, nzt-1
184!
185!--             Predetermine flag to mask topography and wall-bounded grid points.
186                flag       = MERGE( 1.0_wp, 0.0_wp,                            &
187                                    BTEST( wall_flags_0(k,j,i),   3 ) ) 
188                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
189                                    BTEST( wall_flags_0(k,j,i+1), 3 ) )
190                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
191                                    BTEST( wall_flags_0(k,j,i-1), 3 ) )
192                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
193                                    BTEST( wall_flags_0(k,j-1,i), 3 ) )
194                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
195                                    BTEST( wall_flags_0(k,j+1,i), 3 ) )
196!
197!--             Interpolate eddy diffusivities on staggered gridpoints
198                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +               &
199                                   km(k+1,j,i) + km(k+1,j,i+1) )
200                kmxm = 0.25_wp * ( km(k,j,i)   + km(k,j,i-1)   +               &
201                                   km(k+1,j,i) + km(k+1,j,i-1) )
202                kmyp = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
203                                   km(k,j+1,i) + km(k+1,j+1,i) )
204                kmym = 0.25_wp * ( km(k,j,i)   + km(k+1,j,i)   +               &
205                                   km(k,j-1,i) + km(k+1,j-1,i) )
206
207                tend(k,j,i) = tend(k,j,i)                                      &
208                       + ( mask_east *  kmxp * (                               &
209                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
210                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
211                                               )                               &
212                         - mask_west * kmxm *  (                               &
213                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
214                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
215                                               )                               &
216                         ) * ddx                                 * flag        &
217                       + ( mask_north * kmyp * (                               &
218                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
219                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
220                                               )                               &
221                         - mask_south * kmym * (                               &
222                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
223                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
224                                               )                               &
225                         ) * ddy                                 * flag        &
226                       + 2.0_wp * (                                            &
227                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) )   * ddzw(k+1) &
228                                     * rho_air(k+1)                            &
229                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
230                                     * rho_air(k)                              &
231                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
232             ENDDO
233
234!
235!--          Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
236!--          surfaces. Note, in the the flat case, loops won't be entered as
237!--          start_index > end_index. Furtermore, note, no vertical natural surfaces
238!--          so far.           
239!--          Default-type surfaces
240             DO  l = 0, 1
241                surf_s = surf_def_v(l)%start_index(j,i)
242                surf_e = surf_def_v(l)%end_index(j,i)
243                DO  m = surf_s, surf_e
244                   k           = surf_def_v(l)%k(m)
245                   tend(k,j,i) = tend(k,j,i) +                                 &
246                                     surf_def_v(l)%mom_flux_w(m) * ddy
247                ENDDO   
248             ENDDO
249!
250!--          Natural-type surfaces
251             DO  l = 0, 1
252                surf_s = surf_lsm_v(l)%start_index(j,i)
253                surf_e = surf_lsm_v(l)%end_index(j,i)
254                DO  m = surf_s, surf_e
255                   k           = surf_lsm_v(l)%k(m)
256                   tend(k,j,i) = tend(k,j,i) +                                 &
257                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
258                ENDDO   
259             ENDDO
260!
261!--          Urban-type surfaces
262             DO  l = 0, 1
263                surf_s = surf_usm_v(l)%start_index(j,i)
264                surf_e = surf_usm_v(l)%end_index(j,i)
265                DO  m = surf_s, surf_e
266                   k           = surf_usm_v(l)%k(m)
267                   tend(k,j,i) = tend(k,j,i) +                                 &
268                                     surf_usm_v(l)%mom_flux_w(m) * ddy
269                ENDDO   
270             ENDDO
271!
272!--          Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
273!--          surface.
274!--          Default-type surfaces
275             DO  l = 2, 3
276                surf_s = surf_def_v(l)%start_index(j,i)
277                surf_e = surf_def_v(l)%end_index(j,i)
278                DO  m = surf_s, surf_e
279                   k           = surf_def_v(l)%k(m)
280                   tend(k,j,i) = tend(k,j,i) +                                 &
281                                     surf_def_v(l)%mom_flux_w(m) * ddx
282                ENDDO   
283             ENDDO
284!
285!--          Natural-type surfaces
286             DO  l = 2, 3
287                surf_s = surf_lsm_v(l)%start_index(j,i)
288                surf_e = surf_lsm_v(l)%end_index(j,i)
289                DO  m = surf_s, surf_e
290                   k           = surf_lsm_v(l)%k(m)
291                   tend(k,j,i) = tend(k,j,i) +                                 &
292                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
293                ENDDO   
294             ENDDO
295!
296!--          Urban-type surfaces
297             DO  l = 2, 3
298                surf_s = surf_usm_v(l)%start_index(j,i)
299                surf_e = surf_usm_v(l)%end_index(j,i)
300                DO  m = surf_s, surf_e
301                   k           = surf_usm_v(l)%k(m)
302                   tend(k,j,i) = tend(k,j,i) +                                 &
303                                     surf_usm_v(l)%mom_flux_w(m) * ddx
304                ENDDO   
305             ENDDO
306
307          ENDDO
308       ENDDO
309
310    END SUBROUTINE diffusion_w
311
312
313!------------------------------------------------------------------------------!
314! Description:
315! ------------
316!> Call for grid point i,j
317!------------------------------------------------------------------------------!
318    SUBROUTINE diffusion_w_ij( i, j )
319
320       USE arrays_3d,                                                          &         
321           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
322           
323       USE grid_variables,                                                     &     
324           ONLY :  ddx, ddy
325           
326       USE indices,                                                            &           
327           ONLY :  nzb, nzt, wall_flags_0
328           
329       USE kinds
330
331       USE surface_mod,                                                        &
332           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
333
334       IMPLICIT NONE
335
336
337       INTEGER(iwp) ::  i             !< running index x direction
338       INTEGER(iwp) ::  j             !< running index y direction
339       INTEGER(iwp) ::  k             !< running index z direction
340       INTEGER(iwp) ::  l             !< running index of surface type, south- or north-facing wall
341       INTEGER(iwp) ::  m             !< running index surface elements
342       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
343       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
344       
345       REAL(wp) ::  flag              !< flag to mask topography grid points
346       REAL(wp) ::  kmxm              !< diffusion coefficient on leftward side of the w-gridbox - interpolated onto xu-y grid
347       REAL(wp) ::  kmxp              !< diffusion coefficient on rightward side of the w-gridbox - interpolated onto xu-y grid
348       REAL(wp) ::  kmym              !< diffusion coefficient on southward side of the w-gridbox - interpolated onto x-yv grid
349       REAL(wp) ::  kmyp              !< diffusion coefficient on northward side of the w-gridbox - interpolated onto x-yv grid
350       REAL(wp) ::  mask_west         !< flag to mask vertical wall west of the grid point
351       REAL(wp) ::  mask_east         !< flag to mask vertical wall east of the grid point
352       REAL(wp) ::  mask_south        !< flag to mask vertical wall south of the grid point
353       REAL(wp) ::  mask_north        !< flag to mask vertical wall north of the grid point
354
355
356       DO  k = nzb+1, nzt-1
357!
358!--       Predetermine flag to mask topography and wall-bounded grid points.
359          flag       = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i),   3 ) ) 
360          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 3 ) )
361          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 3 ) )
362          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 3 ) )
363          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 3 ) )
364!
365!--       Interpolate eddy diffusivities on staggered gridpoints
366          kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k+1,j,i)+km(k+1,j,i+1) )
367          kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k+1,j,i)+km(k+1,j,i-1) )
368          kmyp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j+1,i)+km(k+1,j+1,i) )
369          kmym = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
370
371          tend(k,j,i) = tend(k,j,i)                                            &
372                       + ( mask_east *  kmxp * (                               &
373                                   ( w(k,j,i+1)   - w(k,j,i)   ) * ddx         &
374                                 + ( u(k+1,j,i+1) - u(k,j,i+1) ) * ddzu(k+1)   &
375                                               )                               &
376                         - mask_west * kmxm *  (                               &
377                                   ( w(k,j,i)     - w(k,j,i-1) ) * ddx         &
378                                 + ( u(k+1,j,i)   - u(k,j,i)   ) * ddzu(k+1)   &
379                                               )                               &
380                         ) * ddx                                 * flag        &
381                       + ( mask_north * kmyp * (                               &
382                                   ( w(k,j+1,i)   - w(k,j,i)   ) * ddy         &
383                                 + ( v(k+1,j+1,i) - v(k,j+1,i) ) * ddzu(k+1)   &
384                                               )                               &
385                         - mask_south * kmym * (                               &
386                                   ( w(k,j,i)     - w(k,j-1,i) ) * ddy         &
387                                 + ( v(k+1,j,i)   - v(k,j,i)   ) * ddzu(k+1)   &
388                                               )                               &
389                         ) * ddy                                 * flag        &
390                       + 2.0_wp * (                                            &
391                         km(k+1,j,i) * ( w(k+1,j,i) - w(k,j,i) ) * ddzw(k+1)   &
392                                     * rho_air(k+1)                            &
393                       - km(k,j,i)   * ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)   &
394                                     * rho_air(k)                              &
395                                  ) * ddzu(k+1) * drho_air_zw(k) * flag
396       ENDDO
397!
398!--    Add horizontal momentum flux v'w' at north- (l=0) and south-facing (l=1)
399!--    surfaces. Note, in the the flat case, loops won't be entered as
400!--    start_index > end_index. Furtermore, note, no vertical natural surfaces
401!--    so far.           
402!--    Default-type surfaces
403       DO  l = 0, 1
404          surf_s = surf_def_v(l)%start_index(j,i)
405          surf_e = surf_def_v(l)%end_index(j,i)
406          DO  m = surf_s, surf_e
407             k           = surf_def_v(l)%k(m)
408             tend(k,j,i) = tend(k,j,i) +                                       &
409                                     surf_def_v(l)%mom_flux_w(m) * ddy
410          ENDDO   
411       ENDDO
412!
413!--    Natural-type surfaces
414       DO  l = 0, 1
415          surf_s = surf_lsm_v(l)%start_index(j,i)
416          surf_e = surf_lsm_v(l)%end_index(j,i)
417          DO  m = surf_s, surf_e
418             k           = surf_lsm_v(l)%k(m)
419             tend(k,j,i) = tend(k,j,i) +                                       &
420                                     surf_lsm_v(l)%mom_flux_w(m) * ddy
421          ENDDO   
422       ENDDO
423!
424!--    Urban-type surfaces
425       DO  l = 0, 1
426          surf_s = surf_usm_v(l)%start_index(j,i)
427          surf_e = surf_usm_v(l)%end_index(j,i)
428          DO  m = surf_s, surf_e
429             k           = surf_usm_v(l)%k(m)
430             tend(k,j,i) = tend(k,j,i) +                                       &
431                                     surf_usm_v(l)%mom_flux_w(m) * ddy
432          ENDDO   
433       ENDDO
434!
435!--    Add horizontal momentum flux u'w' at east- (l=2) and west-facing (l=3)
436!--    surfaces.
437!--    Default-type surfaces
438       DO  l = 2, 3
439          surf_s = surf_def_v(l)%start_index(j,i)
440          surf_e = surf_def_v(l)%end_index(j,i)
441          DO  m = surf_s, surf_e
442             k           = surf_def_v(l)%k(m)
443             tend(k,j,i) = tend(k,j,i) +                                       &
444                                     surf_def_v(l)%mom_flux_w(m) * ddx
445          ENDDO   
446       ENDDO
447!
448!--    Natural-type surfaces
449       DO  l = 2, 3
450          surf_s = surf_lsm_v(l)%start_index(j,i)
451          surf_e = surf_lsm_v(l)%end_index(j,i)
452          DO  m = surf_s, surf_e
453             k           = surf_lsm_v(l)%k(m)
454             tend(k,j,i) = tend(k,j,i) +                                       &
455                                     surf_lsm_v(l)%mom_flux_w(m) * ddx
456          ENDDO   
457       ENDDO
458!
459!--    Urban-type surfaces
460       DO  l = 2, 3
461          surf_s = surf_usm_v(l)%start_index(j,i)
462          surf_e = surf_usm_v(l)%end_index(j,i)
463          DO  m = surf_s, surf_e
464             k           = surf_usm_v(l)%k(m)
465             tend(k,j,i) = tend(k,j,i) +                                       &
466                                     surf_usm_v(l)%mom_flux_w(m) * ddx
467          ENDDO   
468       ENDDO
469
470
471    END SUBROUTINE diffusion_w_ij
472
473 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.