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

Last change on this file since 2978 was 2759, checked in by suehring, 6 years ago

Major bugfix in horizontal diffusion of all scalar quantities at vertical surfaces; further bugfix: density is considered for vertical fluxes of passive scalar and salinity

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