Home

Context Navigation

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

Last change on this file since 4464 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 33.1 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	!
[4360]	17	! Copyright 1997-2020 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 4360 2020-01-07 11:25:50Z Giersch $
[4346]	27	! Introduction of wall_flags_total_0, which currently sets bits based on static
	28	! topography information used in wall_flags_static_0
[4329]	29	!
[4346]	30	! 4329 2019-12-10 15:46:36Z motisi
	31	! Renamed wall_flags_0 to wall_flags_total_0
	32	!
[4329]	33	! 4182 2019-08-22 15:20:23Z scharf
[4182]	34	! Corrected "Former revisions" section
	35	!
	36	! 3927 2019-04-23 13:24:29Z raasch
[3927]	37	! pointer attribute removed from scalar 3d-array for performance reasons
	38	!
	39	! 3761 2019-02-25 15:31:42Z raasch
[3761]	40	! unused variables removed
	41	!
	42	! 3655 2019-01-07 16:51:22Z knoop
[3636]	43	! nopointer option removed
[1321]	44	!
[4182]	45	! Revision 1.1 2000/04/13 14:54:02 schroeter
	46	! Initial revision
	47	!
	48	!
[1]	49	! Description:
	50	! ------------
[1682]	51	!> Diffusion term of scalar quantities (temperature and water content)
[1]	52	!------------------------------------------------------------------------------!
[1682]	53	MODULE diffusion_s_mod
	54
[1]	55
	56	PRIVATE
[2118]	57	PUBLIC diffusion_s
[1]	58
	59	INTERFACE diffusion_s
	60	MODULE PROCEDURE diffusion_s
	61	MODULE PROCEDURE diffusion_s_ij
	62	END INTERFACE diffusion_s
	63
	64	CONTAINS
	65
	66
	67	!------------------------------------------------------------------------------!
[1682]	68	! Description:
	69	! ------------
	70	!> Call for all grid points
[1]	71	!------------------------------------------------------------------------------!
[2232]	72	SUBROUTINE diffusion_s( s, s_flux_def_h_up, s_flux_def_h_down, &
	73	s_flux_t, &
	74	s_flux_lsm_h_up, s_flux_usm_h_up, &
	75	s_flux_def_v_north, s_flux_def_v_south, &
	76	s_flux_def_v_east, s_flux_def_v_west, &
	77	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	78	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	79	s_flux_usm_v_north, s_flux_usm_v_south, &
	80	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	81
[1320]	82	USE arrays_3d, &
[2037]	83	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[1320]	84
	85	USE control_parameters, &
	86	ONLY: use_surface_fluxes, use_top_fluxes
	87
	88	USE grid_variables, &
[2759]	89	ONLY: ddx, ddx2, ddy, ddy2
[1320]	90
	91	USE indices, &
[4346]	92	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
	93	wall_flags_total_0
[1320]	94
	95	USE kinds
[1]	96
[2232]	97	USE surface_mod, &
	98	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	99	surf_usm_v
	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
	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
	113	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	114	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	115	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	116	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	117
	118	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	119	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical 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
	121	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	122	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	123	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing 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
	125	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical natural-type surfaces
	126	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical natural-type surfaces
	127	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-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
	130	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	131	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	132	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-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
[4346]	166	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[2232]	167	!
	168	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
	169	!-- former s_outer 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
[1320]	175	tend(k,j,i) = tend(k,j,i) &
[1340]	176	+ 0.5_wp * ( &
[2232]	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 &
[1340]	182	+ 0.5_wp * ( &
[2232]	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	!
[2232]	191	!-- Apply prescribed horizontal wall heatflux where necessary. First,
	192	!-- determine start and end index for respective (j,i)-index. Please
	193	!-- note, in the flat case following loop will not be entered, as
	194	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
	195	!-- so far.
	196	!-- First, for default-type surfaces
	197	!-- North-facing vertical default-type surfaces
	198	surf_s = surf_def_v(0)%start_index(j,i)
	199	surf_e = surf_def_v(0)%end_index(j,i)
	200	DO m = surf_s, surf_e
	201	k = surf_def_v(0)%k(m)
[2759]	202	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	203	ENDDO
	204	!
	205	!-- South-facing vertical default-type surfaces
	206	surf_s = surf_def_v(1)%start_index(j,i)
	207	surf_e = surf_def_v(1)%end_index(j,i)
	208	DO m = surf_s, surf_e
	209	k = surf_def_v(1)%k(m)
[2759]	210	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	211	ENDDO
	212	!
	213	!-- East-facing vertical default-type surfaces
	214	surf_s = surf_def_v(2)%start_index(j,i)
	215	surf_e = surf_def_v(2)%end_index(j,i)
	216	DO m = surf_s, surf_e
	217	k = surf_def_v(2)%k(m)
[2759]	218	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	219	ENDDO
	220	!
	221	!-- West-facing vertical default-type surfaces
	222	surf_s = surf_def_v(3)%start_index(j,i)
	223	surf_e = surf_def_v(3)%end_index(j,i)
	224	DO m = surf_s, surf_e
	225	k = surf_def_v(3)%k(m)
[2759]	226	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	227	ENDDO
	228	!
	229	!-- Now, for natural-type surfaces.
	230	!-- North-facing
	231	surf_s = surf_lsm_v(0)%start_index(j,i)
	232	surf_e = surf_lsm_v(0)%end_index(j,i)
	233	DO m = surf_s, surf_e
	234	k = surf_lsm_v(0)%k(m)
[2759]	235	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	236	ENDDO
	237	!
	238	!-- South-facing
	239	surf_s = surf_lsm_v(1)%start_index(j,i)
	240	surf_e = surf_lsm_v(1)%end_index(j,i)
	241	DO m = surf_s, surf_e
	242	k = surf_lsm_v(1)%k(m)
[2759]	243	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	244	ENDDO
	245	!
	246	!-- East-facing
	247	surf_s = surf_lsm_v(2)%start_index(j,i)
	248	surf_e = surf_lsm_v(2)%end_index(j,i)
	249	DO m = surf_s, surf_e
	250	k = surf_lsm_v(2)%k(m)
[2759]	251	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	252	ENDDO
	253	!
	254	!-- West-facing
	255	surf_s = surf_lsm_v(3)%start_index(j,i)
	256	surf_e = surf_lsm_v(3)%end_index(j,i)
	257	DO m = surf_s, surf_e
	258	k = surf_lsm_v(3)%k(m)
[2759]	259	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	260	ENDDO
	261	!
	262	!-- Now, for urban-type surfaces.
	263	!-- North-facing
	264	surf_s = surf_usm_v(0)%start_index(j,i)
	265	surf_e = surf_usm_v(0)%end_index(j,i)
	266	DO m = surf_s, surf_e
	267	k = surf_usm_v(0)%k(m)
[2759]	268	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	269	ENDDO
	270	!
	271	!-- South-facing
	272	surf_s = surf_usm_v(1)%start_index(j,i)
	273	surf_e = surf_usm_v(1)%end_index(j,i)
	274	DO m = surf_s, surf_e
	275	k = surf_usm_v(1)%k(m)
[2759]	276	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	277	ENDDO
	278	!
	279	!-- East-facing
	280	surf_s = surf_usm_v(2)%start_index(j,i)
	281	surf_e = surf_usm_v(2)%end_index(j,i)
	282	DO m = surf_s, surf_e
	283	k = surf_usm_v(2)%k(m)
[2759]	284	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	285	ENDDO
	286	!
	287	!-- West-facing
	288	surf_s = surf_usm_v(3)%start_index(j,i)
	289	surf_e = surf_usm_v(3)%end_index(j,i)
	290	DO m = surf_s, surf_e
	291	k = surf_usm_v(3)%k(m)
[2759]	292	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	293	ENDDO
[1]	294
	295	!
	296	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	297	!-- prescribed or computed at bottom and/or top, index k starts/ends at
[2232]	298	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
	299	!-- is given.
	300	DO k = nzb+1, nzt
	301	!
	302	!-- Determine flags to mask topography below and above. Flag 0 is
	303	!-- used to mask topography in general, and flag 8 implies
	304	!-- information about use_surface_fluxes. Flag 9 is used to control
	305	!-- flux at model top.
	306	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
[4346]	307	BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
[2232]	308	mask_top = MERGE( 1.0_wp, 0.0_wp, &
[4346]	309	BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
[2232]	310	MERGE( 1.0_wp, 0.0_wp, &
[4346]	311	BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
[2232]	312	flag = MERGE( 1.0_wp, 0.0_wp, &
[4346]	313	BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	314
[1320]	315	tend(k,j,i) = tend(k,j,i) &
[1340]	316	+ 0.5_wp * ( &
[2232]	317	( kh(k,j,i) + kh(k+1,j,i) ) * &
	318	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
[2037]	319	* rho_air_zw(k) &
[2232]	320	* mask_top &
	321	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	322	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
[2037]	323	* rho_air_zw(k-1) &
[2232]	324	* mask_bottom &
	325	) * ddzw(k) * drho_air(k) &
	326	* flag
[1]	327	ENDDO
	328
	329	!
[2232]	330	!-- Vertical diffusion at horizontal walls.
[1]	331	IF ( use_surface_fluxes ) THEN
[2232]	332	!
	333	!-- Default-type surfaces, upward-facing
	334	surf_s = surf_def_h(0)%start_index(j,i)
	335	surf_e = surf_def_h(0)%end_index(j,i)
	336	DO m = surf_s, surf_e
[1]	337
[2232]	338	k = surf_def_h(0)%k(m)
	339	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
	340	* ddzw(k) * drho_air(k)
[1]	341
[2232]	342	ENDDO
	343	!
	344	!-- Default-type surfaces, downward-facing
	345	surf_s = surf_def_h(1)%start_index(j,i)
	346	surf_e = surf_def_h(1)%end_index(j,i)
	347	DO m = surf_s, surf_e
[1]	348
[2232]	349	k = surf_def_h(1)%k(m)
	350	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
	351	* ddzw(k) * drho_air(k)
[1]	352
[2232]	353	ENDDO
[19]	354	!
[2232]	355	!-- Natural-type surfaces, upward-facing
	356	surf_s = surf_lsm_h%start_index(j,i)
	357	surf_e = surf_lsm_h%end_index(j,i)
	358	DO m = surf_s, surf_e
[19]	359
[2232]	360	k = surf_lsm_h%k(m)
	361	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
	362	* ddzw(k) * drho_air(k)
[19]	363
[2232]	364	ENDDO
	365	!
	366	!-- Urban-type surfaces, upward-facing
	367	surf_s = surf_usm_h%start_index(j,i)
	368	surf_e = surf_usm_h%end_index(j,i)
	369	DO m = surf_s, surf_e
[19]	370
[2232]	371	k = surf_usm_h%k(m)
	372	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
	373	* ddzw(k) * drho_air(k)
	374
	375	ENDDO
	376
[19]	377	ENDIF
[2232]	378	!
	379	!-- Vertical diffusion at the last computational gridpoint along z-direction
	380	IF ( use_top_fluxes ) THEN
	381	surf_s = surf_def_h(2)%start_index(j,i)
	382	surf_e = surf_def_h(2)%end_index(j,i)
	383	DO m = surf_s, surf_e
[19]	384
[2232]	385	k = surf_def_h(2)%k(m)
	386	tend(k,j,i) = tend(k,j,i) &
	387	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
	388	ENDDO
	389	ENDIF
	390
[1]	391	ENDDO
	392	ENDDO
	393
	394	END SUBROUTINE diffusion_s
	395
	396	!------------------------------------------------------------------------------!
[1682]	397	! Description:
	398	! ------------
	399	!> Call for grid point i,j
[1]	400	!------------------------------------------------------------------------------!
[2232]	401	SUBROUTINE diffusion_s_ij( i, j, s, &
	402	s_flux_def_h_up, s_flux_def_h_down, &
	403	s_flux_t, &
	404	s_flux_lsm_h_up, s_flux_usm_h_up, &
	405	s_flux_def_v_north, s_flux_def_v_south, &
	406	s_flux_def_v_east, s_flux_def_v_west, &
	407	s_flux_lsm_v_north, s_flux_lsm_v_south, &
	408	s_flux_lsm_v_east, s_flux_lsm_v_west, &
	409	s_flux_usm_v_north, s_flux_usm_v_south, &
	410	s_flux_usm_v_east, s_flux_usm_v_west )
[1]	411
[1320]	412	USE arrays_3d, &
[2037]	413	ONLY: ddzu, ddzw, kh, tend, drho_air, rho_air_zw
[1320]	414
	415	USE control_parameters, &
	416	ONLY: use_surface_fluxes, use_top_fluxes
	417
	418	USE grid_variables, &
[2759]	419	ONLY: ddx, ddx2, ddy, ddy2
[1320]	420
	421	USE indices, &
[4346]	422	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzt, wall_flags_total_0
[1320]	423
	424	USE kinds
[1]	425
[2232]	426	USE surface_mod, &
	427	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
	428	surf_usm_v
	429
[1]	430	IMPLICIT NONE
	431
[2232]	432	INTEGER(iwp) :: i !< running index x direction
	433	INTEGER(iwp) :: j !< running index y direction
	434	INTEGER(iwp) :: k !< running index z direction
	435	INTEGER(iwp) :: m !< running index surface elements
	436	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	437	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
	438
	439	REAL(wp) :: flag !< flag to mask topography grid points
	440	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	441	REAL(wp) :: mask_east !< flag to mask vertical surface east of the grid point
	442	REAL(wp) :: mask_north !< flag to mask vertical surface north of the grid point
	443	REAL(wp) :: mask_south !< flag to mask vertical surface south of the grid point
	444	REAL(wp) :: mask_west !< flag to mask vertical surface west of the grid point
	445	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	446
	447	REAL(wp), DIMENSION(1:surf_def_v(0)%ns) :: s_flux_def_v_north !< flux at north-facing vertical default-type surfaces
	448	REAL(wp), DIMENSION(1:surf_def_v(1)%ns) :: s_flux_def_v_south !< flux at south-facing vertical default-type surfaces
	449	REAL(wp), DIMENSION(1:surf_def_v(2)%ns) :: s_flux_def_v_east !< flux at east-facing vertical default-type surfaces
	450	REAL(wp), DIMENSION(1:surf_def_v(3)%ns) :: s_flux_def_v_west !< flux at west-facing vertical default-type surfaces
	451	REAL(wp), DIMENSION(1:surf_def_h(0)%ns) :: s_flux_def_h_up !< flux at horizontal upward-facing default-type surfaces
	452	REAL(wp), DIMENSION(1:surf_def_h(1)%ns) :: s_flux_def_h_down !< flux at horizontal donwward-facing default-type surfaces
	453	REAL(wp), DIMENSION(1:surf_lsm_h%ns) :: s_flux_lsm_h_up !< flux at horizontal upward-facing natural-type surfaces
	454	REAL(wp), DIMENSION(1:surf_lsm_v(0)%ns) :: s_flux_lsm_v_north !< flux at north-facing vertical urban-type surfaces
	455	REAL(wp), DIMENSION(1:surf_lsm_v(1)%ns) :: s_flux_lsm_v_south !< flux at south-facing vertical urban-type surfaces
	456	REAL(wp), DIMENSION(1:surf_lsm_v(2)%ns) :: s_flux_lsm_v_east !< flux at east-facing vertical urban-type surfaces
	457	REAL(wp), DIMENSION(1:surf_lsm_v(3)%ns) :: s_flux_lsm_v_west !< flux at west-facing vertical urban-type surfaces
	458	REAL(wp), DIMENSION(1:surf_usm_h%ns) :: s_flux_usm_h_up !< flux at horizontal upward-facing urban-type surfaces
	459	REAL(wp), DIMENSION(1:surf_usm_v(0)%ns) :: s_flux_usm_v_north !< flux at north-facing vertical urban-type surfaces
	460	REAL(wp), DIMENSION(1:surf_usm_v(1)%ns) :: s_flux_usm_v_south !< flux at south-facing vertical urban-type surfaces
	461	REAL(wp), DIMENSION(1:surf_usm_v(2)%ns) :: s_flux_usm_v_east !< flux at east-facing vertical urban-type surfaces
	462	REAL(wp), DIMENSION(1:surf_usm_v(3)%ns) :: s_flux_usm_v_west !< flux at west-facing vertical urban-type surfaces
	463	REAL(wp), DIMENSION(1:surf_def_h(2)%ns) :: s_flux_t !< flux at model top
[3636]	464
[3927]	465	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s !< treated scalar
[1]	466
	467	!
	468	!-- Compute horizontal diffusion
[2232]	469	DO k = nzb+1, nzt
	470	!
	471	!-- Predetermine flag to mask topography and wall-bounded grid points
[4346]	472	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[2232]	473	!
	474	!-- Predetermine flag to mask wall-bounded grid points, equivalent to
	475	!-- former s_outer array
[4346]	476	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 0 ) )
	477	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 0 ) )
	478	mask_south = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j-1,i), 0 ) )
	479	mask_north = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j+1,i), 0 ) )
[2232]	480	!
	481	!-- Finally, determine flag to mask both topography itself as well
	482	!-- as wall-bounded grid points, which will be treated further below
[1]	483
[1320]	484	tend(k,j,i) = tend(k,j,i) &
[1340]	485	+ 0.5_wp * ( &
[2232]	486	mask_east * ( kh(k,j,i) + kh(k,j,i+1) ) &
	487	* ( s(k,j,i+1) - s(k,j,i) ) &
	488	- mask_west * ( kh(k,j,i) + kh(k,j,i-1) ) &
	489	* ( s(k,j,i) - s(k,j,i-1) ) &
	490	) * ddx2 * flag &
[1340]	491	+ 0.5_wp * ( &
[2232]	492	mask_north * ( kh(k,j,i) + kh(k,j+1,i) ) &
	493	* ( s(k,j+1,i) - s(k,j,i) ) &
	494	- mask_south * ( kh(k,j,i) + kh(k,j-1,i) ) &
	495	* ( s(k,j,i) - s(k,j-1,i) ) &
	496	) * ddy2 * flag
[1]	497	ENDDO
	498
	499	!
[2232]	500	!-- Apply prescribed horizontal wall heatflux where necessary. First,
	501	!-- determine start and end index for respective (j,i)-index. Please
	502	!-- note, in the flat case following loops will not be entered, as
	503	!-- surf_s=1 and surf_e=0. Furtermore, note, no vertical natural surfaces
	504	!-- so far.
	505	!-- First, for default-type surfaces
	506	!-- North-facing vertical default-type surfaces
	507	surf_s = surf_def_v(0)%start_index(j,i)
	508	surf_e = surf_def_v(0)%end_index(j,i)
	509	DO m = surf_s, surf_e
	510	k = surf_def_v(0)%k(m)
[2759]	511	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_north(m) * ddy
[2232]	512	ENDDO
	513	!
	514	!-- South-facing vertical default-type surfaces
	515	surf_s = surf_def_v(1)%start_index(j,i)
	516	surf_e = surf_def_v(1)%end_index(j,i)
	517	DO m = surf_s, surf_e
	518	k = surf_def_v(1)%k(m)
[2759]	519	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_south(m) * ddy
[2232]	520	ENDDO
	521	!
	522	!-- East-facing vertical default-type surfaces
	523	surf_s = surf_def_v(2)%start_index(j,i)
	524	surf_e = surf_def_v(2)%end_index(j,i)
	525	DO m = surf_s, surf_e
	526	k = surf_def_v(2)%k(m)
[2759]	527	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_east(m) * ddx
[2232]	528	ENDDO
	529	!
	530	!-- West-facing vertical default-type surfaces
	531	surf_s = surf_def_v(3)%start_index(j,i)
	532	surf_e = surf_def_v(3)%end_index(j,i)
	533	DO m = surf_s, surf_e
	534	k = surf_def_v(3)%k(m)
[2759]	535	tend(k,j,i) = tend(k,j,i) + s_flux_def_v_west(m) * ddx
[2232]	536	ENDDO
	537	!
	538	!-- Now, for natural-type surfaces
	539	!-- North-facing
	540	surf_s = surf_lsm_v(0)%start_index(j,i)
	541	surf_e = surf_lsm_v(0)%end_index(j,i)
	542	DO m = surf_s, surf_e
	543	k = surf_lsm_v(0)%k(m)
[2759]	544	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_north(m) * ddy
[2232]	545	ENDDO
	546	!
	547	!-- South-facing
	548	surf_s = surf_lsm_v(1)%start_index(j,i)
	549	surf_e = surf_lsm_v(1)%end_index(j,i)
	550	DO m = surf_s, surf_e
	551	k = surf_lsm_v(1)%k(m)
[2759]	552	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_south(m) * ddy
[2232]	553	ENDDO
	554	!
	555	!-- East-facing
	556	surf_s = surf_lsm_v(2)%start_index(j,i)
	557	surf_e = surf_lsm_v(2)%end_index(j,i)
	558	DO m = surf_s, surf_e
	559	k = surf_lsm_v(2)%k(m)
[2759]	560	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_east(m) * ddx
[2232]	561	ENDDO
	562	!
	563	!-- West-facing
	564	surf_s = surf_lsm_v(3)%start_index(j,i)
	565	surf_e = surf_lsm_v(3)%end_index(j,i)
	566	DO m = surf_s, surf_e
	567	k = surf_lsm_v(3)%k(m)
[2759]	568	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_v_west(m) * ddx
[2232]	569	ENDDO
	570	!
	571	!-- Now, for urban-type surfaces
	572	!-- North-facing
	573	surf_s = surf_usm_v(0)%start_index(j,i)
	574	surf_e = surf_usm_v(0)%end_index(j,i)
	575	DO m = surf_s, surf_e
	576	k = surf_usm_v(0)%k(m)
[2759]	577	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_north(m) * ddy
[2232]	578	ENDDO
	579	!
	580	!-- South-facing
	581	surf_s = surf_usm_v(1)%start_index(j,i)
	582	surf_e = surf_usm_v(1)%end_index(j,i)
	583	DO m = surf_s, surf_e
	584	k = surf_usm_v(1)%k(m)
[2759]	585	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_south(m) * ddy
[2232]	586	ENDDO
	587	!
	588	!-- East-facing
	589	surf_s = surf_usm_v(2)%start_index(j,i)
	590	surf_e = surf_usm_v(2)%end_index(j,i)
	591	DO m = surf_s, surf_e
	592	k = surf_usm_v(2)%k(m)
[2759]	593	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_east(m) * ddx
[2232]	594	ENDDO
	595	!
	596	!-- West-facing
	597	surf_s = surf_usm_v(3)%start_index(j,i)
	598	surf_e = surf_usm_v(3)%end_index(j,i)
	599	DO m = surf_s, surf_e
	600	k = surf_usm_v(3)%k(m)
[2759]	601	tend(k,j,i) = tend(k,j,i) + s_flux_usm_v_west(m) * ddx
[2232]	602	ENDDO
[1]	603
	604
	605	!
	606	!-- Compute vertical diffusion. In case that surface fluxes have been
[19]	607	!-- prescribed or computed at bottom and/or top, index k starts/ends at
[2232]	608	!-- nzb+2 or nzt-1, respectively. Model top is also mask if top flux
	609	!-- is given.
	610	DO k = nzb+1, nzt
	611	!
	612	!-- Determine flags to mask topography below and above. Flag 0 is
	613	!-- used to mask topography in general, and flag 8 implies
	614	!-- information about use_surface_fluxes. Flag 9 is used to control
	615	!-- flux at model top.
	616	mask_bottom = MERGE( 1.0_wp, 0.0_wp, &
[4346]	617	BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
[2232]	618	mask_top = MERGE( 1.0_wp, 0.0_wp, &
[4346]	619	BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
[2232]	620	MERGE( 1.0_wp, 0.0_wp, &
[4346]	621	BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
[2232]	622	flag = MERGE( 1.0_wp, 0.0_wp, &
[4346]	623	BTEST( wall_flags_total_0(k,j,i), 0 ) )
[1]	624
[1320]	625	tend(k,j,i) = tend(k,j,i) &
[1340]	626	+ 0.5_wp * ( &
[2232]	627	( kh(k,j,i) + kh(k+1,j,i) ) * &
	628	( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
[2037]	629	* rho_air_zw(k) &
[2232]	630	* mask_top &
	631	- ( kh(k,j,i) + kh(k-1,j,i) ) * &
	632	( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
[2037]	633	* rho_air_zw(k-1) &
[2232]	634	* mask_bottom &
	635	) * ddzw(k) * drho_air(k) &
	636	* flag
[1]	637	ENDDO
	638
	639	!
[2232]	640	!-- Vertical diffusion at horizontal walls.
	641	!-- TO DO: Adjust for downward facing walls and mask already in main loop
[1]	642	IF ( use_surface_fluxes ) THEN
[2232]	643	!
	644	!-- Default-type surfaces, upward-facing
	645	surf_s = surf_def_h(0)%start_index(j,i)
	646	surf_e = surf_def_h(0)%end_index(j,i)
	647	DO m = surf_s, surf_e
[1]	648
[2232]	649	k = surf_def_h(0)%k(m)
[1]	650
[2232]	651	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_up(m) &
	652	* ddzw(k) * drho_air(k)
	653	ENDDO
	654	!
	655	!-- Default-type surfaces, downward-facing
	656	surf_s = surf_def_h(1)%start_index(j,i)
	657	surf_e = surf_def_h(1)%end_index(j,i)
	658	DO m = surf_s, surf_e
[1]	659
[2232]	660	k = surf_def_h(1)%k(m)
[1]	661
[2232]	662	tend(k,j,i) = tend(k,j,i) + s_flux_def_h_down(m) &
	663	* ddzw(k) * drho_air(k)
	664	ENDDO
[19]	665	!
[2232]	666	!-- Natural-type surfaces, upward-facing
	667	surf_s = surf_lsm_h%start_index(j,i)
	668	surf_e = surf_lsm_h%end_index(j,i)
	669	DO m = surf_s, surf_e
	670	k = surf_lsm_h%k(m)
	671
	672	tend(k,j,i) = tend(k,j,i) + s_flux_lsm_h_up(m) &
	673	* ddzw(k) * drho_air(k)
	674	ENDDO
	675	!
	676	!-- Urban-type surfaces, upward-facing
	677	surf_s = surf_usm_h%start_index(j,i)
	678	surf_e = surf_usm_h%end_index(j,i)
	679	DO m = surf_s, surf_e
	680	k = surf_usm_h%k(m)
	681
	682	tend(k,j,i) = tend(k,j,i) + s_flux_usm_h_up(m) &
	683	* ddzw(k) * drho_air(k)
	684	ENDDO
	685	ENDIF
	686	!
[19]	687	!-- Vertical diffusion at the last computational gridpoint along z-direction
	688	IF ( use_top_fluxes ) THEN
[2232]	689	surf_s = surf_def_h(2)%start_index(j,i)
	690	surf_e = surf_def_h(2)%end_index(j,i)
	691	DO m = surf_s, surf_e
[19]	692
[2232]	693	k = surf_def_h(2)%k(m)
	694	tend(k,j,i) = tend(k,j,i) &
	695	+ ( - s_flux_t(m) ) * ddzw(k) * drho_air(k)
	696	ENDDO
[19]	697	ENDIF
	698
[1]	699	END SUBROUTINE diffusion_s_ij
	700
	701	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |