Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |