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

Last change on this file since 4329 was 4329, checked in by motisi, 4 years ago

Renamed wall_flags_0 to wall_flags_static_0

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