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

Last change on this file since 4360 was 4360, checked in by suehring, 16 months ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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