Home

Context Navigation

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

Last change on this file since 4326 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 32.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |