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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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