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

Last change on this file since 3655 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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