Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |