Home

Context Navigation

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

Last change on this file since 4069 was 3927, checked in by raasch, 6 years ago
pointer attribute removed from scalar 3d-array for performance reasons
Property svn:keywords set to `Id`
File size: 34.7 KB

Rev	Line
[1873]	1	!> @file diffusion_s.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[2000]	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.
[1036]	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	!
[3655]	17	! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1001]	21	! ------------------
[1341]	22	!
[2233]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: diffusion_s.f90 3927 2019-04-23 13:24:29Z Giersch $
[3927]	27	! pointer attribute removed from scalar 3d-array for performance reasons
	28	!
	29	! 3761 2019-02-25 15:31:42Z raasch
[3761]	30	! unused variables removed
	31	!
	32	! 3655 2019-01-07 16:51:22Z knoop
[3636]	33	! nopointer option removed
	34	!
	35	! 3634 2018-12-18 12:31:28Z knoop
[3634]	36	! OpenACC port for SPEC
	37	!
	38	! 3547 2018-11-21 13:21:24Z suehring
[3547]	39	! variables documented
	40	!
	41	! 2759 2018-01-17 16:24:59Z suehring
[2759]	42	! Major bugfix, horizontal diffusion at vertical surfaces corrected.
	43	!
	44	! 2718 2018-01-02 08:49:38Z maronga
[2716]	45	! Corrected "Former revisions" section
	46	!
	47	! 2696 2017-12-14 17:12:51Z kanani
	48	! Change in file header (GPL part)
[1321]	49	!
[2716]	50	! 2233 2017-05-30 18:08:54Z suehring
	51	!
[2233]	52	! 2232 2017-05-30 17:47:52Z suehring
	53	! Adjustments to new topography and surface concept
	54	!
[2119]	55	! 2118 2017-01-17 16:38:49Z raasch
	56	! OpenACC version of subroutine removed
	57	!
[2038]	58	! 2037 2016-10-26 11:15:40Z knoop
	59	! Anelastic approximation implemented
	60	!
[2001]	61	! 2000 2016-08-20 18:09:15Z knoop
	62	! Forced header and separation lines into 80 columns
	63	!
[1874]	64	! 1873 2016-04-18 14:50:06Z maronga
	65	! Module renamed (removed _mod)
[2118]	66	!
[1851]	67	! 1850 2016-04-08 13:29:27Z maronga
	68	! Module renamed
	69	!
[1692]	70	! 1691 2015-10-26 16:17:44Z maronga
	71	! Formatting corrections.
	72	!
[1683]	73	! 1682 2015-10-07 23:56:08Z knoop
	74	! Code annotations made doxygen readable
	75	!
[1375]	76	! 1374 2014-04-25 12:55:07Z raasch
	77	! missing variables added to ONLY list
	78	!
[1341]	79	! 1340 2014-03-25 19:45:13Z kanani
	80	! REAL constants defined as wp-kind
	81	!
[1321]	82	! 1320 2014-03-20 08:40:49Z raasch
[1320]	83	! ONLY-attribute added to USE-statements,
	84	! kind-parameters added to all INTEGER and REAL declaration statements,
	85	! kinds are defined in new module kinds,
	86	! revision history before 2012 removed,
	87	! comment fields (!:) to be used for variable explanations added to
	88	! all variable declaration statements
[1321]	89	!
[1258]	90	! 1257 2013-11-08 15:18:40Z raasch
	91	! openacc loop and loop vector clauses removed
	92	!
[1132]	93	! 1128 2013-04-12 06:19:32Z raasch
	94	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	95	! j_north
	96	!
[1093]	97	! 1092 2013-02-02 11:24:22Z raasch
	98	! unused variables removed
	99	!
[1037]	100	! 1036 2012-10-22 13:43:42Z raasch
	101	! code put under GPL (PALM 3.9)
	102	!
[1017]	103	! 1015 2012-09-27 09:23:24Z raasch
	104	! accelerator version (*_acc) added
	105	!
[1011]	106	! 1010 2012-09-20 07:59:54Z raasch
	107	! cpp switch __nopointer added for pointer free version
	108	!
[1002]	109	! 1001 2012-09-13 14:08:46Z raasch
	110	! some arrays comunicated by module instead of parameter list
	111	!
[1]	112	! Revision 1.1 2000/04/13 14:54:02 schroeter
	113	! Initial revision
	114	!
	115	!
	116	! Description:
	117	! ------------
[1682]	118	!> Diffusion term of scalar quantities (temperature and water content)
[1]	119	!------------------------------------------------------------------------------!
[1682]	120	MODULE diffusion_s_mod
	121
[1]	122
	123	PRIVATE
[2118]	124	PUBLIC diffusion_s
[1]	125
	126	INTERFACE diffusion_s
	127	MODULE PROCEDURE diffusion_s
	128	MODULE PROCEDURE diffusion_s_ij
	129	END INTERFACE diffusion_s
	130
	131	CONTAINS
	132
	133
	134	!------------------------------------------------------------------------------!
[1682]	135	! Description:
	136	! ------------
	137	!> Call for all grid points
[1]	138	!------------------------------------------------------------------------------!
[2232]	139	SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, &
	140	s_flux_t, &
	141	s_flux_lsm_h_up, s_flux_usm_h_up, &
	142	s_flux_def_v_north, s_flux_def_v_south, &
	143	s_flux_def_v_east, s_flux_def_v_west, &
	144	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	145	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	146	s_flux_usm_v_north, s_flux_usm_v_south, &
	147	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	148
[1320]	149	USE arrays_3d, &
[2037]	150	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[1320]	151
	152	USE control_parameters, &
	153	ONLY: use_surface_fluxes, use_top_fluxes
	154
	155	USE grid_variables, &
[2759]	156	ONLY: ddx, ddx2, ddy, ddy2
[1320]	157
	158	USE indices, &
[3927]	159	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_0
[1320]	160
	161	USE kinds
[1]	162
[2232]	163	USE surface_mod, &
	164	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	165	surf_usm_v
	166
[1]	167	IMPLICIT NONE
	168
[2232]	169	INTEGER(iwp) :: i !< running index x direction
	170	INTEGER(iwp) :: j !< running index y direction
	171	INTEGER(iwp) :: k !< running index z direction
	172	INTEGER(iwp) :: m !< running index surface elements
	173	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	174	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	175
	176	REAL(wp) :: flag !< flag to mask topography grid points
	177	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	178	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
	179	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	180	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	181	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	182	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	183
	184	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	185	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
	186	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
	187	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	188	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	189	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
	190	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
	191	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
	192	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
	193	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical natural-type surfaces
	194	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical natural-type surfaces
	195	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
	196	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	197	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	198	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
	199	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
	200	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
[3636]	201
[3927]	202	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
[1]	203
[3636]	204
[3634]	205	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
	206	!$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
	207	!$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
	208	!$ACC PRESENT(wall_flags_0, kh) &
	209	!$ACC PRESENT(s) &
	210	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
	211	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
	212	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
	213	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
	214	!$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
	215	!$ACC PRESENT(s_flux_t) &
	216	!$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
	217	!$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
	218	!$ACC PRESENT(s_flux_lsm_h_up) &
	219	!$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
	220	!$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
	221	!$ACC PRESENT(s_flux_usm_h_up) &
	222	!$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
	223	!$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
	224	!$ACC PRESENT(tend)
[1]	225	DO i = nxl, nxr
	226	DO j = nys,nyn
	227	!
	228	!-- Compute horizontal diffusion
[2232]	229	DO k = nzb+1, nzt
	230	!
	231	!-- Predetermine flag to mask topography and wall-bounded grid points
	232	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	233	!
	234	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
	235	!-- former s_outer array
	236	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
	237	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
	238	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
	239	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
[1]	240
[1320]	241	tend(k,j,i) = tend(k,j,i) &
[1340]	242	+ 0.5_wp * ( &
[2232]	243	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
	244	* ( s(k,j,i+1) - s(k,j,i) ) &
	245	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
	246	* ( s(k,j,i) - s(k,j,i-1) ) &
	247	) * ddx2 * flag &
[1340]	248	+ 0.5_wp * ( &
[2232]	249	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
	250	* ( s(k,j+1,i) - s(k,j,i) ) &
	251	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
	252	* ( s(k,j,i) - s(k,j-1,i) ) &
	253	) * ddy2 * flag
[1]	254	ENDDO
	255
	256	!
[2232]	257	!-- Apply prescribed horizontal wall heatflux where necessary. First,
	258	!-- determine start and end index for respective (j,i)-index. Please
	259	!-- note, in the flat case following loop will not be entered, as
	260	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
	261	!-- so far.
	262	!-- First, for default-type surfaces
	263	!-- North-facing vertical default-type surfaces
	264	surf_s = surf_def_v(0)%start_index(j,i)
	265	surf_e = surf_def_v(0)%end_index(j,i)
	266	DO m = surf_s, surf_e
	267	k = surf_def_v(0)%k(m)
[2759]	268	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	269	ENDDO
	270	!
	271	!-- South-facing vertical default-type surfaces
	272	surf_s = surf_def_v(1)%start_index(j,i)
	273	surf_e = surf_def_v(1)%end_index(j,i)
	274	DO m = surf_s, surf_e
	275	k = surf_def_v(1)%k(m)
[2759]	276	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	277	ENDDO
	278	!
	279	!-- East-facing vertical default-type surfaces
	280	surf_s = surf_def_v(2)%start_index(j,i)
	281	surf_e = surf_def_v(2)%end_index(j,i)
	282	DO m = surf_s, surf_e
	283	k = surf_def_v(2)%k(m)
[2759]	284	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	285	ENDDO
	286	!
	287	!-- West-facing vertical default-type surfaces
	288	surf_s = surf_def_v(3)%start_index(j,i)
	289	surf_e = surf_def_v(3)%end_index(j,i)
	290	DO m = surf_s, surf_e
	291	k = surf_def_v(3)%k(m)
[2759]	292	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	293	ENDDO
	294	!
	295	!-- Now, for natural-type surfaces.
	296	!-- North-facing
	297	surf_s = surf_lsm_v(0)%start_index(j,i)
	298	surf_e = surf_lsm_v(0)%end_index(j,i)
	299	DO m = surf_s, surf_e
	300	k = surf_lsm_v(0)%k(m)
[2759]	301	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	302	ENDDO
	303	!
	304	!-- South-facing
	305	surf_s = surf_lsm_v(1)%start_index(j,i)
	306	surf_e = surf_lsm_v(1)%end_index(j,i)
	307	DO m = surf_s, surf_e
	308	k = surf_lsm_v(1)%k(m)
[2759]	309	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	310	ENDDO
	311	!
	312	!-- East-facing
	313	surf_s = surf_lsm_v(2)%start_index(j,i)
	314	surf_e = surf_lsm_v(2)%end_index(j,i)
	315	DO m = surf_s, surf_e
	316	k = surf_lsm_v(2)%k(m)
[2759]	317	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	318	ENDDO
	319	!
	320	!-- West-facing
	321	surf_s = surf_lsm_v(3)%start_index(j,i)
	322	surf_e = surf_lsm_v(3)%end_index(j,i)
	323	DO m = surf_s, surf_e
	324	k = surf_lsm_v(3)%k(m)
[2759]	325	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	326	ENDDO
	327	!
	328	!-- Now, for urban-type surfaces.
	329	!-- North-facing
	330	surf_s = surf_usm_v(0)%start_index(j,i)
	331	surf_e = surf_usm_v(0)%end_index(j,i)
	332	DO m = surf_s, surf_e
	333	k = surf_usm_v(0)%k(m)
[2759]	334	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	335	ENDDO
	336	!
	337	!-- South-facing
	338	surf_s = surf_usm_v(1)%start_index(j,i)
	339	surf_e = surf_usm_v(1)%end_index(j,i)
	340	DO m = surf_s, surf_e
	341	k = surf_usm_v(1)%k(m)
[2759]	342	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	343	ENDDO
	344	!
	345	!-- East-facing
	346	surf_s = surf_usm_v(2)%start_index(j,i)
	347	surf_e = surf_usm_v(2)%end_index(j,i)
	348	DO m = surf_s, surf_e
	349	k = surf_usm_v(2)%k(m)
[2759]	350	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	351	ENDDO
	352	!
	353	!-- West-facing
	354	surf_s = surf_usm_v(3)%start_index(j,i)
	355	surf_e = surf_usm_v(3)%end_index(j,i)
	356	DO m = surf_s, surf_e
	357	k = surf_usm_v(3)%k(m)
[2759]	358	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	359	ENDDO
[1]	360
	361	!
	362	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	363	!-- prescribed or computed at bottom and/or top, index k starts/ends at
[2232]	364	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
	365	!-- is given.
	366	DO k = nzb+1, nzt
	367	!
	368	!-- Determine flags to mask topography below and above. Flag 0 is
	369	!-- used to mask topography in general, and flag 8 implies
	370	!-- information about use_surface_fluxes. Flag 9 is used to control
	371	!-- flux at model top.
	372	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
	373	BTEST( wall_flags_0(k-1,j,i), 8 ) )
	374	mask_top = MERGE( 1.0_wp, 0.0_wp, &
	375	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
	376	MERGE( 1.0_wp, 0.0_wp, &
	377	BTEST( wall_flags_0(k+1,j,i), 9 ) )
	378	flag = MERGE( 1.0_wp, 0.0_wp, &
	379	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	380
[1320]	381	tend(k,j,i) = tend(k,j,i) &
[1340]	382	+ 0.5_wp * ( &
[2232]	383	( kh(k,j,i) + kh(k+1,j,i) ) * &
	384	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
[2037]	385	* rho_air_zw(k) &
[2232]	386	* mask_top &
	387	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	388	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
[2037]	389	* rho_air_zw(k-1) &
[2232]	390	* mask_bottom &
	391	) * ddzw(k) * drho_air(k) &
	392	* flag
[1]	393	ENDDO
	394
	395	!
[2232]	396	!-- Vertical diffusion at horizontal walls.
[1]	397	IF ( use_surface_fluxes ) THEN
[2232]	398	!
	399	!-- Default-type surfaces, upward-facing
	400	surf_s = surf_def_h(0)%start_index(j,i)
	401	surf_e = surf_def_h(0)%end_index(j,i)
	402	DO m = surf_s, surf_e
[1]	403
[2232]	404	k = surf_def_h(0)%k(m)
	405	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
	406	* ddzw(k) * drho_air(k)
[1]	407
[2232]	408	ENDDO
	409	!
	410	!-- Default-type surfaces, downward-facing
	411	surf_s = surf_def_h(1)%start_index(j,i)
	412	surf_e = surf_def_h(1)%end_index(j,i)
	413	DO m = surf_s, surf_e
[1]	414
[2232]	415	k = surf_def_h(1)%k(m)
	416	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
	417	* ddzw(k) * drho_air(k)
[1]	418
[2232]	419	ENDDO
[19]	420	!
[2232]	421	!-- Natural-type surfaces, upward-facing
	422	surf_s = surf_lsm_h%start_index(j,i)
	423	surf_e = surf_lsm_h%end_index(j,i)
	424	DO m = surf_s, surf_e
[19]	425
[2232]	426	k = surf_lsm_h%k(m)
	427	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
	428	* ddzw(k) * drho_air(k)
[19]	429
[2232]	430	ENDDO
	431	!
	432	!-- Urban-type surfaces, upward-facing
	433	surf_s = surf_usm_h%start_index(j,i)
	434	surf_e = surf_usm_h%end_index(j,i)
	435	DO m = surf_s, surf_e
[19]	436
[2232]	437	k = surf_usm_h%k(m)
	438	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
	439	* ddzw(k) * drho_air(k)
	440
	441	ENDDO
	442
[19]	443	ENDIF
[2232]	444	!
	445	!-- Vertical diffusion at the last computational gridpoint along z-direction
	446	IF ( use_top_fluxes ) THEN
	447	surf_s = surf_def_h(2)%start_index(j,i)
	448	surf_e = surf_def_h(2)%end_index(j,i)
	449	DO m = surf_s, surf_e
[19]	450
[2232]	451	k = surf_def_h(2)%k(m)
	452	tend(k,j,i) = tend(k,j,i) &
	453	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
	454	ENDDO
	455	ENDIF
	456
[1]	457	ENDDO
	458	ENDDO
	459
	460	END SUBROUTINE diffusion_s
	461
	462	!------------------------------------------------------------------------------!
[1682]	463	! Description:
	464	! ------------
	465	!> Call for grid point i,j
[1]	466	!------------------------------------------------------------------------------!
[2232]	467	SUBROUTINE diffusion_s_ij( i, j, s, &
	468	s_flux_def_h_up, s_flux_def_h_down, &
	469	s_flux_t, &
	470	s_flux_lsm_h_up, s_flux_usm_h_up, &
	471	s_flux_def_v_north, s_flux_def_v_south, &
	472	s_flux_def_v_east, s_flux_def_v_west, &
	473	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	474	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	475	s_flux_usm_v_north, s_flux_usm_v_south, &
	476	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	477
[1320]	478	USE arrays_3d, &
[2037]	479	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[1320]	480
	481	USE control_parameters, &
	482	ONLY: use_surface_fluxes, use_top_fluxes
	483
	484	USE grid_variables, &
[2759]	485	ONLY: ddx, ddx2, ddy, ddy2
[1320]	486
	487	USE indices, &
[3927]	488	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_0
[1320]	489
	490	USE kinds
[1]	491
[2232]	492	USE surface_mod, &
	493	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	494	surf_usm_v
	495
[1]	496	IMPLICIT NONE
	497
[2232]	498	INTEGER(iwp) :: i !< running index x direction
	499	INTEGER(iwp) :: j !< running index y direction
	500	INTEGER(iwp) :: k !< running index z direction
	501	INTEGER(iwp) :: m !< running index surface elements
	502	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	503	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	504
	505	REAL(wp) :: flag !< flag to mask topography grid points
	506	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	507	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
	508	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	509	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	510	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	511	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	512
	513	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	514	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
	515	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
	516	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	517	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	518	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
	519	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
	520	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
	521	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
	522	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical urban-type surfaces
	523	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical urban-type surfaces
	524	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
	525	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	526	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	527	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
	528	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
	529	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
[3636]	530
[3927]	531	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
[1]	532
	533	!
	534	!-- Compute horizontal diffusion
[2232]	535	DO k = nzb+1, nzt
	536	!
	537	!-- Predetermine flag to mask topography and wall-bounded grid points
	538	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
	539	!
	540	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
	541	!-- former s_outer array
	542	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i-1), 0 ) )
	543	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i+1), 0 ) )
	544	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j-1,i), 0 ) )
	545	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j+1,i), 0 ) )
	546	!
	547	!-- Finally, determine flag to mask both topography itself as well
	548	!-- as wall-bounded grid points, which will be treated further below
[1]	549
[1320]	550	tend(k,j,i) = tend(k,j,i) &
[1340]	551	+ 0.5_wp * ( &
[2232]	552	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
	553	* ( s(k,j,i+1) - s(k,j,i) ) &
	554	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
	555	* ( s(k,j,i) - s(k,j,i-1) ) &
	556	) * ddx2 * flag &
[1340]	557	+ 0.5_wp * ( &
[2232]	558	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
	559	* ( s(k,j+1,i) - s(k,j,i) ) &
	560	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
	561	* ( s(k,j,i) - s(k,j-1,i) ) &
	562	) * ddy2 * flag
[1]	563	ENDDO
	564
	565	!
[2232]	566	!-- Apply prescribed horizontal wall heatflux where necessary. First,
	567	!-- determine start and end index for respective (j,i)-index. Please
	568	!-- note, in the flat case following loops will not be entered, as
	569	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
	570	!-- so far.
	571	!-- First, for default-type surfaces
	572	!-- North-facing vertical default-type surfaces
	573	surf_s = surf_def_v(0)%start_index(j,i)
	574	surf_e = surf_def_v(0)%end_index(j,i)
	575	DO m = surf_s, surf_e
	576	k = surf_def_v(0)%k(m)
[2759]	577	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	578	ENDDO
	579	!
	580	!-- South-facing vertical default-type surfaces
	581	surf_s = surf_def_v(1)%start_index(j,i)
	582	surf_e = surf_def_v(1)%end_index(j,i)
	583	DO m = surf_s, surf_e
	584	k = surf_def_v(1)%k(m)
[2759]	585	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	586	ENDDO
	587	!
	588	!-- East-facing vertical default-type surfaces
	589	surf_s = surf_def_v(2)%start_index(j,i)
	590	surf_e = surf_def_v(2)%end_index(j,i)
	591	DO m = surf_s, surf_e
	592	k = surf_def_v(2)%k(m)
[2759]	593	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	594	ENDDO
	595	!
	596	!-- West-facing vertical default-type surfaces
	597	surf_s = surf_def_v(3)%start_index(j,i)
	598	surf_e = surf_def_v(3)%end_index(j,i)
	599	DO m = surf_s, surf_e
	600	k = surf_def_v(3)%k(m)
[2759]	601	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	602	ENDDO
	603	!
	604	!-- Now, for natural-type surfaces
	605	!-- North-facing
	606	surf_s = surf_lsm_v(0)%start_index(j,i)
	607	surf_e = surf_lsm_v(0)%end_index(j,i)
	608	DO m = surf_s, surf_e
	609	k = surf_lsm_v(0)%k(m)
[2759]	610	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	611	ENDDO
	612	!
	613	!-- South-facing
	614	surf_s = surf_lsm_v(1)%start_index(j,i)
	615	surf_e = surf_lsm_v(1)%end_index(j,i)
	616	DO m = surf_s, surf_e
	617	k = surf_lsm_v(1)%k(m)
[2759]	618	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	619	ENDDO
	620	!
	621	!-- East-facing
	622	surf_s = surf_lsm_v(2)%start_index(j,i)
	623	surf_e = surf_lsm_v(2)%end_index(j,i)
	624	DO m = surf_s, surf_e
	625	k = surf_lsm_v(2)%k(m)
[2759]	626	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	627	ENDDO
	628	!
	629	!-- West-facing
	630	surf_s = surf_lsm_v(3)%start_index(j,i)
	631	surf_e = surf_lsm_v(3)%end_index(j,i)
	632	DO m = surf_s, surf_e
	633	k = surf_lsm_v(3)%k(m)
[2759]	634	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	635	ENDDO
	636	!
	637	!-- Now, for urban-type surfaces
	638	!-- North-facing
	639	surf_s = surf_usm_v(0)%start_index(j,i)
	640	surf_e = surf_usm_v(0)%end_index(j,i)
	641	DO m = surf_s, surf_e
	642	k = surf_usm_v(0)%k(m)
[2759]	643	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	644	ENDDO
	645	!
	646	!-- South-facing
	647	surf_s = surf_usm_v(1)%start_index(j,i)
	648	surf_e = surf_usm_v(1)%end_index(j,i)
	649	DO m = surf_s, surf_e
	650	k = surf_usm_v(1)%k(m)
[2759]	651	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	652	ENDDO
	653	!
	654	!-- East-facing
	655	surf_s = surf_usm_v(2)%start_index(j,i)
	656	surf_e = surf_usm_v(2)%end_index(j,i)
	657	DO m = surf_s, surf_e
	658	k = surf_usm_v(2)%k(m)
[2759]	659	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	660	ENDDO
	661	!
	662	!-- West-facing
	663	surf_s = surf_usm_v(3)%start_index(j,i)
	664	surf_e = surf_usm_v(3)%end_index(j,i)
	665	DO m = surf_s, surf_e
	666	k = surf_usm_v(3)%k(m)
[2759]	667	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	668	ENDDO
[1]	669
	670
	671	!
	672	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	673	!-- prescribed or computed at bottom and/or top, index k starts/ends at
[2232]	674	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
	675	!-- is given.
	676	DO k = nzb+1, nzt
	677	!
	678	!-- Determine flags to mask topography below and above. Flag 0 is
	679	!-- used to mask topography in general, and flag 8 implies
	680	!-- information about use_surface_fluxes. Flag 9 is used to control
	681	!-- flux at model top.
	682	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
	683	BTEST( wall_flags_0(k-1,j,i), 8 ) )
	684	mask_top = MERGE( 1.0_wp, 0.0_wp, &
	685	BTEST( wall_flags_0(k+1,j,i), 8 ) ) * &
	686	MERGE( 1.0_wp, 0.0_wp, &
	687	BTEST( wall_flags_0(k+1,j,i), 9 ) )
	688	flag = MERGE( 1.0_wp, 0.0_wp, &
	689	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	690
[1320]	691	tend(k,j,i) = tend(k,j,i) &
[1340]	692	+ 0.5_wp * ( &
[2232]	693	( kh(k,j,i) + kh(k+1,j,i) ) * &
	694	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
[2037]	695	* rho_air_zw(k) &
[2232]	696	* mask_top &
	697	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	698	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
[2037]	699	* rho_air_zw(k-1) &
[2232]	700	* mask_bottom &
	701	) * ddzw(k) * drho_air(k) &
	702	* flag
[1]	703	ENDDO
	704
	705	!
[2232]	706	!-- Vertical diffusion at horizontal walls.
	707	!-- TO DO: Adjust for downward facing walls and mask already in main loop
[1]	708	IF ( use_surface_fluxes ) THEN
[2232]	709	!
	710	!-- Default-type surfaces, upward-facing
	711	surf_s = surf_def_h(0)%start_index(j,i)
	712	surf_e = surf_def_h(0)%end_index(j,i)
	713	DO m = surf_s, surf_e
[1]	714
[2232]	715	k = surf_def_h(0)%k(m)
[1]	716
[2232]	717	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
	718	* ddzw(k) * drho_air(k)
	719	ENDDO
	720	!
	721	!-- Default-type surfaces, downward-facing
	722	surf_s = surf_def_h(1)%start_index(j,i)
	723	surf_e = surf_def_h(1)%end_index(j,i)
	724	DO m = surf_s, surf_e
[1]	725
[2232]	726	k = surf_def_h(1)%k(m)
[1]	727
[2232]	728	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
	729	* ddzw(k) * drho_air(k)
	730	ENDDO
[19]	731	!
[2232]	732	!-- Natural-type surfaces, upward-facing
	733	surf_s = surf_lsm_h%start_index(j,i)
	734	surf_e = surf_lsm_h%end_index(j,i)
	735	DO m = surf_s, surf_e
	736	k = surf_lsm_h%k(m)
	737
	738	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
	739	* ddzw(k) * drho_air(k)
	740	ENDDO
	741	!
	742	!-- Urban-type surfaces, upward-facing
	743	surf_s = surf_usm_h%start_index(j,i)
	744	surf_e = surf_usm_h%end_index(j,i)
	745	DO m = surf_s, surf_e
	746	k = surf_usm_h%k(m)
	747
	748	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
	749	* ddzw(k) * drho_air(k)
	750	ENDDO
	751	ENDIF
	752	!
[19]	753	!-- Vertical diffusion at the last computational gridpoint along z-direction
	754	IF ( use_top_fluxes ) THEN
[2232]	755	surf_s = surf_def_h(2)%start_index(j,i)
	756	surf_e = surf_def_h(2)%end_index(j,i)
	757	DO m = surf_s, surf_e
[19]	758
[2232]	759	k = surf_def_h(2)%k(m)
	760	tend(k,j,i) = tend(k,j,i) &
	761	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
	762	ENDDO
[19]	763	ENDIF
	764
[1]	765	END SUBROUTINE diffusion_s_ij
	766
	767	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |