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

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

OpenACC port for SPEC

  • 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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: diffusion_s.f90 3634 2018-12-18 12:31:28Z knoop $
27! OpenACC port for SPEC
28!
29! 3547 2018-11-21 13:21:24Z suehring
30! variables documented
31!
32! 2759 2018-01-17 16:24:59Z suehring
33! Major bugfix, horizontal diffusion at vertical surfaces corrected.
34!
35! 2718 2018-01-02 08:49:38Z maronga
36! Corrected "Former revisions" section
37!
38! 2696 2017-12-14 17:12:51Z kanani
39! Change in file header (GPL part)
40!
41! 2233 2017-05-30 18:08:54Z suehring
42!
43! 2232 2017-05-30 17:47:52Z suehring
44! Adjustments to new topography and surface concept
45!
46! 2118 2017-01-17 16:38:49Z raasch
47! OpenACC version of subroutine removed
48!
49! 2037 2016-10-26 11:15:40Z knoop
50! Anelastic approximation implemented
51!
52! 2000 2016-08-20 18:09:15Z knoop
53! Forced header and separation lines into 80 columns
54!
55! 1873 2016-04-18 14:50:06Z maronga
56! Module renamed (removed _mod)
57!
58! 1850 2016-04-08 13:29:27Z maronga
59! Module renamed
60!
61! 1691 2015-10-26 16:17:44Z maronga
62! Formatting corrections.
63!
64! 1682 2015-10-07 23:56:08Z knoop
65! Code annotations made doxygen readable
66!
67! 1374 2014-04-25 12:55:07Z raasch
68! missing variables added to ONLY list
69!
70! 1340 2014-03-25 19:45:13Z kanani
71! REAL constants defined as wp-kind
72!
73! 1320 2014-03-20 08:40:49Z raasch
74! ONLY-attribute added to USE-statements,
75! kind-parameters added to all INTEGER and REAL declaration statements,
76! kinds are defined in new module kinds,
77! revision history before 2012 removed,
78! comment fields (!:) to be used for variable explanations added to
79! all variable declaration statements
80!
81! 1257 2013-11-08 15:18:40Z raasch
82! openacc loop and loop vector clauses removed
83!
84! 1128 2013-04-12 06:19:32Z raasch
85! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
86! j_north
87!
88! 1092 2013-02-02 11:24:22Z raasch
89! unused variables removed
90!
91! 1036 2012-10-22 13:43:42Z raasch
92! code put under GPL (PALM 3.9)
93!
94! 1015 2012-09-27 09:23:24Z raasch
95! accelerator version (*_acc) added
96!
97! 1010 2012-09-20 07:59:54Z raasch
98! cpp switch __nopointer added for pointer free version
99!
100! 1001 2012-09-13 14:08:46Z raasch
101! some arrays comunicated by module instead of parameter list
102!
103! Revision 1.1  2000/04/13 14:54:02  schroeter
104! Initial revision
105!
106!
107! Description:
108! ------------
109!> Diffusion term of scalar quantities (temperature and water content)
110!------------------------------------------------------------------------------!
111 MODULE diffusion_s_mod
112 
113
114    PRIVATE
115    PUBLIC diffusion_s
116
117    INTERFACE diffusion_s
118       MODULE PROCEDURE diffusion_s
119       MODULE PROCEDURE diffusion_s_ij
120    END INTERFACE diffusion_s
121
122 CONTAINS
123
124
125!------------------------------------------------------------------------------!
126! Description:
127! ------------
128!> Call for all grid points
129!------------------------------------------------------------------------------!
130    SUBROUTINE diffusion_s( s, s_flux_def_h_up,    s_flux_def_h_down,          &
131                               s_flux_t,                                       &
132                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
133                               s_flux_def_v_north, s_flux_def_v_south,         &
134                               s_flux_def_v_east,  s_flux_def_v_west,          &
135                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
136                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
137                               s_flux_usm_v_north, s_flux_usm_v_south,         &
138                               s_flux_usm_v_east,  s_flux_usm_v_west )
139
140       USE arrays_3d,                                                          &
141           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
142       
143       USE control_parameters,                                                 & 
144           ONLY: use_surface_fluxes, use_top_fluxes
145       
146       USE grid_variables,                                                     &
147           ONLY:  ddx, ddx2, ddy, ddy2
148       
149       USE indices,                                                            &
150           ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb,             &
151                  nzt, wall_flags_0
152       
153       USE kinds
154
155       USE surface_mod,                                                        &
156           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
157                   surf_usm_v
158
159       IMPLICIT NONE
160
161       INTEGER(iwp) ::  i             !< running index x direction
162       INTEGER(iwp) ::  j             !< running index y direction
163       INTEGER(iwp) ::  k             !< running index z direction
164       INTEGER(iwp) ::  m             !< running index surface elements
165       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
166       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
167
168       REAL(wp) ::  flag              !< flag to mask topography grid points
169       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
170       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
171       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
172       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
173       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
174       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
175
176       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
177       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
178       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
179       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
180       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
181       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
182       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
183       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
184       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
185       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical natural-type surfaces
186       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical natural-type surfaces
187       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
188       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
189       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
190       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
191       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
192       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
193#if defined( __nopointer )
194       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s  !< treated scalar
195#else
196       REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !< treated scalar
197#endif
198
199       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
200       !$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
201       !$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
202       !$ACC PRESENT(wall_flags_0, kh) &
203       !$ACC PRESENT(s) &
204       !$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
205       !$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
206       !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
207       !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
208       !$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
209       !$ACC PRESENT(s_flux_t) &
210       !$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
211       !$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
212       !$ACC PRESENT(s_flux_lsm_h_up) &
213       !$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
214       !$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
215       !$ACC PRESENT(s_flux_usm_h_up) &
216       !$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
217       !$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
218       !$ACC PRESENT(tend)
219       DO  i = nxl, nxr
220          DO  j = nys,nyn
221!
222!--          Compute horizontal diffusion
223             DO  k = nzb+1, nzt
224!
225!--             Predetermine flag to mask topography and wall-bounded grid points
226                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 
227!
228!--             Predetermine flag to mask wall-bounded grid points, equivalent to
229!--             former s_outer array
230                mask_west  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
231                mask_east  = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
232                mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
233                mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
234
235                tend(k,j,i) = tend(k,j,i)                                      &
236                                          + 0.5_wp * (                         &
237                        mask_east  * ( kh(k,j,i) + kh(k,j,i+1) )               &
238                                   * ( s(k,j,i+1) - s(k,j,i)   )               &
239                      - mask_west  * ( kh(k,j,i) + kh(k,j,i-1) )               &
240                                   * ( s(k,j,i)   - s(k,j,i-1) )               &
241                                                     ) * ddx2 * flag           &
242                                          + 0.5_wp * (                         &
243                        mask_north * ( kh(k,j,i) + kh(k,j+1,i) )               &
244                                   * ( s(k,j+1,i) - s(k,j,i)   )               &
245                      - mask_south * ( kh(k,j,i) + kh(k,j-1,i) )               &
246                                   * ( s(k,j,i)   - s(k,j-1,i) )               &
247                                                     ) * ddy2 * flag
248             ENDDO
249
250!
251!--          Apply prescribed horizontal wall heatflux where necessary. First,
252!--          determine start and end index for respective (j,i)-index. Please
253!--          note, in the flat case following loop will not be entered, as
254!--          surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
255!--          so far.
256!--          First, for default-type surfaces
257!--          North-facing vertical default-type surfaces
258             surf_s = surf_def_v(0)%start_index(j,i)
259             surf_e = surf_def_v(0)%end_index(j,i)
260             DO  m = surf_s, surf_e
261                k           = surf_def_v(0)%k(m)
262                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
263             ENDDO
264!
265!--          South-facing vertical default-type surfaces
266             surf_s = surf_def_v(1)%start_index(j,i)
267             surf_e = surf_def_v(1)%end_index(j,i)
268             DO  m = surf_s, surf_e
269                k           = surf_def_v(1)%k(m)
270                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
271             ENDDO
272!
273!--          East-facing vertical default-type surfaces
274             surf_s = surf_def_v(2)%start_index(j,i)
275             surf_e = surf_def_v(2)%end_index(j,i)
276             DO  m = surf_s, surf_e
277                k           = surf_def_v(2)%k(m)
278                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
279             ENDDO
280!
281!--          West-facing vertical default-type surfaces
282             surf_s = surf_def_v(3)%start_index(j,i)
283             surf_e = surf_def_v(3)%end_index(j,i)
284             DO  m = surf_s, surf_e
285                k           = surf_def_v(3)%k(m)
286                tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
287             ENDDO
288!
289!--          Now, for natural-type surfaces.
290!--          North-facing
291             surf_s = surf_lsm_v(0)%start_index(j,i)
292             surf_e = surf_lsm_v(0)%end_index(j,i)
293             DO  m = surf_s, surf_e
294                k           = surf_lsm_v(0)%k(m)
295                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
296             ENDDO
297!
298!--          South-facing
299             surf_s = surf_lsm_v(1)%start_index(j,i)
300             surf_e = surf_lsm_v(1)%end_index(j,i)
301             DO  m = surf_s, surf_e
302                k           = surf_lsm_v(1)%k(m)
303                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
304             ENDDO
305!
306!--          East-facing
307             surf_s = surf_lsm_v(2)%start_index(j,i)
308             surf_e = surf_lsm_v(2)%end_index(j,i)
309             DO  m = surf_s, surf_e
310                k           = surf_lsm_v(2)%k(m)
311                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
312             ENDDO
313!
314!--          West-facing
315             surf_s = surf_lsm_v(3)%start_index(j,i)
316             surf_e = surf_lsm_v(3)%end_index(j,i)
317             DO  m = surf_s, surf_e
318                k           = surf_lsm_v(3)%k(m)
319                tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
320             ENDDO
321!
322!--          Now, for urban-type surfaces.
323!--          North-facing
324             surf_s = surf_usm_v(0)%start_index(j,i)
325             surf_e = surf_usm_v(0)%end_index(j,i)
326             DO  m = surf_s, surf_e
327                k           = surf_usm_v(0)%k(m)
328                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
329             ENDDO
330!
331!--          South-facing
332             surf_s = surf_usm_v(1)%start_index(j,i)
333             surf_e = surf_usm_v(1)%end_index(j,i)
334             DO  m = surf_s, surf_e
335                k           = surf_usm_v(1)%k(m)
336                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
337             ENDDO
338!
339!--          East-facing
340             surf_s = surf_usm_v(2)%start_index(j,i)
341             surf_e = surf_usm_v(2)%end_index(j,i)
342             DO  m = surf_s, surf_e
343                k           = surf_usm_v(2)%k(m)
344                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
345             ENDDO
346!
347!--          West-facing
348             surf_s = surf_usm_v(3)%start_index(j,i)
349             surf_e = surf_usm_v(3)%end_index(j,i)
350             DO  m = surf_s, surf_e
351                k           = surf_usm_v(3)%k(m)
352                tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
353             ENDDO
354
355!
356!--          Compute vertical diffusion. In case that surface fluxes have been
357!--          prescribed or computed at bottom and/or top, index k starts/ends at
358!--          nzb+2 or nzt-1, respectively. Model top is also mask if top flux
359!--          is given.
360             DO  k = nzb+1, nzt
361!
362!--             Determine flags to mask topography below and above. Flag 0 is
363!--             used to mask topography in general, and flag 8 implies
364!--             information about use_surface_fluxes. Flag 9 is used to control
365!--             flux at model top.
366                mask_bottom = MERGE( 1.0_wp, 0.0_wp,                           &
367                                     BTEST( wall_flags_0(k-1,j,i), 8 ) ) 
368                mask_top    = MERGE( 1.0_wp, 0.0_wp,                           &
369                                     BTEST( wall_flags_0(k+1,j,i), 8 ) ) *     &
370                              MERGE( 1.0_wp, 0.0_wp,                           &
371                                     BTEST( wall_flags_0(k+1,j,i), 9 ) ) 
372                flag        = MERGE( 1.0_wp, 0.0_wp,                           &
373                                     BTEST( wall_flags_0(k,j,i), 0 ) )
374
375                tend(k,j,i) = tend(k,j,i)                                      &
376                                       + 0.5_wp * (                            &
377                                      ( kh(k,j,i) + kh(k+1,j,i) ) *            &
378                                          ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1)  &
379                                                            * rho_air_zw(k)    &
380                                                            * mask_top         &
381                                    - ( kh(k,j,i) + kh(k-1,j,i) ) *            &
382                                          ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k)    &
383                                                            * rho_air_zw(k-1)  &
384                                                            * mask_bottom      &
385                                                  ) * ddzw(k) * drho_air(k)    &
386                                                              * flag
387             ENDDO
388
389!
390!--          Vertical diffusion at horizontal walls.
391             IF ( use_surface_fluxes )  THEN
392!
393!--             Default-type surfaces, upward-facing               
394                surf_s = surf_def_h(0)%start_index(j,i)
395                surf_e = surf_def_h(0)%end_index(j,i)
396                DO  m = surf_s, surf_e
397
398                   k   = surf_def_h(0)%k(m)
399                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m)              &
400                                       * ddzw(k) * drho_air(k)
401
402                ENDDO
403!
404!--             Default-type surfaces, downward-facing               
405                surf_s = surf_def_h(1)%start_index(j,i)
406                surf_e = surf_def_h(1)%end_index(j,i)
407                DO  m = surf_s, surf_e
408
409                   k   = surf_def_h(1)%k(m)
410                   tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m)            &
411                                       * ddzw(k) * drho_air(k)
412
413                ENDDO
414!
415!--             Natural-type surfaces, upward-facing 
416                surf_s = surf_lsm_h%start_index(j,i)
417                surf_e = surf_lsm_h%end_index(j,i)
418                DO  m = surf_s, surf_e
419
420                   k   = surf_lsm_h%k(m)
421                   tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m)              &
422                                       * ddzw(k) * drho_air(k)
423
424                ENDDO
425!
426!--             Urban-type surfaces, upward-facing     
427                surf_s = surf_usm_h%start_index(j,i)
428                surf_e = surf_usm_h%end_index(j,i)
429                DO  m = surf_s, surf_e
430
431                   k   = surf_usm_h%k(m)
432                   tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m)              &
433                                       * ddzw(k) * drho_air(k)
434
435                ENDDO
436
437             ENDIF
438!
439!--          Vertical diffusion at the last computational gridpoint along z-direction
440             IF ( use_top_fluxes )  THEN
441                surf_s = surf_def_h(2)%start_index(j,i)
442                surf_e = surf_def_h(2)%end_index(j,i)
443                DO  m = surf_s, surf_e
444
445                   k   = surf_def_h(2)%k(m)
446                   tend(k,j,i) = tend(k,j,i)                                   &
447                           + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
448                ENDDO
449             ENDIF
450
451          ENDDO
452       ENDDO
453
454    END SUBROUTINE diffusion_s
455
456!------------------------------------------------------------------------------!
457! Description:
458! ------------
459!> Call for grid point i,j
460!------------------------------------------------------------------------------!
461    SUBROUTINE diffusion_s_ij( i, j, s,                                        &
462                               s_flux_def_h_up,    s_flux_def_h_down,          &
463                               s_flux_t,                                       &
464                               s_flux_lsm_h_up,    s_flux_usm_h_up,            &
465                               s_flux_def_v_north, s_flux_def_v_south,         &
466                               s_flux_def_v_east,  s_flux_def_v_west,          &
467                               s_flux_lsm_v_north, s_flux_lsm_v_south,         &
468                               s_flux_lsm_v_east,  s_flux_lsm_v_west,          &
469                               s_flux_usm_v_north, s_flux_usm_v_south,         &
470                               s_flux_usm_v_east,  s_flux_usm_v_west )       
471
472       USE arrays_3d,                                                          &
473           ONLY:  ddzu, ddzw, kh, tend, drho_air, rho_air_zw
474           
475       USE control_parameters,                                                 & 
476           ONLY: use_surface_fluxes, use_top_fluxes
477       
478       USE grid_variables,                                                     &
479           ONLY:  ddx, ddx2, ddy, ddy2
480       
481       USE indices,                                                            &
482           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
483       
484       USE kinds
485
486       USE surface_mod,                                                        &
487           ONLY :  surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
488                   surf_usm_v
489
490       IMPLICIT NONE
491
492       INTEGER(iwp) ::  i             !< running index x direction
493       INTEGER(iwp) ::  j             !< running index y direction
494       INTEGER(iwp) ::  k             !< running index z direction
495       INTEGER(iwp) ::  m             !< running index surface elements
496       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
497       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
498
499       REAL(wp) ::  flag              !< flag to mask topography grid points
500       REAL(wp) ::  mask_bottom       !< flag to mask vertical upward-facing surface     
501       REAL(wp) ::  mask_east         !< flag to mask vertical surface east of the grid point
502       REAL(wp) ::  mask_north        !< flag to mask vertical surface north of the grid point
503       REAL(wp) ::  mask_south        !< flag to mask vertical surface south of the grid point
504       REAL(wp) ::  mask_west         !< flag to mask vertical surface west of the grid point
505       REAL(wp) ::  mask_top          !< flag to mask vertical downward-facing surface 
506
507       REAL(wp), DIMENSION(1:surf_def_v(0)%ns) ::  s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
508       REAL(wp), DIMENSION(1:surf_def_v(1)%ns) ::  s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
509       REAL(wp), DIMENSION(1:surf_def_v(2)%ns) ::  s_flux_def_v_east  !< flux at east-facing vertical default-type surfaces
510       REAL(wp), DIMENSION(1:surf_def_v(3)%ns) ::  s_flux_def_v_west  !< flux at west-facing vertical default-type surfaces
511       REAL(wp), DIMENSION(1:surf_def_h(0)%ns) ::  s_flux_def_h_up    !< flux at horizontal upward-facing default-type surfaces
512       REAL(wp), DIMENSION(1:surf_def_h(1)%ns) ::  s_flux_def_h_down  !< flux at horizontal donwward-facing default-type surfaces
513       REAL(wp), DIMENSION(1:surf_lsm_h%ns)    ::  s_flux_lsm_h_up    !< flux at horizontal upward-facing natural-type surfaces
514       REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) ::  s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
515       REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) ::  s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
516       REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) ::  s_flux_lsm_v_east  !< flux at east-facing vertical urban-type surfaces
517       REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) ::  s_flux_lsm_v_west  !< flux at west-facing vertical urban-type surfaces
518       REAL(wp), DIMENSION(1:surf_usm_h%ns)    ::  s_flux_usm_h_up    !< flux at horizontal upward-facing urban-type surfaces
519       REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) ::  s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
520       REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) ::  s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
521       REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) ::  s_flux_usm_v_east  !< flux at east-facing vertical urban-type surfaces
522       REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) ::  s_flux_usm_v_west  !< flux at west-facing vertical urban-type surfaces
523       REAL(wp), DIMENSION(1:surf_def_h(2)%ns) ::  s_flux_t           !< flux at model top
524#if defined( __nopointer )
525       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  s !< treated scalar
526#else
527       REAL(wp), DIMENSION(:,:,:), POINTER ::  s  !< treated scalar
528#endif
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.