Home

Context Navigation

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

Last change on this file since 2608 was 2233, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 33.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |