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

Last change on this file since 3854 was 3761, checked in by raasch, 5 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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