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

Last change on this file since 4132 was 3927, checked in by raasch, 5 years ago

pointer attribute removed from scalar 3d-array for performance reasons

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