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

Last change on this file since 4313 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
RevLine 
[1873]1!> @file diffusion_w.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1341]22!
[2233]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: diffusion_w.f90 4182 2019-08-22 15:20:23Z suehring $
[4182]27! Corrected "Former revisions" section
28!
29! 3655 2019-01-07 16:51:22Z knoop
[3634]30! OpenACC port for SPEC
[1321]31!
[4182]32! Revision 1.1  1997/09/12 06:24:11  raasch
33! Initial revision
34!
35!
[1]36! Description:
37! ------------
[1682]38!> Diffusion term of the w-component
[1]39!------------------------------------------------------------------------------!
[1682]40 MODULE diffusion_w_mod
41 
[1]42
43    PRIVATE
[2118]44    PUBLIC diffusion_w
[1]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!------------------------------------------------------------------------------!
[1682]55! Description:
56! ------------
57!> Call for all grid points
[1]58!------------------------------------------------------------------------------!
[1001]59    SUBROUTINE diffusion_w
[1]60
[1320]61       USE arrays_3d,                                                          &         
[2037]62           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]63           
64       USE grid_variables,                                                     &     
[2232]65           ONLY :  ddx, ddy
[1320]66           
67       USE indices,                                                            &           
[2232]68           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[1320]69           
70       USE kinds
[1]71
[2232]72       USE surface_mod,                                                        &
73           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
74
[1]75       IMPLICIT NONE
76
[2232]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
[1320]84       
[2232]85       REAL(wp) ::  flag              !< flag to mask topography grid points
[3547]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
[2232]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
[1001]94
[1]95
96
[3634]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)
[1]107       DO  i = nxl, nxr
108          DO  j = nys, nyn
[2232]109             DO  k = nzb+1, nzt-1
[1]110!
[2232]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!
[1]123!--             Interpolate eddy diffusivities on staggered gridpoints
[3547]124                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +               &
[2232]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) )
[1]132
133                tend(k,j,i) = tend(k,j,i)                                      &
[2232]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
[1]158             ENDDO
159
160!
[2232]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
[1]175!
[2232]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
[1]232
233          ENDDO
234       ENDDO
235
236    END SUBROUTINE diffusion_w
237
238
239!------------------------------------------------------------------------------!
[1682]240! Description:
241! ------------
242!> Call for grid point i,j
[1]243!------------------------------------------------------------------------------!
[1001]244    SUBROUTINE diffusion_w_ij( i, j )
[1]245
[1320]246       USE arrays_3d,                                                          &         
[2037]247           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]248           
249       USE grid_variables,                                                     &     
[2232]250           ONLY :  ddx, ddy
[1320]251           
252       USE indices,                                                            &           
[3241]253           ONLY :  nzb, nzt, wall_flags_0
[1320]254           
255       USE kinds
[1]256
[2232]257       USE surface_mod,                                                        &
258           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
259
[1]260       IMPLICIT NONE
261
[2232]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
[1320]270       
[2232]271       REAL(wp) ::  flag              !< flag to mask topography grid points
[3547]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
[2232]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
[1]280
281
[2232]282       DO  k = nzb+1, nzt-1
[1]283!
[2232]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!
[1]291!--       Interpolate eddy diffusivities on staggered gridpoints
[1322]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) )
[1]296
297          tend(k,j,i) = tend(k,j,i)                                            &
[2232]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
[1]322       ENDDO
323!
[2232]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
[51]338!
[2232]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
[1]349!
[2232]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
[1]395
396
397    END SUBROUTINE diffusion_w_ij
398
[1323]399 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.