source: palm/trunk/SOURCE/diffusion_s.f90 @ 4648

Last change on this file since 4648 was 4583, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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