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

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