Home

Context Navigation

diffusion_v.f90 @ 4896

Last change on this file since 4896 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: 25.3 KB

Rev	Line
[1873]	1	!> @file diffusion_v.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:
[1]	20	! -----------------
[4671]	21	!
	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: diffusion_v.f90 4828 2021-01-05 11:21:41Z raasch $
[4674]	26	! Update ACC directives for downward facing USM and LSM surfaces
	27	!
	28	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	29	! Implementation of downward facing USM and LSM surfaces
	30	!
	31	! 4583 2020-06-29 12:36:47Z raasch
[4583]	32	! file re-formatted to follow the PALM coding standard
	33	!
	34	! 4360 2020-01-07 11:25:50Z suehring
	35	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	36	! information used in wall_flags_static_0
	37	!
[4346]	38	! 4329 2019-12-10 15:46:36Z motisi
[4329]	39	! Renamed wall_flags_0 to wall_flags_static_0
[4583]	40	!
[4329]	41	! 4182 2019-08-22 15:20:23Z scharf
[4182]	42	! Corrected "Former revisions" section
[4583]	43	!
[4182]	44	! 3655 2019-01-07 16:51:22Z knoop
[3634]	45	! OpenACC port for SPEC
[2716]	46	!
[4182]	47	! Revision 1.1 1997/09/12 06:24:01 raasch
	48	! Initial revision
	49	!
	50	!
[1]	51	! Description:
	52	! ------------
[1682]	53	!> Diffusion term of the v-component
[4583]	54	!--------------------------------------------------------------------------------------------------!
[1682]	55	MODULE diffusion_v_mod
[1]	56
[4583]	57
[1]	58	PRIVATE
[2118]	59	PUBLIC diffusion_v
[1]	60
	61	INTERFACE diffusion_v
	62	MODULE PROCEDURE diffusion_v
	63	MODULE PROCEDURE diffusion_v_ij
	64	END INTERFACE diffusion_v
	65
	66	CONTAINS
	67
	68
[4583]	69	!--------------------------------------------------------------------------------------------------!
[1682]	70	! Description:
	71	! ------------
	72	!> Call for all grid points
[4583]	73	!--------------------------------------------------------------------------------------------------!
[1001]	74	SUBROUTINE diffusion_v
[1]	75
[4583]	76	USE arrays_3d, &
	77	ONLY: ddzu, ddzw, drho_air, km, rho_air_zw, tend, u, v, w
	78
	79	USE control_parameters, &
	80	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	81
	82	USE grid_variables, &
[2232]	83	ONLY: ddx, ddy, ddy2
[4583]	84
	85	USE indices, &
[4346]	86	ONLY: nxl, nxr, nyn, nysv, nzb, nzt, wall_flags_total_0
[4583]	87
[1320]	88	USE kinds
[1]	89
[4583]	90	USE surface_mod, &
	91	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	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) :: l !< running index of surface type, south- or north-facing wall
	99	INTEGER(iwp) :: m !< running index surface elements
	100	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	101	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1001]	102
[2232]	103	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	104	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
	105	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
	106	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto yv-zw grid
	107	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto yv-zw grid
[4583]	108	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	109	REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point
	110	REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point
	111	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
[1]	112
[3634]	113	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
	114	!$ACC PRIVATE(surf_e, surf_s, flag, kmxm, kmxp, kmzm, kmzp) &
	115	!$ACC PRIVATE(mask_bottom, mask_east, mask_west, mask_top) &
[4346]	116	!$ACC PRESENT(wall_flags_total_0, km) &
[3634]	117	!$ACC PRESENT(u, v, w) &
	118	!$ACC PRESENT(ddzu, ddzw, drho_air, rho_air_zw) &
	119	!$ACC PRESENT(surf_def_h(0:2), surf_def_v(2:3)) &
[4674]	120	!$ACC PRESENT(surf_lsm_h(0:1), surf_lsm_v(2:3)) &
	121	!$ACC PRESENT(surf_usm_h(0:1), surf_usm_v(2:3)) &
[3634]	122	!$ACC PRESENT(tend)
[1]	123	DO i = nxl, nxr
[106]	124	DO j = nysv, nyn
[1]	125	!
	126	!-- Compute horizontal diffusion
[2232]	127	DO k = nzb+1, nzt
	128
[1]	129	!
[4583]	130	!-- Predetermine flag to mask topography and wall-bounded grid points.
	131	!-- It is sufficient to masked only east- and west-facing surfaces, which need special
	132	!-- treatment for the v-component.
	133	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[4346]	134	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
	135	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
[2232]	136	!
[1]	137	!-- Interpolate eddy diffusivities on staggered gridpoints
[2232]	138	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
	139	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	140
[4583]	141	tend(k,j,i) = tend(k,j,i) + ( &
	142	mask_east * kmxp * ( &
	143	( v(k,j,i+1) - v(k,j,i) ) * ddx &
	144	+ ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	145	) &
	146	- mask_west * kmxm * ( &
	147	( v(k,j,i) - v(k,j,i-1) ) * ddx &
	148	+ ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	149	) &
	150	) * ddx * flag &
	151	+ 2.0_wp * ( &
	152	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	153	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
	154	) * ddy2 * flag
[2232]	155
[1]	156	ENDDO
	157
	158	!
[4583]	159	!-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note,
	160	!-- in the the flat case, loops won't be entered as start_index > end_index. Furtermore,
	161	!-- note, no vertical natural surfaces so far.
[2232]	162	!-- Default-type surfaces
	163	DO l = 2, 3
	164	surf_s = surf_def_v(l)%start_index(j,i)
	165	surf_e = surf_def_v(l)%end_index(j,i)
	166	DO m = surf_s, surf_e
	167	k = surf_def_v(l)%k(m)
[4583]	168	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
	169	ENDDO
[2232]	170	ENDDO
[1]	171	!
[2232]	172	!-- Natural-type surfaces
	173	DO l = 2, 3
	174	surf_s = surf_lsm_v(l)%start_index(j,i)
	175	surf_e = surf_lsm_v(l)%end_index(j,i)
	176	DO m = surf_s, surf_e
	177	k = surf_lsm_v(l)%k(m)
[4583]	178	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
	179	ENDDO
[2232]	180	ENDDO
[1]	181	!
[2232]	182	!-- Urban-type surfaces
	183	DO l = 2, 3
	184	surf_s = surf_usm_v(l)%start_index(j,i)
	185	surf_e = surf_usm_v(l)%end_index(j,i)
	186	DO m = surf_s, surf_e
	187	k = surf_usm_v(l)%k(m)
[4583]	188	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
	189	ENDDO
[2232]	190	ENDDO
	191	!
[4583]	192	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid
	193	!-- diffusive fluxes are masked (flag 10) within this loop, and added further below, else,
	194	!-- simple gradient approach is applied. Model top is also mask if top-momentum flux is
	195	!-- given.
[2232]	196	DO k = nzb+1, nzt
	197	!
[4583]	198	!-- Determine flags to mask topography below and above. Flag 2 is used to mask
	199	!-- topography in general, while flag 8 implies also information about
	200	!-- use_surface_fluxes. Flag 9 is used to control momentum flux at model top.
	201	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	202	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	203	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	204	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2232]	205	!
[1]	206	!-- Interpolate eddy diffusivities on staggered gridpoints
[4583]	207	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	208	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]	209
[4583]	210	tend(k,j,i) = tend(k,j,i) &
	211	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	212	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	213	& ) * rho_air_zw(k) * mask_top &
	214	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	215	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	216	& ) * rho_air_zw(k-1) * mask_bottom &
	217	& ) * ddzw(k) * drho_air(k) * flag
[1]	218	ENDDO
	219
	220	!
[4583]	221	!-- Vertical diffusion at the first grid point above the surface, if the momentum flux at
	222	!-- the bottom is given by the Prandtl law or if it is prescribed by the user.
	223	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
	224	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the
	225	!-- values of the momentum flux becomes too large in this case.
[1]	226	IF ( use_surface_fluxes ) THEN
	227	!
[2232]	228	!-- Default-type surfaces, upward-facing
	229	surf_s = surf_def_h(0)%start_index(j,i)
	230	surf_e = surf_def_h(0)%end_index(j,i)
	231	DO m = surf_s, surf_e
	232	k = surf_def_h(0)%k(m)
[1]	233
[4583]	234	tend(k,j,i) = tend(k,j,i) &
	235	+ ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	236	ENDDO
	237	!
	238	!-- Default-type surfaces, dowward-facing
	239	surf_s = surf_def_h(1)%start_index(j,i)
	240	surf_e = surf_def_h(1)%end_index(j,i)
	241	DO m = surf_s, surf_e
	242	k = surf_def_h(1)%k(m)
[1]	243
[4583]	244	tend(k,j,i) = tend(k,j,i) &
	245	+ ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]	246	ENDDO
[102]	247	!
[2232]	248	!-- Natural-type surfaces, upward-facing
[4671]	249	surf_s = surf_lsm_h(0)%start_index(j,i)
	250	surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]	251	DO m = surf_s, surf_e
[4671]	252	k = surf_lsm_h(0)%k(m)
[2232]	253
[4583]	254	tend(k,j,i) = tend(k,j,i) &
[4671]	255	+ ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	256
	257	ENDDO
[102]	258	!
[4671]	259	!-- Natural-type surfaces, downward-facing
	260	surf_s = surf_lsm_h(1)%start_index(j,i)
	261	surf_e = surf_lsm_h(1)%end_index(j,i)
	262	DO m = surf_s, surf_e
	263	k = surf_lsm_h(1)%k(m)
	264
	265	tend(k,j,i) = tend(k,j,i) &
	266	+ ( - surf_lsm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
	267
	268	ENDDO
	269	!
[2232]	270	!-- Urban-type surfaces, upward-facing
[4671]	271	surf_s = surf_usm_h(0)%start_index(j,i)
	272	surf_e = surf_usm_h(0)%end_index(j,i)
[2232]	273	DO m = surf_s, surf_e
[4671]	274	k = surf_usm_h(0)%k(m)
[102]	275
[4583]	276	tend(k,j,i) = tend(k,j,i) &
[4671]	277	+ ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	278
	279	ENDDO
[4671]	280	!
	281	!-- Urban-type surfaces, downward-facing
	282	surf_s = surf_usm_h(1)%start_index(j,i)
	283	surf_e = surf_usm_h(1)%end_index(j,i)
	284	DO m = surf_s, surf_e
	285	k = surf_usm_h(1)%k(m)
	286
	287	tend(k,j,i) = tend(k,j,i) &
	288	+ ( - surf_usm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
	289
	290	ENDDO
[102]	291	ENDIF
[2232]	292	!
	293	!-- Add momentum flux at model top
[2638]	294	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[2232]	295	surf_s = surf_def_h(2)%start_index(j,i)
	296	surf_e = surf_def_h(2)%end_index(j,i)
	297	DO m = surf_s, surf_e
[102]	298
[2232]	299	k = surf_def_h(2)%k(m)
	300
[4583]	301	tend(k,j,i) = tend(k,j,i) &
	302	+ ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]	303	ENDDO
	304	ENDIF
	305
[1]	306	ENDDO
	307	ENDDO
	308
	309	END SUBROUTINE diffusion_v
	310
	311
[4583]	312	!--------------------------------------------------------------------------------------------------!
[1682]	313	! Description:
	314	! ------------
	315	!> Call for grid point i,j
[4583]	316	!--------------------------------------------------------------------------------------------------!
[1001]	317	SUBROUTINE diffusion_v_ij( i, j )
[1]	318
[4583]	319	USE arrays_3d, &
	320	ONLY: ddzu, ddzw, drho_air, km, tend, u, v, w, rho_air_zw
	321
	322	USE control_parameters, &
	323	ONLY: constant_top_momentumflux, use_surface_fluxes, use_top_fluxes
	324
	325	USE grid_variables, &
[2232]	326	ONLY: ddx, ddy, ddy2
[4583]	327
	328	USE indices, &
[4346]	329	ONLY: nzb, nzt, wall_flags_total_0
[4583]	330
[1320]	331	USE kinds
[1]	332
[4583]	333	USE surface_mod, &
	334	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v
[2232]	335
[1]	336	IMPLICIT NONE
	337
	338
[2232]	339	INTEGER(iwp) :: i !< running index x direction
	340	INTEGER(iwp) :: j !< running index y direction
	341	INTEGER(iwp) :: k !< running index z direction
	342	INTEGER(iwp) :: l !< running index of surface type, south- or north-facing wall
	343	INTEGER(iwp) :: m !< running index surface elements
	344	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
	345	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
[1001]	346
[2232]	347	REAL(wp) :: flag !< flag to mask topography grid points
[3547]	348	REAL(wp) :: kmxm !< diffusion coefficient on leftward side of the v-gridbox - interpolated onto xu-yv grid
	349	REAL(wp) :: kmxp !< diffusion coefficient on rightward side of the v-gridbox - interpolated onto xu-yv grid
	350	REAL(wp) :: kmzm !< diffusion coefficient on bottom of the gridbox - interpolated onto xu-zw grid
	351	REAL(wp) :: kmzp !< diffusion coefficient on top of the gridbox - interpolated onto xu-zw grid
[4583]	352	REAL(wp) :: mask_bottom !< flag to mask vertical upward-facing surface
	353	REAL(wp) :: mask_east !< flag to mask vertical surface south of the grid point
	354	REAL(wp) :: mask_west !< flag to mask vertical surface north of the grid point
[2232]	355	REAL(wp) :: mask_top !< flag to mask vertical downward-facing surface
	356
[1]	357	!
	358	!-- Compute horizontal diffusion
[2232]	359	DO k = nzb+1, nzt
[1]	360	!
[4583]	361	!-- Predetermine flag to mask topography and wall-bounded grid points.
	362	!-- It is sufficient to masked only east- and west-facing surfaces, which need special
	363	!-- treatment for the v-component.
	364	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[4346]	365	mask_east = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i+1), 2 ) )
	366	mask_west = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i-1), 2 ) )
[2232]	367	!
[1]	368	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	369	kmxp = 0.25_wp * ( km(k,j,i)+km(k,j,i+1)+km(k,j-1,i)+km(k,j-1,i+1) )
	370	kmxm = 0.25_wp * ( km(k,j,i)+km(k,j,i-1)+km(k,j-1,i)+km(k,j-1,i-1) )
[1]	371
[4583]	372	tend(k,j,i) = tend(k,j,i) + ( &
	373	mask_east * kmxp * ( &
	374	( v(k,j,i+1) - v(k,j,i) ) * ddx &
	375	+ ( u(k,j,i+1) - u(k,j-1,i+1) ) * ddy &
	376	) &
	377	- mask_west * kmxm * ( &
	378	( v(k,j,i) - v(k,j,i-1) ) * ddx &
	379	+ ( u(k,j,i) - u(k,j-1,i) ) * ddy &
	380	) &
	381	) * ddx * flag &
	382	+ 2.0_wp * ( &
	383	km(k,j,i) * ( v(k,j+1,i) - v(k,j,i) ) &
	384	- km(k,j-1,i) * ( v(k,j,i) - v(k,j-1,i) ) &
[2232]	385	) * ddy2 * flag
[1]	386	ENDDO
	387
	388	!
[4583]	389	!-- Add horizontal momentum flux v'u' at east- (l=2) and west-facing (l=3) surfaces. Note, in the
	390	!-- the flat case, loops won't be entered as start_index > end_index. Furtermore, note, no
	391	!-- vertical natural surfaces so far.
[2232]	392	!-- Default-type surfaces
	393	DO l = 2, 3
	394	surf_s = surf_def_v(l)%start_index(j,i)
	395	surf_e = surf_def_v(l)%end_index(j,i)
	396	DO m = surf_s, surf_e
	397	k = surf_def_v(l)%k(m)
	398	tend(k,j,i) = tend(k,j,i) + surf_def_v(l)%mom_flux_uv(m) * ddx
[4583]	399	ENDDO
[2232]	400	ENDDO
[51]	401	!
[2232]	402	!-- Natural-type surfaces
	403	DO l = 2, 3
	404	surf_s = surf_lsm_v(l)%start_index(j,i)
	405	surf_e = surf_lsm_v(l)%end_index(j,i)
	406	DO m = surf_s, surf_e
	407	k = surf_lsm_v(l)%k(m)
	408	tend(k,j,i) = tend(k,j,i) + surf_lsm_v(l)%mom_flux_uv(m) * ddx
[4583]	409	ENDDO
[2232]	410	ENDDO
[1]	411	!
[2232]	412	!-- Urban-type surfaces
	413	DO l = 2, 3
	414	surf_s = surf_usm_v(l)%start_index(j,i)
	415	surf_e = surf_usm_v(l)%end_index(j,i)
	416	DO m = surf_s, surf_e
	417	k = surf_usm_v(l)%k(m)
	418	tend(k,j,i) = tend(k,j,i) + surf_usm_v(l)%mom_flux_uv(m) * ddx
[4583]	419	ENDDO
[2232]	420	ENDDO
[1]	421	!
[4583]	422	!-- Compute vertical diffusion. In case of simulating a surface layer, respective grid diffusive
	423	!-- fluxes are masked (flag 8) within this loop, and added further below, else, simple gradient
	424	!-- approach is applied. Model top is also mask if top-momentum flux is given.
[2232]	425	DO k = nzb+1, nzt
	426	!
[4583]	427	!-- Determine flags to mask topography below and above. Flag 2 is used to mask topography in
	428	!-- general, while flag 10 implies also information about use_surface_fluxes. Flag 9 is used
	429	!-- to control momentum flux at model top.
	430	mask_bottom = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k-1,j,i), 8 ) )
	431	mask_top = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 8 ) ) * &
	432	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k+1,j,i), 9 ) )
	433	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[2232]	434	!
[1]	435	!-- Interpolate eddy diffusivities on staggered gridpoints
[1340]	436	kmzp = 0.25_wp * ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
	437	kmzm = 0.25_wp * ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
[1]	438
[4583]	439	tend(k,j,i) = tend(k,j,i) &
	440	& + ( kmzp * ( ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	441	& + ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	442	& ) * rho_air_zw(k) * mask_top &
	443	& - kmzm * ( ( v(k,j,i) - v(k-1,j,i) ) * ddzu(k) &
	444	& + ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy &
	445	& ) * rho_air_zw(k-1) * mask_bottom &
	446	& ) * ddzw(k) * drho_air(k) * flag
[1]	447	ENDDO
	448
	449	!
[4583]	450	!-- Vertical diffusion at the first grid point above the surface, if the momentum flux at the
	451	!-- bottom is given by the Prandtl law or if it is prescribed by the user.
	452	!-- Difference quotient of the momentum flux is not formed over half of the grid spacing
	453	!-- (2.0*ddzw(k)) any more, since the comparison with other (LES) models showed that the values
	454	!-- of the momentum flux becomes too large in this case.
[1]	455	IF ( use_surface_fluxes ) THEN
	456	!
[2232]	457	!-- Default-type surfaces, upward-facing
	458	surf_s = surf_def_h(0)%start_index(j,i)
	459	surf_e = surf_def_h(0)%end_index(j,i)
	460	DO m = surf_s, surf_e
	461	k = surf_def_h(0)%k(m)
[1]	462
[4583]	463	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_def_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	464	ENDDO
	465	!
	466	!-- Default-type surfaces, dowward-facing
	467	surf_s = surf_def_h(1)%start_index(j,i)
	468	surf_e = surf_def_h(1)%end_index(j,i)
	469	DO m = surf_s, surf_e
	470	k = surf_def_h(1)%k(m)
[1]	471
[4583]	472	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]	473	ENDDO
[102]	474	!
[2232]	475	!-- Natural-type surfaces, upward-facing
[4671]	476	surf_s = surf_lsm_h(0)%start_index(j,i)
	477	surf_e = surf_lsm_h(0)%end_index(j,i)
[2232]	478	DO m = surf_s, surf_e
[4671]	479	k = surf_lsm_h(0)%k(m)
[2232]	480
[4671]	481	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_lsm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	482
	483	ENDDO
[102]	484	!
[4671]	485	!-- Natural-type surfaces, downward-facing
	486	surf_s = surf_lsm_h(1)%start_index(j,i)
	487	surf_e = surf_lsm_h(1)%end_index(j,i)
	488	DO m = surf_s, surf_e
	489	k = surf_lsm_h(1)%k(m)
	490
	491	tend(k,j,i) = tend(k,j,i) + ( - surf_lsm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
	492
	493	ENDDO
	494	!
[2232]	495	!-- Urban-type surfaces, upward-facing
[4671]	496	surf_s = surf_usm_h(0)%start_index(j,i)
	497	surf_e = surf_usm_h(0)%end_index(j,i)
[2232]	498	DO m = surf_s, surf_e
[4671]	499	k = surf_usm_h(0)%k(m)
[102]	500
[4671]	501	tend(k,j,i) = tend(k,j,i) + ( - ( - surf_usm_h(0)%vsws(m) ) ) * ddzw(k) * drho_air(k)
[2232]	502
	503	ENDDO
[4671]	504	!
	505	!-- Urban-type surfaces, downward-facing
	506	surf_s = surf_usm_h(1)%start_index(j,i)
	507	surf_e = surf_usm_h(1)%end_index(j,i)
	508	DO m = surf_s, surf_e
	509	k = surf_usm_h(1)%k(m)
	510
	511	tend(k,j,i) = tend(k,j,i) + ( - surf_usm_h(1)%vsws(m) ) * ddzw(k) * drho_air(k)
	512
	513	ENDDO
[102]	514	ENDIF
[2232]	515	!
	516	!-- Add momentum flux at model top
[2638]	517	IF ( use_top_fluxes .AND. constant_top_momentumflux ) THEN
[2232]	518	surf_s = surf_def_h(2)%start_index(j,i)
	519	surf_e = surf_def_h(2)%end_index(j,i)
	520	DO m = surf_s, surf_e
[102]	521
[2232]	522	k = surf_def_h(2)%k(m)
	523
[4583]	524	tend(k,j,i) = tend(k,j,i) + ( - surf_def_h(2)%vsws(m) ) * ddzw(k) * drho_air(k)
[2232]	525	ENDDO
	526	ENDIF
	527
[1]	528	END SUBROUTINE diffusion_v_ij
	529
[1321]	530	END MODULE diffusion_v_mod

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

Context Navigation

source: palm/trunk/SOURCE/diffusion_v.f90 @ 4896

Download in other formats: