Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 32.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |