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

Last change on this file since 2707 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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