Home

Context Navigation

diffusion_s.f90 @ 4867

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

Context Navigation

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

Download in other formats: