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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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