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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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