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

Last change on this file since 3159 was 2759, checked in by suehring, 7 years ago

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

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