Home

Context Navigation

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

Last change on this file since 4598 was 4583, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 33.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |