Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |