Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |