Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |