Home

Context Navigation

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

Last change on this file since 4596 was 4583, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 33.4 KB

Rev	Line
[1873]	1	!> @file diffusion_s.f90
[4583]	2	!--------------------------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	4	!
[4583]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[1036]	8	!
[4583]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[1036]	12	!
[4583]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[1036]	15	!
[4360]	16	! Copyright 1997-2020 Leibniz Universitaet Hannover
[4583]	17	!--------------------------------------------------------------------------------------------------!
[1036]	18	!
[484]	19	! Current revisions:
[1001]	20	! ------------------
[1341]	21	!
[2233]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: diffusion_s.f90 4583 2020-06-29 12:36:47Z suehring $
[4583]	26	! file re-formatted to follow the PALM coding standard
	27	!
	28	! 4360 2020-01-07 11:25:50Z suehring
	29	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	30	! information used in wall_flags_static_0
	31	!
[4346]	32	! 4329 2019-12-10 15:46:36Z motisi
	33	! Renamed wall_flags_0 to wall_flags_total_0
[4583]	34	!
[4329]	35	! 4182 2019-08-22 15:20:23Z scharf
[4182]	36	! Corrected "Former revisions" section
[4583]	37	!
[4182]	38	! 3927 2019-04-23 13:24:29Z raasch
[3927]	39	! pointer attribute removed from scalar 3d-array for performance reasons
[4583]	40	!
[3927]	41	! 3761 2019-02-25 15:31:42Z raasch
[3761]	42	! unused variables removed
[4583]	43	!
[3761]	44	! 3655 2019-01-07 16:51:22Z knoop
[3636]	45	! nopointer option removed
[1321]	46	!
[4182]	47	! Revision 1.1 2000/04/13 14:54:02 schroeter
	48	! Initial revision
	49	!
	50	!
[1]	51	! Description:
	52	! ------------
[1682]	53	!> Diffusion term of scalar quantities (temperature and water content)
[4583]	54	!--------------------------------------------------------------------------------------------------!
[1682]	55	MODULE diffusion_s_mod
[1]	56
[4583]	57
[1]	58	PRIVATE
[2118]	59	PUBLIC diffusion_s
[1]	60
	61	INTERFACE diffusion_s
	62	MODULE PROCEDURE diffusion_s
	63	MODULE PROCEDURE diffusion_s_ij
	64	END INTERFACE diffusion_s
	65
	66	CONTAINS
	67
	68
[4583]	69	!--------------------------------------------------------------------------------------------------!
[1682]	70	! Description:
	71	! ------------
	72	!> Call for all grid points
[4583]	73	!--------------------------------------------------------------------------------------------------!
	74	SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, &
	75	s_flux_t, &
	76	s_flux_lsm_h_up, s_flux_usm_h_up, &
	77	s_flux_def_v_north, s_flux_def_v_south, &
	78	s_flux_def_v_east, s_flux_def_v_west, &
	79	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	80	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	81	s_flux_usm_v_north, s_flux_usm_v_south, &
[2232]	82	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	83
[4583]	84	USE arrays_3d, &
[2037]	85	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[4583]	86
	87	USE control_parameters, &
[1320]	88	ONLY: use_surface_fluxes, use_top_fluxes
[4583]	89
	90	USE grid_variables, &
[2759]	91	ONLY: ddx, ddx2, ddy, ddy2
[4583]	92
	93	USE indices, &
	94	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, wall_flags_total_0
	95
[1320]	96	USE kinds
[1]	97
[4583]	98	USE surface_mod, &
	99	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	100
[1]	101	IMPLICIT NONE
	102
[2232]	103	INTEGER(iwp) :: i !< running index x direction
	104	INTEGER(iwp) :: j !< running index y direction
	105	INTEGER(iwp) :: k !< running index z direction
	106	INTEGER(iwp) :: m !< running index surface elements
	107	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	108	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	109
	110	REAL(wp) :: flag !< flag to mask topography grid points
[4583]	111	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	112	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
[2232]	113	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
[4583]	114	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	115	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
[2232]	116	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	117
[4583]	118	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
	119	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	120	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
[2232]	121	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	122	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
	123	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	124	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
[4583]	125	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical natural-type surfaces
[2232]	126	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
	127	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
	128	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical natural-type surfaces
	129	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
[4583]	130	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
[2232]	131	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	132	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	133	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
	134	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
[3636]	135
[3927]	136	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
[1]	137
[3636]	138
[3634]	139	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
	140	!$ACC PRIVATE(surf_e, surf_s, flag, mask_top, mask_bottom) &
	141	!$ACC PRIVATE(mask_north, mask_south, mask_west, mask_east) &
[4346]	142	!$ACC PRESENT(wall_flags_total_0, kh) &
[3634]	143	!$ACC PRESENT(s) &
	144	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
	145	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(0:3)) &
	146	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
	147	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) &
	148	!$ACC PRESENT(s_flux_def_h_up, s_flux_def_h_down) &
	149	!$ACC PRESENT(s_flux_t) &
	150	!$ACC PRESENT(s_flux_def_v_north, s_flux_def_v_south) &
	151	!$ACC PRESENT(s_flux_def_v_east, s_flux_def_v_west) &
	152	!$ACC PRESENT(s_flux_lsm_h_up) &
	153	!$ACC PRESENT(s_flux_lsm_v_north, s_flux_lsm_v_south) &
	154	!$ACC PRESENT(s_flux_lsm_v_east, s_flux_lsm_v_west) &
	155	!$ACC PRESENT(s_flux_usm_h_up) &
	156	!$ACC PRESENT(s_flux_usm_v_north, s_flux_usm_v_south) &
	157	!$ACC PRESENT(s_flux_usm_v_east, s_flux_usm_v_west) &
	158	!$ACC PRESENT(tend)
[1]	159	DO i = nxl, nxr
	160	DO j = nys,nyn
	161	!
	162	!-- Compute horizontal diffusion
[2232]	163	DO k = nzb+1, nzt
	164	!
	165	!-- Predetermine flag to mask topography and wall-bounded grid points
[4583]	166	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[2232]	167	!
[4583]	168	!-- Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer
	169	!-- array
[4346]	170	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
	171	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
	172	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 0 ) )
	173	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
[1]	174
[4583]	175	tend(k,j,i) = tend(k,j,i) &
	176	+ 0.5_wp * ( &
	177	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
	178	* ( s(k,j,i+1) - s(k,j,i) ) &
	179	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
	180	* ( s(k,j,i) - s(k,j,i-1) ) &
	181	) * ddx2 * flag &
	182	+ 0.5_wp * ( &
	183	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
	184	* ( s(k,j+1,i) - s(k,j,i) ) &
	185	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
	186	* ( s(k,j,i) - s(k,j-1,i) ) &
	187	) * ddy2 * flag
[1]	188	ENDDO
	189
	190	!
[4583]	191	!-- Apply prescribed horizontal wall heatflux where necessary. First, determine start and
	192	!-- end index for respective (j,i)-index. Please note, in the flat case following loop will
	193	!-- not be entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural
	194	!-- surfaces so far.
	195	!-- First, for default-type surfaces.
[2232]	196	!-- North-facing vertical default-type surfaces
	197	surf_s = surf_def_v(0)%start_index(j,i)
	198	surf_e = surf_def_v(0)%end_index(j,i)
	199	DO m = surf_s, surf_e
	200	k = surf_def_v(0)%k(m)
[2759]	201	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	202	ENDDO
	203	!
	204	!-- South-facing vertical default-type surfaces
	205	surf_s = surf_def_v(1)%start_index(j,i)
	206	surf_e = surf_def_v(1)%end_index(j,i)
	207	DO m = surf_s, surf_e
	208	k = surf_def_v(1)%k(m)
[2759]	209	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	210	ENDDO
	211	!
	212	!-- East-facing vertical default-type surfaces
	213	surf_s = surf_def_v(2)%start_index(j,i)
	214	surf_e = surf_def_v(2)%end_index(j,i)
	215	DO m = surf_s, surf_e
	216	k = surf_def_v(2)%k(m)
[2759]	217	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	218	ENDDO
	219	!
	220	!-- West-facing vertical default-type surfaces
	221	surf_s = surf_def_v(3)%start_index(j,i)
	222	surf_e = surf_def_v(3)%end_index(j,i)
	223	DO m = surf_s, surf_e
	224	k = surf_def_v(3)%k(m)
[2759]	225	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	226	ENDDO
	227	!
	228	!-- Now, for natural-type surfaces.
	229	!-- North-facing
	230	surf_s = surf_lsm_v(0)%start_index(j,i)
	231	surf_e = surf_lsm_v(0)%end_index(j,i)
	232	DO m = surf_s, surf_e
	233	k = surf_lsm_v(0)%k(m)
[2759]	234	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	235	ENDDO
	236	!
	237	!-- South-facing
	238	surf_s = surf_lsm_v(1)%start_index(j,i)
	239	surf_e = surf_lsm_v(1)%end_index(j,i)
	240	DO m = surf_s, surf_e
	241	k = surf_lsm_v(1)%k(m)
[2759]	242	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	243	ENDDO
	244	!
	245	!-- East-facing
	246	surf_s = surf_lsm_v(2)%start_index(j,i)
	247	surf_e = surf_lsm_v(2)%end_index(j,i)
	248	DO m = surf_s, surf_e
	249	k = surf_lsm_v(2)%k(m)
[2759]	250	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	251	ENDDO
	252	!
	253	!-- West-facing
	254	surf_s = surf_lsm_v(3)%start_index(j,i)
	255	surf_e = surf_lsm_v(3)%end_index(j,i)
	256	DO m = surf_s, surf_e
	257	k = surf_lsm_v(3)%k(m)
[2759]	258	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	259	ENDDO
	260	!
	261	!-- Now, for urban-type surfaces.
	262	!-- North-facing
	263	surf_s = surf_usm_v(0)%start_index(j,i)
	264	surf_e = surf_usm_v(0)%end_index(j,i)
	265	DO m = surf_s, surf_e
	266	k = surf_usm_v(0)%k(m)
[2759]	267	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	268	ENDDO
	269	!
	270	!-- South-facing
	271	surf_s = surf_usm_v(1)%start_index(j,i)
	272	surf_e = surf_usm_v(1)%end_index(j,i)
	273	DO m = surf_s, surf_e
	274	k = surf_usm_v(1)%k(m)
[2759]	275	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	276	ENDDO
	277	!
	278	!-- East-facing
	279	surf_s = surf_usm_v(2)%start_index(j,i)
	280	surf_e = surf_usm_v(2)%end_index(j,i)
	281	DO m = surf_s, surf_e
	282	k = surf_usm_v(2)%k(m)
[2759]	283	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	284	ENDDO
	285	!
	286	!-- West-facing
	287	surf_s = surf_usm_v(3)%start_index(j,i)
	288	surf_e = surf_usm_v(3)%end_index(j,i)
	289	DO m = surf_s, surf_e
	290	k = surf_usm_v(3)%k(m)
[2759]	291	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	292	ENDDO
[1]	293
	294	!
[4583]	295	!-- Compute vertical diffusion. In case that surface fluxes have been prescribed or
	296	!-- computed at bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively.
	297	!-- Model top is also mask if top flux is given.
[2232]	298	DO k = nzb+1, nzt
	299	!
[4583]	300	!-- Determine flags to mask topography below and above. Flag 0 is used to mask
	301	!-- topography in general, and flag 8 implies information about use_surface_fluxes.
	302	!-- Flag 9 is used to control flux at model top.
	303	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	304	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	305	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	306	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	307
[4583]	308	tend(k,j,i) = tend(k,j,i) &
	309	+ 0.5_wp * ( &
	310	( kh(k,j,i) + kh(k+1,j,i) ) * &
	311	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	312	* rho_air_zw(k) &
	313	* mask_top &
	314	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	315	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	316	* rho_air_zw(k-1) &
	317	* mask_bottom &
	318	) * ddzw(k) * drho_air(k) &
[2232]	319	* flag
[1]	320	ENDDO
	321
	322	!
[2232]	323	!-- Vertical diffusion at horizontal walls.
[1]	324	IF ( use_surface_fluxes ) THEN
[2232]	325	!
[4583]	326	!-- Default-type surfaces, upward-facing
[2232]	327	surf_s = surf_def_h(0)%start_index(j,i)
	328	surf_e = surf_def_h(0)%end_index(j,i)
	329	DO m = surf_s, surf_e
[1]	330
[2232]	331	k = surf_def_h(0)%k(m)
[4583]	332	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
[1]	333
[2232]	334	ENDDO
	335	!
[4583]	336	!-- Default-type surfaces, downward-facing
[2232]	337	surf_s = surf_def_h(1)%start_index(j,i)
	338	surf_e = surf_def_h(1)%end_index(j,i)
	339	DO m = surf_s, surf_e
[1]	340
[2232]	341	k = surf_def_h(1)%k(m)
[4583]	342	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
[1]	343
[2232]	344	ENDDO
[19]	345	!
[4583]	346	!-- Natural-type surfaces, upward-facing
[2232]	347	surf_s = surf_lsm_h%start_index(j,i)
	348	surf_e = surf_lsm_h%end_index(j,i)
	349	DO m = surf_s, surf_e
[19]	350
[2232]	351	k = surf_lsm_h%k(m)
[4583]	352	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
[19]	353
[2232]	354	ENDDO
	355	!
[4583]	356	!-- Urban-type surfaces, upward-facing
[2232]	357	surf_s = surf_usm_h%start_index(j,i)
	358	surf_e = surf_usm_h%end_index(j,i)
	359	DO m = surf_s, surf_e
[19]	360
[2232]	361	k = surf_usm_h%k(m)
[4583]	362	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
[2232]	363
	364	ENDDO
	365
[19]	366	ENDIF
[2232]	367	!
	368	!-- Vertical diffusion at the last computational gridpoint along z-direction
	369	IF ( use_top_fluxes ) THEN
	370	surf_s = surf_def_h(2)%start_index(j,i)
	371	surf_e = surf_def_h(2)%end_index(j,i)
	372	DO m = surf_s, surf_e
[19]	373
[2232]	374	k = surf_def_h(2)%k(m)
[4583]	375	tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
[2232]	376	ENDDO
	377	ENDIF
	378
[1]	379	ENDDO
	380	ENDDO
	381
	382	END SUBROUTINE diffusion_s
	383
[4583]	384	!--------------------------------------------------------------------------------------------------!
[1682]	385	! Description:
	386	! ------------
	387	!> Call for grid point i,j
[4583]	388	!--------------------------------------------------------------------------------------------------!
	389	SUBROUTINE diffusion_s_ij( i, j, s, &
	390	s_flux_def_h_up, s_flux_def_h_down, &
	391	s_flux_t, &
	392	s_flux_lsm_h_up, s_flux_usm_h_up, &
	393	s_flux_def_v_north, s_flux_def_v_south, &
	394	s_flux_def_v_east, s_flux_def_v_west, &
	395	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	396	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	397	s_flux_usm_v_north, s_flux_usm_v_south, &
	398	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	399
[4583]	400	USE arrays_3d, &
[2037]	401	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[4583]	402
	403	USE control_parameters, &
[1320]	404	ONLY: use_surface_fluxes, use_top_fluxes
[4583]	405
	406	USE grid_variables, &
[2759]	407	ONLY: ddx, ddx2, ddy, ddy2
[4583]	408
	409	USE indices, &
[4346]	410	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_total_0
[4583]	411
[1320]	412	USE kinds
[1]	413
[4583]	414	USE surface_mod, &
	415	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	416
[1]	417	IMPLICIT NONE
	418
[2232]	419	INTEGER(iwp) :: i !< running index x direction
	420	INTEGER(iwp) :: j !< running index y direction
	421	INTEGER(iwp) :: k !< running index z direction
	422	INTEGER(iwp) :: m !< running index surface elements
	423	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	424	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	425
	426	REAL(wp) :: flag !< flag to mask topography grid points
[4583]	427	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	428	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
[2232]	429	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
[4583]	430	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	431	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
[2232]	432	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	433
[4583]	434	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
	435	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	436	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
[2232]	437	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	438	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
	439	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	440	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
[4583]	441	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical urban-type surfaces
[2232]	442	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
	443	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
	444	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical urban-type surfaces
	445	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
[4583]	446	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
[2232]	447	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	448	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	449	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
	450	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
[3636]	451
[3927]	452	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
[1]	453
[4583]	454
[1]	455	!
	456	!-- Compute horizontal diffusion
[2232]	457	DO k = nzb+1, nzt
	458	!
	459	!-- Predetermine flag to mask topography and wall-bounded grid points
[4583]	460	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[2232]	461	!
[4583]	462	!-- Predetermine flag to mask wall-bounded grid points, equivalent to former s_outer array
[4346]	463	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
	464	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
	465	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 0 ) )
	466	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
[2232]	467	!
[4583]	468	!-- Finally, determine flag to mask both topography itself as well as wall-bounded grid
	469	!-- points, which will be treated further below
[1]	470
[4583]	471	tend(k,j,i) = tend(k,j,i) &
	472	+ 0.5_wp * ( &
	473	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
	474	* ( s(k,j,i+1) - s(k,j,i) ) &
	475	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
	476	* ( s(k,j,i) - s(k,j,i-1) ) &
	477	) * ddx2 * flag &
	478	+ 0.5_wp * ( &
	479	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
	480	* ( s(k,j+1,i) - s(k,j,i) ) &
	481	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
	482	* ( s(k,j,i) - s(k,j-1,i) ) &
[2232]	483	) * ddy2 * flag
[1]	484	ENDDO
	485
	486	!
[4583]	487	!-- Apply prescribed horizontal wall heatflux where necessary. First, determine start and end
	488	!-- index for respective (j,i)-index. Please note, in the flat case following loops will not be
	489	!-- entered, as surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces so far.
	490	!-- First, for default-type surfaces.
[2232]	491	!-- North-facing vertical default-type surfaces
	492	surf_s = surf_def_v(0)%start_index(j,i)
	493	surf_e = surf_def_v(0)%end_index(j,i)
	494	DO m = surf_s, surf_e
	495	k = surf_def_v(0)%k(m)
[2759]	496	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	497	ENDDO
	498	!
	499	!-- South-facing vertical default-type surfaces
	500	surf_s = surf_def_v(1)%start_index(j,i)
	501	surf_e = surf_def_v(1)%end_index(j,i)
	502	DO m = surf_s, surf_e
	503	k = surf_def_v(1)%k(m)
[2759]	504	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	505	ENDDO
	506	!
	507	!-- East-facing vertical default-type surfaces
	508	surf_s = surf_def_v(2)%start_index(j,i)
	509	surf_e = surf_def_v(2)%end_index(j,i)
	510	DO m = surf_s, surf_e
	511	k = surf_def_v(2)%k(m)
[2759]	512	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	513	ENDDO
	514	!
	515	!-- West-facing vertical default-type surfaces
	516	surf_s = surf_def_v(3)%start_index(j,i)
	517	surf_e = surf_def_v(3)%end_index(j,i)
	518	DO m = surf_s, surf_e
	519	k = surf_def_v(3)%k(m)
[2759]	520	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	521	ENDDO
	522	!
	523	!-- Now, for natural-type surfaces
	524	!-- North-facing
	525	surf_s = surf_lsm_v(0)%start_index(j,i)
	526	surf_e = surf_lsm_v(0)%end_index(j,i)
	527	DO m = surf_s, surf_e
	528	k = surf_lsm_v(0)%k(m)
[2759]	529	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	530	ENDDO
	531	!
	532	!-- South-facing
	533	surf_s = surf_lsm_v(1)%start_index(j,i)
	534	surf_e = surf_lsm_v(1)%end_index(j,i)
	535	DO m = surf_s, surf_e
	536	k = surf_lsm_v(1)%k(m)
[2759]	537	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	538	ENDDO
	539	!
	540	!-- East-facing
	541	surf_s = surf_lsm_v(2)%start_index(j,i)
	542	surf_e = surf_lsm_v(2)%end_index(j,i)
	543	DO m = surf_s, surf_e
	544	k = surf_lsm_v(2)%k(m)
[2759]	545	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	546	ENDDO
	547	!
	548	!-- West-facing
	549	surf_s = surf_lsm_v(3)%start_index(j,i)
	550	surf_e = surf_lsm_v(3)%end_index(j,i)
	551	DO m = surf_s, surf_e
	552	k = surf_lsm_v(3)%k(m)
[2759]	553	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	554	ENDDO
	555	!
	556	!-- Now, for urban-type surfaces
	557	!-- North-facing
	558	surf_s = surf_usm_v(0)%start_index(j,i)
	559	surf_e = surf_usm_v(0)%end_index(j,i)
	560	DO m = surf_s, surf_e
	561	k = surf_usm_v(0)%k(m)
[2759]	562	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	563	ENDDO
	564	!
	565	!-- South-facing
	566	surf_s = surf_usm_v(1)%start_index(j,i)
	567	surf_e = surf_usm_v(1)%end_index(j,i)
	568	DO m = surf_s, surf_e
	569	k = surf_usm_v(1)%k(m)
[2759]	570	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	571	ENDDO
	572	!
	573	!-- East-facing
	574	surf_s = surf_usm_v(2)%start_index(j,i)
	575	surf_e = surf_usm_v(2)%end_index(j,i)
	576	DO m = surf_s, surf_e
	577	k = surf_usm_v(2)%k(m)
[2759]	578	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	579	ENDDO
	580	!
	581	!-- West-facing
	582	surf_s = surf_usm_v(3)%start_index(j,i)
	583	surf_e = surf_usm_v(3)%end_index(j,i)
	584	DO m = surf_s, surf_e
	585	k = surf_usm_v(3)%k(m)
[2759]	586	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	587	ENDDO
[1]	588
	589
	590	!
[4583]	591	!-- Compute vertical diffusion. In case that surface fluxes have been prescribed or computed at
	592	!-- bottom and/or top, index k starts/ends at nzb+2 or nzt-1, respectively. Model top is also
	593	!-- mask if top flux is given.
[2232]	594	DO k = nzb+1, nzt
	595	!
[4583]	596	!-- Determine flags to mask topography below and above. Flag 0 is used to mask topography in
	597	!-- general, and flag 8 implies information about use_surface_fluxes. Flag 9 is used to
	598	!-- control flux at model top.
	599	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	600	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	601	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	602	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	603
[4583]	604	tend(k,j,i) = tend(k,j,i) &
	605	+ 0.5_wp * ( &
	606	( kh(k,j,i) + kh(k+1,j,i) ) * &
	607	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
	608	* rho_air_zw(k) &
	609	* mask_top &
	610	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	611	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
	612	* rho_air_zw(k-1) &
	613	* mask_bottom &
	614	) * ddzw(k) * drho_air(k) &
[2232]	615	* flag
[1]	616	ENDDO
	617
	618	!
[2232]	619	!-- Vertical diffusion at horizontal walls.
	620	!-- TO DO: Adjust for downward facing walls and mask already in main loop
[1]	621	IF ( use_surface_fluxes ) THEN
[2232]	622	!
	623	!-- Default-type surfaces, upward-facing
	624	surf_s = surf_def_h(0)%start_index(j,i)
	625	surf_e = surf_def_h(0)%end_index(j,i)
	626	DO m = surf_s, surf_e
[1]	627
[2232]	628	k = surf_def_h(0)%k(m)
[1]	629
[4583]	630	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) * ddzw(k) * drho_air(k)
[2232]	631	ENDDO
	632	!
	633	!-- Default-type surfaces, downward-facing
	634	surf_s = surf_def_h(1)%start_index(j,i)
	635	surf_e = surf_def_h(1)%end_index(j,i)
	636	DO m = surf_s, surf_e
[1]	637
[2232]	638	k = surf_def_h(1)%k(m)
[1]	639
[4583]	640	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) * ddzw(k) * drho_air(k)
[2232]	641	ENDDO
[19]	642	!
[2232]	643	!-- Natural-type surfaces, upward-facing
	644	surf_s = surf_lsm_h%start_index(j,i)
	645	surf_e = surf_lsm_h%end_index(j,i)
	646	DO m = surf_s, surf_e
	647	k = surf_lsm_h%k(m)
	648
[4583]	649	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) * ddzw(k) * drho_air(k)
[2232]	650	ENDDO
	651	!
	652	!-- Urban-type surfaces, upward-facing
	653	surf_s = surf_usm_h%start_index(j,i)
	654	surf_e = surf_usm_h%end_index(j,i)
	655	DO m = surf_s, surf_e
	656	k = surf_usm_h%k(m)
	657
[4583]	658	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) * ddzw(k) * drho_air(k)
[2232]	659	ENDDO
	660	ENDIF
	661	!
[19]	662	!-- Vertical diffusion at the last computational gridpoint along z-direction
	663	IF ( use_top_fluxes ) THEN
[2232]	664	surf_s = surf_def_h(2)%start_index(j,i)
	665	surf_e = surf_def_h(2)%end_index(j,i)
	666	DO m = surf_s, surf_e
[19]	667
[2232]	668	k = surf_def_h(2)%k(m)
[4583]	669	tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
[2232]	670	ENDDO
[19]	671	ENDIF
	672
[1]	673	END SUBROUTINE diffusion_s_ij
	674
	675	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |