Home

Context Navigation

diffusion_s.f90 @ 4867

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

Context Navigation

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

Download in other formats: