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

Last change on this file since 3567 was 3547, checked in by suehring, 6 years ago

variable description added some routines

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