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

Last change on this file since 4458 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
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!
[4360]17! Copyright 1997-2020 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 4360 2020-01-07 11:25:50Z raasch $
[4346]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
[4329]31! Renamed wall_flags_0 to wall_flags_static_0
32!
33! 4182 2019-08-22 15:20:23Z scharf
[4182]34! Corrected "Former revisions" section
35!
36! 3655 2019-01-07 16:51:22Z knoop
[3634]37! OpenACC port for SPEC
[1321]38!
[4182]39! Revision 1.1  1997/09/12 06:24:11  raasch
40! Initial revision
41!
42!
[1]43! Description:
44! ------------
[1682]45!> Diffusion term of the w-component
[1]46!------------------------------------------------------------------------------!
[1682]47 MODULE diffusion_w_mod
48 
[1]49
50    PRIVATE
[2118]51    PUBLIC diffusion_w
[1]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!------------------------------------------------------------------------------!
[1682]62! Description:
63! ------------
64!> Call for all grid points
[1]65!------------------------------------------------------------------------------!
[1001]66    SUBROUTINE diffusion_w
[1]67
[1320]68       USE arrays_3d,                                                          &         
[2037]69           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]70           
71       USE grid_variables,                                                     &     
[2232]72           ONLY :  ddx, ddy
[1320]73           
74       USE indices,                                                            &           
[4346]75           ONLY :  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
[1320]76           
77       USE kinds
[1]78
[2232]79       USE surface_mod,                                                        &
80           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
81
[1]82       IMPLICIT NONE
83
[2232]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
[1320]91       
[2232]92       REAL(wp) ::  flag              !< flag to mask topography grid points
[3547]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
[2232]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
[1001]101
[1]102
103
[3634]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) &
[4346]107       !$ACC PRESENT(wall_flags_total_0, km) &
[3634]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)
[1]114       DO  i = nxl, nxr
115          DO  j = nys, nyn
[2232]116             DO  k = nzb+1, nzt-1
[1]117!
[2232]118!--             Predetermine flag to mask topography and wall-bounded grid points.
119                flag       = MERGE( 1.0_wp, 0.0_wp,                            &
[4346]120                                    BTEST( wall_flags_total_0(k,j,i),   3 ) ) 
[2232]121                mask_east  = MERGE( 1.0_wp, 0.0_wp,                            &
[4346]122                                    BTEST( wall_flags_total_0(k,j,i+1), 3 ) )
[2232]123                mask_west  = MERGE( 1.0_wp, 0.0_wp,                            &
[4346]124                                    BTEST( wall_flags_total_0(k,j,i-1), 3 ) )
[2232]125                mask_south = MERGE( 1.0_wp, 0.0_wp,                            &
[4346]126                                    BTEST( wall_flags_total_0(k,j-1,i), 3 ) )
[2232]127                mask_north = MERGE( 1.0_wp, 0.0_wp,                            &
[4346]128                                    BTEST( wall_flags_total_0(k,j+1,i), 3 ) )
[2232]129!
[1]130!--             Interpolate eddy diffusivities on staggered gridpoints
[3547]131                kmxp = 0.25_wp * ( km(k,j,i)   + km(k,j,i+1)   +               &
[2232]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) )
[1]139
140                tend(k,j,i) = tend(k,j,i)                                      &
[2232]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
[1]165             ENDDO
166
167!
[2232]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
[1]182!
[2232]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
[1]239
240          ENDDO
241       ENDDO
242
243    END SUBROUTINE diffusion_w
244
245
246!------------------------------------------------------------------------------!
[1682]247! Description:
248! ------------
249!> Call for grid point i,j
[1]250!------------------------------------------------------------------------------!
[1001]251    SUBROUTINE diffusion_w_ij( i, j )
[1]252
[1320]253       USE arrays_3d,                                                          &         
[2037]254           ONLY :  ddzu, ddzw, km, tend, u, v, w, drho_air_zw, rho_air
[1320]255           
256       USE grid_variables,                                                     &     
[2232]257           ONLY :  ddx, ddy
[1320]258           
259       USE indices,                                                            &           
[4346]260           ONLY :  nzb, nzt, wall_flags_total_0
[1320]261           
262       USE kinds
[1]263
[2232]264       USE surface_mod,                                                        &
265           ONLY :  surf_def_v, surf_lsm_v, surf_usm_v
266
[1]267       IMPLICIT NONE
268
[2232]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
[1320]277       
[2232]278       REAL(wp) ::  flag              !< flag to mask topography grid points
[3547]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
[2232]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
[1]287
288
[2232]289       DO  k = nzb+1, nzt-1
[1]290!
[2232]291!--       Predetermine flag to mask topography and wall-bounded grid points.
[4346]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 ) )
[2232]297!
[1]298!--       Interpolate eddy diffusivities on staggered gridpoints
[1322]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) )
[1]303
304          tend(k,j,i) = tend(k,j,i)                                            &
[2232]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
[1]329       ENDDO
330!
[2232]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
[51]345!
[2232]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
[1]356!
[2232]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
[1]402
403
404    END SUBROUTINE diffusion_w_ij
405
[1323]406 END MODULE diffusion_w_mod
Note: See TracBrowser for help on using the repository browser.