Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |