source: palm/trunk/SOURCE/diffusion_s.f90

Last change on this file was 4828, checked in by Giersch, 7 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 35.6 KB
Line 
1!> @file diffusion_s.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 terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: diffusion_s.f90 4828 2021-01-05 11:21:41Z banzhafs $
26! Implementation of downward facing USM and LSM surfaces
27!
28! 4583 2020-06-29 12:36:47Z raasch
29! file re-formatted to follow the PALM coding standard
30!
31! 4360 2020-01-07 11:25:50Z suehring
32! Introduction of wall_flags_total_0, which currently sets bits based on static topography
33! information used in wall_flags_static_0
34!
35! 4329 2019-12-10 15:46:36Z motisi
36! Renamed wall_flags_0 to wall_flags_total_0
37!
38! 4182 2019-08-22 15:20:23Z scharf
39! Corrected "Former revisions" section
40!
41! 3927 2019-04-23 13:24:29Z raasch
42! pointer attribute removed from scalar 3d-array for performance reasons
43!
44! 3761 2019-02-25 15:31:42Z raasch
45! unused variables removed
46!
47! 3655 2019-01-07 16:51:22Z knoop
48! nopointer option removed
49!
50! Revision 1.1  2000/04/13 14:54:02  schroeter
51! Initial revision
52!
53!
54! Description:
55! ------------
56!> Diffusion term of scalar quantities (temperature and water content)
57!--------------------------------------------------------------------------------------------------!
58 MODULE diffusion_s_mod
59
60
61    PRIVATE
62    PUBLIC diffusion_s
63
64    INTERFACE diffusion_s
65       MODULE PROCEDURE diffusion_s
66       MODULE PROCEDURE diffusion_s_ij
67    END INTERFACE diffusion_s
68
69 CONTAINS
70
71
72!--------------------------------------------------------------------------------------------------!
73! Description:
74! ------------
75!> Call for all grid points
76!--------------------------------------------------------------------------------------------------!
77    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,                              &
78                               s_flux_t,                                                           &
79                               s_flux_lsm_h_up,    s_flux_lsm_h_down,                              &
80                               s_flux_usm_h_up,    s_flux_usm_h_down,                              &
81                               s_flux_def_v_north, s_flux_def_v_south,                             &
82                               s_flux_def_v_east,  s_flux_def_v_west,                              &
83                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
84                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
85                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
86                               s_flux_usm_v_east,  s_flux_usm_v_west )
87
88       USE arrays_3d,                                                                              &
89           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
90
91       USE control_parameters,                                                                     &
92           ONLY: use_surface_fluxes, use_top_fluxes
93
94       USE grid_variables,                                                                         &
95           ONLY:  ddx, ddx2, ddy, ddy2
96
97       USE indices,                                                                                &
98           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
99
100       USE kinds
101
102       USE surface_mod,                                                                            &
103           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
104
105       IMPLICIT NONE
106
107       INTEGER(iwp) ::  i             !< running index x direction
108       INTEGER(iwp) ::  j             !< running index y direction
109       INTEGER(iwp) ::  k             !< running index z direction
110       INTEGER(iwp) ::  m             !< running index surface elements
111       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
112       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
113
114       REAL(wp) ::  flag              !< flag to mask topography grid points
115       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
116       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
117       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
118       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
119       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
120       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
121
122       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
123       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
124       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
125       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
126       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
127       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
128       REAL(wp), DIMENSION(1:surf_lsm_h(0)%ns) ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
129       REAL(wp), DIMENSION(1:surf_lsm_h(1)%ns) ::  s_flux_lsm_h_down  !< flux at horizontal downward-facing natural-type surfaces
130       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
131       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
132       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
133       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
134       REAL(wp), DIMENSION(1:surf_usm_h(0)%ns) ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
135       REAL(wp), DIMENSION(1:surf_usm_h(1)%ns) ::  s_flux_usm_h_down  !< flux at horizontal downward-facing urban-type surfaces
136       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
137       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
138       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
139       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
140       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
141
142       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< treated scalar
143
144
145       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
146       !$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
147       !$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
148       !$ACC PRESENT(wall_flags_total_0, kh) &
149       !$ACC PRESENT(s) &
150       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
151       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
152       !$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(0:3)) &
153       !$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(0:3)) &
154       !$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
155       !$ACC PRESENT(s_flux_t) &
156       !$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
157       !$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
158       !$ACC PRESENT(s_flux_lsm_h_up, s_flux_lsm_h_down) &
159       !$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
160       !$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
161       !$ACC PRESENT(s_flux_usm_h_up, s_flux_usm_h_down) &
162       !$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
163       !$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
164       !$ACC PRESENT(tend)
165       DO  i = nxl, nxr
166          DO  j = nys,nyn
167!
168!--          Compute horizontal diffusion
169             DO  k = nzb+1, nzt
170!
171!--             Predetermine flag to mask topography and wall-bounded grid points
172                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
173!
174!--             Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer
175!--             array
176                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
177                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
178                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 0 ) )
179                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
180
181                tend(k,j,i) = tend(k,j,i)                                                          &
182                                                  + 0.5_wp * (                                     &
183                                mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                           &
184                                           * ( s(k,j,i+1) - s(k,j,i)   )                           &
185                              - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                           &
186                                           * ( s(k,j,i)   - s(k,j,i-1) )                           &
187                                                             ) * ddx2 * flag                       &
188                                                  + 0.5_wp * (                                     &
189                                mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                           &
190                                           * ( s(k,j+1,i) - s(k,j,i)   )                           &
191                              - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                           &
192                                           * ( s(k,j,i)   - s(k,j-1,i) )                           &
193                                                             ) * ddy2 * flag
194             ENDDO
195
196!
197!--          Apply prescribed horizontal wall heatflux where necessary. First, determine start and
198!--          end index for respective (j,i)-index. Please note, in the flat case following loop will
199!--          not be entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural
200!--          surfaces so far.
201!--          First, for default-type surfaces.
202!--          North-facing vertical default-type surfaces
203             surf_s = surf_def_v(0)%start_index(j,i)
204             surf_e = surf_def_v(0)%end_index(j,i)
205             DO  m = surf_s, surf_e
206                k           = surf_def_v(0)%k(m)
207                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
208             ENDDO
209!
210!--          South-facing vertical default-type surfaces
211             surf_s = surf_def_v(1)%start_index(j,i)
212             surf_e = surf_def_v(1)%end_index(j,i)
213             DO  m = surf_s, surf_e
214                k           = surf_def_v(1)%k(m)
215                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
216             ENDDO
217!
218!--          East-facing vertical default-type surfaces
219             surf_s = surf_def_v(2)%start_index(j,i)
220             surf_e = surf_def_v(2)%end_index(j,i)
221             DO  m = surf_s, surf_e
222                k           = surf_def_v(2)%k(m)
223                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
224             ENDDO
225!
226!--          West-facing vertical default-type surfaces
227             surf_s = surf_def_v(3)%start_index(j,i)
228             surf_e = surf_def_v(3)%end_index(j,i)
229             DO  m = surf_s, surf_e
230                k           = surf_def_v(3)%k(m)
231                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
232             ENDDO
233!
234!--          Now, for natural-type surfaces.
235!--          North-facing
236             surf_s = surf_lsm_v(0)%start_index(j,i)
237             surf_e = surf_lsm_v(0)%end_index(j,i)
238             DO  m = surf_s, surf_e
239                k           = surf_lsm_v(0)%k(m)
240                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
241             ENDDO
242!
243!--          South-facing
244             surf_s = surf_lsm_v(1)%start_index(j,i)
245             surf_e = surf_lsm_v(1)%end_index(j,i)
246             DO  m = surf_s, surf_e
247                k           = surf_lsm_v(1)%k(m)
248                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
249             ENDDO
250!
251!--          East-facing
252             surf_s = surf_lsm_v(2)%start_index(j,i)
253             surf_e = surf_lsm_v(2)%end_index(j,i)
254             DO  m = surf_s, surf_e
255                k           = surf_lsm_v(2)%k(m)
256                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
257             ENDDO
258!
259!--          West-facing
260             surf_s = surf_lsm_v(3)%start_index(j,i)
261             surf_e = surf_lsm_v(3)%end_index(j,i)
262             DO  m = surf_s, surf_e
263                k           = surf_lsm_v(3)%k(m)
264                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
265             ENDDO
266!
267!--          Now, for urban-type surfaces.
268!--          North-facing
269             surf_s = surf_usm_v(0)%start_index(j,i)
270             surf_e = surf_usm_v(0)%end_index(j,i)
271             DO  m = surf_s, surf_e
272                k           = surf_usm_v(0)%k(m)
273                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
274             ENDDO
275!
276!--          South-facing
277             surf_s = surf_usm_v(1)%start_index(j,i)
278             surf_e = surf_usm_v(1)%end_index(j,i)
279             DO  m = surf_s, surf_e
280                k           = surf_usm_v(1)%k(m)
281                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
282             ENDDO
283!
284!--          East-facing
285             surf_s = surf_usm_v(2)%start_index(j,i)
286             surf_e = surf_usm_v(2)%end_index(j,i)
287             DO  m = surf_s, surf_e
288                k           = surf_usm_v(2)%k(m)
289                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
290             ENDDO
291!
292!--          West-facing
293             surf_s = surf_usm_v(3)%start_index(j,i)
294             surf_e = surf_usm_v(3)%end_index(j,i)
295             DO  m = surf_s, surf_e
296                k           = surf_usm_v(3)%k(m)
297                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
298             ENDDO
299
300!
301!--          Compute vertical diffusion. In case that surface fluxes have been prescribed or
302!--          computed at bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively.
303!--          Model top is also mask if top flux is given.
304             DO  k = nzb+1, nzt
305!
306!--             Determine flags to mask topography below and above. Flag 0 is used to mask
307!--             topography in general, and flag 8 implies information about use_surface_fluxes.
308!--             Flag 9 is used to control flux at model top.
309                mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
310                mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) *   &
311                              MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
312                flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
313
314                tend(k,j,i) = tend(k,j,i)                                                          &
315                                       + 0.5_wp * (                                                &
316                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
317                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
318                                                            * rho_air_zw(k)                        &
319                                                            * mask_top                             &
320                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
321                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
322                                                            * rho_air_zw(k-1)                      &
323                                                            * mask_bottom                          &
324                                                  ) * ddzw(k) * drho_air(k)                        &
325                                                              * flag
326             ENDDO
327
328!
329!--          Vertical diffusion at horizontal walls.
330             IF ( use_surface_fluxes )  THEN
331!
332!--             Default-type surfaces, upward-facing
333                surf_s = surf_def_h(0)%start_index(j,i)
334                surf_e = surf_def_h(0)%end_index(j,i)
335                DO  m = surf_s, surf_e
336                   k   = surf_def_h(0)%k(m)
337                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
338                ENDDO
339!
340!--             Default-type surfaces, downward-facing
341                surf_s = surf_def_h(1)%start_index(j,i)
342                surf_e = surf_def_h(1)%end_index(j,i)
343                DO  m = surf_s, surf_e
344                   k   = surf_def_h(1)%k(m)
345                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
346                ENDDO
347!
348!--             Natural-type surfaces, upward-facing
349                surf_s = surf_lsm_h(0)%start_index(j,i)
350                surf_e = surf_lsm_h(0)%end_index(j,i)
351                DO  m = surf_s, surf_e
352                   k   = surf_lsm_h(0)%k(m)
353                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
354                ENDDO
355!
356!--             Natural-type surfaces, downward-facing
357                surf_s = surf_lsm_h(1)%start_index(j,i)
358                surf_e = surf_lsm_h(1)%end_index(j,i)
359                DO  m = surf_s, surf_e
360                   k   = surf_lsm_h(1)%k(m)
361                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_down(m) * ddzw(k) * drho_air(k)
362                ENDDO
363!
364!--             Urban-type surfaces, upward-facing
365                surf_s = surf_usm_h(0)%start_index(j,i)
366                surf_e = surf_usm_h(0)%end_index(j,i)
367                DO  m = surf_s, surf_e
368                   k   = surf_usm_h(0)%k(m)
369                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
370                ENDDO
371!
372!--             Urban-type surfaces, downward-facing
373                surf_s = surf_usm_h(1)%start_index(j,i)
374                surf_e = surf_usm_h(1)%end_index(j,i)
375                DO  m = surf_s, surf_e
376                   k   = surf_usm_h(1)%k(m)
377                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_down(m) * ddzw(k) * drho_air(k)
378                ENDDO
379
380             ENDIF
381!
382!--          Vertical diffusion at the last computational gridpoint along z-direction
383             IF ( use_top_fluxes )  THEN
384                surf_s = surf_def_h(2)%start_index(j,i)
385                surf_e = surf_def_h(2)%end_index(j,i)
386                DO  m = surf_s, surf_e
387
388                   k   = surf_def_h(2)%k(m)
389                   tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
390                ENDDO
391             ENDIF
392
393          ENDDO
394       ENDDO
395
396    END SUBROUTINE diffusion_s
397
398!--------------------------------------------------------------------------------------------------!
399! Description:
400! ------------
401!> Call for grid point i,j
402!--------------------------------------------------------------------------------------------------!
403    SUBROUTINE diffusion_s_ij( i, j, s,                                                            &
404                               s_flux_def_h_up,    s_flux_def_h_down,                              &
405                               s_flux_t,                                                           &
406                               s_flux_lsm_h_up,    s_flux_lsm_h_down,                              &
407                               s_flux_usm_h_up,    s_flux_usm_h_down,                              &
408                               s_flux_def_v_north, s_flux_def_v_south,                             &
409                               s_flux_def_v_east,  s_flux_def_v_west,                              &
410                               s_flux_lsm_v_north, s_flux_lsm_v_south,                             &
411                               s_flux_lsm_v_east,  s_flux_lsm_v_west,                              &
412                               s_flux_usm_v_north, s_flux_usm_v_south,                             &
413                               s_flux_usm_v_east,  s_flux_usm_v_west )
414
415       USE arrays_3d,                                                                              &
416           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
417
418       USE control_parameters,                                                                     &
419           ONLY: use_surface_fluxes, use_top_fluxes
420
421       USE grid_variables,                                                                         &
422           ONLY:  ddx, ddx2, ddy, ddy2
423
424       USE indices,                                                                                &
425           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_total_0
426
427       USE kinds
428
429       USE surface_mod,                                                                            &
430           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
431
432       IMPLICIT NONE
433
434       INTEGER(iwp) ::  i             !< running index x direction
435       INTEGER(iwp) ::  j             !< running index y direction
436       INTEGER(iwp) ::  k             !< running index z direction
437       INTEGER(iwp) ::  m             !< running index surface elements
438       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
439       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
440
441       REAL(wp) ::  flag              !< flag to mask topography grid points
442       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface
443       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
444       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
445       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
446       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface
447       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
448
449       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
450       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
451       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
452       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
453       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
454       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
455       REAL(wp), DIMENSION(1:surf_lsm_h(0)%ns) ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
456       REAL(wp), DIMENSION(1:surf_lsm_h(1)%ns) ::  s_flux_lsm_h_down    !< flux at horizontal downward-facing natural-type surfaces
457       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
458       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
459       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
460       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
461       REAL(wp), DIMENSION(1:surf_usm_h(0)%ns) ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
462       REAL(wp), DIMENSION(1:surf_usm_h(1)%ns) ::  s_flux_usm_h_down  !< flux at horizontal downward-facing urban-type surfaces
463       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
464       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
465       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
466       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
467       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
468
469       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< treated scalar
470
471
472!
473!--    Compute horizontal diffusion
474       DO  k = nzb+1, nzt
475!
476!--       Predetermine flag to mask topography and wall-bounded grid points
477          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
478!
479!--       Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer array
480          mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
481          mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
482          mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 0 ) )
483          mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
484!
485!--       Finally, determine flag to mask both topography itself as well as wall-bounded grid
486!--       points, which will be treated further below
487
488          tend(k,j,i) = tend(k,j,i)                                                                &
489                                          + 0.5_wp * (                                             &
490                            mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )                               &
491                                       * ( s(k,j,i+1) - s(k,j,i)   )                               &
492                          - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )                               &
493                                       * ( s(k,j,i)   - s(k,j,i-1) )                               &
494                                                     ) * ddx2 * flag                               &
495                                          + 0.5_wp * (                                             &
496                            mask_north * ( kh(k,j,i) + kh(k,j+1,i) )                               &
497                                       * ( s(k,j+1,i) - s(k,j,i)   )                               &
498                          - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )                               &
499                                       * ( s(k,j,i)  - s(k,j-1,i)  )                               &
500                                                     ) * ddy2 * flag
501       ENDDO
502
503!
504!--    Apply prescribed horizontal wall heatflux where necessary. First, determine start and end
505!--    index for respective (j,i)-index. Please note, in the flat case following loops will not be
506!--    entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces so far.
507!--    First, for default-type surfaces.
508!--    North-facing vertical default-type surfaces
509       surf_s = surf_def_v(0)%start_index(j,i)
510       surf_e = surf_def_v(0)%end_index(j,i)
511       DO  m = surf_s, surf_e
512          k           = surf_def_v(0)%k(m)
513          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
514       ENDDO
515!
516!--    South-facing vertical default-type surfaces
517       surf_s = surf_def_v(1)%start_index(j,i)
518       surf_e = surf_def_v(1)%end_index(j,i)
519       DO  m = surf_s, surf_e
520          k           = surf_def_v(1)%k(m)
521          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
522       ENDDO
523!
524!--    East-facing vertical default-type surfaces
525       surf_s = surf_def_v(2)%start_index(j,i)
526       surf_e = surf_def_v(2)%end_index(j,i)
527       DO  m = surf_s, surf_e
528          k           = surf_def_v(2)%k(m)
529          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
530       ENDDO
531!
532!--    West-facing vertical default-type surfaces
533       surf_s = surf_def_v(3)%start_index(j,i)
534       surf_e = surf_def_v(3)%end_index(j,i)
535       DO  m = surf_s, surf_e
536          k           = surf_def_v(3)%k(m)
537          tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
538       ENDDO
539!
540!--    Now, for natural-type surfaces
541!--    North-facing
542       surf_s = surf_lsm_v(0)%start_index(j,i)
543       surf_e = surf_lsm_v(0)%end_index(j,i)
544       DO  m = surf_s, surf_e
545          k           = surf_lsm_v(0)%k(m)
546          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
547       ENDDO
548!
549!--    South-facing
550       surf_s = surf_lsm_v(1)%start_index(j,i)
551       surf_e = surf_lsm_v(1)%end_index(j,i)
552       DO  m = surf_s, surf_e
553          k           = surf_lsm_v(1)%k(m)
554          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
555       ENDDO
556!
557!--    East-facing
558       surf_s = surf_lsm_v(2)%start_index(j,i)
559       surf_e = surf_lsm_v(2)%end_index(j,i)
560       DO  m = surf_s, surf_e
561          k           = surf_lsm_v(2)%k(m)
562          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
563       ENDDO
564!
565!--    West-facing
566       surf_s = surf_lsm_v(3)%start_index(j,i)
567       surf_e = surf_lsm_v(3)%end_index(j,i)
568       DO  m = surf_s, surf_e
569          k           = surf_lsm_v(3)%k(m)
570          tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
571       ENDDO
572!
573!--    Now, for urban-type surfaces
574!--    North-facing
575       surf_s = surf_usm_v(0)%start_index(j,i)
576       surf_e = surf_usm_v(0)%end_index(j,i)
577       DO  m = surf_s, surf_e
578          k           = surf_usm_v(0)%k(m)
579          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
580       ENDDO
581!
582!--    South-facing
583       surf_s = surf_usm_v(1)%start_index(j,i)
584       surf_e = surf_usm_v(1)%end_index(j,i)
585       DO  m = surf_s, surf_e
586          k           = surf_usm_v(1)%k(m)
587          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
588       ENDDO
589!
590!--    East-facing
591       surf_s = surf_usm_v(2)%start_index(j,i)
592       surf_e = surf_usm_v(2)%end_index(j,i)
593       DO  m = surf_s, surf_e
594          k           = surf_usm_v(2)%k(m)
595          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
596       ENDDO
597!
598!--    West-facing
599       surf_s = surf_usm_v(3)%start_index(j,i)
600       surf_e = surf_usm_v(3)%end_index(j,i)
601       DO  m = surf_s, surf_e
602          k           = surf_usm_v(3)%k(m)
603          tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
604       ENDDO
605
606
607!
608!--    Compute vertical diffusion. In case that surface fluxes have been prescribed or computed at
609!--    bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively. Model top is also
610!--    mask if top flux is given.
611       DO  k = nzb+1, nzt
612!
613!--       Determine flags to mask topography below and above. Flag 0 is used to mask topography in
614!--       general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
615!--       control flux at model top.
616          mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
617          mask_top    = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) )  *        &
618                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
619          flag        = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
620
621          tend(k,j,i) = tend(k,j,i)                                                                &
622                                       + 0.5_wp * (                                                &
623                                      ( kh(k,j,i) + kh(k+1,j,i) ) *                                &
624                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)                      &
625                                                            * rho_air_zw(k)                        &
626                                                            * mask_top                             &
627                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *                                &
628                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)                        &
629                                                            * rho_air_zw(k-1)                      &
630                                                            * mask_bottom                          &
631                                                  ) * ddzw(k) * drho_air(k)                        &
632                                                              * flag
633       ENDDO
634
635!
636!--    Vertical diffusion at horizontal walls.
637!--    TO DO: Adjust for downward facing walls and mask already in main loop
638       IF ( use_surface_fluxes )  THEN
639!
640!--       Default-type surfaces, upward-facing
641          surf_s = surf_def_h(0)%start_index(j,i)
642          surf_e = surf_def_h(0)%end_index(j,i)
643          DO  m = surf_s, surf_e
644             k   = surf_def_h(0)%k(m)
645             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
646          ENDDO
647!
648!--       Default-type surfaces, downward-facing
649          surf_s = surf_def_h(1)%start_index(j,i)
650          surf_e = surf_def_h(1)%end_index(j,i)
651          DO  m = surf_s, surf_e
652             k   = surf_def_h(1)%k(m)
653             tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
654          ENDDO
655!
656!--       Natural-type surfaces, upward-facing
657          surf_s = surf_lsm_h(0)%start_index(j,i)
658          surf_e = surf_lsm_h(0)%end_index(j,i)
659          DO  m = surf_s, surf_e
660             k   = surf_lsm_h(0)%k(m)
661             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
662          ENDDO
663!
664!--       Natural-type surfaces, downward-facing
665          surf_s = surf_lsm_h(1)%start_index(j,i)
666          surf_e = surf_lsm_h(1)%end_index(j,i)
667          DO  m = surf_s, surf_e
668             k   = surf_lsm_h(1)%k(m)
669             tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_down(m) * ddzw(k) * drho_air(k)
670          ENDDO
671!
672!--       Urban-type surfaces, upward-facing
673          surf_s = surf_usm_h(0)%start_index(j,i)
674          surf_e = surf_usm_h(0)%end_index(j,i)
675          DO  m = surf_s, surf_e
676             k   = surf_usm_h(0)%k(m)
677             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
678          ENDDO
679!
680!--       Urban-type surfaces, upward-facing
681          surf_s = surf_usm_h(1)%start_index(j,i)
682          surf_e = surf_usm_h(1)%end_index(j,i)
683          DO  m = surf_s, surf_e
684             k   = surf_usm_h(1)%k(m)
685             tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_down(m) * ddzw(k) * drho_air(k)
686          ENDDO
687       ENDIF
688!
689!--    Vertical diffusion at the last computational gridpoint along z-direction
690       IF ( use_top_fluxes )  THEN
691          surf_s = surf_def_h(2)%start_index(j,i)
692          surf_e = surf_def_h(2)%end_index(j,i)
693          DO  m = surf_s, surf_e
694             k   = surf_def_h(2)%k(m)
695             tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
696          ENDDO
697       ENDIF
698
699    END SUBROUTINE diffusion_s_ij
700
701 END MODULE diffusion_s_mod
Note: See TracBrowser for help on using the repository browser.