Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |